2.6 Hyperelasticity

Constitutive Restrictions

Now, we need to make constitutive assumptions on the dependent variables which appears in Clausius-Duhem inequality, eq.(2.4-9), namely, \(\psi\), \(\eta\), \(\boldsymbol{\sigma}\) and \(\mathbf{q}\). These functions could depend on a number of state variables, such as strain, rate of deformation, temperature, temperature gradient, etc. For example, the pressure in a gas is generally assumed to depend on temperature and density, so it would follow that \(\boldsymbol{\sigma}\) for a gas should exhibit the same dependencies.

For solid materials we have found that that the mass density is dependent on the strain (eq.(2.5-35), \(\rho=\rho_{r}\left(\det\mathbf{C}\right)^{-1/2}=\rho_{r}/J\), where \(\rho_{r}\) is a constant); so long as a dependence on the strain is provided, an explicit dependency on mass density is not required.

For the purpose of modeling elastic solids under isothermal conditions, we would like to restrict our constitutive models to be independent of the temperature or its gradient. For an elastic solid, the rate of load application does not have an effect; therefore, there should be no dependency on the rate of deformation \(\mathbf{D}\). Furthermore, for an elastic solid, upon removal of the loading, deformations should disappear completely. This implies that there should be no dependence on the history (or path) of loading, only on the current state of loading.

Given these restrictions, we assume that the non-observable functions of state \(\psi\), \(\eta\), \(\boldsymbol{\sigma}\) and \(\mathbf{q}\) depend only on the observable solid matrix strain, such as \(\mathbf{E}\). Based on the principle of equipresence, all functions of state depend only on \(\mathbf{E}\). In particular, since \(\psi\) is the only function of state which is differentiated with respect to time in the Clausius-Duhem inequality, eq.(2.4-9), we express it explicitly as

\[ \begin{equation} \psi=\psi\left(\mathbf{E}\right)\,.\label{eq95-1} \end{equation} \]

Given this dependence, we can use the chain rule of differentiation to get

\[ \begin{equation} \dot{\psi}=\frac{d\psi}{d\mathbf{E}}:\dot{\mathbf{E}}=\frac{d\psi}{d\mathbf{E}}:\mathbf{F}^{T}\cdot\mathbf{D}\cdot\mathbf{F}\label{eq96-1} \end{equation} \]

where we have made use of eq.(2.5-33). Making use of the identity

\[ \begin{equation} \frac{d\psi}{d\mathbf{E}}:\mathbf{F}^{T}\cdot\mathbf{D}\cdot\mathbf{F}=\mathbf{F}\cdot\frac{d\psi}{d\mathbf{E}}\cdot\mathbf{F}^{T}:\mathbf{D}\label{eq100-1} \end{equation} \]

we find that the entropy inequality reduces to

\[ \begin{equation} -\rho\eta\dot{\theta}-\frac{1}{\theta}\mathbf{q}\cdot\grad\theta+\left(\boldsymbol{\sigma}-\rho\mathbf{F}\cdot\frac{d\psi}{d\mathbf{E}}\cdot\mathbf{F}^{T}\right):\mathbf{D}\geqslant0\label{eq101-1} \end{equation} \]

This inequality must hold for arbitrary processes, i.e., for arbitrary changes in \(\mathbf{D}\), \(\dot{\theta}\) and \(\grad\theta\) (all of which are independent, observable variables of state), under our self-imposed constraint that \(\psi\), \(\eta\), \(\boldsymbol{\sigma}\) and \(\mathbf{q}\) only depend on the state variable \(\mathbf{E}\). For example, since the temperature \(\theta\) in the leftmost term of eq.\eqref{eq101-1} may either increase or decrease at various times in any process, the sign of \(\dot{\theta}\) may be variably positive or negative. However, the sign of the product \(-\rho\eta\dot{\theta}\) may not be set strictly positive since \(\rho\eta\) does not depend on \(\dot{\theta}\). The same argument applies to the remaining terms in the entropy inequality, each of which may vary independently of other terms. Therefore, the only way that the entropy inequality may be satisfied for arbitrary processes, given our choice of state variables, is to have each term be independently equal to zero. Thus, the entropy inequality requires that

\[ \begin{equation} \eta=0\label{eq102-1} \end{equation} \]

\[ \begin{equation} \boldsymbol{\sigma}=\rho\mathbf{F}\cdot\frac{d\psi}{d\mathbf{E}}\cdot\mathbf{F}^{T}\label{eq103-1} \end{equation} \]

\[ \begin{equation} \mathbf{q}=\mathbf{0}\label{eq105-1} \end{equation} \]

Equation \eqref{eq102-1} shows that the entropy must be zero for our choice of state variable, whereas Eq.\eqref{eq103-1} is the fundamental constraint we seek for formulating constitutive relations in elasticity problems. The constraint of eq.\eqref{eq105-1} implies that isothermal processes must be adiabatic, as there cannot be heat flow in such problems.

Let the Helmholtz free energy per unit reference volume (energy density) be denoted by

\[ \begin{equation} \Psi_{r}\left(\mathbf{E}\right)\equiv\rho_{r}\psi\left(\mathbf{E}\right)\label{eq106-1} \end{equation} \]

Then, using eq.(2.5-35) and recalling that \(\rho_{r}\) is constant, eq.\eqref{eq103-1} reduces to

\[ \begin{equation} \boldsymbol{\sigma}=J^{-1}\mathbf{F}\cdot\frac{d\Psi_{r}}{d\mathbf{E}}\cdot\mathbf{F}^{T}\,.\label{eq:hyperelastic-stress-E} \end{equation} \]

This formula represents the constraint imposed by the second law of thermodynamics on the formulation of a constitutive relation for an elastic solid. For historical reasons, this formula is described as the hyperelasticity constraint. The method described for obtaining this formula was first proposed by Coleman and Noll in 1963 ¹. Prior to that formulation, other constitutive relations had been proposed for elastic solids (such as hypoelasticity), though it is now believed that the only thermodynamically valid formulation is the one given here.

In hyperelasticity theory, since the free energy density \(\Psi_{r}\) only depends on the strain, we may also call it the strain energy density. A common alternative symbol for the strain energy density is \(W\left(\mathbf{E}\right)\). Since it is common to include multiple state variables in a general formulation of \(\Psi_{r}\) for a material, the derivative \(d\Psi_{r}/d\mathbf{E}\) appearing in eq.\eqref{eq103-1} is often written as a partial derivative, \(\partial\Psi_{r}/\partial\mathbf{E}\). It is also common to use \(\mathbf{C}\) instead of \(\mathbf{E}\) as a suitable alternative measure of strain under finite deformation and rotation. In that case, based on eq.(2.5-21), the expression of eq.\eqref{eq103-1} takes the form

\[ \begin{equation} \boldsymbol{\sigma}=2J^{-1}\mathbf{F}\cdot\frac{\partial\Psi_{r}}{\partial\mathbf{C}}\cdot\mathbf{F}^{T}\,.\label{eq:hyperelastic-stress-C} \end{equation} \]

Other Stress Tensors

The reaction force \(\mathbf{f}\) exerted on a surface \(S\) whose outward normal is \(\mathbf{n}\) is given by

\[ \mathbf{f}=\int_{S}\mathbf{t}\,da=\int_{S}\boldsymbol{\sigma}\cdot\mathbf{n}\,da \]

where \(da\) is an elemental area of \(S\) in the current configuration. Here, we used the relation between the traction vector \(\mathbf{t}\) and the Cauchy stress \(\boldsymbol{\sigma}\), given in eq.(2.3-1). Note that the directed area \(d\mathbf{a}=\mathbf{n}\,da\) in this integral may be related to \(d\mathbf{A}\) in the reference configuration, using Nanson's formula in eq.(2.5-46) to produce

\[ \mathbf{f}=\int_{S_{r}}\boldsymbol{\sigma}\cdot J\mathbf{F}^{-T}\cdot d\mathbf{A}=\int_{S_{r}}J\boldsymbol{\sigma}\cdot\mathbf{F}^{-T}\cdot\mathbf{n}_{r}\,dA\,. \]

Here, \(\mathbf{n}_{r}\) is the unit outward normal on \(S_{r}\) in the reference configuration. Equating these two forms of \(\mathbf{f}\) and using the divergence theorem produces

\[ \mathbf{f}=\int_{V}\divg\boldsymbol{\sigma}\,dV=\int_{V_{r}}\Divg\left(J\boldsymbol{\sigma}\cdot\mathbf{F}^{-T}\right)dV_{r}\,, \]

where the divergence operator in the material frame is denoted by \(\Divg\left(\cdot\right)\). The expression in the argument of this divergence operator is known as the first Piola-Kirchhoff stress tensor,

\[ \begin{equation} \mathbf{P}=J\boldsymbol{\sigma}\cdot\mathbf{F}^{-T}\,.\label{eq50} \end{equation} \]

Note that \(\mathbf{P}\), like \(\mathbf{F}\), is not symmetric. Also, like \(\mathbf{F}\), \(\mathbf{P}\) is known as a two-point tensor, meaning it is neither a material nor a spatial tensor. This is demonstrated more easily by expressing \(\boldsymbol{\sigma}\) in terms of its contravariant components \(\sigma^{ij}\),

\[ \begin{equation} \boldsymbol{\sigma}=\sigma^{ij}\mathbf{g}_{i}\otimes\mathbf{g}_{j}\label{eq:Cauchy-cov-cont} \end{equation} \]

and noting that

\[ \begin{equation} \mathbf{F}^{-1}=\mathbf{G}_{i}\otimes\mathbf{g}^{i}\label{eq:Finv-cov-cont} \end{equation} \]

so that the component form of \(\mathbf{P}\) in eq.\eqref{eq50} reduces to

\[ \mathbf{P}=J\sigma^{ij}\left(\mathbf{g}_{i}\otimes\mathbf{g}_{j}\right)\cdot\left(\mathbf{g}^{k}\otimes\mathbf{G}_{k}\right)=J\sigma^{ij}\mathbf{g}_{i}\otimes\mathbf{G}_{j}\,, \]

confirming that it is a two-point tensor.

Two other stress measures are often used in continuum mechanics. Examining eq.\eqref{eq:hyperelastic-stress-E} or eq.\eqref{eq:hyperelastic-stress-C}, we may define the second Piola-Kirchhoff (PK2) stress tensor as

\[ \begin{equation} \mathbf{S}=J\mathbf{F}^{-1}\cdot\boldsymbol{\sigma}\cdot\mathbf{F}^{-T}\,,\label{eq51} \end{equation} \]

where, for a hyperelastic material,

\[ \begin{equation} \mathbf{S}=\frac{\partial\Psi_{r}}{\partial\mathbf{E}}=2\frac{\partial\Psi_{r}}{\partial\mathbf{C}}\,.\label{eq:PK2-hyperelasticity} \end{equation} \]

Substituting the component form of \(\boldsymbol{\sigma}\) in eq.\eqref{eq:Cauchy-cov-cont} into eq.\eqref{eq51}, and using \eqref{eq:Finv-cov-cont}, we find that

\[ \begin{equation} \mathbf{S}=J\sigma^{ij}\mathbf{G}_{i}\otimes\mathbf{G}_{j}\,,\label{eq:PK2-cov-cont} \end{equation} \]

showing that \(\mathbf{S}\) is a material tensor. This component form suggests that may also define the Kirchhoff stress tensor as the spatial version of \(\mathbf{S}\),

\[ \begin{equation} \boldsymbol{\tau}=J\boldsymbol{\sigma}=J\sigma^{ij}\mathbf{g}_{i}\otimes\mathbf{g}_{j}\,.\label{eq49} \end{equation} \]

Since \(\boldsymbol{\sigma}\) is symmetric, it follows that \(\boldsymbol{\tau}\) is symmetric.

The inverse relations for these various stresses are are

\[ \begin{equation} \boldsymbol{\mathbf{\sigma}}=\frac{1}{J}\boldsymbol{\tau}\,,\quad\boldsymbol{\sigma}=\frac{1}{J}\mathbf{P}\cdot\mathbf{F}^{T}\,,\quad\boldsymbol{\sigma}=\frac{1}{J}\mathbf{F}\cdot\mathbf{S}\cdot\mathbf{F}^{T}\,.\label{eq52} \end{equation} \]

Directional Derivative of the Stress

The directional derivative of the 2\(^{\mathrm{nd}}\) PK stress tensor needs to be calculated for the linearization of the finite element equations. For a hyperelastic material, as expressed in \eqref{eq:PK2-hyperelasticity}, a linear relationship between the directional derivative of \(\mathbf{S}\) and the linearized strain \(D\mathbf{E}\left[\mathbf{u}\right]\) can be obtained:

\[ \begin{equation} D\mathbf{S}\left[\mathbf{u}\right]=\boldsymbol{\mathbb{C}}:D\mathbf{E}\left[\mathbf{u}\right]\,.\label{eq54} \end{equation} \]

Here, \(\boldsymbol{\mathbb{C}}\) is a fourth-order tensor known as the material elasticity tensor. Its Cartesian components and tensorial definition are given by

\[ \begin{equation} \begin{aligned}\mathbb{C}_{IJKL} & =\frac{\partial S_{IJ}}{\partial E_{KL}}=\frac{4\partial^{2}\Psi_{r}}{\partial C_{IJ}\partial C_{KL}}\\ \boldsymbol{\mathbb{C}} & =\frac{\partial\mathbf{S}}{\partial\mathbf{E}}=\frac{4\partial^{2}\Psi_{r}}{\partial\mathbf{C}\partial\mathbf{C}} \end{aligned} \,.\label{eq55} \end{equation} \]

The spatial equivalent – the spatial elasticity tensor \(\boldsymbol{\mathcal{C}}\) – can be obtained from

\[ \begin{equation} \begin{aligned}\mathcal{C}_{ijkl} & =\frac{1}{J}F_{iI}F_{jJ}F_{kK}F_{lL}\mathbb{C}_{IJKL}\\ \boldsymbol{\mathcal{C}} & =\frac{1}{J}\left(\mathbf{F}\oslash\mathbf{F}\right):\boldsymbol{\mathbb{C}}:\left(\mathbf{F}^{T}\oslash\mathbf{F}^{T}\right) \end{aligned}.\label{eq56} \end{equation} \]

Isotropic Hyperelasticity

The hyperelastic constitutive equations discussed so far are unrestricted in their application. Isotropic material symmetry is defined by requiring the constitutive behavior to be independent of any material axes. Consequently, \(\Psi_{r}\) must only be a function of the invariants of \(\mathbf{C}\) (or \(\mathbf{E}\)),

\[ \begin{equation} \Psi_{r}=\Psi_{r}\left(I_{1},I_{2},I_{3}\right)\,,\label{eq:isotropic-SED} \end{equation} \]

where the invariants of \(\mathbf{C}\) are defined here as,

\[ \begin{equation} I_{1}=\tr\mathbf{C}=\mathbf{I}:\mathbf{C},\,I_{2}=\frac{1}{2}\left[\left(\tr\mathbf{C}\right)^{2}-\tr\mathbf{C}^{2}\right],\,I_{3}=\det\mathbf{C}=J^{2}\,.\label{eq:C-invariants} \end{equation} \]

As a result of the isotropic restriction, the second Piola-Kirchhoff stress tensor can be written as,

\[ \begin{equation} \mathbf{S}=2\frac{\partial\Psi}{\partial\mathbf{C}}=2\frac{\partial\Psi}{\partial I_{1}}\frac{\partial I_{1}}{\partial\mathbf{C}}+2\frac{\partial\Psi}{\partial I_{2}}\frac{\partial I_{2}}{\partial\mathbf{C}}+2\frac{\partial\Psi}{\partial I_{3}}\frac{\partial I_{3}}{\partial\mathbf{C}}\,.\label{eq:PK2-isotropic} \end{equation} \]

The second order tensors formed by the derivatives of the invariants with respect to \(\mathbf{C}\) can be evaluated as

\[ \begin{equation} \frac{\partial I_{1}}{\partial\mathbf{C}}=\mathbf{I}\,,\,\frac{\partial I_{2}}{\partial\mathbf{C}}=I_{1}\mathbf{I}-\mathbf{C}\,,\,\frac{\partial I_{3}}{\partial\mathbf{C}}=I_{3}\mathbf{C}^{-1}\,.\label{eq:C-invariant-deriv} \end{equation} \]

Introducing expressions \eqref{eq:C-invariant-deriv} into equation \eqref{eq:PK2-isotropic} enables the second Piola-Kirchhoff stress to be evaluated as,

\[ \begin{equation} \mathbf{S}=2\left(\Psi_{1}+I_{1}\Psi_{2}+I_{2}\Psi_{3})\right)\mathbf{I}-2\left(\Psi_{2}+I_{1}\Psi_{3}\right)\mathbf{C}+\Psi_{3}\mathbf{C}^{2}\,,\label{eq:PK2-isotropic-final} \end{equation} \]

where \(\Psi_{1}=\partial\Psi/\partial I_{1}\), \(\Psi_{2}=\partial\Psi/\partial I_{2}\), and \(\Psi_{3}=\partial\Psi/\partial I_{3}\).

The Cauchy stresses can now be obtained from the second Piola-Kirchhoff stresses by using \eqref{eq52}:

\[ \begin{equation} \boldsymbol{\sigma}=J^{-1}\left(2\left(\Psi_{1}+I_{1}\Psi_{2}+I_{2}\Psi_{3})\right)\mathbf{b}-2\left(\Psi_{2}+I_{1}\Psi_{3}\right)\mathbf{b}^{2}+\Psi_{3}\mathbf{b}^{3}\right)\,,\label{eq:Cauchy-isotropic} \end{equation} \]

where

\[ \begin{equation} \mathbf{b}=\mathbf{F}\cdot\mathbf{F}^{T}\,.\label{eq:left-Cauchy-Green} \end{equation} \]

is the left Cauchy-Green tensor. Note that in this equation \(\Psi_{1}\), \(\Psi_{2}\), and \(\Psi_{3}\) still involve derivatives with respect to the invariants of \(\mathbf{C}\). However, since the invariants of \(\mathbf{b}\) are identical to those of \(\mathbf{C}\), the quantities\(\Psi_{1}\), \(\Psi_{2}\) and \(\Psi_{3}\) may also be considered to be the derivatives with respect to the invariants of \(\mathbf{b}\).

Isotropic Elasticity in Principal Directions

For isotropic materials, the principal directions of the strain and stress tensors are the same. Let the eigenvalues of \(\mathbf{C}\) be denoted by \(\lambda_{i}^{2}\) (\(i=1,2,3)\), then the strain energy density may be given as a function of these eigenvalues, \(\Psi\left(\lambda_{1}^{2},\lambda_{2}^{2},\lambda_{3}^{2}\right)\). To derive the expression for the stress, recognize that

\[ \begin{equation} \frac{\partial\lambda_{i}^{2}}{\partial\mathbf{C}}=\mathbf{N}_{i}\otimes\mathbf{N}_{i}\equiv\mathbf{A}_{i}\,,\label{eq68} \end{equation} \]

where the \(\mathbf{N}_{i}\) are the eigenvectors of \(\mathbf{C}\). It follows that the second Piola-Kirchhoff stress may be represented as

\[ \begin{equation} \mathbf{S}=\sum\limits_{i=1}^{3}S_{i}\mathbf{A}_{i}\,,\label{eq69} \end{equation} \]

where

\[ \begin{equation} S_{i}=2\frac{\partial\Psi}{\partial\lambda_{i}^{2}}\,.\label{eq70} \end{equation} \]

To evaluate the material elasticity tensor, recognize that

\[ \begin{equation} \frac{\partial\mathbf{A}_{i}}{\partial\mathbf{C}}=\frac{1}{\lambda_{i}^{2}-\lambda_{j}^{2}}\left(\mathbf{A}_{i}\odot\mathbf{A}_{j}+\mathbf{A}_{j}\odot\mathbf{A}_{i}\right)+\frac{1}{\lambda_{i}^{2}-\lambda_{k}^{2}}\left(\mathbf{A}_{i}\odot\mathbf{A}_{k}+\mathbf{A}_{k}\odot\mathbf{A}_{i}\right)\,,\label{eq71} \end{equation} \]

where \(i,j,k\) form a permutation over \(1,2,3\). Then it can be shown that the material elasticity tensor is given by

\[ \begin{equation} \begin{aligned}\boldsymbol{\mathbb{C}} & =\sum\limits_{i=1}^{3}4\frac{\partial^{2}\Psi}{\partial\lambda_{i}^{2}\partial\lambda_{i}^{2}}\mathbf{A}_{i}\otimes\mathbf{A}_{i}\\ & +\sum\limits_{i=1}^{3}\sum\limits_{j=i+1}^{3}4\frac{\partial^{2}\Psi}{\partial\lambda_{i}^{2}\partial\lambda_{j}^{2}}\left(\mathbf{A}_{i}\otimes\mathbf{A}_{j}+\mathbf{A}_{j}\otimes\mathbf{A}_{i}\right)\\ & +\sum\limits_{i=1}^{3}\sum\limits_{j=i+1}^{3}2\frac{S_{i}-S_{j}}{\lambda_{i}^{2}-\lambda_{j}^{2}}\left(\mathbf{A}_{i}\odot\mathbf{A}_{j}+\mathbf{A}_{j}\odot\mathbf{A}_{i}\right) \end{aligned} \,.\label{eq72} \end{equation} \]

When eigenvalues coincide, L'Hospital's rule may be used to evaluate the coefficient in the last term,

\[ \begin{equation} \lim\limits_{\lambda_{j}^{2}\to\lambda_{i}^{2}}2\frac{S_{i}-S_{j}}{\lambda_{i}^{2}-\lambda_{j}^{2}}=4\left(\frac{\partial^{2}\Psi}{\partial\lambda_{j}^{2}\partial\lambda_{j}^{2}}-\frac{\partial^{2}\Psi}{\partial\lambda_{i}^{2}\partial\lambda_{j}^{2}}\right)\,.\label{eq73} \end{equation} \]

The double summations in eq.\eqref{eq72} are arranged such that the summands represent fourth-order tensors with major and minor symmetries.

In the spatial frame, the Cauchy stress is given by

\[ \begin{equation} \boldsymbol{\sigma}=\sum\limits_{i=1}^{3}\sigma_{i}\mathbf{a}_{i}\,,\label{eq:IEPD-stress} \end{equation} \]

where

\[ \begin{equation} \mathbf{a}_{i}=\mathbf{n}_{i}\otimes\mathbf{n}_{i}\,,\label{eq:IEPD-tensor-basis} \end{equation} \]

and \(\mathbf{n}_{i}=\left(\mathbf{F}\cdot\mathbf{N}_{i}\right)/\lambda_{i}\) are the eigenvectors of \(\mathbf{b}\). The principal normal stresses are

\[ \begin{equation} \sigma_{i}=\frac{\lambda_{i}}{J}\frac{\partial\Psi}{\partial\lambda_{i}}\,.\label{eq:IEPD-principal-stress} \end{equation} \]

The spatial elasticity tensor is given by

\[ \begin{equation} \begin{aligned}\boldsymbol{\mathcal{C}} & =\sum\limits_{i=1}^{3}\left(J^{-1}\lambda_{i}^{2}\frac{\partial^{2}\Psi}{\partial\lambda_{i}^{2}}-\sigma_{i}\right)\mathbf{a}_{i}\otimes\mathbf{a}_{i}\\ & +\sum\limits_{i=1}^{3}\sum\limits_{j=i+1}^{3}J^{-1}\lambda_{i}\lambda_{j}\frac{\partial^{2}\Psi}{\partial\lambda_{i}\partial\lambda_{j}}\left(\mathbf{a}_{i}\otimes\mathbf{a}_{j}+\mathbf{a}_{j}\otimes\mathbf{a}_{i}\right)\\ & +\sum\limits_{i=1}^{3}\sum\limits_{j=i+1}^{3}2\frac{\lambda_{j}^{2}\sigma_{i}-\lambda_{i}^{2}\sigma_{j}}{\lambda_{i}^{2}-\lambda_{j}^{2}}\left(\mathbf{a}_{i}\odot\mathbf{a}_{j}+\mathbf{a}_{j}\odot\mathbf{a}_{i}\right) \end{aligned} \,.\label{eq:IEPD-elasticity} \end{equation} \]

Transversely Isotropic Hyperelasticity

Transverse isotropy can be introduced by adding a vector field representing the material preferred direction explicitly into the strain energy ². We require that the strain energy depends on a unit vector field \(\mathbf{A}\), which describes the local fiber direction in the undeformed configuration. When the material undergoes deformation, the vector \(\mathbf{A}\left(\mathbf{X}\right)\) may be described by a unit vector field \(\mathbf{a}\left(\boldsymbol{\varphi}\left(\mathbf{X}\right)\right)\). In general, the fibers will also undergo length change. The fiber stretch, \(\lambda\), can be determined in terms of the deformation gradient and the fiber direction in the undeformed configuration,

\[ \begin{equation} \lambda\mathbf{a}=\mathbf{F}\cdot\mathbf{A}\,.\label{eq95} \end{equation} \]

Also, since \(\mathbf{a}\) is a unit vector,

\[ \begin{equation} \lambda^{2}=\mathbf{A}\cdot\mathbf{C}\cdot\mathbf{A}\,.\label{eq96} \end{equation} \]

The strain energy function for a transversely isotropic material, \(\Psi\left(\mathbf{C},\mathbf{A},\mathbf{X}\right)\) is an isotropic function of \(\mathbf{C}\) and \(\mathbf{A}\otimes\mathbf{A}\). It can be shown ³ that the following set of invariants are sufficient to describe the material fully:

\[ \begin{equation} I_{1}=\tr\mathbf{C}\,,\quad I_{2}=\frac{1}{2}\left[\left(\tr\mathbf{C}\right)^{2}-\tr\mathbf{C}^{2}\right]\,,\quad I_{3}=\det C=J^{2}\,,\label{eq97} \end{equation} \]

\[ \begin{equation} I_{4}=\mathbf{A}\cdot\mathbf{C}\cdot\mathbf{A}\,,\quad I_{5}=\mathbf{A}\cdot\mathbf{C}^{2}\cdot\mathbf{A}\,.\label{eq98} \end{equation} \]

The strain energy function can be written in terms of these invariants such that

\[ \begin{equation} \Psi\left(\mathbf{C},\mathbf{A},\mathbf{X}\right)=\Psi\left(I_{1}\left(\mathbf{C}\right),I_{2}\left(\mathbf{C}\right),I_{3}\left(\mathbf{C}\right),I_{4}\left(\mathbf{C},\mathbf{A}\right),I_{5}\left(\mathbf{C},\mathbf{A}\right)\right)\,.\label{eq99} \end{equation} \]

The second Piola-Kirchhoff can now be obtained in the standard manner:

\[ \begin{equation} \mathbf{S}=2\frac{\partial\Psi}{\partial\mathbf{C}}=2\sum\limits_{i=1}^{5}\frac{\partial\Psi}{\partial I_{i}}\frac{\partial I_{i}}{\partial\mathbf{C}}\,.\label{eq100} \end{equation} \]

In the transversely isotropic constitutive models described in Chapter 5 it is further assumed that the strain energy function can be split into the following terms:

\[ \begin{equation} \Psi\left(\mathbf{C},\mathbf{A}\right)=\Psi_{1}\left(I_{1},I_{2},I_{3}\right)+\Psi_{2}\left(I_{4}\right)+\Psi_{3}\left(I_{1},I_{2},I_{3}I_{4}\right)\,.\label{eq101} \end{equation} \]

The strain energy function \(\Psi_{1}\) represents the material response of the isotropic ground substance matrix, \(\Psi_{2}\) represents the contribution from the fiber family (e.g. collagen), and \(\Psi_{3}\) is the contribution from interactions between the fibers and matrix. The form \eqref{eq101} generalizes many constitutive equations that have been successfully used in the past to describe biological soft tissues e.g. ⁴⁵⁶. While this relation represents a large simplification when compared to the general case, it also embodies almost all of the material behavior that one would expect from transversely isotropic, large deformation matrix-fiber composites.

Incompressibility

A material is incompressible when its mass density remains constant, \(\rho=\rho_{r}\), or equivalently according to eq.(2.5-34), its relative volume satisfies

\[ \begin{equation} J=1\label{eq:incompressibility-constraint} \end{equation} \]

everywhere throughout the material. No real material is truly incompressible, since true incompressibility implies that the speed of sound in that medium is infinite, which is not physically possible. Thus, one should view this assumption as an idealization of the behavior of real materials. We may replace eq.\eqref{eq:incompressibility-constraint} with the equivalent form \(\dot{J}=0\), and make use of eq.(2.5-38) in the form

\[ \begin{equation} \divg\mathbf{v}=\mathbf{I}:\mathbf{D}=0\,,\label{eq:incompressibility-redux} \end{equation} \]

where we recognized that the divergence of the velocity is equal to the trace of the velocity gradient \(\mathbf{L}\), hence, the trace of its symmetric part \(\mathbf{D}\).

Since \(J=\left(\det\mathbf{C}\right)^{1/2}\), the incompressibility constraint of eq.\eqref{eq:incompressibility-constraint} implies that the symmetric strain tensor \(\mathbf{E}=\frac{1}{2}\left(\mathbf{C}-\mathbf{I}\right)\) used as a state variable in the derivation of the hyperelasticity relation (see Section Constitutive Restrictions) does not truly have six independent components. Therefore, we cannot claim that \(\mathbf{E}\) can vary arbitrarily when enforcing the entropy inequality. Instead, we introduce the incompressibility constraint of eq.\eqref{eq:incompressibility-redux} into the Clausius-Duhem inequality, eq.(2.4-9), using the method of Lagrange multipliers,

\[ -\rho\left(\dot{\psi}+\eta\dot{\theta}\right)+\left(\boldsymbol{\sigma}+\lambda\mathbf{I}\right):\mathbf{D}-\frac{1}{\theta}\mathbf{q}\cdot\grad\theta\geqslant0\,, \]

where \(\lambda\) is the Lagrange multiplier. Now, we may proceed with the constitutive assumption that \(\psi=\psi\left(\mathbf{E}\right)\), as done in Section Constitutive Restrictions, until we conclude that the stress \(\boldsymbol{\sigma}\) is constrained to have the form

\[ \begin{equation} \boldsymbol{\sigma}=p\mathbf{I}+2J^{-1}\mathbf{F}\cdot\frac{\partial\Psi_{r}}{\partial\mathbf{C}}\cdot\mathbf{F}^{T}\,,\quad p\equiv-\lambda\,.\label{eq:Cauchy-hyperelastic-incomp} \end{equation} \]

This expression shows that the incompressibility constraint requires us to have a parameter \(p\) (a pressure) that accounts for that part of the strain energy density that would normally be stored in the material as a result of volumetric strain. Analytically, we would solve for this additional scalar parameter \(p\) using the constraint of eq.\eqref{eq:incompressibility-constraint}. Finally, substituting eq.\eqref{eq:Cauchy-hyperelastic-incomp} into eq.\eqref{eq51}, the PK2 stress in an incompressible material takes the form

\[ \begin{equation} \mathbf{S}=Jp\mathbf{C}^{-1}+2\frac{\partial\Psi_{r}}{\partial\mathbf{C}}\,.\label{eq:PK2-hyperelastic-incomp} \end{equation} \]

Nearly-Incompressible Hyperelasticity

When dealing with incompressible and nearly incompressible materials, it proves useful to separate the volumetric and the deviatoric (distortional) components of the deformation gradient \(\mathbf{F}\). Such a separation must ensure that the deviatoric part of the deformation gradient, namely \(\mathbf{\tilde{F}}\), does not produce any change in volume. Noting that the determinant of the deformation gradient gives the volume ratio, the determinant of \(\mathbf{\tilde{F}}\) must therefore satisfy,

\[ \begin{equation} \det\mathbf{\tilde{F}}=1\,.\label{eq38} \end{equation} \]

This condition can be achieved by choosing \(\mathbf{\tilde{F}}\) as,

\[ \begin{equation} \mathbf{\tilde{F}}=J^{-1/3}\mathbf{F}\,.\label{eq39} \end{equation} \]

The right Cauchy-Green tensor can also be split into volumetric and deviatoric components. With the use of,\eqref{eq39} the deviatoric right Cauchy-Green tensor is given by

\[ \begin{equation} \mathbf{\tilde{C}}=\tilde{\mathbf{F}}^{T}\cdot\tilde{\mathbf{F}}=J^{-2/3}\mathbf{C}\,.\label{eq43-1} \end{equation} \]

The process of defining constitutive equations in the case of nearly incompressible hyperelasticity is simplified by uncoupling the strain energy density \(\Psi_{r}\left(\mathbf{C}\right)\) into a volumetric energy component \(U\left(J\right)\) and a distortional component \(\tilde{\Psi}\left(\mathbf{C}\right)\):

\[ \begin{equation} \Psi\left(\mathbf{C}\right)=\tilde{\Psi}\left(\mathbf{C}\right)+U\left(J\right)\,.\label{eq:UC-SED} \end{equation} \]

The second Piola-Kirchhoff tensor for a material defined by \eqref{eq:UC-SED} is obtained in the standard manner with the help of equation \eqref{eq:PK2-isotropic}.

\[ \begin{equation} \begin{aligned}\mathbf{S} & =2\frac{\partial\Psi}{\partial\mathbf{C}}\\ & =2\frac{\partial\tilde{\Psi}}{\partial\mathbf{C}}+2\frac{dU}{dJ}\frac{\partial J}{\partial C}\\ & =2\frac{\partial\tilde{\Psi}}{\partial\mathbf{C}}+pJ\mathbf{C}^{-1} \end{aligned} \,,\label{eq:UC-PK2-redux} \end{equation} \]

where the pressure \(p\) is defined as

\[ \begin{equation} p=\frac{dU}{dJ}\,.\label{eq:UC-p} \end{equation} \]

An example for \(U\) that will be used later in the definition of the constitutive models is

\[ \begin{equation} U\left(J\right)=\frac{1}{2}\kappa\left(\ln J\right)^{2}\,.\label{eq:UC-U-exa} \end{equation} \]

The parameter \(\kappa\) will be used later as a penalty factor that will enforce the (nearly-) incompressible constraint. However, \(\kappa\) can represent a true material coefficient, namely the bulk modulus, for a compressible material that happens to have a hyperelastic strain energy function in the form of \eqref{eq:UC-SED}. In the case where the dilatational energy is given by \eqref{eq:UC-U-exa}, the pressure is

\[ \begin{equation} p=\kappa\frac{\ln J}{J}\,.\label{eq:UC-p-exa} \end{equation} \]

Equation \eqref{eq:UC-PK2-redux} can be further developed by applying the chain rule to the first term:

\[ \begin{equation} \mathbf{S}=pJ\mathbf{C}^{-1}+J^{-2/3}\dev\tilde{\mathbf{S}}\,,\label{eq:UC-PK2-redux2} \end{equation} \]

where the fictitious second Piola-Kirchoff tensor ⁷ is defined by,

\[ \begin{equation} \tilde{\mathbf{S}}=2\frac{\partial\tilde{\Psi}}{\partial\tilde{\mathbf{C}}}\,,\label{eq:UC-PK2-tilde} \end{equation} \]

and \(\Dev\) is the deviator operator in the reference frame:

\[ \begin{equation} \Dev\left(\cdot\right)=\left(\cdot\right)-\frac{1}{3}\left(\left(\cdot\right):\mathbf{C}\right)\mathbf{C}^{-1}.\label{eq:UC-Dev} \end{equation} \]

The Cauchy stress can then be obtained from eq.\eqref{eq52}\(_{\mathrm{3}}\):

\[ \begin{equation} \boldsymbol{\sigma}=p\mathbf{I}+\dev\tilde{\boldsymbol{\sigma}}\,,\label{eq:UC-Cauchy-stress} \end{equation} \]

where

\[ \begin{equation} \tilde{\boldsymbol{\sigma}}=\frac{2}{J}\mathbf{\tilde{F}}\cdot\frac{\partial\tilde{\Psi}}{\partial\mathbf{\tilde{C}}}\cdot\mathbf{\tilde{F}}^{T}\,,\label{eq:UC-Cauchy-stress-tilde} \end{equation} \]

and

\[ \begin{equation} \dev\left(\cdot\right)=\left(\cdot\right)-\frac{1}{3}\left(\left(\cdot\right):\mathbf{I}\right)\mathbf{I}\label{eq:UC-dev} \end{equation} \]

is the deviator operator in the spatial frame.

The following expression will be useful in the following development.

\[ \begin{equation} \begin{aligned}\frac{d\tilde{C}_{IJ}}{dC_{KL}} & =J^{-2/3}\left(\frac{1}{2}\left(\delta_{IK}\delta_{JL}+\delta_{IL}\delta_{JK}\right)-\frac{1}{3}\tilde{C}_{IJ}\tilde{C}_{KL}^{-1}\right)\\ \frac{\partial\tilde{\mathbf{C}}}{\partial\mathbf{C}} & =J^{-2/3}\left(\mathbf{I}\odot\mathbf{I}-\frac{1}{3}\tilde{\mathbf{C}}\otimes\tilde{\mathbf{C}}^{-1}\right) \end{aligned} \,.\label{eq89} \end{equation} \]

Notice that the contraction with a symmetric tensor \(\mathbf{A}\) results in,

\[ \begin{equation} \begin{aligned}\frac{d\tilde{C}_{IJ}}{dC_{KL}}A_{IJ} & =J^{-2/3}\Dev A_{KL}\\ \mathbf{A}:\frac{\partial\tilde{\mathbf{C}}}{\partial\mathbf{C}} & =J^{-2/3}\Dev\mathbf{A} \end{aligned} \,.\label{eq90} \end{equation} \]

The elasticity tensor, defined in \eqref{eq55}, takes on the following form.

\[ \begin{equation} \begin{aligned}\mathbb{C}_{IJKL} & =\left(J^{2}\frac{dp}{dJ}+Jp\right)C_{IJ}^{-1}C_{KL}^{-1}-2pJ\mathcal{I}_{IJKL}\\ & -\frac{2}{3}J^{-2/3}\left(\Dev\tilde{S}_{IJ}C_{KL}^{-1}+C_{IJ}^{-1}\Dev\tilde{S}_{KL}\right)\\ & +\frac{2}{3}\tilde{S}_{RS}\tilde{C}_{RS}\left(\mathcal{I}_{IJKL}-\frac{1}{3}C_{KL}^{-1}C_{IJ}^{-1}\right)+J^{-4/3}\hat{\mathbb{C}}_{IJKL}\\ \boldsymbol{\mathbb{C}} & =\left(J^{2}\frac{dp}{dJ}+Jp\right)\mathbf{C}^{-1}\otimes\mathbf{C}^{-1}-2pJ\boldsymbol{\mathcal{I}}\\ & -\frac{2}{3}J^{-2/3}\left(\Dev\tilde{\mathbf{S}}\otimes\mathbf{C}^{-1}+\mathbf{C}^{-1}\otimes\Dev\tilde{\mathbf{S}}\right)\\ & +\frac{2}{3}\left(\tilde{\mathbf{S}}:\tilde{\mathbf{C}}\right)\left(\boldsymbol{\mathcal{I}}-\frac{1}{3}\mathbf{C}^{-1}\otimes\mathbf{C}^{-1}\right)+J^{-4/3}\hat{\boldsymbol{\mathbb{C}}} \end{aligned} \,,\label{eq91} \end{equation} \]

where

\[ \begin{equation} \begin{aligned}\hat{\mathcal{\mathbb{C}}}_{IJKL} & =\tilde{\mathcal{\mathbb{C}}}_{IJKL}-\frac{1}{3}\left(\tilde{\mathcal{\mathbb{C}}}_{IJRS}\tilde{C}_{RS}\tilde{C}_{KL}^{-1}+\tilde{\mathcal{\mathbb{C}}}_{RSKL}\tilde{C}_{RS}\tilde{C}_{IJ}^{-1}\right)+\frac{1}{9}\tilde{C}_{IJ}^{-1}\tilde{C}_{RS}\tilde{\mathbb{C}}_{RSMN}\tilde{C}_{MN}\tilde{C}_{KL}^{-1}\\ \hat{\mathcal{\boldsymbol{\mathbb{C}}}} & =\tilde{\mathcal{\boldsymbol{\mathbb{C}}}}-\frac{1}{3}\left(\tilde{\mathcal{\boldsymbol{\mathbb{C}}}}:\tilde{\mathbf{C}}\otimes\tilde{\mathbf{C}}^{-1}+\tilde{\mathbf{C}}^{-1}\otimes\tilde{\mathbf{C}}:\tilde{\boldsymbol{\mathcal{\mathbb{C}}}}\right)+\frac{1}{9}\tilde{\mathbf{C}}^{-1}\otimes\tilde{\mathbf{C}}:\tilde{\boldsymbol{\mathbb{C}}}:\tilde{\mathbf{C}}\otimes\tilde{\mathbf{C}}^{-1} \end{aligned} \,.\label{eq92} \end{equation} \]

The spatial elasticity tensor follows from

\[ \begin{equation} \begin{aligned}\boldsymbol{\mathcal{C}} & =\left(J\frac{dp}{dJ}+p\right)\mathbf{I}\otimes\mathbf{I}-2p\mathbf{I}\odot\mathbf{I}\\ & -\frac{2}{3}\left(\dev\tilde{\boldsymbol{\sigma}}\otimes\mathbf{I}+\mathbf{I}\otimes\dev\tilde{\boldsymbol{\sigma}}\right)\\ & +\frac{2}{3}\left(\tilde{\boldsymbol{\sigma}}:\mathbf{I}\right)\left(\mathbf{I}\odot\mathbf{I}-\frac{1}{3}\mathbf{I}\otimes\mathbf{I}\right)+\hat{\mathcal{\boldsymbol{C}}} \end{aligned} \,,\label{eq:UC-spatial-elasticity} \end{equation} \]

where

\[ \begin{equation} \begin{aligned}\hat{\mathcal{C}}_{ijkl} & =\frac{1}{J}\tilde{F}_{iI}\tilde{F}_{jJ}\tilde{F}_{kK}\tilde{F}_{lL}\hat{\mathcal{\mathbb{C}}}_{IJKL}\\ \hat{\boldsymbol{\mathcal{C}}} & =\frac{1}{J}\left(\tilde{\mathbf{F}}\oslash\tilde{\mathbf{F}}\right):\hat{\boldsymbol{\mathbb{C}}}:\left(\tilde{\mathbf{F}}^{T}\oslash\tilde{\mathbf{F}}^{T}\right)\\ & =\tilde{\boldsymbol{\mathcal{C}}}-\frac{1}{3}\left(\mathbf{I}\otimes\mathbf{I}:\tilde{\boldsymbol{\mathcal{C}}}+\tilde{\boldsymbol{\mathcal{C}}}:\mathbf{I}\otimes\mathbf{I}\right)+\frac{1}{9}\left(\mathbf{I}:\tilde{\boldsymbol{\mathcal{C}}}:\mathbf{I}\right)\mathbf{I}\otimes\mathbf{I} \end{aligned} \,.\label{eq:UC-spatial-elasticity-hat} \end{equation} \]

Tension-Bearing Fiber Materials

In biomechanics we often find it convenient to model fibrous or fibrillar materials using one-dimensional fibers that can only sustain tension. Typically, such fibers are represented with a strain energy density function

\[ \begin{equation} \Psi_{r}\left(\mathbf{C}\right)=H\left(I_{n}-1\right)\Psi_{n}\left(I_{n}\right)\,,\quad I_{n}=\mathbf{n}_{r}\cdot\mathbf{C}\cdot\mathbf{n}_{r}\label{eq:fiber-model} \end{equation} \]

where \(\mathbf{n}_{r}\) is the unit vector along the fiber in its reference configuration, and \(I_{n}\) is the square of the stretch ratio along the fiber. The Heaviside unit step function \(H\left(I_{n}-1\right)\) in eq.\eqref{eq:fiber-model} ensures that the fiber contributes strain energy only when it is under tension (\(I_{n}>1\)); thus, the constitutive model \(\Psi_{n}\left(I_{n}\right)\) for the tensile response of the fiber must reduce to zero when \(I_{n}=1\).

Using the hyperelasticity relations presented above, the Cauchy stress in this fiber material can be evaluated as

\[ \begin{equation} \boldsymbol{\sigma}=H\left(I_{n}-1\right)2J^{-1}\frac{d\Psi_{n}}{dI_{n}}\left(\mathbf{F}\cdot\mathbf{n}_{r}\right)\otimes\left(\mathbf{F}\cdot\mathbf{n}_{r}\right)\,,\label{eq:fiber-stress} \end{equation} \]

and the spatial elasticity tensor is

\[ \begin{equation} \mathcal{C}=4J^{-1}\frac{d^{2}\Psi_{n}}{dI_{n}^{2}}\left(\mathbf{F}\cdot\mathbf{n}_{r}\right)\otimes\left(\mathbf{F}\cdot\mathbf{n}_{r}\right)\otimes\left(\mathbf{F}\cdot\mathbf{n}_{r}\right)\otimes\left(\mathbf{F}\cdot\mathbf{n}_{r}\right)\,.\label{eq:fiber-elasticity} \end{equation} \]

If we denote the unit vector along the fiber in the current configuration as \(\mathbf{n}\equiv I_{n}^{-1/2}\mathbf{F}\cdot\mathbf{n}_{r}\), the above expressions may be rewritten as \(\boldsymbol{\sigma}=H\left(I_{n}-1\right)2J^{-1}I_{n}\Psi_{n}^{\prime}\left(I_{n}\right)\mathbf{n}\otimes\mathbf{n}\) and \(\mathcal{C}=4J^{-1}I_{n}^{2}\Psi_{n}^{\prime\prime}\left(I_{n}\right)\mathbf{n}\otimes\mathbf{n}\otimes\mathbf{n}\otimes\mathbf{n}\). As explained in ⁸, these fiber models must be combined with a ground matrix in order to produce a stable material response. In FEBio this can be done by using a constrained mixture of solid constituents (for example, see Section Simple Solid Mixtures).

In the classical fiber mechanics literature it was suggested that uncoupled fiber formulations could also be implemented, whereby

\[ \begin{equation} \Psi_{r}\left(\mathbf{C}\right)=H\left(\tilde{I}_{n}-1\right)\tilde{\Psi}_{n}\left(\tilde{I}_{n}\right)+U\left(J\right)\,,\quad\tilde{I}_{n}=\mathbf{n}_{r}\cdot\tilde{\mathbf{C}}\cdot\mathbf{n}_{r}\,,\label{eq:fiber-model-UC} \end{equation} \]

and the Cauchy stress \(\boldsymbol{\sigma}\) is given by eq.\eqref{eq:UC-Cauchy-stress} where

\[ \begin{equation} \tilde{\boldsymbol{\sigma}}=H\left(\tilde{I}_{n}-1\right)2J^{-1}\frac{d\tilde{\Psi}_{n}}{d\tilde{I}_{n}}\left(\tilde{\mathbf{F}}\cdot\mathbf{n}_{r}\right)\otimes\left(\tilde{\mathbf{F}}\cdot\mathbf{n}_{r}\right)\,.\label{eq:fiber-stress-UC} \end{equation} \]

Uncoupled fiber formulations of this kind are available in FEBio. More recently however, several studies have demonstrated that using \(H\left(\tilde{I}_{n}-1\right)\) to detect whether a fiber is in tension or compression is non-physical, since \(\tilde{I}_{n}\ne I_{n}\) and \(I_{n}\) is the sole true measure of the tensile stretch in the fiber ⁹¹⁰¹¹. Therefore, uncoupled fiber formulations have fallen out of favor, even though FEBio still allows users to employ these for historical reasons. It is now recommended to use the standard (unconstrained) fiber models, also available in FEBio, with the formulation given in Eqs.\eqref{eq:fiber-model}-\eqref{eq:fiber-elasticity}, even when the ground matrix uses an uncoupled formulation.

---

1. Coleman, Bernard D; Noll, Walter. "The thermodynamics of elastic materials with heat conduction and viscosity." Arch Ration Mech An, vol. 13, pp. 167--178 (1963). ↩

2. Weiss, J.A.; Maker, B.N.; Govindjee, S.. "Finite element implementation of incompressible, transversely isotropic hyperelasticity." Computer Methods in Applications of Mechanics and Engineering, vol. 135, pp. 107-128 (1996). ↩

3. Spencer, Anthony James Merril. "Continuum Theory of the Mechanics of Fibre-Reinforced Composites." Springer-Verlag (1984). ↩

4. Horowitz, A.; Sheinman, I.; Lanir, Y.; Perl, M.; Sideman, S.. "Nonlinear Incompressilbe Finite Element for Simulating Loading of Cardiac Tissue- part I: two Dimensional Formulation for Thin Myocardial Strips." Journal of Biomechanical Engineering, Transactions of the ASME, vol. 110, pp. 57-61 (1988). ↩

5. Humphrey, J. D.; Strumpf, R. K.; Yin, F. C. P.. "Determination of a constitutive relation for passive myocardium. I. A new functional form." Journal of Biomechanical Engineering, Transactions of the ASME, vol. 112, pp. 333-339 (1990). ↩

6. Humphrey, J. D.; Yin, F. C. P.. "On constitutive Relations and Finite Deformations of Passive Cardiac Tissue: I. A Pseudostrain-Energy Function." Journal of Biomechanical Engineering, Transactions of the ASME, vol. 109, pp. 298-304 (1987). ↩

7. Holzapfel, Gerhard A. "Nonlinear solid mechanics: a continuum approach for engineering." Wiley (2000). ↩

8. Ateshian, G. A.. "Anisotropy of fibrous tissues in relation to the distribution of tensed and buckled fibers." J Biomech Eng, vol. 129, pp. 240-9 (2007). ↩

9. Sansour, Carlo. "On the physical assumptions underlying the volumetric-isochoric split and the case of anisotropy." European Journal of Mechanics-A/Solids, vol. 27, pp. 28--39 (2008). ↩

10. Helfenstein, J; Jabareen, M; Mazza, Edoardo; Govindjee, S. "On non-physical response in models for fiber-reinforced hyperelastic materials." International Journal of Solids and Structures, vol. 47, pp. 2056--2061 (2010). ↩

11. G{\"u}ltekin, Osman; Dal, H{\"u}sn{\"u}; Holzapfel, Gerhard A. "On the quasi-incompressible finite element analysis of anisotropic hyperelastic materials." Computational mechanics, vol. 63, pp. 443--453 (2019). ↩