# Neo-Hookean

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 {cite}`holzapfel2000nonlinear` to distinguish between the current and initial configurations.
As explained in the {ref}`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 $\trace \left( \bm{E} \right), \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} + 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$

$$
  \begin{aligned}
    \psi \left(\bm{E} \right) &= \firstlame V(J) - \secondlame \log J + \frac \secondlame 2 \left( \trace \bm{C} - 3 \right) \\
                             &= \frac{\firstlame}{4} \left( J^2 - 1 -2 \log J \right) - \secondlame \log J + \secondlame \trace \bm{E},
  \end{aligned}
$$ (neo-hookean-energy)

where $J = \lvert \bm{F} \rvert = \sqrt{\lvert \bm{C} \rvert}$ is the determinant of deformation (i.e., volume change {eq}`volume-transformation-and-J`) and $\firstlame$ and $\secondlame$ are the Lamé parameters in the infinitesimal strain limit.

:::{dropdown} 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 {cite}`holzapfel2000nonlinear`.
To require this condition, the volumetric part of the strain energy function has to satisfy the convexity condition

$$
  \frac{\partial^2 \psi_{\text{vol}}}{\partial J^2} \geq 0.
$$

One of the commonly used form for $\psi_{\text{vol}}$ is

$$
  \psi_{\text{vol}}(J) = \firstlame V(J) = \frac{\firstlame}{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 {eq}`neo-hookean-energy` and the non-convex form are plotted vs time increment.

```{altair-plot}
:hide-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.X("increment"),
   alt.Y("energy", scale=alt.Scale(type="sqrt")),
   alt.Color("model"),
   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 {eq}`strain-energy-grad` for the model {eq}`neo-hookean-energy`, we make use of {eq}`invariants-derivative`, {eq}`derivative-J-wrt-C` and {eq}`derivative-wrt-E` to derive

$$
  \bm{S} = \firstlame J V' \bm{C}^{-1} + \secondlame \left( \bm{I} - \bm{C}^{-1} \right) = \frac{\firstlame}{2} \left( J^2 - 1 \right)\bm{C}^{-1} + \secondlame \left( \bm{I} - \bm{C}^{-1} \right),
$$ (neo-hookean-stress)

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

:::{note}
One can linearize {eq}`neo-hookean-stress` around $\bm{E} = 0$, for which $\bm{C} = \bm{I} + 2 \bm{E} \to \bm{I}$ and $J \to 1 + \trace \bm{E}$, therefore {eq}`neo-hookean-stress` reduces to

$$
\bm{S} = \firstlame \trace \bm{E} \bm{I} + 2 \secondlame \bm{E},
$$ (eq-st-venant-kirchoff)

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

$$
\kappa_f(J) = \frac{1}{\lvert \log J \rvert}, \quad \kappa_g(J) = \frac{\lvert J \rvert}{\lvert J-1 \rvert}
$$ (condition-number-logJ)

where in the case of nearly incompressible materials, {eq}`condition-number-logJ` 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

$$
\mathtt{J_{-1}} = J-1 = u_{0,0} + u_{1,1} + u_{0,0} u_{1,1} - u_{0,1} u_{1,0},
$$ (JM1)

and if $\log J$ appears in a stress we use `log1p` as

$$
\log J = \mathtt{log1p} \left( \mathtt{J_{-1}} \right),
$$ (log1p)

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} - \bm{C}^{-1}$ term.
We use an equivalent form of {eq}`neo-hookean-stress` for Neo-Hookean materials as

$$
\bm{S} = \frac{\firstlame}{2} \mathtt{J_{-1}} \left(\mathtt{J_{-1}} + 2 \right) \bm{C}^{-1} + 2 \secondlame \bm{C}^{-1} \bm{E},
$$ (neo-hookean-stress-stable)

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 {eq}`eq-green-lagrange-strain`.
For example, if $u_{i,j} \sim 10^{-8}$, then naive computation of $\bm{I} - \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 {eq}`neo-hookean-stress` using {eq}`S-to-tau`

$$
  \bm{\tau} = \bm{F}\bm{S}\bm{F}^T = \firstlame J V' \bm{I} + \secondlame \left( \bm{b} - \bm{I} \right) = \frac{\firstlame}{2} \left( J^2 - 1 \right) \bm{I} + \secondlame \left( \bm{b} - \bm{I} \right),
$$ (neo-hookean-tau)

One can arrive at the same relation {eq}`neo-hookean-tau` by expressing {eq}`neo-hookean-energy` in terms of current configuration invariant $\trace \bm{e}$ and use equation {eq}`strain-energy-grad-current`.

:::{tip}
In our simulation we use the stable version of Kirchhoff stress tensor {eq}`neo-hookean-tau` as

$$
\bm{\tau} = \frac{\firstlame}{2} \mathtt{J_{-1}} \left(\mathtt{J_{-1}} + 2 \right) \bm{I} + 2 \secondlame \bm{e}.
$$ (neo-hookean-tau-stable)

where $\mathtt{J_{-1}}$ is computed by {eq}`JM1`.
:::