Hyperelastic Mixed Models¶

Isotropic models¶

Neo-Hookean isochoric-split¶

Constitutive theory¶

The neo-Hookean decoupled strain energy function is in terms of $J$ and modified invariant $\mathbb{\bar{I}}_1$

(154)¶\[ \begin{aligned} \psi \left(\bm{C} \right) &= \psi_{\text{vol}}(J) + \psi_{\text{iso}}(\bar{\bm{C}})\\ &=\bulk \, V(J) + \frac \secondlame 2 \left( \mathbb{\bar{I}_1} - 3 \right). \end{aligned} \]

where $\bulk$ is the bulk modulus of the material and $V(J)$ is a strictly convex function.

One approach to obtain mixed formulation is to consider the hydrostatic pressure $-\frac{\trace \bm \sigma}{3}$ as a variable, where in terms of energy is

(155)¶\[ \begin{aligned} p_{\text{hyd}} = - \frac{\trace \bm \sigma}{3} &= -\frac{1}{3J} \trace (\bm F \bm S \bm{F}^T) = -\frac{1}{3J} \trace (\bm S \bm F^T \bm F) = -\frac{1}{3J} \trace (\bm S \bm C) \\ &= -\frac{1}{3J} \trace \left( \frac{\partial \psi_{\text{vol}}}{\partial J} ~ J \bm{I} + \sum_{i=1}^3 \frac{\partial \psi_{\text{iso}}}{\partial \mathbb{\bar{I}}_i} \frac{\partial \mathbb{\bar{I}}_i}{\partial \bm E} \bm C \right) \\ &= - \frac{\partial \psi_{\text{vol}}}{\partial J} = -\bulk \frac{\partial V(J)}{\partial J} = -\bulk V', \end{aligned} \]

where we have used $\trace (\bm A \bm B) = \trace (\bm B \bm A)$ and

\[ \begin{aligned} \trace \left( \frac{\partial \mathbb{\bar{I}}_1}{\partial \bm E} \bm C \right) &= 2 J^{-2/3} \trace\left( \bm C -\frac{1}{3} \mathbb{I}_1 \bm{I} \right) = 0, \\ \trace \left( \frac{\partial \mathbb{\bar{I}}_2}{\partial \bm E} \bm C \right) &= 2 J^{-4/3} \trace \left( \mathbb{I}_1 \bm C - \bm C^2 - \frac{2}{3} \mathbb{I}_2 \bm{I} \right) = 0, \\ \trace \left( \frac{\partial \mathbb{\bar{I}}_3}{\partial \bm E} \bm C \right) &= 0. \end{aligned} \]

However, we define a pressure-like variable similar to mixed linear elasticity (150) in terms of bulk and primal bulk (151) moduli

(156)¶\[ p = - (\bulk - \bulk_p) \frac{\partial V}{\partial J} = p_{\text{hyd}} + \bulk_p \frac{\partial V}{\partial J}, \]

and considered it as independent field variable in our mixed hyperelastic formulation.

The second Piola-Kirchoff tensor is computed via (144)

(157)¶\[ \begin{aligned} \bm{S} = \frac{\partial \psi}{\partial \bm{E}} &= \underbrace{\frac{\partial \psi_{\text{vol}}}{\partial J}}_{-p_{\text{hyd}}} \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}_{\text{vol}} + \bm{S}_{\text{iso}}\\ &= \left(\bulk_p J \, V' - p \, J \right) \bm{C}^{-1} + \secondlame J^{-2/3}\left( \bm{I} - \frac{1}{3}\mathbb{{I}_1}\bm{C}^{-1} \right). \end{aligned} \]

where the invariants are defined in (68) and (143) and we have used (146), (145), and (147).

Tip

In our simulation we use the stable version of (157) as

(158)¶\[ \bm{S} = \left(\bulk_p J \, V' - p \, J \right) \bm{C}^{-1} + 2 \secondlame J^{-2/3} \bm{C}^{-1} \bm{E}_{\text{dev}}, \]

where $\bm{E}_{\text{dev}} = \bm{E} - \frac{1}{3}\trace \bm{E} \, \bm{I} $ is the deviatoric part of Green-Lagrange strain tensor.

The isochoric neo-Hookean stress relation can be represented in current configuration by pushing forward (157) using (25)

(159)¶\[ \bm{\tau} = \bm{F}\bm{S}\bm{F}^T = \left(\bulk_p J \, V' - p \, J \right) \bm{I} + \secondlame J^{-2/3}\left( \bm{b} - \frac{1}{3}\mathbb{{I}_1} \bm{I} \right). \]

Tip

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

(160)¶\[ \bm{\tau} = \left(\bulk_p J \, V' - p \, J \right) \bm{I} + 2 \secondlame J^{-2/3} \bm{e}_{\text{dev}} \]

where $\bm{e}_{\text{dev}} = \bm{e} - \frac{1}{3}\trace \bm{e} \, \bm{I} $ is the deviatoric part of Green-Euler strain tensor

Note

We can solve the isochoric model with single field displacement if the material is not incompressible which leads to the second Piola-Kirchoff tensor

(161)¶\[ \begin{aligned} \bm{S} = \frac{\partial \psi}{\partial \bm{E}} &= \frac{\partial \psi_{\text{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}_{\text{vol}} + \bm{S}_{\text{iso}}\\ &= \bulk J \, V' \, \bm{C}^{-1} + \secondlame J^{-2/3}\left( \bm{I} - \frac{1}{3}\mathbb{{I}_1}\bm{C}^{-1} \right). \end{aligned} \]

and the kirchoff stress

(162)¶\[ \bm{\tau} = \bm{F}\bm{S}\bm{F}^T = \bulk J \, V' \, \bm{I} + \secondlame J^{-2/3}\left( \bm{b} - \frac{1}{3}\mathbb{{I}_1} \bm{I} \right). \]

Note that the stable form of the above stresses can be derived similar to (158), and (160).

Strong and weak formulations: initial configuration¶

We can state the mixed 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 and pressure-like variables $(\bm u, p) \in \mathcal{V} \times \mathcal{Q} $, such that:

(163)¶\[ \begin{aligned} -\nabla_X \cdot \bm P - \rho_0 \bm g &= 0. \qquad \text{in $\Omega_0$} \\ -\frac{\partial V}{\partial J} - \frac{p}{\bulk - \bulk_p} &= 0. \qquad \text{in $\Omega$} \\ \bm{u} &= \bar{\bm u}, \quad \text{on $\partial \Omega^{D}_0$} \\ \bm P \cdot \bm N &= \bar{\bm t}. \qquad \text{on $\partial \Omega^{N}_0$} \end{aligned} \]

and by multiplying the strong form (163) by test functions $(\bm{v}, q)$ and integrate by parts, we can obtain the weak form for finite-strain incompressible hyperelasticity as: find $(\bm{u}, p) \in \mathcal{V} \times \mathcal{Q} \subset H^1 \left( \Omega_0 \right) \times L^2 \left( \Omega_0 \right)$ such that

(164)¶\[ \begin{aligned} R_u(\bm u, p) &\coloneqq \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 \bar{\bm t}} \, dA = 0, \quad \forall \bm{v} \in \mathcal{V}, \\ R_p(\bm u, p) &\coloneqq \int_{\Omega_0} q \cdot L J \, dV = 0, \quad \forall q\in \mathcal{Q}, \end{aligned} \]

where $L = -V' - \frac{p}{\bulk - \bulk_p} $. 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.

Linearization: initial configuration¶

As explained in Linearization: initial configuration, we need the Jacobian form of the (164) which is: find $(\diff \bm{u}, \diff p) \in \mathcal{V} \times \mathcal{Q}$ such that

(165)¶\[ \begin{aligned} \int_{\Omega_0} \nabla_X \bm{v} \tcolon \diff \bm{P} dV &= -R_u, \quad \forall \bm{v} \in \mathcal{V}, \\ \int_{\Omega_0} q \cdot \left( \diff L J + L \diff J \right)dV &= -R_p, \quad \forall q \in \mathcal{Q} \end{aligned} \]

where

(166)¶\[ \begin{aligned} \diff \bm{P} &= \diff \bm{F}\, \bm{S} + \bm{F} \diff \bm{S}, \\ \diff L J + L \diff J &= J \frac{\partial L}{\partial \bm{E}} \tcolon \diff \bm{E} + J \frac{\partial L}{\partial p} \diff p + L \frac{\partial J}{\partial \bm{E}} \tcolon \diff \bm{E} \end{aligned} \]

with $\diff \bm{F}$ and $\diff \bm E$ defined in (97). The linearization of the second equation of (166) is

(167)¶\[ \begin{aligned} \diff L J + L \diff J &= \left( -V^{''} \diff J - \frac{\diff p}{\bulk - \bulk_p} \right) J + L J (\bm{C}^{-1} \tcolon \diff \bm{E}), \\ &= \left( - J^2 \, V^{''} - J \, V' - J \frac{p}{\bulk - \bulk_p} \right) \bm{C}^{-1} \tcolon \diff \bm{E} - J \frac{\diff p}{\bulk - \bulk_p} \end{aligned} \]

The linearization of the second Piola-Kirchhoff stress tensor, $\diff \bm{S}$, depends upon the material model and derived below for different hyperelastic materials.

Note

For the fully incompressible case, we can set $\nu=0.5$ which leads to infinite bulk modulus. Therefore, all terms in residual (164) and jacobian (165) with coefficient $\frac{1}{\bulk -\bulk_p}$ goes to zero in the numerical implementation. In face we have

\[ L = -V', \\ \diff L J + L \diff J = \left( - J^2 \, V^{''} - J \, V' \right) \bm{C}^{-1} \tcolon \diff \bm{E}, \]

which leads to a saddle point system $\begin{bmatrix} \bm A & \bm B \\ \bm C & \bm 0 \end{bmatrix}$.

Note

If we use single field isochoric model (161) the linearization of the volumetric stress becomes

(172)¶\[ \diff \bm{S}_{\text{vol}} = \left( \bulk J^2 V^{''} + \bulk J \, V'\right) (\bm{C}^{-1} \tcolon \diff \bm{E}) \bm{C}^{-1} + \bulk J \, V' \diff \bm{C}^{-1}. \]

$\diff \bm{S}_{\text{vol}}^p = 0$ and we only need to solve the first equation of (164).

Strong and weak formulations: current configuration¶

Similar to what we have shown in (101), we can write the first equation of (164) in the current configuration (see (103) ) which yields to the Jacobian form

(173)¶\[ \begin{aligned} \int_{\Omega_0} \nabla_x \bm{v} \tcolon \left( \diff \bm{\tau} - \bm{\tau} \left( \nabla_x \diff \bm{u} \right)^T \right) dV &= -r_u, \quad \forall \bm{v} \in \mathcal{V}, \\ \int_{\Omega_0} q \cdot \left( \diff L J + L \diff J \right)dV &= -r_p, \quad \forall q \in \mathcal{Q} \end{aligned} \]

where the second equation is similar to (165) and we can use (109)

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

to compute the linearization of the first equation in current configuration for different material models.

Linearization: current configuration¶

Based on the split (168), we can derive

(174)¶\[ \bm{F} \diff \bm{S}_{\text{vol}}^u \bm{F}^{T} = \left( \bulk_p J^2 V^{''} + \bulk_p J \, V' - p J \right) \trace (\diff \bm \epsilon) \bm{I} - 2 \left(\bulk_p J \, V' - p \, J \right) \diff \bm \epsilon. \]

(175)¶\[ \bm{F} \diff \bm{S}_{\text{vol}}^p \bm{F}^{T} = -J \diff p \bm{I}, \]

where $\diff \bm \epsilon$ is defined in (107) and we have used (113) and (114). The isochoric part can be simplified as

(176)¶\[ \begin{aligned} \bm{F} \diff \bm{S}_{\text{iso}} \bm{F}^{T} &= -\frac{2}{3} \mu J^{-2/3} \left( 2 \trace (\diff \bm \epsilon) \bm{e}_{\text{dev}} + \trace(\diff \bm e) \bm{I} - \mathbb{{I}_1} \diff \bm \epsilon \right), \end{aligned} \]

where $\mathbb{{I}_1} = \trace(\bm{b})$.

Note

If we use single field isochoric model the linearization of the volumetric stress becomes

(177)¶\[ \bm{F} \diff \bm{S}_{\text{vol}} \bm{F}^{T} = \left( \bulk J^2 V^{''} + \bulk J \, V'\right) \trace (\diff \bm \epsilon) \bm{I} -2 \bulk J \, V' \diff \bm \epsilon. \]

$\diff \bm{S}_{\text{vol}}^p = 0$ and we only need to solve the first equation of (173).

Command-line interface¶

To enable the hyperelastic isochoric-split neo-Hookean model, use the model option -model elasticity-isochoric-neo-hookean-current or -model elasticity-isochoric-neo-hookean-initial and set the material parameters listed in Isochoric-split neo-Hookean model options. Any parameter without a default option is required.

Table 22 Isochoric-split neo-Hookean model options¶
Option	Description
`-model elasticity-isochoric-neo-hookean-current \| elasticity-isochoric-neo-hookean-initial`	Required to enable the isochoric-split neo-Hookean model, `\|` denotes an option for the user.
`-E [real]`	Young’s modulus, $E > 0$
`-nu [real]`	Poisson’s ratio, $\nu \leq 0.5$.

An example using the isochoric-split neo-Hookean model can be run via

$ ./bin/ratel-quasistatic -options_file examples/ymls/ex02-quasistatic-elasticity-isochoric-neo-hookean-current-face-forces.yml

To enable the hyperelastic mixed neo-Hookean model, use the model option -model elasticity-mixed-neo-hookean-current or -model elasticity-mixed-neo-hookean-initial and set the material parameters listed in Mixed neo-Hookean model options. Any parameter without a default option is required.

Table 23 Mixed neo-Hookean model options¶
Option	Description	Default value
`-model elasticity-mixed-neo-hookean-current \| elasticity-mixed-neo-hookean-initial`	Required to enable the mixed neo-Hookean model, `\|` denotes an option for the user.
`-E [real]`	Young’s modulus, $E > 0$
`-nu [real]`	Poisson’s ratio, $\nu \leq 0.5$.
`-nu_primal [real]`	Primal part of the Poisson’s ratio for Jacobian operator, $-1 \leq \nu_p < \nu \leq 0.5$. This parameter may improve solver performance and should not affect the solution.	`-1.0`
`-nu_primal_pc [real]`	Primal part of the Poisson’s ratio for preconditioner operator, $-1 \leq \nu_{pc} < \nu \leq 0.5$. This parameter may improve solver performance and should not affect the solution.	`-nu_primal` value

An example using the mixed neo-Hookean model can be run via

$ ./bin/ratel-quasistatic -options_file examples/ymls/ex02-quasistatic-elasticity-mixed-neo-hookean-current-pcjacobi.yml

Mooney-Rivlin isochoric-split¶

Constitutive theory¶

The Mooney-Rivlin decoupled strain energy function is in terms of $J$ and modified invariants $\mathbb{\bar{I}}_1, \mathbb{\bar{I}}_2$ (143)

(178)¶\[ \begin{aligned} \psi \left(\bm{C} \right) &= \psi_{\text{vol}}(J) + \psi_{\text{iso}}(\bar{\bm{C}})\\ &=\bulk \, V(J) + \frac{\secondlame_1}{2} \left( \mathbb{\bar{I}_1} - 3 \right) + \frac{\secondlame_2}{2} \left( \mathbb{\bar{I}_2} - 3 \right). \end{aligned} \]

where $\bulk$ is the bulk modulus of the material and $V(J)$ is a strictly convex function.

The second Piola-Kirchoff tensor is computed via (144)

(179)¶\[ \begin{aligned} \bm{S} = \frac{\partial \psi}{\partial \bm{E}} &= \frac{\partial \psi_{\text{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}} + \frac{\partial \psi}{\partial \mathbb{\bar{I}_2}} \frac{\partial \mathbb{\bar{I}_2}}{\partial \bm{E}} = \bm{S}_{\text{vol}} + \bm{S}_{\text{iso}}\\ &= \left(\bulk_p J \, V' - p \, J \right) \bm{C}^{-1} + \secondlame_1 J^{-2/3}\left( \bm{I} - \frac{1}{3}\mathbb{{I}_1}\bm{C}^{-1} \right) + \secondlame_2 J^{-4/3}\left( \mathbb{{I}_1} \bm{I} - \bm{C} - \frac{2}{3}\mathbb{{I}_2}\bm{C}^{-1} \right). \end{aligned} \]

where the invariants are defined in (68) and (143) and we have used (146), (145), and (147).

The pressure-like variable $p$ in (179) is defined similar to (156).

Tip

Similar to neo-Hookean isochoric model (158) we can derive a stable form of (179) as

(180)¶\[ \begin{aligned} \bm{S} &= \left(\bulk_p J \, V' - p \, J \right) \bm{C}^{-1} + \left( 2 \secondlame_1 J^{-2/3} + 4 \secondlame_2 J^{-4/3} \right) \bm{C}^{-1} \bm{E}_{\text{dev}} \\ &+ 2 \secondlame_2 J^{-4/3} \left( \mathbb{{I}_1}(\bm{E})\bm{I} - \bm{E} \right) - \frac{4}{3} \secondlame_2 J^{-4/3} \left( \mathbb{{I}_1}(\bm{E}) + 2 \mathbb{{I}_2}(\bm{E})\right) \bm{C}^{-1}, \end{aligned} \]

where $\bm{E}_{\text{dev}} = \bm{E} - \frac{1}{3}\trace \bm{E} \, \bm{I} $ is the deviatoric part of Green-Lagrange strain tensor and $\mathbb{{I}_1}(\bm{E}), \mathbb{{I}_2}(\bm{E})$ are first and second invariants of the tensor $\bm{E}$.

The isochoric Mooney-Rivlin stress relation can be represented in current configuration by pushing forward (179) using (25)

(181)¶\[ \bm{\tau} = \bm{F}\bm{S}\bm{F}^T = \left(\bulk_p J \, V' - p \, J \right) \bm{I} + \secondlame_1 J^{-2/3}\left( \bm{b} - \frac{1}{3}\mathbb{{I}_1} \bm{I} \right) + \secondlame_2 J^{-4/3}\left( \mathbb{{I}_1} \bm{b} - \bm{b}^2 - \frac{2}{3}\mathbb{{I}_2} \bm{I} \right) \]

Tip

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

(182)¶\[ \begin{aligned} \bm{\tau} &= \left(\bulk_p J \, V' - p \, J \right) \bm{I} + \left( 2 \secondlame_1 J^{-2/3} + 4 \secondlame_2 J^{-4/3} \right) \bm{e}_{\text{dev}} \\ &+ 2 \secondlame_2 J^{-4/3} \left( \mathbb{{I}_1}(\bm{e})\bm{I} - \bm{e} \right) \bm{b} - \frac{4}{3} \secondlame_2 J^{-4/3} \left( \mathbb{{I}_1}(\bm{e}) + 2 \mathbb{{I}_2}(\bm{e})\right) \bm{I}, \end{aligned} \]

where $\bm{e}_{\text{dev}} = \bm{e} - \frac{1}{3}\trace \bm{e} \, \bm{I} $ is the deviatoric part of Green-Euler strain tensor and $\mathbb{{I}_1}(\bm{e}), \mathbb{{I}_2}(\bm{e})$ are first and second invariants of the tensor $\bm{e}$.

Note

We can solve the isochoric model with single field displacement if the material is not incompressible which leads to the second Piola-Kirchoff tensor

(183)¶\[ \begin{aligned} \bm{S} = \frac{\partial \psi}{\partial \bm{E}} &= \frac{\partial \psi_{\text{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}} + \frac{\partial \psi}{\partial \mathbb{\bar{I}_2}} \frac{\partial \mathbb{\bar{I}_2}}{\partial \bm{E}} = \bm{S}_{\text{vol}} + \bm{S}_{\text{iso}}\\ &= \bulk J \, V' \,\bm{C}^{-1} + \secondlame_1 J^{-2/3}\left( \bm{I} - \frac{1}{3}\mathbb{{I}_1}\bm{C}^{-1} \right) + \secondlame_2 J^{-4/3}\left( \mathbb{{I}_1} \bm{I} - \bm{C} - \frac{2}{3}\mathbb{{I}_2}\bm{C}^{-1} \right). \end{aligned} \]

and the kirchoff stress

(184)¶\[ \bm{\tau} = \bulk J \, V' \, \bm{I} + \secondlame_1 J^{-2/3}\left( \bm{b} - \frac{1}{3}\mathbb{{I}_1} \bm{I} \right) + \secondlame_2 J^{-4/3}\left( \mathbb{{I}_1} \bm{b} - \bm{b}^2 - \frac{2}{3}\mathbb{{I}_2} \bm{I} \right) \]

Note that the stable form of the above stresses can be derived similar to (180), and (182).

The strong and weak formulation for initial configuration and current configuration is described in Strong and weak formulations: initial configuration and Strong and weak formulations: current 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 general process for Newton linearization in the initial configuration is described for the neo-Hookean model Linearization: initial configuration.

The general process for Newton linearization in the current configuration is described for the neo-Hookean model Linearization: current configuration.

Command-line interface¶

To enable the hyperelastic isochoric-split Mooney-Rivlin model, use the model option -model elasticity-isochoric-mooney-rivlin-current or -model elasticity-isochoric-mooney-rivlin-initial and set the material parameters listed in Isochoric-split Mooney-Rivlin model options. Any parameter without a default option is required.

Table 24 Isochoric-split Mooney-Rivlin model options¶
Option	Description
`-model elasticity-isochoric-mooney-rivlin-current \| elasticity-isochoric-mooney-rivlin-initial`	Required to enable the isochoric-split Mooney-Rivlin model, `\|` denotes an option for the user.
`-E [real]`	Young’s modulus, $E > 0$
`-nu [real]`	Poisson’s ratio, $\nu \leq 0.5$.
`-mu_1 [real]`	Mooney-Rivlin material constant, $\mu_1 > 0$
`-mu_2 [real]`	Mooney-Rivlin material constant, $\mu_2 > 0$

An example using the isochoric-split Mooney-Rivlin model can be run via

$ ./bin/ratel-quasistatic -options_file examples/ymls/ex02-quasistatic-elasticity-isochoric-mooney-rivlin-current-face-forces.yml

Ogden isochoric-split¶

Constitutive theory¶

The strain energy of neo-Hookean and Mooney-Rivlin are written in terms of invariants of $\bm{C}$. The decoupled Ogden strain energy function is in terms of modified principal stretch (142)

(190)¶\[ \begin{aligned} \psi \left(\bar{\lambda}_1, \bar{\lambda}_2, \bar{\lambda}_3, J \right) &= \psi_{\text{vol}}(J) + \psi_{\text{iso}}\left( \bar{\lambda}_1, \bar{\lambda}_2, \bar{\lambda}_3 \right)\\ &=\bulk \, V(J) + \sum_{i=1}^3 \omega(\bar{\lambda}_i). \end{aligned} \]

where

(191)¶\[ \omega(\bar{\lambda}_i) = \sum_{j=1}^N \frac{m_j}{\alpha_j} \left(\bar{\lambda}_i^{\alpha_j} -1 \right) \]

with the consistency condition

(192)¶\[ 2 \secondlame = \sum_{j=1}^N m_j \alpha_j \quad \text{with} \quad m_j \alpha_j > 0 \]

obtained from comparison with the linear theory.

Tip

For $N=1, \alpha_1 = 2, m_1 = \secondlame$ the Ogden model (190) simplifies to

\[ \psi \left( \bar{\lambda}_1, \bar{\lambda}_2, \bar{\lambda}_3, J \right) = \frac{\bulk}{4} \left( J^2 - 1 -2 \log J \right) + \frac{\secondlame}{2} \left(\bar{\lambda}_1^2 + \bar{\lambda}_2^2 + \bar{\lambda}_3^2 - 3\right) \]

which is the same as the decoupled neo-Hookean energy function (154). Also for $J=1, N=2, \alpha_1 = 2, \alpha_2 = -2, m_1 = \secondlame_1, m_2 = -\secondlame_2$ and by using $J=1=\lambda_1 \lambda_2 \lambda_3 = \bar{\lambda}_1 \bar{\lambda}_2 \bar{\lambda}_3$ constrain, we can simplify (190) as

\[ \begin{aligned} \psi \left( \bar{\lambda}_1, \bar{\lambda}_2, \bar{\lambda}_3, J=1 \right) &= \frac{m_1}{2} \left(\bar{\lambda}_1^2 + \bar{\lambda}_2^2 + \bar{\lambda}_3^2 - 3\right) - \frac{m_2}{2} \left(\bar{\lambda}_1^{-2} + \bar{\lambda}_2^{-2} + \bar{\lambda}_3^{-2} - 3\right) \\ &= \frac{\secondlame_1}{2} \left(\mathbb{\bar{I}}_1 - 3\right) + \frac{\secondlame_2}{2} \left(\mathbb{\bar{I}}_2 - 3\right) \end{aligned} \]

which is the Mooney-Rivlin isochoric model.

We differentiate $\psi$ as in the isochoric neo-Hookean case (157) to yield the second Piola-Kirchoff tensor,

(193)¶\[ \begin{aligned} \bm{S} = \frac{\partial \psi}{\partial \bm{E}} &= \frac{\partial \psi_{\text{vol}}}{\partial J} \frac{\partial J}{\partial \bm{E}} + \frac{\partial \psi_{\text{iso}}}{\partial \bm{E}} =\bm{S}_{\text{vol}} + \bm{S}_{\text{iso}}\\ &= \left(\bulk_p J \, V' - p \, J \right) \bm{C}^{-1} + \sum_{i=1}^3 S_i^{\text{iso}} \hat{\bm{N}_i} \hat{\bm{N}_i}^T. \end{aligned} \]

where the pressure-like variable $p$ is defined via (156) and we have used (72) and (75). By employing $\frac{\partial J}{\partial {\lambda}_i} = J {\lambda}_i^{-1}$ and relation (142)

(194)¶\[ \begin{aligned} S_i^{\text{iso}} = \frac{1}{\lambda_i}\frac{\partial \psi_{\text{iso}}}{\partial \lambda_i} =\frac{1}{\lambda_i}\frac{\partial \bar{\lambda}_k}{\partial \lambda_i} \frac{\partial \psi_{\text{iso}}}{\partial \bar{\lambda}_k} &=\frac{J^{-1/3}}{\lambda_i} \left(\delta_{ik} -\frac{1}{3}\lambda_i^{-1}\lambda_k \right) \frac{\partial \psi_{\text{iso}}}{\partial \bar{\lambda}_k} \\ &=\frac{J^{-1/3}}{\lambda_i} \left(\delta_{ik} -\frac{1}{3}\bar{\lambda}_i^{-1}\bar{\lambda}_k \right) \frac{\partial \psi_{\text{iso}}}{\partial \bar{\lambda}_k}. \end{aligned} \]

in which

(195)¶\[ \frac{\partial \psi_{\text{iso}}}{\partial \bar{\lambda}_k} = \frac{\partial \omega (\bar{\lambda}_k)}{ \partial \bar{\lambda}_k} = \sum_{j=1}^N m_j \bar{\lambda}_k^{\alpha_j - 1} \]

Tip

To derive an equivalent numerically stable form of the isochoric part of (193) we rewrite the $\psi_{\text{iso}}$ by substituting $\bar{\lambda}_i = J^{-1/3}\lambda_i$ as

(196)¶\[ \psi_{\text{iso}}({\lambda}_i) = \sum_{j=1}^N \frac{m_j}{\alpha_j} \left[\left({\lambda}_1^{\alpha_j} + {\lambda}_2^{\alpha_j} + {\lambda}_3^{\alpha_j} \right)J^{-\alpha_j/3} -3 \right] \]

then

(197)¶\[ \begin{aligned} S_1^{\text{iso}} = \frac{1}{\lambda_1}\frac{\partial \psi_{\text{iso}}}{\partial \lambda_1} &= \frac{1}{\lambda_1} \sum_{j=1}^N \frac{m_j}{\alpha_j} \left[ \left(\alpha_j \lambda_1^{\alpha_j -1} \right)J^{-\alpha_j/3} - \frac{\alpha_j}{3} J^{-\alpha_j/3} \lambda_1^{-1} \left( {\lambda}_1^{\alpha_j} + {\lambda}_2^{\alpha_j} + {\lambda}_3^{\alpha_j} \right)\right] \\ &= \frac{1}{\lambda_1} \sum_{j=1}^N m_j \left[ \lambda_1^{\alpha_j -1} - \frac{1}{3} \lambda_1^{-1} \left( {\lambda}_1^{\alpha_j} + {\lambda}_2^{\alpha_j} + {\lambda}_3^{\alpha_j} \right)\right] J^{-\alpha_j/3} \\ &= \frac{1}{\lambda_1^2} \sum_{j=1}^N \frac{m_j}{3} \left[ 2\lambda_1^{\alpha_j} - {\lambda}_2^{\alpha_j} - {\lambda}_3^{\alpha_j} \right] J^{-\alpha_j/3} \\ &= \frac{1}{\lambda_1^2} \sum_{j=1}^N \frac{m_j}{3} \left[ 2\left(\lambda_1^{\alpha_j} -1 \right) - \left({\lambda}_2^{\alpha_j} -1 \right) - \left({\lambda}_3^{\alpha_j} -1 \right)\right] J^{-\alpha_j/3} \\ &=\frac{1}{\lambda_1^2} \sum_{j=1}^N \frac{m_j}{3} \left[ 2(e^{\alpha_j \ell_1}-1) - (e^{\alpha_j \ell_2}-1) - (e^{\alpha_j \ell_3}-1) \right] J^{-\alpha_j/3} \\ &=\frac{1}{1 + 2\lambda_1^E} \sum_{j=1}^N \frac{m_j}{3} \left[ 2\operatorname{\tt expm1}(\alpha_j \ell_1) - \operatorname{\tt expm1}(\alpha_j \ell_2) - \operatorname{\tt expm1}(\alpha_j \ell_3) \right] J^{-\alpha_j/3} \end{aligned} \]

where the eigen pair $(\lambda_i, \hat{\bm{N}_i})$ are computed by solving the eigen problem of Green-Lagrange strain tensor $\bm{E} \hat{\bm{N}_i} = \lambda_i^E \hat{\bm{N}_i}$ with $\lambda_i = \sqrt{1 + 2 \lambda_i^E}$ and and $\ell_i = \log \lambda_i = \frac 1 2 \operatorname{\tt log1p}(2\lambda_i^E)$. Following the above approach we have

(198)¶\[ \begin{aligned} S_2^{\text{iso}} = \frac{1}{\lambda_2}\frac{\partial \psi_{\text{iso}}}{\partial \lambda_2} =\frac{1}{1 + 2\lambda_2^E} \sum_{j=1}^N \frac{m_j}{3} \left[ -\operatorname{\tt expm1}(\alpha_j \ell_1) +2 \operatorname{\tt expm1}(\alpha_j \ell_2) - \operatorname{\tt expm1}(\alpha_j \ell_3) \right] J^{-\alpha_j/3} \\ S_3^{\text{iso}} = \frac{1}{\lambda_3}\frac{\partial \psi_{\text{iso}}}{\partial \lambda_3} =\frac{1}{1 + 2\lambda_3^E} \sum_{j=1}^N \frac{m_j}{3} \left[ -\operatorname{\tt expm1}(\alpha_j \ell_1) - \operatorname{\tt expm1}(\alpha_j \ell_2) + 2 \operatorname{\tt expm1}(\alpha_j \ell_3) \right] J^{-\alpha_j/3} \end{aligned} \]

substituting the new definition of $S_i^{\text{iso}}$ (197), (198) into (193) gives the stable form of the second Piola-Kirchhoff stress for Ogden model.

The Kirchhoff stress tensor $\bm{\tau}$ for Ogden model is given by

(199)¶\[ \bm{\tau} = \bm{F}\bm{S}\bm{F}^T = \bm{\tau}_{\text{vol}} + \bm{\tau}_{\text{iso}} = \left(\bulk_p J \, V' - p \, J \right) \bm{I} + \sum_{i=1}^3 \tau_i^{\text{iso}} \hat{\bm{n}_i} \hat{\bm{n}_i}^T. \]

where we have used $\bm{F} \hat{\bm{N}_i} = \lambda_i \hat{\bm{n}_i}$ which $\hat{\bm{n}_i}$ are eigenvectors of $\bm{b}$ and $\tau_i^{\text{iso}} = \lambda_i^2 S_i^{\text{iso}}$.

Tip

Similar to the initial configuration, we can compute the $\bm{\tau}_{\text{iso}}$ in stable form by using the stable coefficients

(200)¶\[ \begin{aligned} \tau_1^{\text{iso}} &= \lambda_1 \frac{\partial \psi_{\text{iso}}}{\partial \lambda_1} = \sum_{j=1}^N \frac{m_j}{3} \left[ 2\operatorname{\tt expm1}(\alpha_j \ell_1) - \operatorname{\tt expm1}(\alpha_j \ell_2) - \operatorname{\tt expm1}(\alpha_j \ell_3) \right] J^{-\alpha_j/3} \\ \tau_2^{\text{iso}} &= \lambda_2 \frac{\partial \psi_{\text{iso}}}{\partial \lambda_2} = \sum_{j=1}^N \frac{m_j}{3} \left[ -\operatorname{\tt expm1}(\alpha_j \ell_1) +2 \operatorname{\tt expm1}(\alpha_j \ell_2) - \operatorname{\tt expm1}(\alpha_j \ell_3) \right] J^{-\alpha_j/3} \\ \tau_3^{\text{iso}} &= \lambda_3 \frac{\partial \psi_{\text{iso}}}{\partial \lambda_3} = \sum_{j=1}^N \frac{m_j}{3} \left[ -\operatorname{\tt expm1}(\alpha_j \ell_1) - \operatorname{\tt expm1}(\alpha_j \ell_2) + 2 \operatorname{\tt expm1}(\alpha_j \ell_3) \right] J^{-\alpha_j/3} \end{aligned} \]

and computing the eigen pair $(\lambda_i, \hat{\bm{n}_i})$ of Green-Euler strain tensor $\bm{e} \hat{\bm{n}_i} = \lambda_i^e \hat{\bm{e}_i}$ with $\lambda_i = \sqrt{1 + 2 \lambda_i^e}$ and $\ell_i = \log \lambda_i = \frac 1 2 \operatorname{\tt log1p}(2\lambda_i^e)$.

Note

It should be noted that the Ogden model can be implemented with single field displacement as described in neo-Hookean isochoric model (161) which in that case we have

(201)¶\[ \begin{aligned} \bm{S} = \frac{\partial \psi}{\partial \bm{E}} &= \frac{\partial \psi_{\text{vol}}}{\partial J} \frac{\partial J}{\partial \bm{E}} + \frac{\partial \psi_{\text{iso}}}{\partial \bm{E}} =\bm{S}_{\text{vol}} + \bm{S}_{\text{iso}}\\ &= \bulk J \, V' \, \bm{C}^{-1} + \sum_{i=1}^3 S_i^{\text{iso}} \hat{\bm{N}_i} \hat{\bm{N}_i}^T. \end{aligned} \]

Linearization¶

The general process for Newton linearization in the initial configuration is described for the neo-Hookean model Linearization: initial configuration.

The general process for Newton linearization in the current configuration is described for the neo-Hookean model Linearization: current configuration.

Command-line interface¶

To enable the hyperelastic Ogden model, use the model option -model elasticity-isochoric-ogden-current or -model elasticity-isochoric-ogden-initial and set the material parameters listed in Isochoric-split Ogden model options. Any parameter without a default option is required.

Table 25 Isochoric-split Ogden model options¶
Option	Description
`-model elasticity-isochoric-ogden-current \| elasticity-isochoric-ogden-initial`	Required to enable the isochoric-split Ogden model, `\|` denotes an option for the user.
`-E [real]`	Young’s modulus, $E > 0$
`-nu [real]`	Poisson’s ratio, $\nu \leq 0.5$.
`-m [m_1, m_2, m_3]`	Ogden material constant
`-alpha [alpha_1, alpha_2, alpha_3]`	Ogden material constant, $2 \secondlame = \sum_{j=1}^N m_j \alpha_j,$ with $m_j \alpha_j > 0$

An example using the isochoric-split Ogden model can be run via

$ ./bin/ratel-quasistatic -options_file examples/ymls/ex02-quasistatic-elasticity-isochoric-ogden-current-face-forces.yml

To enable the hyperelastic mixed Ogden model, use the model option -model elasticity-mixed-ogden-current or -model elasticity-mixed-ogden-initial and set the material parameters listed in Mixed Ogden model options. Any parameter without a default option is required.

Table 26 Mixed Ogden model options¶
Option	Description	Default value
`-model elasticity-mixed-ogden-current \| elasticity-mixed-ogden-initial`	Required to enable the mixed Ogden model, `\|` denotes an option for the user.
`-E [real]`	Young’s modulus, $E > 0$
`-nu [real]`	Poisson’s ratio, $\nu \leq 0.5$.
`-nu_primal [real]`	Primal part of the Poisson’s ratio for Jacobian operator, $-1 \leq \nu_p < \nu \leq 0.5$. This parameter may improve solver performance and should not affect the solution.	`-1.0`
`-nu_primal_pc [real]`	Primal part of the Poisson’s ratio for preconditioner operator, $-1 \leq \nu_{pc} < \nu \leq 0.5$. This parameter may improve solver performance and should not affect the solution.	`-nu_primal` value
`-m [m_1, m_2, m_3]`	Ogden material constant
`-alpha [alpha_1, alpha_2, alpha_3]`	Ogden material constant, $2 \secondlame = \sum_{j=1}^N m_j \alpha_j,$ with $m_j \alpha_j > 0$

An example using the mixed Ogden model can be run via

$ ./bin/ratel-quasistatic -options_file examples/ymls/ex02-quasistatic-elasticity-mixed-ogden-initial-pcjacobi.yml

Anisotropic models¶

Holzapfel-Gasser-Ogden (HGO) Model¶

Fiber distributions¶

The HGO model is a micromechanical model which assumes that a hyperelastic groundmatrix material, an , is reinforced by one or more bundles of fibers. Each bundle of fibers $\alpha$ is assumed to be represented by a continuous density function of fiber directions $\rho(\bm{A}_{\alpha})$ around a preferred unit direction vector $\bm{A}_{\alpha}$, defined in the initial configuration. While the possible choices for the density function $\rho$ are quite expansive — any general orthotropic fiber distribution is permissible [GOH05] — it is more convenient to only consider fiber bundles with transversely isotropic distributions. This simplifying assumption results in a compact representation of all possible fiber distribution via a generalized structure tensor $\bm{H}_{\alpha}$ which depends only on the preferred direction $\bm{A}_{\alpha}$ and a fiber dispersion parameter $\kappa \in [0, 1/3]$, given by

\[ \bm{H}_{\alpha} = \kappa \bm{1} + (1 - 3\kappa) (\bm{A}_{\alpha} \otimes \bm{A}_{\alpha}). \]

These distributions range from a Dirac delta function at $\kappa=0$, where all fibers are exactly oriented along $\bm{A}_\alpha$, to a sphere at $\kappa=1/3$, where the fibers are arranged isotropically with no preferred direction. The distributions for $\kappa \in (0, 1/3)$ are analogous to normalized $\pi$-periodic von Mises distributions such that the angle between fiber directions and the preferred direction $\bm{A}_{\alpha}$ falls within increasingly steep bell curves as $\kappa \to 0$, see [GOH05].

Note that while each fiber distribution is transversely isotropic, there is no theoretical limit on the number of fiber distributions. As such, the full composite material may be transversely isotropic, orthotropic, or generally anisotropic. Technically, different fiber distributions could have different dispersion parameters or other material parameters, although in our implementation (as well as in Abaqus and [GOH05]), all distributions are assumed to share the same material and dispersion parameters.

Strain energy function¶

The strain energy function $\psi$ is assumed to be an additive decomposition of the isotropic groundmatrix material $\psi_g$ and the anisotropic strain energy contributions of each set of fibers $\bar{\psi}_{\text{aniso},\alpha}$. As suggested by the notation, these anisotropic contributions are assumed to be isochoric in the original implementation by [GOH05]

Note

This assumption is shown to be problematic for nearly-incompressible and compressible groundmatrix materials under pure dilational or applied hydrostatic stress in [NGD+14, VDMO13]. Correct behavior can be restored by using the full strain tensor $\bm{C}$, rather than its isochoric part, when computing the anisotropic strain energy. We retain the isochoric implementation in this work in order to compare against existing work.

The total strain energy function is then given by

\[ \psi(\bm{C}, \bm{A}_{\alpha}) = \psi_g (\bm{C}) + \sum_{\alpha=1}^N \bar{\psi}_{\text{aniso},\alpha} (\bar{\bm{C}}, \bm{H}_{\alpha}), \]

The groundmatrix is formulated as an incompressible neo-Hookean material, with the strain energy function given by Perturbed Lagrange-multiplier method as

\[ \psi_g (\bm{C}) = \psi_{\text{vol}}(J) + \bar{\psi}_{\text{iso}}(\bar{\bm{C}}) = \frac{k}{2} (U(J))^2 + \frac{\mu}{2} (\bar{I}_1 - 3), \]

Note that unlike in the notation of Neo-Hookean isochoric-split, we here use an overbar to denote isochoric quantities and subscripts of “iso” and “aniso” to differentiate the isotropic and anisotropic contributions to the strain energy.

The anisotropic contributions to the strain energy are defined in terms of the structure tensors $\bm{H}_{\alpha}$. It is convenient for an efficient numerical implementation to push forward the structure tensors to the current configuration, yielding

\[ \bar{\bm{h}}_{\alpha} = \bar{\bm{F}} \bm{H}_{\alpha} \bar{\bm{F}}^T = \kappa \bar{\bm{b}} + (1 - 3\kappa) (\bar{\bm{a}}_{\alpha} \otimes \bar{\bm{a}}_{\alpha}), \]

where $\bar{\bm{a}}_{\alpha} = \bar{\bm{F}} \bm{A}_{\alpha}$ is the direction of the ${\alpha}$ fiber distribution in the current configuration. The anisotropic strain energy for the fiber distribution ${\alpha}$ is then given by

\[ \bar{\psi}_{\text{aniso},\alpha} = \frac{k_1}{2 k_2} (\exp(k_2 \left\langle \bar{E}_{\alpha} \right\rangle^2) - 1), \]

where $k_1, k_2$ are material parameters, $\left\langle x \right\rangle = \frac{1}{2} (|x| + x)$ is the Macaulay bracket function, and

\[ \bar{E}_{\alpha} = \bm{H}_{\alpha} : \bar{\bm{C}} - 1 = \trace(\bar{\bm{h}}_{\alpha}) - 1 \]

is a strain-like pseudo-invariant of the left Cauchy-Green tensor and the structure tensor. A critical physical phenomenon which this model aims to capture is that for real collagen fibers found in arterial walls, the fibers cannot resist any compressive loads. This is enforced by the Macaulay bracket function in the above, which ensures $\bar{E}_{\alpha}=1$ for compressive loads.

Note

In [GOH05], this is enforced by only including the anisotropic contributions $(1 - 3\kappa) \bm{A}_{\alpha} \otimes \bm{A}_{\alpha}$ to $\bm{H}_{\alpha}$ if $\bar{\bm{C}} : \bm{A}_{\alpha} \otimes \bm{A}_{\alpha} > 1$; that is, if the strain in the direction of the reference fiber direction is positive. The implementation in Abaqus uses the slightly modified constraint $\bar{E}_{\alpha} > 0$. These two constraints are not equivalent. For example, if $\bm{A}_{\alpha} \otimes \bm{A}_{\alpha} : \bar{\bm{C}} \leq 1$, then the approach of [GOH05] yields $\bar{E}_{\alpha} = \bm{H}_{\alpha} : \bar{\bm{C}} - 1 = \kappa \bar{I}_1 - 1$. If $ \kappa \bar{I}_1 < 1$, as is the case when $\kappa = 0$, then $\bar{E}_{\alpha} < 0$ in contradiction with the approach implemented by Abaqus. We adopt the approach of Abaqus in our implementation for purposes of verification against a known implementation.

Stress tensors¶

The second Piola-Kirchhoff stress tensor for hyperelastic materials is known to be

\[ \bm{S} = 2 \frac{\partial \psi}{\partial \bm{C}} = 2 \frac{\partial \psi_{\text{vol}}}{\partial \bm{C}} + 2 \frac{\partial \bar{\psi}}{\partial \bm{C}} + 2 \sum_{\alpha = 1}^N \frac{\partial \bar{\psi}_{\text{aniso},\alpha}}{\partial \bm{C}}. \]

Derivations of the isotropic terms can be found in Neo-Hookean isochoric-split. We focus on the anisotropic stress tensor(s) in the following. Note, the chain rule implies that

\[ \bar{\bm{S}}_{\text{aniso},\alpha} = 2 \frac{\partial \bar{\psi}_{\text{aniso},\alpha}}{\partial \bm{C}} = 2 \frac{\partial \bar{\psi}_{\text{aniso},\alpha}}{\partial \bar{E}_{\alpha}} \frac{\partial \bar{E}_{\alpha}}{\partial \bar{\bm{C}}} : \frac{\partial \bar{\bm{C}}}{\partial \bm{C}}. \]

It is convenient to treat $\bar{\psi}_{\text{aniso},\alpha}$ as a function of $\bar{E}_{\alpha}$ and to apply the chain rule. To this end, the derivatives of $\bar{\psi}_{\text{aniso},\alpha}$ are given by

(208)¶\[ \begin{aligned} \bar{\psi}'_{\text{aniso},\alpha} &= \frac{\diff \bar{\psi}_{\text{aniso},\alpha}}{\diff \bar{E}_{\alpha}} = k_1 \left\langle \bar{E}_{\alpha} \right\rangle \exp\left(k_2 \left\langle \bar{E}_{\alpha} \right\rangle^2\right) \\ {\bar{\psi}''}_{\text{aniso},\alpha} &= \frac{\diff{}^2 \bar{\psi}_{\text{aniso},\alpha}}{\diff \bar{E}_{\alpha}^2} = (k_1 k_2 \left\langle \bar{E}_{\alpha} \right\rangle^2 + k_1 H(\bar{E}_{\alpha})) \exp\left(k_2 \left\langle \bar{E}_{\alpha} \right\rangle^2\right) \end{aligned} \]

where $H(\cdot)$ is the Heavyside function. Further, $ \frac{\partial \bar{E}_{\alpha}}{\partial \bar{\bm{C}}} = \frac{\partial (\bm{H}_{\alpha} : \bar{\bm{C}} - 1)}{\partial \bar{\bm{C}}} = \bm{H}_{\alpha}$. Thus,

(209)¶\[ \bar{\bm{S}}_{\text{aniso},\alpha} = 2 J^{-2/3} \bar{\psi}_{\text{aniso},\alpha}' \bm{H}_{\alpha} : \mathbb{P}^T = 2 J^{-2/3} \bar{\psi}_{\text{aniso},\alpha}' \mathbb{P} : \bm{H}_{\alpha}, \]

where $\mathbb{P}$ is the fourth order projection tensor defined by $\mathbb{P} = \mathbb{I} - \frac{1}{3} \bm{C}^{-1} \otimes \bm{C}$ and $\mathbb{I}$ is the fourth order identity tensor. The Kirchhoff stress tensor in the current configuration is then given by (209) as,

\[ \begin{aligned} \bar{\bm{\tau}}_{\text{aniso},\alpha} = \bm{F} \bar{\bm{S}}_{\text{aniso},\alpha} \bm{F}^T &= 2 J^{-2/3} \bar{\psi}_{\text{aniso},\alpha}' \dev((J^{1/3}\bar{\bm{F}})\bm{H}_{\alpha} (J^{1/3}\bar{\bm{F}}^T)) \\ &= 2 \bar{\psi}_{\text{aniso},\alpha}' \dev(\bar{\bm{h}}_{\alpha}). \end{aligned} \]

For completeness, the first Piola-Kirchhoff stress tensor for the full material is given by

\[ \begin{aligned} \bm{P} &= \bm{F}\bm{S}_{\text{vol}} + \bm{F}\bar{\bm{S}}_{\text{iso}} + \sum_{\alpha=1}^N \bm{F}\bar{\bm{S}}_{\text{aniso},\alpha} \\ &= (k_p U - p) J U' \bm{F}\bm{C}^{-1} + 2\mu \bm{F} J^{-2/3} \bm{C}^{-1} \dev(\bm{E}) + \sum_{\alpha=1}^N 2 J^{-2/3} \bar{\psi}_{\text{aniso},\alpha}' \bm{F} (\mathbb{P} : \bm{H}_{\alpha}) \end{aligned} \]

The strong and weak formulation for current configuration is described in Strong and weak formulations: current configuration. 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 general process for Newton linearization in the initial configuration is described for the neo-Hookean model Linearization: initial configuration.

Command-line interface¶

To enable the HGO model, use the model option -model elasticity-mixed-anisotropic-neo-hookean-current and set the material parameters listed in HGO model options. Any parameter without a default option is required.

Table 27 HGO model options¶
Option	Description	Default Value
`-model elasticity-mixed-anisotropic-neo-hookean-current`	Required to enable the HGO model.
`-E [real]`	Young’s modulus, $E > 0$
`-nu [real]`	Poisson’s ratio, $\nu \leq 0.5$. Note, this model is only validated for $\nu = 0.5$.
`-nu_primal [real]`	Primal part of the Poisson’s ratio for Jacobian operator, $-1 \leq \nu_p < \nu \leq 0.5$. This parameter may improve solver performance and should not affect the solution. See Linear elastic for more details.	`-1.0`
`-nu_primal_pc [real]`	Primal part of the Poisson’s ratio for preconditioner operator, $-1 \leq \nu_p < \nu \leq 0.5$. This parameter may improve solver performance and should not affect the solution. See Linear elastic for more details.	`-nu_primal` value
`-kappa [real]`	Fiber dispersion parameter, $\kappa\in[0,1/3]$. See Fiber distributions.	`0.0`
`-k1 [positive real]`	Stress-like anisotropic parameter, $k_1 > 0$.
`-k2 [positive real]`	Dimensionless exponential coeficient anisotropic parameter, $k_2 > 0$.
`-directions [x1,y1,z1][,x2,y2,z2][,x3,y3,z3]`	Reference directions for transversely anisotropic fiber bundles. Up to 3 directions may be specified. Directions do not need to be normalized or orthogonal, but should not be collinear.