Neo-Hookean Materials#

In the total Lagrangian approach for the Neo-Hookean constitutive model, the model is formulated with respect to the initial configuration. Our notation for the Neo-Hookean constitutive model is inspired by [Hol00] to distinguish between the current and initial configurations. As explained in the Continuum Mechanics section, we denote coordinates in the reference frame by capital letters and the current frame by lowercase letters.

For the constitutive modeling of Neo-Hookean hyperelasticity at finite strain, we will assume without loss of generality that \(\bm{E}\) is diagonal and take its set of eigenvalues as the invariants. It is clear that there can be only three invariants, and there are many alternate choices, such as \(\operatorname{trace} \left( \bm{E} \right), \operatorname{trace}\left( \bm{E}^2 \right), \lvert \bm{E} \rvert\), and combinations thereof. It is common in the literature for invariants to be taken from \(\bm{C} = \bm{I}_3 + 2 \bm{E}\) instead of \(\bm{E}\).

For example, if we take the compressible Neo-Hookean model,

(13)#\[ \begin{aligned} \psi \left(\bm{E} \right) &= \frac{\lambda}{4} \left( J^2 - 1 -2 \log J \right) - \mu \log J + \frac \mu 2 \left( \operatorname{trace} \bm{C} - 3 \right) \\ &= \frac{\lambda}{4} \left( J^2 - 1 -2 \log J \right) - \mu \log J + \mu \operatorname{trace} \bm{E}, \end{aligned} \]

where \(J = \lvert \bm{F} \rvert = \sqrt{\lvert \bm{C} \rvert}\) is the determinant of deformation (i.e., volume change (6)) and \(\lambda\) and \(\mu\) are the Lamé parameters in the infinitesimal strain limit.

To evaluate the gradient of the strain energy density functional (12), we make use of

\[ \frac{\partial J}{\partial \bm{E}} = \frac{\partial \sqrt{\lvert \bm{C} \rvert}}{\partial \bm{E}} = \lvert \bm{C} \rvert^{-1/2} \lvert \bm{C} \rvert \bm{C}^{-1} = J \bm{C}^{-1}, \]

where the factor of \(\frac{1}{2}\) has been absorbed due to \(\bm{C} = \bm{I}_3 + 2 \bm{E}.\) Carrying through the differentiation (12) for the model (13), we arrive at

(14)#\[ \bm{S} = \frac{\lambda}{2} \left( J^2 - 1 \right)\bm{C}^{-1} + \mu \left( \bm{I}_3 - \bm{C}^{-1} \right), \]

which is the second Piola-Kirchoff tensor for Neo-Hookean materials.

Note

One can linearize (14) around \(\bm{E} = 0\), for which \(\bm{C} = \bm{I}_3 + 2 \bm{E} \to \bm{I}_3\) and \(J \to 1 + \operatorname{trace} \bm{E}\), therefore (14) reduces to

(15)#\[ \bm{S} = \frac{\lambda}{2} \operatorname{trace} \bm{E} \left( \operatorname{trace} \bm{E} + 2\right) \bm{I}_3 + 2 \mu \bm{E}, \]

which is the St. Venant-Kirchoff model (constitutive linearization without geometric linearization). This model can be used for geometrically nonlinear mechanics (e.g., snap-through of thin structures), but is inappropriate for large strain.

Tip

In most of the constitutive models we have \(f(J) = \log J\) and \(g(J) = J -1\) functions with the condition numbers

(16)#\[ \kappa_f(J) = \frac{1}{\lvert \log J \rvert}, \quad \kappa_g(J) = \frac{\lvert J \rvert}{\lvert J-1 \rvert} \]

where in the case of nearly incompressible materials, (16) becomes very large and thus ill-conditioned. To compute these functions in a numerically stable way, suppose we have the \(2\times 2\) non-symmetric matrix \(\bm{F} = \left( \begin{smallmatrix} 1 + u_{0,0} & u_{0,1} \\ u_{1,0} & 1 + u_{1,1} \end{smallmatrix} \right)\). Then we compute

(17)#\[ \mathtt{J_{-1}} = J-1 = u_{0,0} + u_{1,1} + u_{0,0} u_{1,1} - u_{0,1} u_{1,0}, \]

and if \(\log J\) appears in a stress we use log1p as

(18)#\[ \log J = \mathtt{log1p} \left( \mathtt{J_{-1}} \right), \]

which gives accurate results even in the limit when the entries \(u_{i,j}\) are very small or \(J\approx 1\).

Another source of error in constitutive models is loss of significance which occurs in \(\bm{I}_3 - \bm{C}^{-1}\) term. We use an equivalent form of (14) for Neo-Hookean materials as

(19)#\[ \bm{S} = \frac{\lambda}{2} \mathtt{J_{-1}} \left(\mathtt{J_{-1}} + 2 \right) \bm{C}^{-1} + 2 \mu \bm{C}^{-1} \bm{E}, \]

which is more numerically stable for small \(\bm{E}\), and thus preferred for computation. Note that the product \(\bm{C}^{-1} \bm{E}\) is also symmetric, and that \(\bm{E}\) should be computed using (9). For example, if \(u_{i,j} \sim 10^{-8}\), then naive computation of \(\bm{I}_3 - \bm{C}^{-1}\) and \(\log J\) will have a relative accuracy of order \(10^{-8}\) in double precision and no correct digits in single precision. When using the stable choices above, these quantities retain full \(\varepsilon_{\text{machine}}\) relative accuracy.

The Neo-Hookean stress relation can be represented in current configuration by pushing forward (14) using (7)

(20)#\[ \bm{\tau} = \bm{F}\bm{S}\bm{F}^T = \frac{\lambda}{2} \left( J^2 - 1 \right) \bm{I}_{3} + \mu \left( \bm{b} - \bm{I}_{3} \right), \]

One can arrive at the same relation (20) by expressing (13) in terms of current configuration invariant \(\operatorname{trace} \bm{e}\) and use the second equation of (12).

Tip

In our simulation we use the stable version of Kirchhoff stress tensor (20) as

(21)#\[ \bm{\tau} = \frac{\lambda}{2} \mathtt{J_{-1}} \left(\mathtt{J_{-1}} + 2 \right) \bm{I}_{3} + 2 \mu \bm{e}. \]

where \(\mathtt{J_{-1}}\) is computed by (17).

Isochoric-split#

For the incompressible and nearly incompressible materials where the motions are constraint by spatial conditions, it is most beneficial to split the deformation locally into a so-called volumetric (dilational) part and an isochoric (distortional) part. In order to be able to define strain energy functions that are separable into volumetric and isochoric parts, we decompose the deformation gradient in a multiplicative sense by

(22)#\[ \bm{F} = (J^{1/3} \bm{I}_3) \bar{\bm{F}} = J^{1/3} \bar{\bm{F}} \]

where \(J^{1/3} \bm{I}_3\) describes the purely volumetric deformation and \(\bar{\bm{F}}\) captures the isochoric or volume-preserving deformation since \(\lvert \bar{\bm{F}} \rvert = \lvert J^{-1/3} \bm{F} \rvert = 1\) . Similar decomposition can be obtained for other tensor such as the right Cauchy-Green tensor \(\bm{C}\) by defining as

(23)#\[ \bar{\bm{C}} = \bar{\bm{F}}^T \bar{\bm{F}} = J^{-2/3} \bm{C} \]

where we call \(\bar{\bm{F}}\) and \(\bar{\bm{C}}\) the modified deformation gradient and the modified right Cauchy-Green tensor, respectively.

Based on our kinematic assumption (23) we can postulate a decoupled strain energy function for incompressible Neo-Hookean model as

(24)#\[ \begin{aligned} \psi \left(\bm{C} \right) &= \psi_{vol}(J) + \psi_{iso}(\bar{\bm{C}})\\ &=\frac{\kappa}{4} \left( J^2 - 1 -2 \log J \right) + \frac \mu 2 \left( \mathbb{\bar{I}_1} - 3 \right). \end{aligned} \]

where \(\mathbb{\bar{I}_1} = \operatorname{trace} \bar{\bm{C}}\) and \(\kappa\) is the bulk modulus of the material.

The second Piola-Kirchoff tensor is computed via

(25)#\[ \begin{aligned} \bm{S} = \frac{\partial \psi}{\partial \bm{E}} &= \frac{\partial \psi_{vol}}{\partial J} \frac{\partial J}{\partial \bm{E}} + \frac{\partial \psi}{\partial \mathbb{\bar{I}_1}} \frac{\partial \mathbb{\bar{I}_1}}{\partial \bm{E}} =\bm{S}_{vol} + \bm{S}_{iso}\\ &= -p J \bm{C}^{-1} + \mu J^{-2/3}\left( \bm{I}_3 - \frac{1}{3}\mathbb{{I}_1}\bm{C}^{-1} \right). \end{aligned} \]

where \(\mathbb{{I}_1} = \operatorname{trace} \bm{C}, \mathbb{\bar{I}_1} = J^{-2/3} \mathbb{{I}_1}, \) and we have used

\[ \begin{aligned} \frac{\partial J}{\partial \bm{E}} &= J \bm{C}^{-1}, & \frac{\partial \mathbb{{I}_1}}{\partial \bm{E}} &= 2 \bm{I}_3. \end{aligned} \]

The pressure \(p\) in (25) is defined by

(26)#\[ p = - \frac{\partial \psi_{vol}}{\partial J} = -\frac{\kappa}{2 J} (J^2 -1). \]

Tip

In our simulation we use the stable version of (25) and (26) as

(27)#\[ \bm{S} = -p J \bm{C}^{-1} + 2 \mu J^{-2/3} \bm{C}^{-1} \bm{E}_{dev}, \]

(28)#\[ p = -\frac{\kappa}{2 J} \mathtt{J_{-1}} \left(\mathtt{J_{-1}} + 2 \right), \]

where \(\bm{E}_{dev} = \bm{E} - \frac{1}{3}\operatorname{trace} \bm{E} \, \bm{I}_{3} \) is the deviatoric part of Green-Lagrange strain tensor and \(\mathtt{J_{-1}}\) is computed by (17).