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.

Convex and non-convex energy comparison

The global existence theory of the solution in elasticity is based on the polyconvexity of the strain energy function [Hol00]. To require this condition, the volumetric part of the strain energy function has to satisfy the convexity condition

\[ \frac{\partial^2 \psi_{vol}}{\partial J^2} \geq 0. \]

One of the commonly used form for \(\psi_{vol}\) is

\[ \psi_{vol}(J) = \frac{\lambda}{2} (\log J)^2 \]

which is valid when \(J \approx 1\). We apply traction and slip boundary conditions on the faces of the block so the deformation is volumetric only. The convex strain energy (13) and the non-convex form are plotted vs time increment.

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.


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.


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).


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).


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). \]


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).