Holzapfel-Gasser-Ogden (HGO) Anisotropic Hyperelastic 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-Volumetric 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-Volumetric 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

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

(122)\[ \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 (122) 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} \]

Linearization

Representation of \(\diff \bar{\bm{S}}_{\text{aniso},\alpha} \) for Holzapfel-Gasser-Odgen models

As shown in [Hol00], \(\diff \bm{S} = \mathbb{C} : \diff \bm{E} = 2 \frac{\partial \bm{S}}{\partial \bm{C}} : \diff \bm{E}\), where \(\diff \bm{E} = \frac{1}{2} (\diff \bm{F}^T \bm{F} + \bm{F}^T \diff \bm{F})\) and \(\mathbb{C} = \mathbb{C}_{\text{vol}} + \bar{\mathbb{C}} + \sum_{\alpha} \bar{\mathbb{C}}_{\text{aniso},\alpha}\) is the fourth-order elasticity tensor. Again, we derive only the anisotropic part of \(\mathbb{C}\), which for each fiber distribution \(\alpha\) is given as

\[ \begin{aligned} \bar{\mathbb{C}}_{\text{aniso},\alpha} = 2 \frac{\partial \bar{\bm{S}}_{\text{aniso},\alpha}}{\partial \bm{C}} &= 4 J^{-4/3} \bar{\psi}_{\text{aniso},\alpha}'' (\mathbb{P}:\bm{H}_{\alpha})\otimes (\mathbb{P} : \bm{H}_{\alpha}) \\ &\quad - \frac{4}{3} J^{-2/3}\bar{\psi}_{\text{aniso},\alpha}'((\mathbb{P} : \bm{H}_{\alpha}) \otimes \bm{C}^{-1} \\ &\quad+ \bm{C}^{-1} \otimes \bm{H}_{\alpha} + (\bm{H}_{\alpha} : \bm{C})\frac{\partial \bm{C}^{-1}}{\partial \bm{C}}) \end{aligned} \]

where \(\bar{\psi}_{\text{aniso},\alpha}'\) and \(\bar{\psi}_{\text{aniso},\alpha}''\) are given in (121). Then, \(\diff \bar{\bm{S}}_{\text{aniso},\alpha}\) can be found as

(123)\[ \begin{aligned} \diff \bar{\bm{S}}_{\text{aniso},\alpha} &= \bar{\mathbb{C}}_{\text{aniso},\alpha} : \diff \bm{E} \\ &= 4 J^{-4/3} \bar{\psi}_{\text{aniso},\alpha}'' ((\mathbb{P}:\bm{H}_{\alpha}) : \diff \bm{E})(\mathbb{P} : \bm{H}_{\alpha}) \\ &\quad - \frac{4}{3} J^{-2/3}\bar{\psi}_{\text{aniso},\alpha}'((\bm{C}^{-1} : \diff \bm{E})(\mathbb{P} : \bm{H}_{\alpha}) \\ &\quad + (\bm{H}_{\alpha} : \diff \bm{E})\bm{C}^{-1} + (\bm{H}_{\alpha} : \bm{C})\bm{C}^{-1} \diff \bm{E} \bm{C}^{-1}). \end{aligned} \]

Representation of \(\bm F \diff \bar{\bm{S}}_{\text{aniso},\alpha} \bm F^T\) for Holzapfel-Gasser-Odgen models

We push forward (123) to get

\[ \begin{aligned} \bm{F}\diff \bar{\bm{S}}_{\text{aniso},\alpha} \bm{F}^T &= 4 \bar{\psi}_{\text{aniso},\alpha}'' (\dev \bar{\bm{h}}_{\text{aniso},\alpha} : \diff \bm{\epsilon})\dev \bar{\bm{h}}_{\text{aniso},\alpha} \\ &\quad- 4/3 \bar{\psi}_{\text{aniso},\alpha}' (\trace \diff \bm{\epsilon}\dev \bar{\bm{h}}_{\text{aniso},\alpha} \\ &\quad- (\bar{E}_{\alpha} + 1) \dev \diff \bm{\epsilon} + (\dev \bar{\bm{h}}_{\text{aniso},\alpha} : \diff \bm{\epsilon}) \bm{1}), \end{aligned} \]

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 11 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 Incompressibility 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 Incompressibility 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.

Example usage of the HGO model

An example using the HGO model can be run via

$ ./bin/ratel-quasistatic examples/ymls/ex02-quasistatic-elasticity-mixed-anisotropic-neo-hookean-current-pcfieldsplit-pmg.yml

The example YML file shows the model options in action, and includes sample displacement curves for two values of \(\kappa\). This example is based on the canonical arterial wall validation problem from [GOH05], which can be reproduced exactly under mesh refinement.

# Ratel problem specification
# Mixed Neo-Hookean hyperelastic material at finite strain

orders: 2,0

dm_plex:
  dim: 3
  box_faces: 5,2,2
  box_upper: 5,1.5,0.25
  simplex: 0

model: elasticity-mixed-anisotropic-neo-hookean-current
E: 0.02292
nu: 0.5
nu_primal: 0.3
nu_primal_pc: 0.3
kappa: 0.226
k1: 0.9966
k2: 524.6
directions: 0.76582002124983765498,-0.64305497047522943672,0, -0.76582002124983765498,-0.64305497047522943672,0

bc:
  slip: 1,3,6
  slip_6_components: 0
  slip_3_components: 1
  slip_1_components: 2
  clamp: 5
  clamp_5_translate: 1.32346761,0,0 # kappa=0.226
  # clamp_5_translate: 2.12120795,0,0 # kappa=0

pc:
  use_amat:
  type: fieldsplit
  fieldsplit:
    type: schur
    schur_fact_type: upper
    schur_precondition: a11
    off_diag_use_amat:

fieldsplit:
  displacement:
    ksp_type: preonly
    pc_type: pmg
  pressure:
    ksp_type: preonly
    pc_type: jacobi

ksp:
  type: gmres
  norm_type: preconditioned
  gmres_modifiedgramschmidt:
  gmres_restart: 200

snes:
  linesearch_type: bt

ts:
  dt: 0.01
  max_time: 0.1

expected_strain_energy: 4.390917228133e-05

surface_force_faces: 5