Hyperelasticity

Ratel solves the momentum balance equations using unstructured high-order finite/spectral element spatial discretizations by matrix-free approach. We present here the notation and mathematical formulation of matrix-free method for a general Dirichlet problem \(\bm{f}(\bm{u}) = 0\): find \(\bm{u} \in \mathcal{V} \subset H^1 \left( \Omega_0 \right)\) such that

(76)\[ \langle\bm{v}, \bm{f}(\bm{u}) \rangle= \int_{\Omega_0}{\bm{v} \cdot \bm{f}_0 (\bm{u}, \nabla \bm{u}) + \nabla \bm{v} \tcolon \bm{f}_1 (\bm{u}, \nabla \bm{u})} \, dV = 0, \quad \forall \bm{v} \in \mathcal{V}, \]

where the operators \(\bm{f}_0\) and \(\bm{f}_1\) contain all possible sources in the problem. In order to solve (76) with Newton-Krylov iterative solvers we need its Jacobian form: find \(\diff \bm{u} \in \mathcal{V} \subset H^1 \left( \Omega_0 \right)\) such that

(77)\[ \langle\bm{v}, \bm{J}(\bm{u}) \diff \bm{u} \rangle= \int_{\Omega_0}{\bm{v} \cdot \diff \bm{f}_0 + \nabla \bm{v} \tcolon \diff \bm{f}_1} \, dV, \]

where the linearization of operators \(\bm{f}_0\) and \(\bm{f}_1\) are

\[ \diff \bm{f}_i = \frac{\partial \bm{f}_i}{\partial \bm{u}} \diff \bm{u} + \frac{\partial \bm{f}_i}{\partial \nabla \bm{u}} \nabla \diff \bm{u}, \quad i=1,2 \]

It should be noted that the gradient in the (76), (77) depends on the configuration system and it could be with respect to initial configuration \(\bm{X}\) i.e., \(\nabla_X \bm{v}\) (Lagrangian approach) or current configuration \(\bm{x}\) i.e., \(\nabla_x \bm{v}\) (Eulerian approach). In the following we specify these operators and their linearization for each problem.

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

(78)\[ \begin{CD} {\overbrace{\bm{S} \left( \bm{E} \right)}^{\text{Finite Strain Hyperelastic}}} @>{\text{constitutive}}>{\text{linearization}}> {\overbrace{\bm{S} = \mathsf{C} \bm{E}}^{\text{St. Venant-Kirchoff}}} \\ @V{\text{geometric}}V{\begin{smallmatrix}\bm{E} \to \bm{\epsilon} \\ \bm{S} \to \bm{\sigma} \end{smallmatrix}} V @V{\begin{smallmatrix}\bm{E} \to \bm{\epsilon} \\ \bm{S} \to \bm{\sigma} \end{smallmatrix}}V{\text{geometric}} V \\ {\underbrace{\bm{\sigma} \left( \bm{\epsilon} \right)}_\text{Small Strain Hyperelastic}} @>{\text{constitutive}}>\text{linearization}> {\underbrace{\bm{\sigma} = \mathsf{C} \bm{\epsilon}}_\text{Linear Elastic}} \end{CD} \]

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 reference 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 reference frame and by small letters the current one.

Weak form

We multiply balance of linear momentum explained in Continuum Mechanics for the static case, by a test function \(\bm{v}\) and integrate by parts to obtain the weak form for finite-strain hyperelasticity: find \(\bm{u} \in \mathcal{V} \subset H^1 \left( \Omega_0 \right)\) such that

(79)\[ \int_{\Omega_0}{\nabla_X \bm{v} \tcolon \bm{P}} \, dV - \int_{\Omega_0}{\bm{v} \cdot \rho_0 \bm{g}} \, dV - \int_{\partial \Omega_0}{\bm{v} \cdot \left( \bm{P} \cdot \hat{\bm{N}} \right)} \, dS = 0, \quad \forall \bm{v} \in \mathcal{V}, \]

(cf. (76) without traction, \(\bm{f}_0 = −\rho_0 \bm{g} \) and \(\bm{f}_1 = \bm{P}\) ) where \(\bm{P} \cdot \hat{\bm{N}}|_{\partial\Omega}\) is replaced by any prescribed force/traction boundary condition written in terms of 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. Discretization of (79) 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.

Newton linearization

To derive a Newton linearization, we need the Jacobian form of (79): find \(\diff \bm{u} \in \mathcal{V}\) such that

\[ \int_{\Omega_0} \nabla_X \bm{v} \tcolon \diff \bm{P} dV = \text{rhs}, \quad \forall \bm{v} \in \mathcal{V}, \]

(cf. (77), \(\diff \bm{f}_0 = 0 \) and \(\diff \bm{f}_1 = \diff \bm{P}\) ) where

(80)\[ \diff \bm{P} = \frac{\partial \bm{P}}{\partial \bm{F}} \tcolon \diff \bm{F} = \diff \bm{F}\, \bm{S} + \bm{F} \underbrace{\frac{\partial \bm{S}}{\partial \bm{E}} \tcolon \diff \bm{E}}_{\diff \bm{S}} \]

with \(\diff \bm{F} = \nabla_X\diff \bm{u}\) and

\[ \diff \bm{E} = \frac{\partial \bm{E}}{\partial \bm{F}} \tcolon \diff \bm{F} = \frac{1}{2} \left(\diff \bm{F}^T \bm{F} + \bm{F}^T \diff \bm{F} \right). \]

The linearization of the second Piola-Kirchhoff stress tensor, \(\diff \bm{S}\), depends upon the material model.

Deriving \(\diff\bm{S}\) for Neo-Hookean material

For the Neo-Hookean model (35),

(81)\[ \diff \bm{S} = \frac{\partial \bm{S}}{\partial \bm{E}} \tcolon \diff \bm{E} = \lambda J^2 \left(\bm{C}^{-1} \tcolon \diff \bm{E} \right) \bm{C}^{-1} + 2 \left( \mu - \frac{\lambda}{2} \left(J^2 - 1 \right) \right) \bm{C}^{-1} \diff \bm{E} \, \bm{C}^{-1}, \]

where we have used

\[ \diff \bm{C}^{-1} = \frac{\partial \bm{C}^{-1}}{\partial \bm{E}} \tcolon \diff \bm{E} = -2 \bm{C}^{-1} \diff \bm{E} \, \bm{C}^{-1} . \]

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 (81) in index notation,

(82)\[ \begin{aligned} \diff \bm{S}_{IJ} &= \frac{\partial \bm{S}_{IJ}}{\partial \bm{E}_{KL}} \diff \bm{E}_{KL} \\ &= \lambda J^2 \left( \bm{C}^{-1}_{KL} \diff \bm{E}_{KL} \right) \bm{C}^{-1}_{IJ} + 2 \left( \mu - \frac{\lambda}{2} \left(J^2 - 1 \right) \right) \bm{C}^{-1}_{IK} \diff \bm{E}_{KL} \bm{C}^{-1}_{LJ} \\ &= \underbrace{\left( \lambda J^2 \bm{C}^{-1}_{IJ} \bm{C}^{-1}_{KL} + 2 \left( \mu - \frac{\lambda}{2} \left(J^2 - 1 \right) \right) \bm{C}^{-1}_{IK} \bm{C}^{-1}_{JL} \right)}_{\mathsf C_{IJKL}} \diff \bm{E}_{KL} , \end{aligned} \]

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

1. recover \(\bm C^{-1}\) and \(J-1\) (either stored at quadrature points or recomputed),

2. proceed with \(3\times 3\) matrix products as in (81) or the second line of (82) to compute \(\diff \bm S\) while avoiding computation or storage of higher order tensors, and

3. conclude by (80), where \(\bm S\) is either stored or recomputed from its definition exactly as in the nonlinear residual evaluation.

Deriving \(\diff\bm{S}\) for Mooney-Rivlin material

Similar to the linearization of the constitutive relationship for Neo-Hookean materials (81), we differentiate (53) using variational notation,

(83)\[ \begin{aligned} \diff \bm{S} &= \lambda J^2 \left(\bm{C}^{-1} \tcolon \diff \bm{E} \right) \bm{C}^{-1} \\ &\quad + 2 \left( \mu_1 + 2\mu_2 - \frac{\lambda}{2} \left(J^2 - 1 \right) \right) \bm{C}^{-1} \diff \bm{E} \bm{C}^{-1} \\ &\quad + 2 \mu_2 \left( \trace \left( \diff \bm{E} \right) \bm{I}_3 - \diff \bm{E} \right) . \end{aligned} \]

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.

Cancellation vs symmetry

Some cancellation is possible (at the expense of symmetry) if we substitute (81) into (80),

(84)\[ \begin{aligned} \diff \bm P &= \diff \bm F\, \bm S + \lambda J^2 \left(\bm C^{-1} : \diff \bm E \right) \bm F^{-T} + 2\left(\mu - \frac{\lambda}{2} \left(J^2 - 1 \right) \right) \bm F^{-T} \diff\bm E \, \bm C^{-1} \\ &= \diff \bm F\, \bm S + \lambda J^2 \left(\bm F^{-T} : \diff \bm F \right) \bm F^{-T} + \left(\mu - \frac{\lambda}{2} \left(J^2 - 1 \right) \right) \bm F^{-T} \left( \bm F^T \diff \bm F + \diff \bm F^T \bm F \right) \bm C^{-1} \\ &= \diff \bm F\, \bm S + \lambda J^2 \left( \bm F^{-T} : \diff \bm F \right) \bm F^{-T} + \left(\mu - \frac{\lambda}{2} \left(J^2 - 1 \right) \right) \Big( \diff \bm F\, \bm C^{-1} + \bm F^{-T} \diff \bm F^T \bm F^{-T} \Big), \end{aligned} \]

where we have exploited \(\bm F \bm C^{-1} = \bm F^{-T}\) and

\[ \begin{aligned} \bm C^{-1} \tcolon \diff \bm E = \bm C_{IJ}^{-1} \diff \bm E_{IJ} &= \frac 1 2 \bm F_{Ik}^{-1} \bm F_{Jk}^{-1} (\bm F_{\ell I} \diff \bm F_{\ell J} + \diff \bm F_{\ell I} \bm F_{\ell J}) \\ &= \frac 1 2 \Big( \delta_{\ell k} \bm F_{Jk}^{-1} \diff \bm F_{\ell J} + \delta_{\ell k} \bm F_{Ik}^{-1} \diff \bm F_{\ell I} \Big) \\ &= \bm F_{Ik}^{-1} \diff \bm F_{kI} = \bm F^{-T} \tcolon \diff \bm F. \end{aligned} \]

We prefer to compute with (81) because (84) 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.

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.

1. Body forces (such as gravity), traction, and contact are more easily defined in the current configuration.

2. Mesh quality in the initial configuration can be very bad for large deformation.

3. 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 (79) to the Jacobian weak form in current configuration (89). One may push forward to the current configuration and then linearize or linearize in initial configuration and then push forward, as summarized below.

(85)\[ \begin{CD} {\overbrace{\nabla_X \bm{v} \tcolon \bm{FS}}^{\text{Initial Residual}}} @>{\text{push forward}}>{}> {\overbrace{\nabla_x \bm{v} \tcolon \bm{\tau}}^{\text{Current Residual}}} \\ @V{\text{linearize}}V{\begin{smallmatrix} \diff \bm{F} = \nabla_X\diff \bm{u} \\ \diff \bm{S} \left( \diff \bm{E} \right) \end{smallmatrix}}V @V{\begin{smallmatrix} \diff\nabla_x\bm{v} = -\nabla_x\bm{v} \nabla_x \diff \bm{u} \\ \diff \bm{\tau} \left( \diff \bm{\epsilon} \right) \end{smallmatrix}}V{\text{linearize}}V \\ {\underbrace{\nabla_X\bm{v}\tcolon \left( \diff \bm{F}\bm{S} + \bm{F}\diff \bm{S} \right)}_\text{Initial Jacobian}} @>{\text{push forward}}>{}> {\underbrace{\nabla_x\bm{v}\tcolon \left( \diff \bm{\tau} -\bm{\tau} \left( \nabla_x \diff \bm{u} \right)^T \right)}_\text{Current Jacobian}} \end{CD} \]

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 (79) 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

\[ \nabla_X \bm{v} \tcolon \bm{P} = \nabla_X \bm{v} \tcolon \bm{\tau} \bm{F}^{-T} = \nabla_X \bm{v} \bm{F}^{-1} \tcolon \bm{\tau} = \nabla_x \bm{v} \tcolon \bm{\tau} \]

therefore, the weak form in terms of \(\bm{\tau}\) and \(\nabla_x\) with integral over \(\Omega_0\) is

(86)\[ \int_{\Omega_0}{\nabla_x \bm{v} \tcolon \bm{\tau}} \, dV - \int_{\Omega_0}{\bm{v} \cdot \rho_0 \bm{g}} \, dV - \int_{\partial \Omega_0}{\bm{v} \cdot \left( \bm{P}\cdot\hat{\bm{N}} \right)} \, dS = 0, \quad \forall \bm{v} \in \mathcal{V}. \]

(cf. (76) without traction, \(\bm{f}_0 = −\rho_0 \bm{g} \) and \(\bm{f}_1 = \bm{\tau}\) )

Linearize in current configuration

To derive a Newton linearization of (86), first we define

(87)\[ \nabla_x \diff \bm{u} = \nabla_X \diff \bm{u} \ \bm{F}^{-1} = \diff \bm{F} \bm{F}^{-1} \]

Then by expanding the directional derivative of \(\nabla_x \bm{v} \tcolon \bm{\tau}\), we arrive at

(88)\[ \diff \ \left(\nabla_x \bm{v} \tcolon \bm{\tau} \right) = \diff \ \left( \nabla_x \bm{v} \right) \tcolon \bm{\tau} + \nabla_x \bm{v} \tcolon \diff \bm{\tau} . \]

The first term of (88) can be written as

\[ \begin{aligned} \diff \ (\nabla_x \bm{v})\tcolon \bm{\tau} &= \diff \ \left( \nabla_X \bm{v} \bm{F}^{-1} \right) \tcolon \bm{\tau} = \left( \underbrace{\nabla_X \left(\diff \bm{v}\right)}_{0}\bm{F}^{-1} + \nabla_X \bm{v}\diff \bm{F}^{-1} \right) \tcolon \bm{\tau}\\ &= \left( -\nabla_X \bm{v} \bm{F}^{-1}\diff \bm{F}\bm{F}^{-1} \right) \tcolon \bm{\tau} = \left( -\nabla_x \bm{v} \diff \bm{F}\bm{F}^{-1} \right) \tcolon \bm{\tau}\\ &= \left( -\nabla_x \bm{v} \nabla_x \diff \bm{u} \right) \tcolon \bm{\tau} = -\nabla_x \bm{v}\tcolon\bm{\tau} \left( \nabla_x \diff \bm{u} \right)^T \,, \end{aligned} \]

where we have used \(\diff \bm{F}^{-1}=-\bm{F}^{-1} \diff \bm{F} \bm{F}^{-1}\) and (87). Using this and (88) in (86) yields the Jacobian form in the current configuration: find \(\diff \bm{u} \in \mathcal{V}\) such that

(89)\[ \int_{\Omega_0} \nabla_x \bm{v} \tcolon \left( \diff \bm{\tau} - \bm{\tau} \left( \nabla_x \diff \bm{u} \right)^T \right) = \text{rhs}. \]

(cf. (77), \(\diff \bm{f}_0 = 0 \) and \(\diff \bm{f}_1 = \diff \bm{\tau} - \bm{\tau} \left( \nabla_x \diff \bm{u} \right)^T \) )

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

\[ \begin{aligned} \diff \bm{b} &= \diff \bm{F} \bm{F}^T + \bm{F} \diff \bm{F}^T \\ &= \nabla_x \diff \bm{u} \ \bm{b} + \bm{b} \ (\nabla_x \diff \bm{u})^T \\ &= (\nabla_x \diff\bm u)(\bm b - \bm I_3) + (\bm b - \bm I_3) (\nabla_x \diff\bm u)^T + 2 \diff\bm\epsilon \end{aligned} \]

where

(90)\[ \diff\bm\epsilon \equiv \frac{1}{2}\Big(\nabla_x \diff\bm{u} + (\nabla_x \diff\bm{u})^T \Big) . \]

and

\[ \begin{aligned} \diff \, (J^2 -1) &= 2 J \frac{\partial J}{\partial b} \tcolon \diff \bm{b} = J^2 \bm{b}^{-1}\tcolon \diff \bm{b}\\ &= J^2 \bm b^{-1} \tcolon \Big(\nabla_x \diff\bm u \ \bm b + \bm b (\nabla_x \diff\bm u)^T \Big) \\ &= 2 J^2 \trace (\nabla_x \diff\bm u) \\ &= 2 J^2 \trace \diff\bm\epsilon . \end{aligned} \]

Substituting into (41) gives

(91)\[ \begin{aligned} \diff \bm{\tau} &= \mu \diff \bm{b} + \lambda J^2 \trace (\diff\bm\epsilon) \bm I_3 \\ &= \underbrace{2 \mu \diff\bm\epsilon + \lambda J^2 \trace \left( \diff\bm\epsilon \right) \bm I_3 - \lambda \left( J^2 -1 \right) \diff\bm\epsilon}_{\bm F \diff\bm S \bm F^T} \\ &\quad + \left( \nabla_x \diff\bm u \right) \underbrace{\Big( \mu \left( \bm b - \bm I_3 \right) + \frac{\lambda}{2} \left( J^2 -1 \right) \bm I_3 \Big)}_{\bm\tau} \\ &\quad + \underbrace{\Big( \mu \left( \bm b - \bm I_3 \right) + \frac{\lambda}{2} \left( J^2 -1 \right) \bm I_3 \Big)}_{\bm\tau} \left( \nabla_x \diff\bm u \right)^T , \end{aligned} \]

where the final expression has been identified according to

(92)\[ \diff\bm\tau = \diff\ \left( \bm F \bm S \bm F^T \right) = \left( \nabla_x \diff\bm u \right) \bm\tau + \bm F \diff\bm S \bm F^T + \bm\tau\left( \nabla_x \diff\bm u \right)^T. \]

Collecting terms, we may thus opt to use either of the two forms

(93)\[ \begin{aligned} \diff \bm{\tau} -\bm{\tau}(\nabla_x \diff\bm{u})^T &= (\nabla_x \diff\bm u)\bm\tau + \bm F \diff\bm S \bm F^T \\ &= (\nabla_x \diff\bm u)\bm\tau + \lambda J^2 \trace(\diff\bm\epsilon) \bm I_3 + 2 \left( \mu - \frac{\lambda}{2} \left( J^2 -1\right) \right) \diff\bm\epsilon, \end{aligned} \]

with the last line showing the especially compact representation available for Neo-Hookean materials.

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

\[ \begin{aligned} \diff \mathbb{I_1} &= \frac{\partial \mathbb{I_1}}{\partial \bm{b}}\tcolon \diff \bm{b} = \bm{I}_3 \tcolon \diff \bm{b} = \trace (\diff \bm{b})\\ \end{aligned} \]

Substituting into (55) gives

(94)\[ \begin{aligned} \diff \bm{\tau} &= \mu_1 \diff \bm{b} + \lambda J^2 \trace (\diff\bm\epsilon) \bm I_3 + \mu_2 \left(\trace (\diff \bm{b}) \bm{b} + \mathbb{I_1}\diff \bm{b} - \bm{b}\diff \bm{b} - \diff \bm{b} \bm{b}\right)\\ &=2 \mu_1 \diff\bm\epsilon + \lambda J^2 \trace (\diff\bm\epsilon) \bm I_3 - \lambda \left( J^2 -1 \right) \diff\bm\epsilon + \left(\nabla_x \diff\bm u \right) \Big( \mu_1 (\bm b - \bm I_3) +\frac{\lambda}{2} \left( J^2 -1 \right) \bm I_3 \Big)\\ &\quad + \Big( \mu_1 (\bm b - \bm I_3) + \frac{\lambda}{2} \left( J^2 -1 \right) \bm I_3 \Big) (\nabla_x \diff\bm u)^T + \mu_2 \trace (\diff \bm{b}) \bm{b} \\ &\quad + \mu_2\mathbb{I_1} \left(\nabla_x \diff\bm u \right) \bm{b} + \mu_2\mathbb{I_1} \bm{b} \left(\nabla_x \diff\bm u \right)^T - \mu_2 \left(\nabla_x \diff\bm u \right) \bm{b}^2 - \mu_2 \bm{b} \left(\nabla_x \diff\bm u \right)^T \bm{b} \\ &\quad - \mu_2 \bm{b} \left(\nabla_x \diff\bm u \right) \bm{b} - \mu_2 \bm{b}^2 \left(\nabla_x \diff\bm u \right)^T + 4\mu_2 \diff\bm\epsilon - 4\mu_2 \diff\bm\epsilon \\ &= \underbrace{\lambda J^2 \trace(\diff\bm \epsilon) \bm I_3 + 2 \left( \mu_1 + 2\mu_2 - \frac{\lambda}{2} \left( J^2 -1 \right) \right) \diff \bm{\epsilon} + \mu_2 \trace \left( \diff \bm{b} \right) \bm{b} - 2\mu_2 \bm{b} \diff \bm{\epsilon} \bm{b}}_{\bm F \diff\bm S \bm F^T} \\ &\quad + (\nabla_x \diff\bm u)\underbrace{\Big( \frac{\lambda}{2} \left( J^2 -1 \right) \bm{I}_{3} + \mu_1 \left( \bm{b} - \bm{I}_3 \right) + \mu_2 \left( \mathbb{I_1} \bm{b} - 2 \bm{I}_{3} - \bm{b}^2 \right) \Big)}_{\bm\tau} \\ &\quad + \underbrace{\Big( \frac{\lambda}{2} \left( J^2 -1 \right) \bm{I}_{3} + \mu_1 \left( \bm{b} - \bm{I}_3 \right) + \mu_2 \left( \mathbb{I_1} \bm{b} - 2 \bm{I}_{3} - \bm{b}^2 \right) \Big)}_{\bm\tau} (\nabla_x \diff\bm u)^T , \end{aligned} \]

where the final expression has been identified according to

\[ \diff\bm\tau = \diff\ (\bm F \bm S \bm F^T) = (\nabla_x \diff\bm u) \bm\tau + \bm F \diff\bm S \bm F^T + \bm\tau(\nabla_x \diff\bm u)^T. \]

Collecting terms, we may thus opt to use either of the two forms

(95)\[ \begin{aligned} \diff \bm{\tau} -\bm{\tau}(\nabla_x \diff\bm{u})^T &= (\nabla_x \diff\bm u)\bm\tau + \bm F \diff\bm S \bm F^T \\ &= (\nabla_x \diff\bm u)\bm\tau + \lambda J^2 \trace(\diff\bm \epsilon) \bm I_3 + 2 \left( \mu_1 + 2\mu_2 - \frac{\lambda}{2} \left( J^2 -1 \right) \right) \diff \bm{\epsilon} \\ &\quad + \mu_2 \trace \left( \diff \bm{b} \right) \bm{b} - 2\mu_2 \bm{b} \diff \bm{\epsilon} \bm{b}.\\ \end{aligned} \]

with the last line showing the especially compact representation available for Mooney-Rivlin materials.

Linearize, then push forward

We can move the derivatives to the current configuration via

\[ \nabla_X \bm{v} \tcolon \diff \bm{P} = \left( \nabla_X \bm{v} \right) \bm{F}^{-1} \tcolon \diff \bm{P} \bm{F}^T = \nabla_x \bm{v} \tcolon \diff \bm{P} \bm{F}^T \]

and expand

(96)\[ \begin{aligned} \diff \bm{P} \bm{F}^T &= \diff \bm{F} \bm{S} \bm{F}^T + \bm{F} \diff \bm{S} \bm{F}^T \\ &= \underbrace{\diff \bm{F} \bm{F}^{-1}}_{\nabla_x \diff \bm{u}} \underbrace{\bm{F} \bm{S} \bm{F}^T}_{\bm{\tau}} + \bm{F} \diff \bm{S} \bm{F}^T . \end{aligned} \]

Representation of \(\bm F \diff\bm S \bm F^T\) for Neo-Hookean materials

Now we push (81) forward via

\[ \begin{aligned} \bm F \diff\bm S \bm F^T &= \lambda J^2 \left( \bm C^{-1} \tcolon \diff\bm E \right) \bm F \bm C^{-1} \bm F^T + 2 \left( \mu - \frac{\lambda}{2} \left( J^2 -1 \right) \right) \bm F \bm C^{-1} \diff\bm E \, \bm C^{-1} \bm F^T \\ &= \lambda J^2 \left( \bm C^{-1} \tcolon \diff\bm E \right) \bm I_3 + 2 \left( \mu - \frac{\lambda}{2} \left( J^2 -1 \right) \right) \bm F^{-T} \diff\bm E \, \bm F^{-1} \\ &= \lambda J^2 \trace(\nabla_x \diff\bm u) \bm I_3 + 2 \left( \mu - \frac{\lambda}{2} \left( J^2 -1 \right) \right) \diff\bm \epsilon \end{aligned} \]

where we have used

(97)\[ \begin{aligned} \bm C^{-1} \tcolon \diff\bm E &= \bm F^{-1} \bm F^{-T} \tcolon \bm F^T \diff\bm F \\ &= \trace(\bm F^{-1} \bm F^{-T} \bm F^T \diff \bm F) \\ &= \trace(\bm F^{-1} \diff\bm F) \\ &= \trace(\diff \bm F \bm F^{-1}) \\ &= \trace(\nabla_x \diff\bm u) \end{aligned} \]

and

(98)\[ \begin{aligned} \bm F^{-T} \diff\bm E \, \bm F^{-1} &= \frac 1 2 \bm F^{-T} \left(\bm F^T \diff\bm F + \diff\bm F^T \bm F\right) \bm F^{-1} \\ &= \frac 1 2 \left(\diff \bm F \bm F^{-1} + \bm F^{-T} \diff\bm F^T\right) \\ &= \frac 1 2 \Big(\nabla_x \diff\bm u + (\nabla_x\diff\bm u)^T \Big) \equiv \diff\bm\epsilon. \end{aligned} \]

Collecting terms, the weak form of the Newton linearization for Neo-Hookean materials in the current configuration is

(99)\[ \int_{\Omega_0} \nabla_x \bm v \tcolon \left( (\nabla_x \diff\bm u) \bm\tau + \lambda J^2 \trace(\diff\bm\epsilon)\bm I_3 + 2 \left(\mu - \frac{\lambda}{2} \left( J^2 -1 \right)\right)\diff \bm\epsilon \right) dV = \text{rhs}, \]

which equivalent to Algorithm 2 of [DPA+20] and requires only derivatives with respect to the current configuration. Note that (93) and (99) have recovered the same representation using different algebraic manipulations.

Representation of \(\bm F \diff\bm S \bm F^T\) for Mooney-Rivlin materials

We push forward \(\diff\bm S\) (83) which yields to

\[ \begin{aligned} \bm F \diff\bm S \bm F^T &= \lambda J^2 (\bm C^{-1} \tcolon \diff\bm E) \bm F \bm C^{-1} \bm F^T + 2 \left( \mu_1 + 2\mu_2 - \frac{\lambda}{2} \left( J^2 -1 \right) \right) \bm F \bm C^{-1} \diff\bm E \, \bm C^{-1} \bm F^T \\ &\quad + 2 \mu_2 \left( \trace \left( \diff \bm{E} \right) \bm{b} - \bm{F} \diff \bm{E} \bm{F}^T \right) \\ &= \lambda J^2 \trace(\diff\bm \epsilon) \bm I_3 + 2 \left( \mu_1 + 2\mu_2 - \frac{\lambda}{2} \left( J^2 -1 \right) \right) \diff \bm{\epsilon} \\ &\quad + 2 \mu_2 \left( \frac{1}{2} \trace \left( \diff \bm{b} \right) \bm{I}_3 - \bm{b} \diff \bm{\epsilon} \right)\bm{b} \end{aligned} \]

where we have used

(100)\[ \begin{aligned} \trace\left(\diff\bm E \right)&= \trace\left(\bm{F}^T \diff\bm{F} \right) = \trace\left(\bm{F} \diff\bm{F}^T \right) \\ &= \trace\left(\diff\bm F \bm{F}^T \right) = \trace\left(\diff \bm e \right) \\ \end{aligned} \]

and

(101)\[ \begin{aligned} \bm{F} \diff\bm{E} \, \bm F^{T} &= \frac{1}{2} \bm{F} \left(\bm{F}^T \diff\bm{F} + \diff\bm{F}^T \bm{F} \right) \bm{F}^{T} \\ &= \frac{1}{2} \left(\bm{F}\bm{F}^T \diff\bm{F}\bm{F}^{-1}\bm{F}\bm{F}^{T} + \bm{F}\bm{F}^{T}\bm{F}^{-T}\diff\bm{F}^T \bm{F}\bm{F}^T \right) \\ &= \frac 1 2 \bm{b}\Big(\nabla_x \diff\bm u + (\nabla_x\diff\bm u)^T \Big)\bm{b} =\bm{b} \, \diff\bm\epsilon \, \bm{b}. \end{aligned} \]

Collecting terms, we arrive at

(102)\[ \begin{aligned} \diff \bm{P} \bm{F}^T &= \left(\nabla_x \diff \bm{u} \right) \bm{\tau} + \bm{F} \diff \bm{S} \bm{F}^T \\ &= \left( \nabla_x \diff \bm{u} \right) \bm{\tau} + \lambda J^2 \trace(\diff\bm \epsilon) \bm I_3 + 2 \left( \mu_1 + 2\mu_2 - \frac{\lambda}{2} \left( J^2 -1 \right) \right) \diff \bm{\epsilon} \\ &\quad + \mu_2 \trace \left( \diff \bm{b} \right) \bm{b} - 2 \mu_2 \bm{b} \diff \bm{\epsilon} \bm{b} \end{aligned} \]

Note that (95) and (102) have recovered the same representation using different algebraic manipulations.

Pressure boundary condition

One of the important load case is the pressure boundary loading which is caused by liquids or gases on the surface of the solid structure. The pressure boundary load depends upon the current state of the deformation and can be considered as a traction vector \(\bar{\bm{t}} = \bm{\sigma} \cdot \hat{\bm{n}} = -p\hat{\bm{n}}\) per unit current surface, acting in the direction of the outward unit normal \(\hat{\bm{n}}\). Therefore, the following surface integral will be add to the weak form (79), or (86)

(103)\[ \int_{\partial \Omega}{\bm{v} \cdot \left( -p \hat{\bm{n}} \right)} \, ds, \quad \forall \bm{v} \in \mathcal{V}. \]

where the normal on current surface is computed by

(104)\[ \hat{\bm{n}} = \frac{ \frac{\partial \bm{x}}{\partial \xi_1} \times \frac{\partial \bm{x}}{\partial \xi_2} }{\lvert \frac{\partial \bm{x}}{\partial \xi_1} \times \frac{\partial \bm{x}}{\partial \xi_2} \rvert} \]

in which \( \xi_1, \xi_2 \in [-1,1]^2 \) are reference coordinate system on the face. If we write the surface area in terms of reference coordinate as \(ds = \lvert \frac{\partial \bm{x}}{\partial \xi_1} \times \frac{\partial \bm{x}}{\partial \xi_2} \rvert d\xi_1 d\xi_2\), the pressure load (103) can be simplified to

(105)\[ -\int_{[-1,1]^2} \bm{v} \cdot p \left(\frac{\partial \bm{x}}{\partial \xi_1} \times \frac{\partial \bm{x}}{\partial \xi_2} \right) d\xi_1 d\xi_2 \]

To achieve second order convergence in our Newton solver, we will need the linearization of (105) which for the constant pressure is given by

(106)\[ -\int_{[-1,1]^2} \bm{v} \cdot p \left( \frac{\partial \diff \bm{u}}{\partial \xi_1} \times \frac{\partial \bm{x}}{\partial \xi_2} + \frac{\partial \bm{x}}{\partial \xi_1} \times \frac{\partial \diff \bm{u}}{\partial \xi_2}\right) d\xi_1 d\xi_2 . \]