Hyperelastic materials¶
Isotropic models¶
Neo-Hookean¶
Constitutive theory¶
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 \(\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\)
where \(J = \lvert \bm{F} \rvert = \sqrt{\lvert \bm{C} \rvert}\) is the determinant of deformation (i.e., volume change (9)) and \(\firstlame\) and \(\secondlame\) 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
One of the commonly used form for \(\psi_{\text{vol}}\) is
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 (83) 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.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 (69) for the model (83), we make use of (70), (146) and (147) to derive
which is the second Piola-Kirchoff tensor for neo-Hookean materials.
Note
One can linearize (84) around \(\bm{E} = 0\), for which \(\bm{C} = \bm{I} + 2 \bm{E} \to \bm{I}\) and \(J \to 1 + \trace \bm{E}\), therefore (84) reduces to
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
where in the case of nearly incompressible materials, (86) 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
and if \(\log J\) appears in a stress we use log1p as
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 (84) for neo-Hookean materials as
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 (17). 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 (84) using (25)
One can arrive at the same relation (90) by expressing (83) in terms of current configuration invariant \(\trace \bm{e}\) and use equation (71).
Strong and weak formulations: initial configuration¶
In the total Lagrangian approach for the neo-Hookean hyperelasticity problem, the discrete equations are formulated with respect to the initial configuration. In this formulation, we solve for displacement \(\bm{u} \left( \bm{X} \right)\) in the initial frame \(\bm{X}\). The notation for elasticity at finite strain is inspired by [Hol00] to distinguish between the current and initial configurations. As explained in the Continuum Mechanics section, we denote by capital letters the initial frame and by small letters the current one.
We can state the hyperelastic strong form as follow: Given body force \(\rho_0 \bm g\), Dirichlet boundary \(\bar{\bm u}\) and applied traction \(\bar{\bm t}\), find the displacement variable \(\bm u \in \mathcal{V}\) such that:
and its weak formulation as:
where \((.)_0\) indicates the initial configuration, \(\bm P = \bm F \bm S\) is the first Piola-Kirchhoff stress tensor (see (22) and (24)), \(\bm N\) is the unit normal vector, and \(_X\) in \(\nabla_X\) indicates that the gradient is calculated with respect to the initial configuration. This equation contains material/constitutive nonlinearities in defining \(\bm{S}(\bm{E})\), as well as geometric nonlinearities through \(\bm{P} = \bm{F}\, \bm{S}\), \(\bm{E}(\bm{F})\), and the body force \(\bm{g}\), which must be pulled back from the current configuration to the initial configuration.
Note
Linear elasticity and small-strain hyperelasticity can both by obtained from the finite-strain hyperelastic formulation by linearization of geometric and constitutive nonlinearities. The effect of these linearizations is sketched in the diagram below, where \(\bm{\sigma}\) and \(\bm{\varepsilon}\) are stress and strain, respectively, in the small strain regime, while \(\bm{S}\) and \(\bm{E}\) are their finite-strain generalizations (second Piola-Kirchoff tensor and Green-Lagrange strain tensor, respectively) defined in the initial configuration, and \(\mathsf{C}\) is a linearized constitutive model.
Linearization: initial configuration¶
Discretization of (93) produces a finite-dimensional system of nonlinear algebraic equations, which we solve using Newton-Raphson methods. One attractive feature of Galerkin discretization is that we can arrive at the same linear system by discretizing the Newton linearization of the continuous form; that is, discretization and differentiation (Newton linearization) commute. In general for solving a nonlinear equation \(R(\bm u) = 0\), we use the first-order (Taylor’s) expansion
where \(\diff\) and \(\diff \bm u\) denote the linearization operator and the increment of the displacement field \(\bm u\). The reminder \(o(\diff \bm u)\) is a small error that tends to zero faster than \(\diff \bm u \to 0\). The second term \(\diff R\) is the linear change in \(R\) due to \(\diff \bm u\) at \(\bm u\). In fact, it is the directional derivative of \(R\) at a given \(\bm u\) (fixed) in the direction of the incremental displacement field.
The linearized form of (93) (assume constant traction and body force) may be stated as: Find \(\diff \bm u \in \mathcal{V}\) such that
where
and the linearization of deformation gradient and Green-Lagrange strain tensors are
The linearization of the second Piola-Kirchhoff stress tensor, \(\diff \bm{S}\), depends upon the material model and derived for different materials in the following.
Deriving \(\diff\bm{S}\) for neo-Hookean material
For the neo-Hookean model (84),
where
and for \(V(J) = \frac{1}{4} \left( J^2 - 1 - 2 \log J \right)\) given in strain energy (83), the derivatives are
The quantity \({\partial \bm{S}} / {\partial \bm{E}}\) is known as the incremental elasticity tensor, and is analogous to the linear elasticity tensor \(\mathsf C\).
\(\diff \bm{S}\) in index notation
It is sometimes useful to express (98) in index notation,
where we have identified the effective elasticity tensor \(\mathsf C = \mathsf C_{IJKL}\). It is generally not desirable to store \(\mathsf C\), but rather to use the earlier expressions so that only \(3 \times 3\) tensors (most of which are symmetric) must be manipulated. That is, given the linearization point \(\bm F\) and solution increment \(\diff \bm F = \nabla_X (\diff \bm u)\) (which we are solving for in the Newton step), we compute \(\diff \bm P\) via
recover \(\bm C^{-1}\) and \(J-1\) (either stored at quadrature points or recomputed),
proceed with \(3\times 3\) matrix products as in (98) or the second line of (99) to compute \(\diff \bm S\) while avoiding computation or storage of higher order tensors, and
conclude by (96), where \(\bm S\) is either stored or recomputed from its definition exactly as in the nonlinear residual evaluation.
Cancellation vs symmetry
Some cancellation is possible (at the expense of symmetry) if we substitute (98) into (96),
where we have exploited \(\bm F \bm C^{-1} = \bm F^{-T}\) and
We prefer to compute with (98) because (100) is more expensive, requiring access to (non-symmetric) \(\bm F^{-1}\) in addition to (symmetric) \(\bm C^{-1} = \bm F^{-1} \bm F^{-T}\), having fewer symmetries to exploit in contractions, and being less numerically stable.
Strong and weak formulations: current configuration¶
In the preceding discussion, all equations have been formulated in the initial configuration. This may feel convenient in that the computational domain is clearly independent of the solution, but there are some advantages to defining the equations in the current configuration.
Body forces (such as gravity), traction, and contact are more easily defined in the current configuration.
Mesh quality in the initial configuration can be very bad for large deformation.
The required storage and numerical representation can be smaller in the current configuration.
Most of the benefit in case 3 can be attained solely by moving the Jacobian representation to the current configuration [DPA+20], though residual evaluation may also be slightly faster in current configuration. There are multiple commuting paths from the nonlinear weak form in initial configuration (93) to the Jacobian weak form in current configuration (106). One may push forward to the current configuration and then linearize or linearize in initial configuration and then push forward, as summarized below.
We will follow both paths for consistency and because both intermediate representations may be useful for implementation.
Push forward, then linearize¶
The first term of (93) can be rewritten in terms of the symmetric Kirchhoff stress tensor \(\bm{\tau} = J \bm{\sigma} = \bm{P} \bm{F}^T = \bm{F} \bm{S} \bm{F}^T\) as
therefore, the weak form in terms of \(\bm{\tau}\) and \(\nabla_x\) with integral over \(\Omega_0\) is
Linearize in current configuration¶
To derive a Newton linearization of (103), first we define
Then by expanding the directional derivative of \(\nabla_x \bm{v} \tcolon \bm{\tau}\), we arrive at
The first term of (105) can be written as
where we have used \(\diff \bm{F}^{-1}=-\bm{F}^{-1} \diff \bm{F} \bm{F}^{-1}\) and (104). Using this and (105) in (103) yields the Jacobian form in the current configuration: find \(\diff \bm{u} \in \mathcal{V}\) such that
Deriving \(\diff\bm\tau\) for neo-Hookean material
To derive a useful expression of \(\diff\bm\tau\) for neo-Hookean materials, we will use the representations
where
and
Substituting into (90) gives
where the final expression has been identified according to
Collecting terms, we may thus opt to use either of the two forms
with the last line showing the especially compact representation available for neo-Hookean materials.
Linearize, then push forward¶
We can move the derivatives to the current configuration via
and expand
Representation of \(\bm F \diff\bm S \bm F^T\) for neo-Hookean materials
Now we push (98) forward via
where we have used
and
Collecting terms, the weak form of the Newton linearization for neo-Hookean materials in the current configuration is
which equivalent to Algorithm 2 of [DPA+20] and requires only derivatives with respect to the current configuration. Note that (110) and (115) have recovered the same representation using different algebraic manipulations.
Command-line interface¶
To enable the hyperelastic neo-Hookean model, use the model option -model elasticity-neo-hookean-current or -model elasticity-neo-hookean-initial and set the material parameters listed in Neo-Hookean model options. Any parameter without a default option is required.
Option |
Description |
|---|---|
|
Required to enable the neo-Hookean model, |
|
Young’s modulus, \(E > 0\) |
|
Poisson’s ratio, \(\nu \leq 0.5\). |
An example using the neo-Hookean model can be run via
$ ./bin/ratel-quasistatic -options_file examples/ymls/ex02-quasistatic-elasticity-neo-hookean-current-face-forces.yml
Mooney-Rivlin¶
Constitutive theory¶
While the neo-Hookean model depends on just two scalar invariants, \(\mathbb{I}_1 = \trace \bm{C} = 3 + 2 \trace \bm{E}\) and \(J\), Mooney-Rivlin models depend on the additional invariant, \(\mathbb{I}_2\) defined in (68). A coupled Mooney-Rivlin strain energy density (cf. neo-Hookean (83)) is [Hol00]
We differentiate \(\psi\) as in the neo-Hookean case (84) to yield the second Piola-Kirchoff tensor,
where we have used (70). This is a common model for vulcanized rubber, with a shear modulus (defined for the small-strain limit) of \(\secondlame_1 + \secondlame_2\) that should be significantly smaller than the first Lamé parameter \(\firstlame\).
Mooney-Rivlin strain energy comparison
We apply traction to a block and plot integrated strain energy \(\psi\) as a function of the loading parameter.
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 = pd.read_csv(source_path("doc/data/NH-strain.csv"))
nh["model"] = "Neo-Hookean"
nh["parameters"] = "E=2.8, nu=0.4"
mr = pd.read_csv(source_path("doc/data/MR-strain.csv"))
mr["model"] = "Mooney-Rivlin; neo-Hookean equivalent"
mr["parameters"] = "mu_1=1, mu_2=0, nu=.4"
mr1 = pd.read_csv(source_path("doc/data/MR-strain1.csv"))
mr1["model"] = "Mooney-Rivlin"
mr1["parameters"] = "mu_1=0.5, mu_2=0.5, nu=.4"
df = pd.concat([nh, mr, mr1])
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()
Tip
Similar to the neo-Hookean materials, the stable form for the Mooney-Rivlin model in initial configuration (117) can be written as
The Kirchhoff stress tensor \(\bm{\tau}\) for Mooney-Rivlin model is given by
Tip
The stable Kirchhoff stress tensor version of (119) is given by
where \(\mathtt{J_{-1}}\) is computed by (87)
The strong and weak formulation for current and initial configurations are described in Strong and weak formulations: current configuration and Strong and weak formulations: initial configuration, respectively. Herein, we have only changed \(\bm{S} = \frac{\partial\psi}{\partial\bm{E}}\), thus only the linearization needs to be updated from the neo-Hookean model.
Linearization: initial configuration¶
General process for Newton linearization in the initial configuration is described for the neo-Hookean model (Linearization: initial configuration).
Deriving \(\diff\bm{S}\) for Mooney-Rivlin material
Similar to the linearization of the constitutive relationship for neo-Hookean materials (98), we differentiate (117) using variational notation,
Note that this agrees with the linearization of the neo-Hookean constitutive relationship {eq}eq-neo-hookean-incremental-stress if \(\mu_1 = \mu, \mu_2 = 0\). Moving from neo-Hookean to Mooney-Rivlin modifies the second term and adds the third.
Linearization: current configuration¶
General process for Newton linearization in the current configuration is described for the neo-Hookean model (Strong and weak formulations: current configuration).
Representation of \(\bm F \diff\bm S \bm F^T\) for Mooney-Rivlin materials
We push forward \(\diff\bm S\) (121) which yields to
where we have used
and
Collecting terms, we arrive at
Note that (126) and (124) have recovered the same representation using different algebraic manipulations.
Deriving \(\diff\bm\tau\) for Mooney-Rivlin material
To derive a useful expression of \(\diff\bm\tau\) for Mooney-Rivlin materials, we will use the representations of \(\diff \bm{b}, \diff (J^2 -1)\) as derived in previous section and
Substituting into (119) gives
where the final expression has been identified according to
Collecting terms, we may thus opt to use either of the two forms
with the last line showing the especially compact representation available for Mooney-Rivlin materials.
Command-line interface¶
To enable the hyperelastic Mooney-Rivlin model, use the model option -model elasticity-mooney-rivlin-current or -model elasticity-mooney-rivlin-initial and set the material parameters listed in Mooney-Rivlin model options. Any parameter without a default option is required.
Option |
Description |
|---|---|
|
Required to enable the Mooney-Rivlin model, |
|
Young’s modulus, \(E > 0\) |
|
Poisson’s ratio, \(\nu \leq 0.5\). |
|
Mooney-Rivlin material constant, \(\mu_1 > 0\) |
|
Mooney-Rivlin material constant, \(\mu_2 > 0\) |
An example using the Mooney-Rivlin model can be run via
$ ./bin/ratel-quasistatic -options_file examples/ymls/ex02-quasistatic-elasticity-mooney-rivlin-initial-slip.yml
Hencky hyperelasticity¶
In Hencky-type hyperelasticity, the strain energy \(\psi\) is expressed as a function of the (Eulerian) logarithmic strain tensor \(\textbf{e}^e\), following the classical quadratic form of small-strain elasticity,
where \(\mathsf C\) is the fourth order elasticity tensor in (80). As shown in (14), \(\textbf{e}^e\) can be equally defined via the (elastic) left stretch tensor \(\textbf{v}^e\) or the (elastic) left Cauchy-Green strain tensor \(\bm{b}^e=(\textbf{v}^e)^2=\bm{F}^e(\bm{F}^e)^T\)
with eigenvalues \(\lambda^b_i=\lambda^2_i\) and eigenvectors \(\hat{\bm{n}}_i\). To ensure numerical stability when \(\lambda^b_i \approx 1\), we compute \(\textbf{e}^e\) using the eigenvalues of the Green-Euler strain tensor \(\lambda^e_i\) as [SGTB24]
The (symmetric) Kirchhoff stress tensor \(\bm{\tau}\) can thus be expressed as
Assuming isotropic elastic response and taking advantage of certain properties of derivatives of isotropic tensor functions, the following identity can be derived [dSNEARJ98]
Substituting (131) into (130) finally gives
which similarly retains the functional format of linear elasticity.
Note
[NEM16] discuss the physical interpretation and appealing properties of the quadratic Hencky free energy, making a case for its use in moderate strain applications. However, at very large strains (principal stretches \(\lambda_i>1\)), the loss of convexity and inconsistencies with experimental observations have led to a preference for alternative free energy formulations that involve additional elastic constants, such as the Ogden material. For this regime, the same authors introduce the exponentiated Hencky free energy that maintain convexity at the cost of two additional elastic constants.
The strong and weak formulation for current and initial configurations are described in Strong and weak formulations: current configuration and Strong and weak formulations: initial configuration, respectively. Herein, we have only changed \(\bm{S} = \frac{\partial\psi}{\partial\bm{E}}\), thus only the linearization needs to be updated from the neo-Hookean model.
Linearization¶
The consistent tangent for Hencky materials can be conveniently expressed as the double contraction of three fourth-order tensors [dSNEARJ98]
with
In (133), \(\mathsf C\) is the linear elasticity tensor in (80). Also note that \(\mathsf N\) is the derivative of tensor-valued function of symmetric tensor, for details see [C98] [BJ16],
Expression for \(\mathsf{B}\)¶
The change of \(\bm{b}\) with \(\bm{F}\) is obtained as,
Expression for \(\mathsf{N}\)¶
The change \(\diff\textbf{e}\) due to a perturbation \(\diff\bm{b}\) can be computed from (128) as
Expression for \(\mathsf{C}\)¶
Following (132), the change of \(\tau\) as a function of \(\textbf{e}^{e}\) increment takes the form
where \(\textbf{e}_{\text{dev}}\) is the deviatoric part of \(\textbf{e}^{e}\).
Command-line interface¶
To enable the hyperelastic Hencky model, use the model option -model elasticity-hencky-current or model elasticity-hencky-initial and set the material parameters listed in Hencky hyperelasticity model options. Any parameter without a default option is required.
Option |
Description |
|---|---|
|
Required to enable the Hencky hyperelasticity model, |
|
Young’s modulus, \(E > 0\) |
|
Poisson’s ratio, \(\nu \leq 0.5\). |
An example using the Hencky hyperelasticity model can be run via
$ ./bin/ratel-quasistatic -options_file examples/ymls/ex02-quasistatic-elasticity-hencky-initial-principal.yml