Polynomial Hyperelasticity

Overview

The polynomial hyperelastic model is a model in which the strain energy is a polynomial function of the invariants of the deformation tensor. The polynomial hyperelastic model

• is a phenomenologial model,
• has an arbitrary number of parameters, though most materials use fewer than 4
• is typically fit to uniaxial extension/compression, equibiaxial extension/compression, or shear experimental data
• usually valid for strains less than 100%

See Also

• User Defined Materials
• Linear Elastic Material
• Mooney-Rivlien Hyperelasticity

Contents

1. Fundamental Equations
2. Model Implementation
3. Model Verification

Fundamental Equations

The strain energy in the polynomial hyperelastic model is defined as

$$ W = \sum_{i,j=0}^N c_{ij}\left(\overline{I}_1 - 3\right)^i \left(\overline{I}_2 - 3\right)^j + + \frac{1}{D_1}\left(J-1\right)^2 $$

where the $c_{ij}$ and $D_1$ are material parameters, the $\overline{I}_i$ are the isochoric invariants of the right Cauchy deformation tensor $\pmb{C} = \pmb{F}^T{\cdot}\pmb{F}$, where $\pmb{F}$ is deformation gradient tensor, and $J$ is the determinant of the deformation gradient.

The Second Piola-Kirchhoff Stress Tensor

$$ \pmb{T} = 2\frac{\partial W}{\partial \pmb{C}} $$

For isotropic hyperelasticity $W=W\left(\overline{I}_1, \overline{I}_2, J\right)$ and

$$ \pmb{T} = 2\left(\frac{\partial W}{\partial \overline{I}_1}\frac{\partial \overline{I}_1}{\partial \pmb{C}} + \frac{\partial W}{\partial \overline{I}_2}\frac{\partial \overline{I}_2}{\partial \pmb{C}} + \frac{\partial W}{\partial J}\frac{\partial J}{\partial \pmb{C}}\right) = \pmb{A} \cdot \pmb{B} $$

where

$$ \pmb{A} = \left[\frac{\partial W}{\partial \overline{I}_1} \quad \frac{\partial W}{\partial \overline{I}_2} \quad \frac{\partial W}{\partial J}\right] $$

and

\begin{align} \pmb{B} &= \left[ \frac{\partial \overline{I}_1}{\partial \pmb{C}} \quad \frac{\partial \overline{I}_2}{\partial \pmb{C}} \quad \frac{\partial J}{\partial \pmb{C}} \right] \\ &= \left[ J^{-2/3}\left(\pmb{I} - \frac{1}{3}I_1 \pmb{C}^{-1} \right) \quad J^{-4/3}\left( I_1\pmb{I} - \pmb{C} - \frac{2}{3}I_2\pmb{C}^{-1} \right)\quad \frac{1}{2}J\pmb{C}^{-1} \right] \end{align}

Elastic stiffness in the material frame is given by

\begin{align} \mathbb{L} &= 4\frac{\partial^2 W}{\partial\pmb{C}\partial\pmb{C}} = 4\frac{\partial}{\partial \pmb{C}}\left( \pmb{A}\cdot\pmb{B}\right) \\ &= 4\left(\frac{\partial \pmb{A}}{\partial \pmb{C}}\cdot\pmb{B} + \pmb{A}\cdot\frac{\partial \pmb{B}}{\partial \pmb{C}}\right) \end{align}

where

\begin{equation} \frac{\partial \pmb{A}}{\partial \pmb{C}} = \pmb{H} \cdot \pmb{B}, \quad \pmb{H} = \begin{bmatrix} \frac{\partial^2 W}{\partial \overline{I}_1 \partial\overline{I}_1} & \frac{\partial^2 W}{\partial \overline{I}_1 \partial \overline{I}_2} & \frac{\partial^2 W}{\partial \overline{I}_1 \partial J} \\ \frac{\partial^2 W}{\partial \overline{I}_2 \partial \overline{I}_1} & \frac{\partial^2 W}{\partial \overline{I}_2 \partial \overline{I}_2} & \frac{\partial^2 W}{\partial \overline{I}_2 \partial J} \\ \frac{\partial^2 W}{\partial J \partial \overline{I}_1} & \frac{\partial^2 W}{\partial J \partial \overline{I}_2} & \frac{\partial^2 W}{\partial J \partial J} \end{bmatrix} \end{equation}

and

\begin{equation} \frac{\partial \pmb{B}}{\partial \pmb{C}} = \begin{Bmatrix} \frac{1}{3}J^{-2/3}\left[ \pmb{I}\pmb{C}^{-1} - \pmb{C}^{-1}\pmb{I} - I_1\left( \pmb{C}^{-1}\odot\pmb{C}{^-1} - \frac{1}{3}\pmb{C}^{-1}\pmb{C}{^-1}\right) \right] \\ \frac{2}{3}J^{-4/3}\left[ \frac{3}{2}\left(\mathbb{I}_1 - \mathbb{I}_2\right) + \left(\pmb{C}^{-1}\pmb{C} + \pmb{C}\pmb{C}^{-1} \right) - I_1\left( \pmb{C}^{-1}\pmb{I} + \pmb{I}\pmb{C}^{-1} \right) - I_2\left( \pmb{C}^{-1}\odot\pmb{C}^{-1} - \frac{2}{3}\pmb{C}^{-1}\pmb{C}^{-1} \right) \right] \\ \frac{1}{4} J\left( \pmb{C}^{-1}\pmb{C}^{-1} - 2\pmb{C}^{-1}\odot\pmb{C}^{-1} \right) \end{Bmatrix} \end{equation}

The operator $\odot$ is defined such that

$$ \mathbb{X}_{ijkl} = \pmb{A} \odot \pmb{B} = (A_{ik} B_{jl} + A_{il} B_{jk}) / 2 $$

Requirements of Objectivity

The constitutive equations for the hyperelastic material evaluate the stress directly in the reference configuration. The components of the stress are identified as the components of the Second Piola-Kirchhoff stress $\pmb{T}$. The push forward of the Second Piola-Kirchhoff stress gives the Cauchy stress $\pmb{\sigma}$ in the spatial configuration, as required by most finite element packages. The push forward of the corresponding material stiffness $\mathbb{L}$ does not, however, correspond to the rate of Cauchy stress but the Truesdell rate of the Cauchy stress. Furthermore, the rate of Cauchy stress is not objective, requiring finite element packages to use other so-called objective stress rates in the incremental solution of the momentum equation. In the following sections, it is demonstrated that the push forward of the material stiffness does not correspond to the rate of Cauchy stress and equations relating the stiffness corresponding to the Jaumann rate of the Kirchhoff stress to the push forward of the material stiffness are developed.

The rate of change of Cauchy stress

\begin{align} \dot{\pmb{\sigma}} &= \frac{d}{dt}\left( \frac{1}{J}\pmb{F}\cdot\pmb{T}\cdot{\pmb{F}}^T \right) \\ &= \frac{1}{J}\left[ -\frac{\dot{J}}{J}\pmb{F}\cdot\pmb{T}\cdot{\pmb{F}}^T + \dot{\pmb{F}}\cdot\pmb{T}\cdot{\pmb{F}}^T + \pmb{F} \cdot\frac{d\pmb{T}}{d\pmb{E}}{:}\dot{\pmb{E}} \cdot{\pmb{F}}^T + \pmb{F}\cdot\pmb{T}\cdot\dot{\pmb{F}}^T \right] \end{align}

With the following identities

\begin{equation} \mathrm{tr}{\pmb{d}} = \frac{\dot{J}}{J} \quad \dot{\pmb{F}} = \pmb{L}\cdot\pmb{F} \quad \dot{\pmb{E}} = \pmb{F}^T\cdot\pmb{d}\cdot\pmb{F} \quad \dot{\pmb{F}}^T = \pmb{F}^T\cdot\pmb{L}^T \end{equation}

\begin{align} \dot{\pmb{\sigma}} &= \frac{1}{J}\left[ -\mathrm{tr}{\pmb{d}}\pmb{\sigma} + \pmb{L}\cdot\pmb{F}\cdot\pmb{T}\cdot{\pmb{F}}^T + \pmb{F}\cdot\mathbb{L}{:}\left( \pmb{F}^T\cdot\pmb{d}\cdot\pmb{F} \right)\cdot\pmb{F}^T + \pmb{F}\cdot\pmb{T}\cdot\pmb{F}^T\cdot\pmb{L}^T \right] \\ &= -\mathrm{tr}{\pmb{d}}\pmb{\sigma} + \pmb{L}\cdot\pmb{\sigma} + \frac{1}{J}\left( \pmb{F}\cdot\pmb{F}\cdot \mathbb{L} \cdot\pmb{F}^T\cdot\pmb{F}^T \right){:}\pmb{d} + \pmb{\sigma}\cdot\pmb{L}^T \\ &= -\mathrm{tr}{\pmb{d}}\pmb{\sigma} + \pmb{L}\cdot\pmb{\sigma} + \mathbb{C}{:}\pmb{d} + \pmb{\sigma}\cdot\pmb{L}^T \end{align}

rearranging

\begin{align} \dot{\pmb{\sigma}} - \pmb{\sigma}\cdot\pmb{L}^T - \pmb{L}\cdot\pmb{\sigma} + \mathrm{tr}{\pmb{d}}\pmb{\sigma} &= \mathbb{C}{:}\pmb{d} \\ \mathring{\pmb{\sigma}} &= \mathbb{C}{:}\pmb{d} \end{align}

where

\begin{align} \mathring{\pmb{\sigma}} = \dot{\pmb{\sigma}} - \pmb{\sigma}\cdot\pmb{L}^T - \pmb{L}\cdot\pmb{\sigma} + \mathrm{tr}{\pmb{d}}\pmb{\sigma} \end{align}

is the Truesdell rate of the Cauchy stress. The Truesdell rate of the Cauchy stress is related to the Lie derivative of the Kirchhoff stress by

\begin{equation} \mathring{\pmb{\sigma}} = \frac{1}{J}\pmb{F}\cdot\dot{\pmb{T}}\pmb{F}^T = \frac{1}{J}\pmb{F}\cdot\left[\frac{d}{dt}\left( J\pmb{F}^{-1}\cdot\pmb{\sigma}\pmb{F}^{-T} \right)\right]\cdot\pmb{F}^T = \frac{1}{J}\pmb{F}\cdot\left[\frac{d}{dt}\left( \pmb{F}^{-1}\cdot\pmb{\tau}\pmb{F}^{-T} \right)\right]\cdot\pmb{F}^{T} = \frac{1}{J}\mathscr{L}\left({\pmb{\tau}}\right) \end{equation}

Thus, the push forward of the material time derivative of $\pmb{T}$ is the Truesdell rate of the Cauchy stress. Correspondingly, the Truesdell rate of the Kirchhoff stress is related to the push forward of the material stiffness $\mathbb{L}$ by

\begin{equation} \mathscr{L}\left({\pmb{\tau}}\right) = J\mathbb{C}{:}\pmb{d} \end{equation}

In many finite element packages, the Jaumann rate of the Kirchhoff stress, and not the Truesdell rate, is required. Thus, the elastic stiffness must correspond to the Jaummann rate.. The Jaumann rate of the Kirchhoff stress is given by

\begin{align} \stackrel{\nabla}{\pmb{\tau}} &= { \dot{\pmb{\tau}} + \pmb{\tau}\cdot\pmb{w} - \pmb{w}\cdot\pmb{\tau} } \\ &= \dot{J}\pmb{\sigma} + J\dot{\pmb{\sigma}} + J\pmb{\sigma}\cdot\pmb{w} - \pmb{w}\cdot J\pmb{\sigma} \\ &= J\left(\mathrm{tr}{\pmb{d}}\pmb{\sigma} + \dot{\pmb{\sigma}} + \pmb{\sigma}\cdot\pmb{w} - \pmb{w}\cdot J\pmb{\sigma}\right) \\ &= J\mathbb{D}{:}\pmb{d} \end{align}

where $\mathbb{D}$ is the stiffness corresponding to the Jaumann rate. Subtracting $\mathring{\pmb{\tau}}$ from $\stackrel{\nabla}{\pmb{\tau}}$, the Jaumann stiffness can be cast in terms of $\mathbb{C}$

\begin{align} \left(\mathbb{D} - \mathbb{C}\right){:}\pmb{d} &= \left(\mathrm{tr}{\pmb{d}}\pmb{\sigma} + \dot{\pmb{\sigma}} + \pmb{\sigma}\cdot\pmb{w} - \pmb{w}\cdot\pmb{\sigma}\right) - \left(\dot{\pmb{\sigma}} - \pmb{\sigma}\cdot{\pmb{L}}^T - \pmb{L}\cdot\pmb{\sigma} + \mathrm{tr}{\pmb{d}}\pmb{\sigma}\right) \\ &= \pmb{\sigma}\cdot\pmb{w} + \pmb{\sigma}\cdot{\pmb{L}}^T + \pmb{L}\cdot\pmb{\sigma} - \pmb{w}\cdot\pmb{\sigma} \\ &= \pmb{\sigma}\cdot\pmb{d} + \pmb{d}\cdot\pmb{\sigma} \\ \end{align}

Using indicial notation, and the fact that $\pmb{d}$ and $\pmb{\sigma}$ are symmetric,

\begin{align} \left(D_{ijkl} - C_{ijkl}\right)d_{kl} &= \sigma_{im}d_{mj} + d_{im}\sigma_{mj} \\ &= \left( \sigma_{ik}\delta_{jl} + \delta_{il}\sigma_{jk} \right)d_{kl} \end{align}

from which

\begin{align} \left( D_{ijkl} - C_{ijkl} - \sigma_{ik}\delta_{jl} - \delta_{il}\sigma_{jk} \right)d_{kl} = 0_{ij} \end{align}

Since the preceding equation must hold for all $d_{kl}$, the stiffness corresponding to the Jaumman rate of the Kirchhoff stress is related to the stiffness corresponding to the Truesdell rate of the Kirchhoff stress as

\begin{align} D_{ijkl} = C_{ijkl} + \sigma_{ik}\delta_{jl} + \delta_{il}\sigma_{jk} \end{align}

Note that the stiffness must be made minor symmetric

\begin{align} D_{ijkl} = C_{ijkl} + \frac{1}{2}\left( \sigma_{ik}\delta_{jl} + \sigma_{il}\delta_{jk} + \delta_{il}\sigma_{jk} + \delta_{ik}\sigma_{jl} \right) \end{align}

In symbolic notation,

$$ \mathbb{D} = \mathbb{C} + \pmb{\sigma}\odot\pmb{I} + \pmb{I}\odot\pmb{\sigma} $$

which is the correct stiffness corresponding to the Jaumann rate of the Kirchhoff stress, in terms of the push forward of the material stiffness. Equations relating $\mathbb{C}$ to stiffnesses corresponding to other objective rates are derived similarly.

Matmodlab Implementation

A polynomial material model is implemented as a standard Matmodlab material in matmodlab2/materials/polyhyper.py as a subclass of Material class. The model defines the following (required) attributes and methods

• name: name by which the model is referenced
• eval: method the updates the material state

Additionally, several helper functions are imported from various locations in Matmodlab:

• matmodlab.materials.tensor

  • symsq: Computes $\pmb{F}^T{\cdot}\pmb{F}$ and returns the result as a 6D first order tensor
  • det, inv: Computes the determinant and inverse of second order tensors (stored as matrix or array)
  • invariants: Computes the invariants of a second order tensor
  • push: Performs the push forward operation (for both second and fourth order tensors)
  • dyad: Computes the dyadic product of two vectors or second order tensors
  • symshuffle: Computes the product $X_{ijkl} = .5 (A_{ik} B_{jl} + A_{il} B_{jk}$)
  • II1, II5: Fourth order identity tensors
  • I6: Identity tensor stored as an array of length 6

The model implementation can be viewed by executing the following cell.

Verification

Uniaxial Stress

For an incompressible isotropic material, uniaxial stress is produced by the following deformation state

$$ [F] = \begin{bmatrix}\lambda && \\ & \frac{1}{\sqrt{\lambda}} & \\ & & \frac{1}{\sqrt{\lambda}} \end{bmatrix} $$

The stress difference $\sigma_{\text{axial}} - \sigma_{\text{lateral}}$ is given by

We now exercise the Mooney-Rivlin material model using Matmodlab