5.2 Compressible Materials

Isotropic Elasticity

The linear elastic material model as described in Section Linear Elasticity is only valid for small strains and small rotations. A first modification to this model to the range of nonlinear deformations is given by the St. Venant-Kirchhoff model ¹, which in FEBio is referred to as isotropic elasticity. This model is objective for large strains and rotations. For the isotropic case it can be derived from the following hyperelastic strain-energy function:

\[ \begin{equation} W=\frac{1}{2}\lambda\left(\tr\mathbf{E}\right)^{2}+\mu\mathbf{E}:\mathbf{E}\,.\label{eq411} \end{equation} \]

The second Piola-Kirchhoff stress can be derived from this:

\[ \begin{equation} \mathbf{S}=\frac{\partial W}{\partial\mathbf{E}}=\lambda\left(\tr\mathbf{E}\right)\mathbf{I}+2\mu\mathbf{E}\,.\label{eq412} \end{equation} \]

Note that these equations are similar to the corresponding equations in the linear elastic case, except that the small strain tensor is replaced by the Green-Lagrange elasticity tensor \(\mathbf{E}\). The material elasticity tensor is then given by,

\[ \begin{equation} \boldsymbol{\mathbb{C}}=\frac{\partial\mathbf{S}}{\partial\mathbf{E}}=\lambda\mathbf{I}\otimes\mathbf{I}+2\mu\mathbf{I}\odot\mathbf{I}\,.\label{eq413} \end{equation} \]

It is important to note that although this model is objective, it should only be used for small strains. For large strains, the response can be somewhat strange if not completely unrealistic. For example, it can be shown that under uni-axial tension the stress becomes infinite and the volume tends to zero for finite strains. Therefore, for large strains it is highly recommended to avoid this material and instead use one of the other non-linear material models described below. The Cauchy stress is

\[ \begin{equation} \boldsymbol{\sigma}=\frac{1}{J}\left(\lambda\tr\mathbf{E}-\mu\right)\mathbf{b}+\frac{\mu}{J}\mathbf{b}^{2}\,,\label{eq414} \end{equation} \]

where \(\tr\mathbf{E}=\left(\tr\mathbf{b}-3\right)/2\), whereas the spatial elasticity tensor is

\[ \begin{equation} \boldsymbol{\mathcal{C}}=\frac{\lambda}{J}\mathbf{b}\otimes\mathbf{b}+\frac{2}{J}\mu\mathbf{b}\odot\mathbf{b}\,.\label{eq415} \end{equation} \]

Orthotropic Elasticity

An extension of the St. Venant-Kirchhoff model ¹ to orthotropic symmetry is provided in FEBio, referred to as orthotropic elasticity. This model is objective for large strains and rotations. It can be derived from the following hyperelastic strain-energy function:

\[ \begin{equation} W=\sum\limits_{a=1}^{3}\mu_{a}\mathbf{A}_{a}^{0}:\mathbf{E}^{2}+\frac{1}{2}\sum\limits_{b=1}^{3}\lambda_{ab}\left(\mathbf{A}_{a}^{0}:\mathbf{E}\right)\left(\mathbf{A}_{b}^{0}:\mathbf{E}\right)\,,\label{eq416} \end{equation} \]

where \(\mathbf{A}_{a}^{0}=\mathbf{a}_{a}^{0}\otimes\mathbf{a}_{a}^{0}\) is the structural tensor corresponding to one of the three mutually orthogonal planes of symmetry whose unit outward normal is \(\mathbf{a}_{a}^{0}\) (\(\mathbf{a}_{a}^{0}\cdot\mathbf{a}_{b}^{0}=\delta_{ab}\)). The material constants are the three shear moduli \(\mu_{a}\) and six moduli \(\lambda_{ab}\), where \(\lambda_{ba}=\lambda_{ab}\). They may be related to the Young's moduli \(E_{a}\), shear moduli \(G_{ab}\) and Poisson's ratios \(\nu_{ab}\) via

\[ \begin{equation} \begin{array}{l} \left[\begin{array}{cccccc} \lambda_{11}+2\mu_{1} & \lambda_{12} & \lambda_{13} & 0 & 0 & 0\\ \lambda_{12} & \lambda_{22}+2\mu_{2} & \lambda_{23} & 0 & 0 & 0\\ \lambda_{13} & \lambda_{23} & \lambda_{33}+2\mu_{3} & 0 & 0 & 0\\ 0 & 0 & 0 & \left(\mu_{1}+\mu_{2}\right)/2 & 0 & 0\\ 0 & 0 & 0 & 0 & \left(\mu_{2}+\mu_{3}\right)/2 & 0\\ 0 & 0 & 0 & 0 & 0 & \left(\mu_{3}+\mu_{1}\right)/2 \end{array}\right]^{-1}=\\ \left[\begin{array}{cccccc} 1/E_{1} & -\nu_{12}/E_{1} & -\nu_{13}/E_{1} & 0 & 0 & 0\\ -\nu_{21}/E_{2} & 1/E_{2} & -\nu_{23}/E_{2} & 0 & 0 & 0\\ -\nu_{31}/E_{3} & -\nu_{32}/E_{3} & 1/E_{3} & 0 & 0 & 0\\ 0 & 0 & 0 & 1/G_{12} & 0 & 0\\ 0 & 0 & 0 & 0 & 1/G_{23} & 0\\ 0 & 0 & 0 & 0 & 0 & 1/G_{31} \end{array}\right]\,. \end{array}\label{eq417} \end{equation} \]

The second Piola-Kirchhoff stress can be derived from this strain energy density function:

\[ \begin{equation} \begin{aligned}\mathbf{S} & =\frac{\partial W}{\partial\mathbf{E}}=\sum\limits_{a=1}^{3}\mu_{a}\left(\mathbf{A}_{a}^{0}\cdot\mathbf{E}+\mathbf{E}\cdot\mathbf{A}_{a}^{0}\right)\\ & +\frac{1}{2}\sum\limits_{b=1}^{3}\lambda_{ab}\left[\left(\mathbf{A}_{a}^{0}:\mathbf{E}\right)\mathbf{A}_{b}^{0}+\left(\mathbf{A}_{b}^{0}:\mathbf{E}\right)\mathbf{A}_{a}^{0}\right]\,. \end{aligned} \label{eq418} \end{equation} \]

Note that these equations are similar to the corresponding equations in the linear orthotropic elastic case, except that the small strain tensor is replaced by the Green-Lagrange elasticity tensor \(\mathbf{E}\). The material elasticity tensor is then given by,

\[ \begin{equation} \boldsymbol{\mathbb{C}}=\frac{\partial\mathbf{S}}{\partial\mathbf{E}}=\sum\limits_{a=1}^{3}\mu_{a}\left(\mathbf{A}_{a}^{0}\odot\mathbf{I}+\mathbf{I}\odot\mathbf{A}_{a}^{0}\right)+\frac{1}{2}\sum\limits_{b=1}^{3}\lambda_{ab}\left(\mathbf{A}_{a}^{0}\otimes\mathbf{A}_{b}^{0}+\mathbf{A}_{b}^{0}\otimes\mathbf{A}_{a}^{0}\right)\,.\label{eq419} \end{equation} \]

It is important to note that although this model is objective, it should only be used for small strains. For large strains, the response can be somewhat strange if not completely unrealistic. For example, it can be shown that under uni-axial tension the stress becomes infinite and the volume tends to zero for finite strains. Therefore, for large strains it is highly recommended to avoid this material and instead use one of the other non-linear material models described below. The Cauchy stress is

\[ \begin{equation} \begin{aligned}\boldsymbol{\sigma} & =\sum\limits_{a=1}^{3}\frac{\mu_{a}}{2J}\left(\mathbf{A}_{a}\cdot\left(\mathbf{b}-\mathbf{I}\right)+\left(\mathbf{b}-\mathbf{I}\right)\cdot\mathbf{A}_{a}\right)\\ & +\frac{1}{2}\sum\limits_{b=1}^{3}\frac{\lambda_{ab}}{2J}\left[\left(\mathbf{A}_{a}:\mathbf{I}-1\right)\mathbf{A}_{b}+\left(\mathbf{A}_{b}:\mathbf{I}-1\right)\mathbf{A}_{a}\right]\,, \end{aligned} \label{eq420} \end{equation} \]

where \(\mathbf{A}_{a}=\mathbf{F}\cdot\mathbf{A}_{a}^{0}\cdot\mathbf{F}^{T}\) and the spatial elasticity tensor is

\[ \begin{equation} \boldsymbol{\mathcal{C}}=\sum\limits_{a=1}^{3}\frac{\mu_{a}}{J}\left(\mathbf{A}_{a}\odot\mathbf{b}+\mathbf{b}\odot\mathbf{A}_{a}\right)+\frac{1}{2}\sum\limits_{b=1}^{3}\frac{\lambda_{ab}}{J}\left({\mathbf{A}_{a}\otimes\mathbf{A}_{b}+\mathbf{A}_{b}\otimes\mathbf{A}_{a}}\right)\,.\label{eq421} \end{equation} \]

Neo-Hookean Hyperelasticity

This is a compressible neo-Hookean material. It is derived from the following hyperelastic strain energy function ¹:

\[ \begin{equation} W=\frac{\mu}{2}\left(I_{1}-3\right)-\mu\ln J+\frac{\lambda}{2}\left(\ln J\right)^{2}\,.\label{eq422} \end{equation} \]

The parameters \(\mu\) and \(\lambda\) are the Lamé parameters from linear elasticity. This model reduces to the isotropic linear elastic model for small strains and rotations.

The Cauchy stress is given by,

\[ \begin{equation} \boldsymbol{\sigma}=\frac{\mu}{J}\left(\mathbf{b}-\mathbf{I}\right)+\frac{\lambda}{J}\left(\ln J\right)\mathbf{I}\,,\label{eq423} \end{equation} \]

and the spatial elasticity tensor is given by

\[ \begin{equation} \boldsymbol{\mathcal{C}}=\frac{\lambda}{J}\mathbf{I}\otimes\mathbf{I}+\frac{2}{J}\left(\mu-\lambda\ln J\right)\mathbf{I}\odot\mathbf{I}\,.\label{eq424} \end{equation} \]

The neo-Hookean material is an extension of Hooke's law for the case of large deformations. It is useable for certain plastics and rubber-like substances. A generalization of this model is the Mooney-Rivlin material, which is often used to describe the elastic response of biological tissue.

In FEBio this constitutive model uses a standard displacement-based element formulation and a "coupled" strain energy, so care must be taken when modeling materials with nearly-incompressible material behavior to avoid element locking.

Natural Neo-Hookean

This is a compressible isotropic neo-Hookean material that uses the natural (Hencky) strain tensor invariants to formulate its strain energy density. These invariants are reviewed in ². The left Hencky strain is evaluated from \(\boldsymbol{\eta}=\ln\mathbf{V}\) where \(\mathbf{V}\) is the left stretch tensor in the polar decomposition of the deformation gradient \(\mathbf{F}=\mathbf{V}\cdot\mathbf{R}\). To evaluate \(\boldsymbol{\eta}\) we first evaluate the left Cauchy-Green tensor \(\mathbf{b}=\mathbf{V}^{2}\) from \(\mathbf{F}\) as in eq.(2.6-28) and get its eigenvalues \(\lambda_{i}^{2}\) and eigenvectors \(\mathbf{n}_{i}\). Then

\[ \begin{equation} \boldsymbol{\eta}=\sum_{i=1}^{3}\left(\ln\lambda_{i}\right)\mathbf{n}_{i}\otimes\mathbf{n}_{i}\,.\label{eq:left-Hencky} \end{equation} \]

The invariants \(K_{i}\) of the natural strain tensor are

\[ \begin{equation} \begin{aligned}K_{1} & =\tr\boldsymbol{\eta}=\ln J & \text{amount of dilatation}\\ K_{2} & =\left|\dev\boldsymbol{\eta}\right|=\sqrt{\dev\boldsymbol{\eta}:\dev\boldsymbol{\eta}} & \text{amount of distortion}\\ K_{3} & =3\sqrt{6}\det\boldsymbol{\Phi} & \text{mode of distortion} \end{aligned} \label{eq:Hencky-isotropic-invariants} \end{equation} \]

where \(J=\det\mathbf{F}\) as usual, and

\[ \begin{equation} \boldsymbol{\Phi}=\frac{1}{K_{2}}\dev\boldsymbol{\eta}\,.\label{eq:Hencky-3rd-basis} \end{equation} \]

It can be shown that

\[ \begin{equation} \boldsymbol{\eta}=\frac{1}{3}K_{1}\mathbf{I}+K_{2}\boldsymbol{\Phi}\,.\label{eq:left-Hencky-redux} \end{equation} \]

Note that \(K_{2}\boldsymbol{\Phi}\to\mathbf{0}\) as \(K_{2}\to0\). It also follows that \(\boldsymbol{\eta}:\boldsymbol{\eta}=\frac{1}{3}K_{1}^{2}+K_{2}^{2}\). As explained in ², \(K_{1}\in\left(-\infty,\infty\right)\) with positive \(K_{1}\) implying expansion and negative \(K_{1}\) implying contraction. Similarly, \(K_{2}\in[0,\infty)\), with \(K_{2}=0\) implying no distortion. Finally, \(K_{3}\in\left[-1,1\right]\) with \(K_{3}=1\) representing uniaxial extension, \(K_{3}=-1\) representing uniaxial contraction and \(K_{3}=0\) representing pure shear.

For the natural neo-Hookean material the strain energy density is

\[ \begin{equation} W=\frac{\kappa}{2}K_{1}^{2}+\mu K_{2}^{2}\label{eq:NNH-sed} \end{equation} \]

where \(\kappa\) is the material's bulk modulus and \(\mu\) is its shear modulus. To evaluate the Cauchy stress \(\boldsymbol{\sigma}\) and spatial elasticity tensor \(\boldsymbol{\mathcal{C}}\), we use the framework of isotropic elasticity in principal directions (Section Isotropic Elasticity in Principal Directions). This requires us to express \(K_{1}\) and \(K_{2}\) in terms of the eigenvalues \(\lambda_{i}\),

\[ \begin{equation} \begin{aligned}K_{1} & =\ln\left(\lambda_{1}\lambda_{2}\lambda_{3}\right)\\ K_{2} & =\frac{1}{3}\sqrt{\left(\ln\frac{\lambda_{1}^{2}}{\lambda_{2}\lambda_{3}}\right)^{2}+\left(\ln\frac{\lambda_{2}^{2}}{\lambda_{3}\lambda_{1}}\right)^{2}+\left(\ln\frac{\lambda_{3}^{2}}{\lambda_{1}\lambda_{2}}\right)^{2}} \end{aligned} \,.\label{eq:NNH-invariants} \end{equation} \]

Now the stress \(\boldsymbol{\sigma}\) is given by eq.(2.6-35) where, based on eq.(2.6-37), the principal normal stresses are evaluated as

\[ \begin{equation} \sigma_{i}=\frac{1}{3J}\left[\left(3\kappa+4\mu\right)\ln\lambda_{i}+\left(3\kappa-2\mu\right)\left(\ln\lambda_{j}+\ln\lambda_{k}\right)\right]\,,\label{eq:NNH-principal-stress} \end{equation} \]

with \(i,j,k\) forming a permutation over \(1,2,3\). Similarly, the spatial elasticity tensor \(\boldsymbol{\mathcal{C}}\) is given by eq.(2.6-38) where we substitute

\[ \begin{equation} J^{-1}\lambda_{i}^{2}\frac{\partial^{2}W}{\partial\lambda_{i}^{2}}-\sigma_{i}=\frac{1}{3J}\left[\left(3\kappa+4\mu\right)\left(1-2\ln\lambda_{i}\right)-2\left(3\kappa-2\mu\right)\left(\ln\lambda_{j}+\ln\lambda_{k}\right)\right]\label{eq:NNH-elasticity-coef1} \end{equation} \]

and

\[ \begin{equation} J^{-1}\lambda_{j}\lambda_{k}\frac{\partial^{2}W}{\partial\lambda_{j}\partial\lambda_{k}}=\frac{3\kappa-2\mu}{3J}\,.\label{eq:NNH-elasticity-coef2} \end{equation} \]

Finally, in the limiting case when pairs of eigenvalues are repeated, we substitute

\[ \begin{equation} \lim_{\lambda_{k}^{2}\to\lambda_{j}^{2}}2\frac{\lambda_{k}^{2}\sigma_{j}-\lambda_{j}^{2}\sigma_{k}}{\lambda_{j}^{2}-\lambda_{k}^{2}}=\frac{2\mu}{J}-\frac{4\left(3\kappa+\mu\right)\ln\lambda_{j}+2\left(3\kappa-2\mu\right)\ln\lambda_{i}}{3J}\,.\label{eq:NNH-elasticity-coef3} \end{equation} \]

Ogden Unconstrained

The Ogden unconstrained material is defined using the following hyperelastic strain energy function:

\[ \begin{equation} W\left(\lambda_{1},\lambda_{2},\lambda_{3}\right)=\frac{1}{2}c_{p}\left(J-1\right)^{2}+\sum\limits_{k=1}^{N}\frac{c_{k}}{m_{k}^{2}}\left(\lambda_{1}^{m_{k}}+\lambda_{2}^{m_{k}}+\lambda_{3}^{m_{k}}-3-m_{k}\ln J\right)\,.\label{eq425} \end{equation} \]

Here, \(\lambda_{i}\) are the principal stretches and \(c_{p}\), \(c_{k}\) and \(m_{k}\) are material parameters.

The Cauchy stress tensor for this material may be obtained using the general formula for isotropic elasticity in principal directions given in (2.6-35), with

\[ \begin{equation} \sigma_{i}=c_{p}\left(J-1\right)+\sum\limits_{k=1}^{N}\frac{1}{J}\frac{c_{k}}{m_{k}}\left(\lambda_{i}^{m_{k}}-1\right)\,.\label{eq426} \end{equation} \]

Similarly, the spatial elasticity tensor is given by

\[ \begin{equation} \begin{aligned}\boldsymbol{\mathcal{C}} & =\sum\limits_{i=1}^{3}\left(c_{p}+\sum\limits_{k=1}^{N}\frac{1}{J}\frac{c_{k}}{m_{k}}\left[\left(m_{k}-2\right)\lambda_{i}^{m_{k}}+2\right]\right)\mathbf{a}_{i}\otimes\mathbf{a}_{i}\\ & +\sum\limits_{i=1}^{3}\sum\limits_{j=i+1}^{3}c_{p}\left(2J-1\right)\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{eq427} \end{equation} \]

where \(\mathbf{a}_{i}=\mathbf{n}_{i}\otimes\mathbf{n}_{i}\) and \(\mathbf{n}_{i}\) are the eigenvectors of \(\mathbf{b}\). In the limit when eigenvalues coincide,

\[ \begin{equation} \lim\limits_{\lambda_{j}\to\lambda_{i}}2\frac{\sigma_{i}\lambda_{j}^{2}-\sigma_{j}\lambda_{i}^{2}}{\lambda_{i}^{2}-\lambda_{j}^{2}}=2c_{p}\left(1-J\right)+\sum\limits_{k=1}^{N}\frac{1}{J}\frac{c_{k}}{m_{k}}\left[2+\left(m_{k}-2\right)\lambda_{i}^{m_{k}}\right]\,.\label{eq428} \end{equation} \]

In the reference configuration the elasticity tensor reduces to

\[ \begin{equation} \left.\boldsymbol{\mathcal{C}}\right|_{\mathbf{b}=\mathbf{I}}=c_{p}\mathbf{I}\otimes\mathbf{I}+\left(\sum\limits_{k=1}^{N}c_{k}\right)\mathbf{I}\odot\mathbf{I}\,,\label{eq429} \end{equation} \]

which has the form of Hooke's law for infinitesimal isotropic elasticity (see Section Linear Elasticity), with equivalent Lamé coefficients \(c_{p}\equiv\lambda\) and \(2\mu\equiv\sum\nolimits_{k=1}^{N}c_{k}\).

Holmes-Mow

The coupled hyperelastic strain-energy function for this material is given by ³,

\[ \begin{equation} \Psi\left(I_{1},I_{2},J\right)=\frac{1}{2}c\left(e^{Q}-1\right)\,,\label{eq430} \end{equation} \]

where \(I_{1}\)and \(I_{2}\)are the first and second invariants of the right Cauchy-Green tensor and \(J\)the jacobian of the deformation. Furthermore,

\[ \begin{equation} \begin{aligned}Q & =\frac{\beta}{\lambda+2\mu}\left[\left(2\mu-\lambda\right)\left(I_{1}-3\right)+\lambda\left(I_{2}-3\right)-\left(\lambda+2\mu\right)\ln J^{2}\right]\,,\\ c & =\frac{\lambda+2\mu}{2\beta}\,, \end{aligned} \label{eq431} \end{equation} \]

and \(\lambda\) and \(\mu\) are the Lamé parameters. The corresponding Cauchy stress tensor is

\[ \begin{equation} \boldsymbol{\sigma}=\frac{1}{2J}e^{Q}\left(\left[2\mu+\lambda\left(I_{1}-1\right)\right]\mathbf{b}-\lambda\mathbf{b}^{2}-\left(\lambda+2\mu\right)\mathbf{I}\right)\,,\label{eq432} \end{equation} \]

and the spatial elasticity tensor is

\[ \begin{equation} \boldsymbol{\mathcal{C}}=\frac{4\beta}{\lambda+2\mu}Je^{-Q}\boldsymbol{\sigma}\otimes\boldsymbol{\sigma}+J^{-1}e^{Q}\left[\lambda\left(\mathbf{b}\otimes\mathbf{b}-\mathbf{b}\odot\mathbf{b}\right)+\left(\lambda+2\mu\right)\mathbf{I}\odot\mathbf{I}\right]\,.\label{eq433} \end{equation} \]

In the limit as \(\beta\to0\) we find that

\[ \begin{equation} \lim_{\beta\to0}\Psi\left(I_{1},I_{2},J\right)=\frac{1}{4}\left[2\mu\left(I_{1}-3\right)+\lambda\left(I_{2}-I_{1}\right)-\left(\lambda+2\mu\right)\ln J^{2}\right]\label{eq:433-1} \end{equation} \]

in which case

\[ \begin{equation} \boldsymbol{\sigma}=\frac{1}{2J}\left(\left(2\mu-\lambda\right)\mathbf{b}+\lambda\left(I_{1}\mathbf{b}-\mathbf{b}^{2}\right)-\left(\lambda+2\mu\right)\left(\ln J\right)\mathbf{I}\right)\label{eq:433-2} \end{equation} \]

Conewise Linear Elasticity

Curnier et al. ⁴ formulated a model for describing bimodular elastic solids exhibiting orthotropic material symmetry. This can be derived from the following hyperelastic strain-energy function:

\[ \begin{equation} W=\sum\limits_{a=1}^{3}\mu_{a}\mathbf{A}_{a}^{0}:\mathbf{E}^{2}+\frac{1}{2}\lambda_{aa}\left[\mathbf{A}_{a}^{0}:\mathbf{E}\right]\left(\mathbf{A}_{a}^{0}:\mathbf{E}\right)+\sum\limits_{b=1,\,b\ne a}^{3}\frac{1}{2}\lambda_{ab}\left(\mathbf{A}_{a}^{0}:\mathbf{E}\right)\left(\mathbf{A}_{b}^{0}:\mathbf{E}\right)\,,\label{eq434} \end{equation} \]

where \(\mathbf{A}_{a}^{0}=\mathbf{a}_{a}^{0}\otimes\mathbf{a}_{a}^{0}\) is the structural tensor corresponding to one of the three mutually orthogonal planes of symmetry whose unit outward normal is \(\mathbf{a}_{a}^{0}\) (\(\mathbf{a}_{a}^{0}\cdot\mathbf{a}_{b}^{0}=\delta_{ab})\). The bimodular response is described by

\[ \begin{equation} \lambda_{aa}\left[\mathbf{A}_{a}^{0}:\mathbf{E}\right]=\begin{cases} \lambda_{+aa} & \mathbf{A}_{a}^{0}:\mathbf{E}\geqslant0\\ \lambda_{-aa} & \mathbf{A}_{a}^{0}:\mathbf{E}<0 \end{cases}\,.\label{eq435} \end{equation} \]

The material constants are the three shear moduli \(\mu_{a}\), three tensile moduli \(\lambda_{+aa}\), three compressive moduli \(\lambda_{-aa}\), and three moduli \(\lambda_{ab}\) (\(b\ne a)\), where \(\lambda_{ba}=\lambda_{ab}\). The second Piola-Kirchhoff stress can be derived from this strain energy density function:

\[ \begin{equation} \begin{aligned}\mathbf{S} & =\frac{\partial W}{\partial\mathbf{E}}=\sum\limits_{a=1}^{3}\mu_{a}\left(\mathbf{A}_{a}^{0}\cdot\mathbf{E}+\mathbf{E}\cdot\mathbf{A}_{a}^{0}\right)\\ & +\lambda_{aa}\left[\mathbf{A}_{a}^{0}:\mathbf{E}\right]\left(\mathbf{A}_{a}^{0}:\mathbf{E}\right)\mathbf{A}_{a}^{0}+\sum\limits_{b=1,\,b\ne a}^{3}\lambda_{ab}\left(\mathbf{A}_{a}^{0}:\mathbf{E}\right)\mathbf{A}_{b}^{0}\,. \end{aligned} \label{eq436} \end{equation} \]

The material elasticity tensor is then given by,

\[ \begin{equation} \begin{aligned} & \boldsymbol{\mathbb{C}}=\frac{\partial\mathbf{S}}{\partial\mathbf{E}}=\sum\limits_{a=1}^{3}\mu_{a}\left(\mathbf{A}_{a}^{0}\odot\mathbf{I}+\mathbf{I}\odot\mathbf{A}_{a}^{0}\right)\\ & +\lambda_{aa}\left[\mathbf{A}_{a}^{0}:\mathbf{E}\right]\mathbf{A}_{a}^{0}\otimes\mathbf{A}_{a}^{0}+\sum\limits_{b=1,\,b\ne a}^{3}\lambda_{ab}\mathbf{A}_{a}^{0}\otimes\mathbf{A}_{b}^{0} \end{aligned} \,.\label{eq437} \end{equation} \]

It is important to note that although this model is objective, it should only be used for small strains. For large strains, the response may be unrealistic. The Cauchy stress is

\[ \begin{equation} \begin{aligned}\boldsymbol{\sigma} & =J^{-1}\left(\sum\limits_{a=1}^{3}\frac{\mu_{a}}{2}\left(\mathbf{A}_{a}\cdot\left(\mathbf{b}-\mathbf{I}\right)+\left(\mathbf{b}-\mathbf{I}\right)\cdot\mathbf{A}_{a}\right)\right.\\ & \left.+\lambda_{aa}\left[K_{a}\right]K_{a}\mathbf{A}_{a}+\sum\limits_{b=1,\,b\ne a}^{3}\lambda_{ab}K_{a}\mathbf{A}_{b}\right)\,, \end{aligned} \label{eq438} \end{equation} \]

where \(\mathbf{A}_{a}=\mathbf{F}\cdot\mathbf{A}_{a}^{0}\cdot\mathbf{F}^{T}\) and \(K_{a}=\frac{1}{2}\left({\mathbf{A}_{a}:\mathbf{I}-1}\right)\). The spatial elasticity tensor is

\[ \begin{equation} \boldsymbol{\mathcal{C}}=J^{-1}\left(\sum\limits_{a=1}^{3}\mu_{a}\left(\mathbf{A}_{a}\odot\mathbf{b}+\mathbf{b}\odot\mathbf{A}_{a}\right)+\lambda_{aa}\left[K_{a}\right]\mathbf{A}_{a}\otimes\mathbf{A}_{a}+\sum\limits_{b=1,\,b\ne a}^{3}\lambda_{ab}\mathbf{A}_{a}\otimes\mathbf{A}_{b}\right)\,.\label{eq439} \end{equation} \]

In the special case of cubic symmetry the number of material constants reduces to four,

\[ \begin{equation} \begin{aligned}\lambda_{+11} & =\lambda_{+22}=\lambda_{+33}\equiv\lambda_{+1}\\ \lambda_{-11} & =\lambda_{-22}=\lambda_{-33}\equiv\lambda_{-1}\\ \lambda_{12} & =\lambda_{23}=\lambda_{31}\equiv\lambda_{2}\\ \mu_{1} & =\mu_{2}=\mu_{3}\equiv\mu \end{aligned} \,.\label{eq440} \end{equation} \]

Donnan Equilibrium Swelling

The swelling pressure is described by the equations for ideal Donnan equilibrium, assuming that the material is porous, with a charged solid matrix, and the external bathing environment consists of a salt solution of monovalent counter-ions. Since osmotic swelling must be resisted by a solid material, this material is not stable on its own. It must be combined with an elastic material that resists the swelling.

The Cauchy stress for this material is the stress from the Donnan equilibrium response ⁵:

\[ \begin{equation} \boldsymbol{\sigma}=-\pi\,\mathbf{I}\,,\label{eq441} \end{equation} \]

where \(\pi\) is the osmotic pressure, given by

\[ \begin{equation} \pi=R\theta\left(\sqrt{\left(c^{F}\right)^{2}+\left(\bar{c}^{\ast}\right)^{2}}-\bar{c}^{\ast}\right),\label{eq442} \end{equation} \]

\(\bar{c}^{*}\) is the bath osmolarity (twice the concentration) and \(c^{F}\) is the fixed charge density in the current configuration, related to the reference configuration via,

\[ \begin{equation} c^{F}=\frac{\varphi_{0}^{w}}{J-1+\varphi_{0}^{w}}c_{0}^{F}\,,\label{eq443} \end{equation} \]

where \(J=\det\mathbf{F}\) is the relative volume, \(R\) is the universal gas constant and \(\theta\) is the absolute temperature.

Note that \(c_{0}^{F}\) may be negative or positive. The gel porosity \(\varphi_{0}^{w}\) is unitless and must be in the range \(0<\varphi_{0}^{w}<1\). The corresponding spatial elasticity tensor is ⁶

\[ \begin{equation} \begin{aligned}\boldsymbol{\mathcal{C}} & =\frac{R\theta J\left(c^{F}\right)^{2}}{\left(J-1+\varphi_{0}^{w}\right)\sqrt{\left(c^{F}\right)^{2}+\left(\bar{c}^{\ast}\right)^{2}}}\mathbf{I}\otimes\mathbf{I}\\ & +R\theta\left[\sqrt{\left(c^{F}\right)^{2}+\left(\bar{c}^{\ast}\right)^{2}}-\bar{c}^{\ast}\right]\left(2\mathbf{I}\odot\mathbf{I}-\mathbf{I}\otimes\mathbf{I}\right)\,. \end{aligned} \label{eq444} \end{equation} \]

Perfect Osmometer Equilibrium Osmotic Pressure

The swelling pressure is described by the equations for a perfect osmometer, assuming that the material is porous, containing an interstitial solution whose solutes cannot be exchanged with the external bathing environment. Similarly, solutes in the external bathing solution cannot be exchanged with the interstitial fluid of the porous material. Therefore, osmotic pressurization occurs when there is an imbalance between the interstitial and bathing solution osmolarities. Since osmotic swelling must be resisted by a solid matrix, this material is not stable on its own. It must be combined with an elastic material that resists the swelling.

The Cauchy stress for this material is the stress from the perfect osmometer equilibrium response ⁷:

\[ \begin{equation} \boldsymbol{\sigma}=-\pi\,\mathbf{I}\,,\label{eq445} \end{equation} \]

where \(\pi\) is the osmotic pressure, given by

\[ \begin{equation} \pi=R\theta\left(\bar{c}-\bar{c}^{\ast}\right)\,.\label{eq446} \end{equation} \]

Here, \(R\) is the universal gas constant and \(\theta\) is the absolute temperature, \(\bar{c}^{\ast}\) is the external bath osmolarity and \(\bar{c}\) is the interstitial fluid osmolarity in the current configuration, related to the reference configuration osmolarity \(\bar{c}_{0}\) via,

\[ \begin{equation} \bar{c}=\frac{\varphi_{0}^{w}}{J-1+\varphi_{0}^{w}}\bar{c}_{0}\,.\label{eq447} \end{equation} \]

Though this material is porous, this is not a full-fledged poroelastic material. The behavior described by this material is strictly valid only after the transient response of interstitial fluid and solute fluxes has subsided. The corresponding spatial elasticity tensor is

\[ \begin{equation} \boldsymbol{\mathcal{C}}=R\theta\left[\frac{J\,\bar{c}}{J-1+\varphi_{0}^{w}}\mathbf{I}\otimes\mathbf{I}+\left(\bar{c}-\bar{c}^{\ast}\right)\left(2\mathbf{I}\odot\mathbf{I}-\mathbf{I}\otimes\mathbf{I}\right)\right]\,.\label{eq448} \end{equation} \]

Large Poisson's Ratio Ligament

This material captures the transversely isotropic behavior of tendon and ligaments while enforcing a large Poisson's ratio. The material utilizes a three part strain energy equation:

\[ W=W_{\text{fiber}}+W_{\text{matrix}}+W_{\text{vol}}\,, \]

where:

\[ \begin{aligned}W_{\text{fiber}} & =\frac{1}{2}\frac{c_{1}}{c_{2}}\left(e^{c_{2}\left(\lambda-1\right)^{2}-1}\right)\,,\\ W_{\text{matrix}} & =\frac{\mu}{2}\left(I_{1}-3\right)-\mu\ln\left(\sqrt{I_{3}}\right)\,,\\ W_{\text{vol}} & =\frac{\kappa}{2}\left(\ln\left(\frac{I_{5}-I_{1}I_{4}+I_{2}}{I_{4}^{2(m-v_{0})}e^{-4m\left(\lambda-1\right)}}\right)\right)^{2}\,. \end{aligned} \]

The transversely isotropic strain energy \(W_{\text{fiber}}\) takes into account the behavior of the collagen fibers. The isotropic strain energy \(W_{\text{matrix}}\) takes into account the mechanical contribution of the extrafibrillar matrix and provides the majority of support when loaded transverse to the fiber direction. The variables \(c_{\mathrm{1}}\), \(c_{2}\) and \(\mu\) are material parameters controlling the stress-strain response of the material.

The volumetric strain energy \(W_{\text{vol}}\) acts as a penalty term which enforces a Poisson's ratio based on user selection of the parameters \(m\) and \(v_{\mathrm{0}}\). The variable \(\kappa\) acts as a penalty parameter. Raising \(\kappa\) will cause the prescribed Poisson's ratio to be enforced. The Poisson's ratio in question is given by the following function:

\[ v_{\text{apparent}}=-\frac{\lambda^{m-v_{0}}e^{-m\left(\lambda-1\right)}-1}{\lambda-1}\,. \]

Porous Neo-Hookean Material

Consider a porous neo-Hookean material with referential porosity \(\varphi_{r}^{w}\). The pores are compressible but the skeleton is intrinsically incompressible. Thus, upon pore closure, the material behavior needs to switch from compressible to incompressible.

In the current configuration, the porosity is given by

\[ \varphi^{w}=\frac{J-1+\varphi_{r}^{w}}{J}\,. \]

We may define a new variable,

\[ \bar{J}\equiv\frac{J-1+\varphi_{r}^{w}}{\varphi_{r}^{w}}=\frac{J-\varphi_{r}^{s}}{1-\varphi_{r}^{s}}\,, \]

which represents the pore volume ratio. It is equal to \(1\) when \(J=1\) and is equal to \(J\) when \(\varphi_{r}^{w}=1\) (or \(\varphi_{r}^{s}=1-\varphi_{r}^{w}=0\)). Now,

\[ \varphi^{w}=\varphi_{r}^{w}\frac{\bar{J}}{J}\,, \]

and

\[ \frac{\partial\bar{J}}{\partial J}=\frac{1}{\varphi_{r}^{w}}\,. \]

Pore closure occurs when \(\varphi^{w}=0\), which corresponds to \(J=\varphi_{r}^{s}\) and \(\bar{J}=0\).

Let us also define a modified deformation gradient,

\[ \bar{\mathbf{F}}=\left(\frac{\bar{J}}{J}\right)^{1/3}\mathbf{F}\,, \]

such that \(\det\bar{\mathbf{F}}=\bar{J}\). Let the corresponding modified right Cauchy-Green tensor be given by

\[ \bar{\mathbf{C}}=\bar{\mathbf{F}}^{T}\cdot\bar{\mathbf{F}}=\left(\frac{\bar{J}}{J}\right)^{2/3}\mathbf{C}\,, \]

so that

\[ \frac{\partial\bar{\mathbf{C}}}{\partial\mathbf{C}}=\left(\frac{\bar{J}}{J}\right)^{2/3}\left(\frac{\varphi_{r}^{s}}{3\left(J-\varphi_{r}^{s}\right)}\mathbf{C}\otimes\mathbf{C}^{-1}+\mathbf{I}\odot\mathbf{I}\right)\,. \]

The constitutive relation for the strain energy density of the compressible porous neo-Hookean material may be given by

\[ \Psi_{r}=\frac{\mu}{2}\left(\bar{I}_{1}-3\right)-\mu\ln\bar{J}\,, \]

where \(\bar{I}_{1}=\text{tr}\bar{\mathbf{C}}\). This relation shows that the material develops an infinite strain energy density as \(\bar{J}\) approaches zero. From this expression, the 2nd Piola-Kirchhoff stress is given by

\[ \mathbf{S}=2\frac{\partial\Psi_{r}}{\partial\mathbf{C}}=\mu\left[\left(\frac{\bar{J}}{J}\right)^{2/3}\mathbf{I}+\frac{1}{J-\varphi_{r}^{s}}\left(\varphi_{r}^{s}\left(\frac{\bar{J}}{J}\right)^{2/3}\frac{I_{1}}{3}-J\right)\mathbf{C}^{-1}\right]\,. \]

When \(\mathbf{C}=\mathbf{I}\) we can verify that \(\mathbf{S}=\mathbf{0}\). The corresponding Cauchy stress is

\[ \boldsymbol{\sigma}=\frac{\mu}{J}\left[\left(\frac{\bar{J}}{J}\right)^{2/3}\mathbf{b}+\frac{1}{J-\varphi_{r}^{s}}\left(\varphi_{r}^{s}\left(\frac{\bar{J}}{J}\right)^{2/3}\frac{I_{1}}{3}-J\right)\mathbf{I}\right]\,, \]

where \(\mathbf{b}\) is the left Cauchy-Green tensor.

The material elasticity tensor is given by

\[ \begin{aligned}\mathbb{C} & =2\frac{\partial\mathbf{S}}{\partial\mathbf{C}}\\ & =\frac{2}{3}g\left(J\right)\left(\mathbf{I}\otimes\mathbf{C}^{-1}+\mathbf{C}^{-1}\otimes\mathbf{I}\right)+\left(J\frac{dg}{dJ}\frac{I_{1}}{3}+J\frac{dh}{dJ}\right)\mathbf{C}^{-1}\otimes\mathbf{C}^{-1}\\ & -2\left[g\left(J\right)\frac{I_{1}}{3}+h\left(J\right)\right]\mathbf{C}^{-1}\odot\mathbf{C}^{-1} \end{aligned} \,, \]

where

\[ \begin{aligned}f\left(J\right) & =\mu\left(\frac{\bar{J}}{J}\right)^{2/3}\\ g\left(J\right) & =\frac{\varphi_{r}^{s}}{J-\varphi_{r}^{s}}f\left(J\right)\\ h\left(J\right) & =-\mu\frac{J}{J-\varphi_{r}^{s}} \end{aligned} \,, \]

and

\[ \begin{aligned}J\frac{dg}{dJ} & =\mu\frac{\left(2\varphi_{r}^{s}-3J\right)\varphi_{r}^{s}}{3\left(J-\varphi_{r}^{s}\right)^{2}}\left(\frac{\bar{J}}{J}\right)^{2/3}\\ J\frac{dh}{dJ} & =\mu\frac{J\varphi_{r}^{s}}{\left(J-\varphi_{r}^{s}\right)^{2}} \end{aligned} \,. \]

Then, the spatial elasticity tensor may be evaluated as

\[ \boldsymbol{\mathcal{C}}=J^{-1}\left[\frac{2}{3}g\left(J\right)\left(\mathbf{b}\otimes\mathbf{I}+\mathbf{I}\otimes\mathbf{b}\right)+\left(J\frac{dg}{dJ}\frac{I_{1}}{3}+J\frac{dh}{dJ}\right)\mathbf{I}\otimes\mathbf{I}-2\left[g\left(J\right)\frac{I_{1}}{3}+h\left(J\right)\right]\mathbf{I}\odot\mathbf{I}\right]\,. \]

In the limit of infinitesimal strains and rotations, when \(\mathbf{b}=\mathbf{I}\) and \(J=1\), we find that

\[ \begin{aligned}f\left(1\right) & =\mu & J\frac{dg}{dJ} & =\mu\frac{\left(2\varphi_{r}^{s}-3\right)\varphi_{r}^{s}}{3\left(1-\varphi_{r}^{s}\right)^{2}}\\ g\left(J\right) & =\mu\frac{\varphi_{r}^{s}}{1-\varphi_{r}^{s}} & J\frac{dh}{dJ} & =\mu\frac{\varphi_{r}^{s}}{\left(1-\varphi_{r}^{s}\right)^{2}}\\ h\left(J\right) & =-\mu\frac{1}{1-\varphi_{r}^{s}}\mu \end{aligned} \]

and

\[ \boldsymbol{\mathcal{C}}=\frac{2\mu}{3}\left(\frac{1}{\left(\varphi_{r}^{w}\right)^{2}}-1\right)\mathbf{I}\otimes\mathbf{I}+2\mu\mathbf{I}\odot\mathbf{I}\,. \]

Thus, by comparison to a standard neo-Hookean material, this porous neo-Hookean material has an effective Young's modulus equal to

\[ E=\frac{3\mu}{1+\frac{1}{2}\left(\varphi_{r}^{w}\right)^{2}}\,, \]

and an effective Poisson's ratio equal to

\[ \nu=\frac{1-\left(\varphi_{r}^{w}\right)^{2}}{2+\left(\varphi_{r}^{w}\right)^{2}}\,. \]

The two material properties that need to be provided are \(E\) and the referential porosity \(\varphi_{r}^{w}\) (or referential solid volume fraction \(\varphi_{r}^{s}=1-\varphi_{r}^{w}\)). Poisson's ratio in the limit of infinitesimal strains is dictated by the porosity according to the above formula. In particular, a highly porous material (\(\varphi_{r}^{w}\to1\)) has an effective (infinitesimal strain) Poisson ratio that approaches zero (\(\nu\to0\)) and \(E\to2\mu\). A low porosity material (\(\varphi_{r}^{w}\to0\)) has \(\nu\to\frac{1}{2}\) and \(E\to3\mu\), which is the expected behavior of an incompressible neo-Hookean solid. Note that setting \(\varphi_{r}^{w}=0\) would not produce good numerical behavior, since the Cauchy stress in an incompressible material would need to be supplemented by a pressure term (a Lagrange multiplier that enforces the incompressibility constraint). Nevertheless, this compressible porous neo-Hookean material behaves well even for values of \(\varphi^{w}\) as low as \(\sim0.015\).

Cell Growth

The cell growth material implements a swelling pressure \(\pi\) such that the Cauchy stress is given by

\[ \boldsymbol{\sigma}=-\pi\mathbf{I}\,, \]

where

\[ \pi=RT\left(\frac{c_{r}}{J-\varphi_{r}^{s}}-c_{e}\right)\,. \]

Here, \(c_{r}\) represents the referential molar concentration of intracellular solutes (moles of solutes per mixture volume in the reference configuration), \(\varphi_{r}^{s}\) is the referential intracellular solid volume fraction, and \(J=\det\mathbf{F}\) is the determinant of the deformation gradient, representing the volume ratio in the current configuration. The extracellular osmolarity is \(c_{e}\). This model assumes that neither intracellular solutes nor extracellular solutes may transport across the cell membrane passively. When a cell divides, it must use active transport mechanisms to bring in membrane-impermeant extracellular solutes inside, some of which are converted into intracellular solid matrix (e.g., cytoskeletal structures). As the intracellular osmolarity increases, water is transported into the cell, thus causing it to swell. The process of cell division is not modeled explicitly in this continuum representation, though the net effect is that cell proliferation leads to an increase in intracellular osmotic pressure, which generally translates into an increase in volume (unless the cell growth is constrained significantly). In a cell growth model, the initial condition (when \(J=1\)) should be selected for \(c_{r}\) and \(\varphi_{r}^{s}\) such that \(\pi=0\), thus

\[ \frac{c_{r}}{1-\varphi_{r}^{s}}=c_{e}\quad\text{initial condition}\,. \]

The spatial elasticity tensor associated with this osmotic pressure is

\[ \boldsymbol{\Pi}=-\left(\pi+J\frac{\partial\pi}{\partial J}\right)\mathbf{I}\otimes\mathbf{I}+2\pi\mathbf{I}\odot\mathbf{I}\,, \]

where

\[ J\frac{\partial\pi}{\partial J}=-RT\frac{J}{\left(J-\varphi_{r}^{s}\right)^{2}}c_{r}\,. \]

This elasticity tensor has the same form as that of an isotropic elastic material whose effective Young's modulus \(E_{Y}\) and Poisson's ratio \(\nu\) are given by

\[ \begin{aligned}E_{Y} & =\pi\frac{\left(\pi+3J\frac{\partial\pi}{\partial J}\right)}{J\frac{\partial\pi}{\partial J}}\\ \nu & =\frac{\pi+J\frac{\partial\pi}{\partial J}}{2J\frac{\partial\pi}{\partial J}} \end{aligned} \,. \]

In the reference configuration, when \(\pi=0\), it follows that \(E_{Y}=0\) and \(\nu=\frac{1}{2}\).

\[ \begin{equation} \tilde{\boldsymbol{\sigma}}=\tilde{\boldsymbol{\sigma}}_{m}+\tilde{\boldsymbol{\sigma}}_{f}\,.\label{eq471-1} \end{equation} \]

Fiber with Exponential Power Law

This material model describes a constitutive model for fibers, where a single fiber family follows an exponential power law strain energy function. The Cauchy stress is given by:

\[ \begin{equation} \boldsymbol{\sigma}=2J^{-1}H\left(I_{n}-I_{0}\right)I_{n}\frac{\partial\Psi_{r}}{\partial I_{n}}\mathbf{n}\otimes\mathbf{n}\,,\label{eq:fepl-stress} \end{equation} \]

and the corresponding spatial elasticity tensor is

\[ \begin{equation} \boldsymbol{\mathcal{C}}=4J^{-1}H\left(I_{n}-I_{0}\right)I_{n}^{2}\frac{\partial^{2}\Psi_{r}}{\partial I_{n}^{2}}\mathbf{n}\otimes\mathbf{n}\otimes\mathbf{n}\otimes\mathbf{n}\,,\label{eq:fepl-elasticity} \end{equation} \]

where \(I_{n}=\lambda_{n}^{2}=\mathbf{n}_{r}\cdot\mathbf{C}\cdot\mathbf{n}_{r}\) is the square of the fiber stretch, \(\mathbf{n}_{r}\) is the fiber orientation in the reference configuration,

\[ \begin{equation} \mathbf{n}_{r}=\sin\varphi\cos\theta\,\mathbf{e}_{1}+\sin\varphi\sin\theta\,\mathbf{e}_{2}+\cos\varphi\,\mathbf{e}_{3}\,,\label{eq:fepl-fiber-nr} \end{equation} \]

and \(\mathbf{n}=\mathbf{F}\cdot\mathbf{n}_{r}/\lambda_{n}\). The function \(H\left(\cdot\right)\) is the unit step function that enforces the tension-only contribution. Thus, the stress and elasticity tensors are non-zero only when \(I_{n}>I_{0}\), where \(I_{0}=\lambda_{0}^{2}\) is the square of the stretch at which the fiber's tensile response engages. By default we may take \(I_{0}=1\), though the actual value of \(I_{0}\) may be set by the user. The fiber strain energy density is given by

\[ \begin{equation} \Psi_{r}=\frac{\xi}{\alpha\beta}\left(\exp\left[\alpha\left(I_{n}-I_{0}\right)^{\beta}\right]-1\right)\,,\label{eq:fepl-fiber-Psi} \end{equation} \]

where \(\xi>0\), \(\alpha\geqslant0\)and \(\beta\geqslant2\). From this expression we get

\[ \begin{equation} \begin{aligned}\frac{\partial\Psi_{r}}{\partial I_{n}} & =\xi\left(I_{n}-I_{0}\right)^{\beta-1}\exp\left[\alpha\left(I_{n}-I_{0}\right)^{\beta}\right]\\ \frac{\partial^{2}\Psi_{r}}{\partial I_{n}^{2}} & =\xi\left(\beta\left(1+\alpha\left(I_{n}-I_{0}\right)^{\beta}\right)-1\right)\left(I_{n}-I_{0}\right)^{\beta-2}\exp\left[\alpha\left(I_{n}-I_{0}\right)^{\beta}\right] \end{aligned} \,.\label{eq:fepl-fiber-dPsi} \end{equation} \]

Note: In the limit when \(\alpha\to0\), this expression produces a power law,

\[ \begin{equation} \lim\limits_{\alpha\to0}\Psi_{r}=\frac{\xi}{\beta}\left(I_{n}-I_{0}\right)^{\beta}\,.\label{eq:fepl-Psi-alpha0} \end{equation} \]

Note: According to \eqref{eq:fepl-elasticity} and \eqref{eq:fepl-fiber-dPsi}, when \(\beta>2\) the fiber modulus is zero at the strain origin (\(I_{n}=I_{0})\). Therefore, use \(\beta>2\) when a smooth transition in the stress is desired from compression to tension.

There is an option to also add a shear modulus \(\mu\) to account for the interaction of a fiber with the ground matrix. This additional contribution does not depend on whether the fiber is in tension. It has a strain energy density

\[ \begin{equation} W=\frac{\mu}{4}\left(K_{n}-1-2\left(I_{n}-1\right)\right)\,,\quad K_{n}=\mathbf{n}_{r}\cdot\mathbf{C}^{2}\cdot\mathbf{n}_{r}\label{eq:fepl-shear-sed} \end{equation} \]

The corresponding stress is

\[ \begin{equation} \begin{aligned}\boldsymbol{\sigma} & =2\frac{I_{n}}{J}\left(\frac{\partial W}{\partial I_{n}}\mathbf{n}\otimes\mathbf{n}+\frac{\partial W}{\partial K_{n}}\left(\mathbf{n}\otimes\mathbf{n}\cdot\mathbf{b}+\mathbf{b}\cdot\mathbf{n}\otimes\mathbf{n}\right)\right)\\ & =\frac{\mu}{2}\frac{I_{n}}{J}\left(\mathbf{n}\otimes\mathbf{n}\cdot\left(\mathbf{b}-\mathbf{I}\right)+\left(\mathbf{b}-\mathbf{I}\right)\cdot\mathbf{n}\otimes\mathbf{n}\right) \end{aligned} \label{eq:fepl-shear-stress} \end{equation} \]

where \(\mathbf{b}=\mathbf{F}\cdot\mathbf{F}^{T}\) is the left Cauchy-Green tensor. The elasticity tensor is

\[ \begin{equation} \begin{aligned}\boldsymbol{\mathcal{C}} & =4\frac{I_{n}^{2}}{J}\left(\frac{\partial^{2}W}{\partial I_{n}^{2}}\mathbf{N}\otimes\mathbf{N}+\frac{\partial^{2}W}{\partial K_{n}^{2}}\left(\mathbf{N}\cdot\mathbf{b}+\mathbf{b}\cdot\mathbf{N}\right)\otimes\left(\mathbf{N}\cdot\mathbf{b}+\mathbf{b}\cdot\mathbf{N}\right)\right)\\ & +4\frac{I_{n}^{2}}{J}\frac{\partial^{2}W}{\partial I_{n}\partial K_{n}}\left(\mathbf{N}\otimes\left(\mathbf{N}\cdot\mathbf{b}+\mathbf{b}\cdot\mathbf{N}\right)+\left(\mathbf{N}\cdot\mathbf{b}+\mathbf{b}\cdot\mathbf{N}\right)\otimes\mathbf{N}\right)\\ & +\frac{4I_{n}}{J}\frac{\partial W}{\partial K_{n}}\left(\mathbf{N}\odot\mathbf{b}+\mathbf{b}\odot\mathbf{N}\right)\\ & =\mu\frac{I_{n}}{J}\left(\mathbf{N}\odot\mathbf{b}+\mathbf{b}\odot\mathbf{N}\right) \end{aligned} \,,\label{eq:fepl-shear-elasticity} \end{equation} \]

where \(\mathbf{N}=\mathbf{n}\otimes\mathbf{n}\).

Fiber with Natural Neo-Hookean Response

This model is an adaptation of the natural neo-Hookean material presented in Section \eqref{subsec:Natural-Neo-Hookean}. Consider that the state of strain in a fiber is given by the unidirectional natural (left Hencky) strain along the fiber,

\[ \begin{equation} \boldsymbol{\eta}=\ln\lambda_{n}\,\mathbf{n}\otimes\mathbf{n}\label{eq:fnh-strain} \end{equation} \]

where \(\lambda_{n}\) is the stretch ratio along the current fiber unit vector \(\mathbf{n}\). For this special state of strain the invariants \(K_{i}\) of the natural strain tensor reduce to \(K_{1}=\ln\lambda_{n}\), \(K_{2}=\sqrt{3/2}\ln\lambda_{n}\), and \(K_{3}=1\). In this case, a natural neo-Hookean fiber response is given by

\[ \begin{equation} \Psi_{r}=\frac{\xi}{2}H\left(\ln\lambda_{n}\right)\left(\ln\lambda_{n}\right)^{2}\,.\label{eq:fnh-sed} \end{equation} \]

The corresponding Cauchy stress is

\[ \begin{equation} \boldsymbol{\sigma}=\xi J^{-1}H\left(\ln\lambda_{n}\right)\ln\lambda_{n}\mathbf{n}\otimes\mathbf{n}\label{eq:fnh-stress} \end{equation} \]

and the elasticity tensor is

\[ \begin{equation} \boldsymbol{\mathcal{C}}=\xi J^{-1}H\left(\ln\lambda_{n}\right)\left(1-2\ln\lambda_{n}\right)\mathbf{N}\otimes\mathbf{N}\label{eq:fnh-elasticity} \end{equation} \]

where \(\xi\) is the fiber modulus. This expression shows that the components of the elasticity tensor become negative when \(1-2\ln\lambda_{n}<0\), or equivalently when \(\lambda_{n}>e^{\frac{1}{2}}\). However, the stress always remains positive.

If we want the tensile response to engage only beyond a threshold stretch ratio \(\lambda_{0}\), we may rewrite the strain energy density as

\[ \begin{equation} \Psi_{r}=\frac{\xi}{2}H\left(\ln\lambda_{n}-\ln\lambda_{0}\right)\left(\ln\lambda_{n}-\ln\lambda_{0}\right)^{2}=\frac{\xi}{2}H\left(\ln\frac{\lambda_{n}}{\lambda_{0}}\right)\left(\ln\frac{\lambda_{n}}{\lambda_{0}}\right)^{2}\,.\label{eq:fnh-sed-redux} \end{equation} \]

Then

\[ \begin{equation} \boldsymbol{\sigma}=\xi J^{-1}H\left(\ln\frac{\lambda_{n}}{\lambda_{0}}\right)\ln\frac{\lambda_{n}}{\lambda_{0}}\mathbf{n}\otimes\mathbf{n}\,,\label{eq:fnh-stress-redux} \end{equation} \]

and

\[ \begin{equation} \boldsymbol{\mathcal{C}}=\xi J^{-1}H\left(\ln\frac{\lambda_{n}}{\lambda_{0}}\right)\left(1-2\ln\frac{\lambda_{n}}{\lambda_{0}}\right)\mathbf{N}\otimes\mathbf{N}\,.\label{eq:fnh-elasticity-redux} \end{equation} \]

HGO Unconstrained

A coupled formulation of the Holzapfel-Gasser-Ogden material is implemented in the “HGO unconstrained” material.

The strain energy density is given by,

\[ \Psi_{r}=\frac{c}{2}\left(I_{1}-3\right)-c\ln J+\frac{k_{1}}{2k_{2}}\sum_{\alpha}\left(\exp\left(k_{2}\left\langle E_{\alpha}\right\rangle ^{2}\right)-1\right)+\frac{K_{0}}{2}\left(\frac{J^{2}-1}{2}-\ln J\right) \]

The fiber strain is

\[ E_{\alpha}=\kappa\left(I_{1}-3\right)+\left(1-3\kappa\right)\left(I_{4\alpha}-1\right) \]

where \(I_{1}=\tr\mathbf{C}\) and \(I_{4\alpha}=\mathbf{a}_{\alpha r}\cdot\mathbf{C}\cdot\mathbf{a}_{\alpha r}\). The Macaulay brackets around \(\left\langle E_{\alpha}\right\rangle\) indicate that this term is zero when \(E_{\alpha}<0\) and equal to \(E_{\alpha}\) when this strain is positive.There are two fiber families along the vectors \(\mathbf{a}_{\alpha r}\) (\(\alpha=1,2\)), lying in the \(\left\{ \mathbf{e}_{1},\mathbf{e}_{2}\right\}\) plane of the local material axes \(\left\{ \mathbf{e}_{1},\mathbf{e}_{2},\mathbf{e}_{3}\right\}\), making an angle \(\pm\gamma\) with respect to \(\mathbf{e}_{1}\). Each fiber family has a dispersion \(\kappa\), where \(0\le\kappa\le\frac{1}{3}\). When \(\kappa=0\) there is no fiber dispersion, implying that all the fibers in that family act along the angle \(\pm\gamma\); the value \(\kappa=\frac{1}{3}\) represents an isotropic fiber dispersion.

The Cauchy stress is given by,

\[ \mathbf{\sigma}=\frac{1}{J}\left[c\left(\mathbf{b}-\mathbf{1}\right)+\frac{k}{2}\left(J^{2}-1\right)\mathbf{1}\right]+\frac{1}{J}\left[{\displaystyle \underset{\alpha}{\sum}2k_{1}\left\langle E_{\alpha}\right\rangle \exp\left(k_{2}\left\langle E_{\alpha}\right\rangle ^{2}\right)\mathbf{h}_{\alpha}}\right] \]

where, \(\mathbf{h}_{\alpha}=\kappa\mathbf{b}+\left(1-3\kappa\right)\mathbf{A}_{\alpha}\), and \(\mathbf{A}_{\alpha}=\mathbf{a}_{\alpha}\otimes\mathbf{a}_{\alpha}\).

The spatial elasticity tensor is given by,

\[ \boldsymbol{c}=\frac{1}{J}\left\{ kJ^{2}\mathbf{1}\otimes\mathbf{1}+\left[2c-k\left(J^{2}-1\right)\right]\mathbf{1}\odot\mathbf{1}\right\} +\frac{1}{J}\left[{\displaystyle \underset{\alpha}{\sum}4k_{1}\left(1+2k_{2}\left\langle E_{\alpha}\right\rangle ^{2}\right)\exp\left(k_{2}\left\langle E_{\alpha}\right\rangle ^{2}\right)\mathbf{h}_{\alpha}\otimes\mathbf{h}_{\alpha}}\right] \]

---

1. Bonet, Javier; Wood, Richard D.. "Nonlinear continuum mechanics for finite element analysis." Cambridge University Press (1997). ↩↩↩

2. Criscione, John C; Humphrey, Jay D; Douglas, Andrew S; Hunter, William C. "An invariant basis for natural strain which yields orthogonal stress response terms in isotropic hyperelasticity." J. Mech. Phys. Solids, vol. 48, pp. 2445--2465 (2000). ↩↩

3. Holmes, M. H.; Mow, V. C.. "The nonlinear characteristics of soft gels and hydrated connective tissues in ultrafiltration." J Biomech, vol. 23, pp. 1145-56 (1990). ↩

4. Curnier, A.; Qi-Chang, He; Zysset, P.. "Conewise linear elastic materials." J Elasticity, vol. 37, pp. 1-38 (1994). ↩

5. Ateshian, G. A.; Rajan, V.; Chahine, N. O.; Canal, C. E.; Hung, C. T.. "Modeling the matrix of articular cartilage using a continuous fiber angular distribution predicts many observed phenomena." J Biomech Eng, vol. 131, pp. 061003 (2009). ↩

6. Azeloglu, Evren U; Albro, Michael B; Thimmappa, Vikrum A; Ateshian, Gerard A; Costa, Kevin D. "Heterogeneous transmural proteoglycan distribution provides a mechanism for regulating residual stresses in the aorta." Am J Physiol Heart Circ Physiol, vol. 294, pp. H1197-205 (2008). ↩

7. Ateshian, G. A.; Likhitpanichkul, M.; Hung, C. T.. "A mixture theory analysis for passive transport in osmotic loading of cells." J Biomech, vol. 39, pp. 464-75 (2006). ↩