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

The compressible Neo-Hookean model is a function of two invariants \(\mathbb{I}_1(\bm{C})\) and \(J\)

(34)#\[ \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 (34) and the non-convex form are plotted vs time increment.

Click to show code
import altair as alt
import pandas as pd
def source_path(rel):
    import os
    return os.path.join(os.path.dirname(os.environ["DOCUTILSCONFIG"]), rel)

nh_convex = pd.read_csv(source_path("doc/data/NH-convex-energy.csv"))
nh_convex["model"] = "Convex Neo-Hookean"
nh_convex["parameters"] = "E=2.8, nu=0.4"

nh_nonconvex = pd.read_csv(source_path("doc/data/NH-nonconvex-energy.csv"))
nh_nonconvex["model"] = "Non-convex Neo-Hookean"
nh_nonconvex["parameters"] = "E=2.8, nu=0.4"

df = pd.concat([nh_convex, nh_nonconvex])
highlight = alt.selection_point(
   on = "mouseover",
   nearest = True,
   fields=["model", "parameters"],
base = alt.Chart(df).encode(
   alt.Y("energy", scale=alt.Scale(type="sqrt")),
   alt.Tooltip(("model", "parameters")),
   opacity=alt.condition(highlight, alt.value(1), alt.value(.5)),
   size=alt.condition(highlight, alt.value(2), alt.value(1)),
base.mark_point().add_params(highlight) + base.mark_line()

To evaluate the gradient of the strain energy density functional (18) for the model (34), we make use of (19), (32) and (33) to derive

(35)#\[ \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 (35) 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 (35) reduces to

(36)#\[ \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

(37)#\[ \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, (37) 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

(38)#\[ \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

(39)#\[ \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 (35) for Neo-Hookean materials as

(40)#\[ \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 (12). 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 (35) using (7)

(41)#\[ \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 (41) by expressing (34) in terms of current configuration invariant \(\operatorname{trace} \bm{e}\) and use equation (20).


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

(42)#\[ \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 (38).