6. Hyperelasticity

We’ve derived various measures of strain based on deformations that are presumably known. We’ve also derived various measures of stress, emphasizing the physical meaning of the stress tensor, for the case when normal and shear stresses on a body are known. To solve problems where the objective is to find deformations (strains) and internal stresses within a structure, knowing, for example, the forces externally applied to the structure, we need to know the stiffnesses of the materials involved. Just as the stiffness of a spring relates force to deformation, as Hooke famously showed in 1678 “ut tensio, sic vis,” meaning,“as the extension, so the force,” the stiffness of a material can be thought of as relating stress to strain. This relationship is called a constitutive relationship.

While many engineering materials, in the elastic regime, can be idealized as linear and only undergo infinitesimal strain, many other types of materials, like rubbery elastomers, behave nonlinearly and reversible (elastic) even when subjected to very large deformations. While in practice, such materials typically display viscoelastic behavior even at low strain rates, at infinitesimally slow strain rates, they approach hyperelastic behavior. Moreover, hyperelastic material models underlie even the most complex models that incorporate viscoelasticity, damage, etc. Thus, hyperelasticity is essential for modeling polymers in FEA. Several specific kinds of rubber models will be considered in examples at the end of this chapter.

A hyperelastic model is a general constitutive relationship, relating stress and strain, for large deformation (large strain), nonlinear, elastic (reversible), isotropic, materials. While different materials can be best characterized by different “strain energy density functions,” and certainly have different material constants (found via experimentation), our general constitutive expression will be derived in closed form, and in 3D.

It may be useful to note that “curve-fitting” techniques that finds engineering stress as a function of engineering strain, based on uniaxial stress-strain experimental data, for example, are becoming more common. Extending such a method to accommodate deformation in three dimensions is a challenge that usually necessitates that we start with a 3D “hyperelastic” energy function. As we will see in the examples at the end of this chapter, “hyperelastic” functions that relate stress and strain, while nonlinear in nature, actually do not require much higher order algebra, yet have been proven to work for a wide variety of hyperelastic materials. We will see that these functions are fairly simple in the sense that they require the determination of only a few experimental constants. While these experimental constants are typically found from uniaxial tests or shear tests, the hyperelastic constitutive relationships are always written in tensor form (3D), using stress and strain tensors that form a work-conjugate pair [Hill]. This enables our hyperelastic constitutive relationships to be implemented into FEA quite naturally.

Hyperelastic functions fall into two major categories: micro-mechanical models and phenomenological models. This will be discussed toward the end of this chapter, and several hyperelastic functions will be presented. In addition, it is important to be able to fully characterize a particular material once a hyperelastic function has been chosen. Experimental data is necessary for this purpose, and several methods for efficiently characterizing materials will be demonstrated. Two methods, which are particularly “cutting edge,” will be presented at the conclusion of this chapter.

This may be an appropriate time to point out that this text will develop a constitutive relationship for two branches of elasticity: hyperelasticity (this chapter) and linear infinitesimal elasticity (final chapter). Given that hyperelasticity encompasses elastic material nonlinearity under large strain, one might ask how we should approach a material that is elastic linear under large strain, or elastic nonlinear under small strain. It turns out that these special cases would be treated just like any nonlinear hyperelastic material – i.e. one would use the basic hyperelastic constitutive relationship that will be derived in this chapter. The constitutive relationship that will be derived in the final chapter, on linear infinitesimal elasticity, will be valid only for materials that are linear and subjected to small deformation (less than 1% strain) [Freed][Bergstrom].

The following will simply be stated, but a proof can be found in Appendix C.1:

(1)   \begin{equation*} \boldsymbol{\sigma}=g(\mathbf{B}) \end{equation*}

note: Eq. 1 is perhaps obvious since \boldsymbol{\sigma^*}=\mathbf{Q} \cdot \boldsymbol{\sigma} \cdot \mathbf{Q}^{T} and \mathbf{B^*}=\mathbf{Q} \cdot \mathbf{B} \cdot \mathbf{Q}^{T} from last chapter (we called this a work-conjugate pair based on derivation from \mathbf{Q} superimposed on \mathbf{dx}).

Now, we know that the strain energy is related to the product of force and deformation.

So, let’s define \phi to be the strain energy density function.
A basic explanation of the concept of strain energy density (strain energy per unit volume) can be found in [Ugural]. A rigorous derivation of strain energy density (Helmholtz free energy) can be found in [Bergstrom].

\boldsymbol{\hat{\sigma}}=\frac{d\phi}{d\mathbf{E}}=2\frac{d\phi}{d\mathbf{C}}, since \mathbf{E}=\frac{1}{2}(\mathbf{C-I})

Remember, \boldsymbol{\hat{\sigma}} and \mathbf{E} (or \mathbf{C}) are a work-conjugate pair and from eq. 2 in Section 4: Alternative Measures of Stress we know:

\boldsymbol{\hat{\sigma}}=(det\mathbf{F})\mathbf{F}^{-1} \cdot \boldsymbol{\sigma} \cdot \mathbf{F}^{-T}=2\frac{d\phi}{d\mathbf{C}}
\longrightarrow \boldsymbol{\sigma}=\frac{2}{det\mathbf{F}}\mathbf{F}\frac{d\phi}{d\mathbf{C}} \cdot \mathbf{F}^{T}

\frac{d\phi}{d\mathbf{C}}=\mathbf{R}^{T} \cdot \frac{d\phi}{d\mathbf{B}} \cdot \mathbf{R} ; The proof is straightforward (can be inferred from Appendix C.1), so long as it is recognized that \phi=\phi(\mathbf{\xi}), where \xi can be replaced by either \mathbf{C} or \mathbf{B}, as they have the same eigenvalues and invariants, for example.

So, \frac{2}{det\mathbf{F}}\mathbf{F}\frac{d\phi}{d\mathbf{C}} \cdot \mathbf{F}^{T}=\frac{2}{det\mathbf{F}}\mathbf{V} \cdot \underbrace{\mathbf{R} \cdot \mathbf{R}^{T}}_{\mathbf{I}} \frac{d\phi}{d\mathbf{B}} \cdot \underbrace{\mathbf{R} \cdot \mathbf{R}^{T}}_{\mathbf{I}} \cdot \mathbf{V}

det\mathbf{F}=det(\mathbf{V \cdot R})=det\mathbf{V}, since det\mathbf{R}=1 (rigid body rotation)
(det\mathbf{V}=det\mathbf{F}=1 if material is incompressible)


(2)   \begin{equation*} \boldsymbol{\sigma}=\frac{2}{det\mathbf{V}}\mathbf{V}\frac{d\phi}{d\mathbf{B}}\cdot \mathbf{V} \end{equation*}

Since \mathbf{B=V}^2 \longrightarrow \boldsymbol{\sigma}=\frac{2}{det(\mathbf{B}^{\frac{1}{2}})}\mathbf{B}\frac{d\phi}{d\mathbf{B}}
(symmetrized form)

(3)   \begin{equation*} \boldsymbol{\sigma}=\frac{1}{det\mathbf{B}^{1/2}}\left(\mathbf{B}\frac{d\phi}{d\mathbf{B}}+\frac{d\phi}{d\mathbf{B}}\mathbf{B}\right) \end{equation*}

Eq. 3 is possible since \mathbf{B}\frac{d\phi}{d\mathbf{B}}=\frac{d\phi}{d\mathbf{B}}\mathbf{B}

For isotropy, it is possible to write \phi=\phi(I_B, II_B, III_B), where I_B, II_B, III_B are each a function of \mathbf{B} [Hughes]. This is called the “representation theorem for invariants” [Holzapfel]. Since we will look at specific strain energy functions, \phi, later on, which happen to be of such a form, we need to apply the chain rule to \frac{d\phi}{d\mathbf{B}}.

(4)   \begin{equation*} \longrightarrow \frac{d\phi}{d\mathbf{B}}=\frac{\partial \phi}{\partial I_B}\frac{d I_B}{d \mathbf{B}}+\frac{\partial \phi}{\partial II_B}\frac{d II_B}{d \mathbf{B}}+ \frac{\partial \phi}{\partial III_B}\frac{d III_B}{d \mathbf{B}} \end{equation*}

We need to determine \frac{d I_B}{d \mathbf{B}}, \frac{d II_B}{d \mathbf{B}}, and \frac{d III_B}{d \mathbf{B}}:

The general expressions for the second order tensor invariants I_B and II_B will just be stated, as it is a math topic:

I_B=tr\mathbf{B} ; II_B=\frac{1}{2}[(tr\mathbf{B})^2-tr(\mathbf{B^2})] (skipped work)

The complete derivations of \frac{d I_B}{d \mathbf{B}}, \frac{d II_B}{d \mathbf{B}}, and \frac{d III_B}{d \mathbf{B}} are given in Appendix C.2.

(5)   \begin{equation*} \frac{dI_B}{d\mathbf{B}}=\mathbf{I} \end{equation*}

(6)   \begin{equation*} \frac{d II_B}{d \mathbf{B}} = I_B\mathbf{I}-\mathbf{B} \end{equation*}

(7)   \begin{equation*} \frac{dIII_B}{d\mathbf{B}}=\mathbf{B}^2-I_B\mathbf{B}+II_B\mathbf{I} \end{equation*}

Eq. 5, Eq. 6, Eq. 7 \longrightarrow Eq. 4

(8)   \begin{equation*} \longrightarrow \frac{d\phi}{d\mathbf{B}}=\frac{\partial \phi}{\partial I_B}\mathbf{I}+\frac{\partial \phi}{\partial II_B}(I_B\mathbf{I}-\mathbf{B})+\frac{\partial \phi}{\partial III_B}(\mathbf{B}^2-I_B\mathbf{B}+II_B\mathbf{I}) \end{equation*}

Recall \boldsymbol{\sigma}=\frac{1}{det\mathbf{B}^{1/2}}(\mathbf{B}\frac{d \phi}{d\mathbf{B}}+\frac{d\phi}{d\mathbf{B}}\mathbf{B})

\frac{1}{det\mathbf{B}^{1/2}}\bigg(\mathbf{B}\bigg[\frac{\partial \phi}{\partial I_B}I_B+\frac{\partial \phi}{\partial II_B}(I_B\mathbf{I}-\mathbf{B})+\frac{\partial \phi}{\partial III_B}(\mathbf{B}^2-I_B\mathbf{B}+II_B\mathbf{I})\bigg]
+\bigg[\frac{\partial \phi}{\partial I_B}I_B+\frac{\partial \phi}{\partial II_B}(I_B\mathbf{I}-\mathbf{B})+\frac{\partial \phi}{\partial III_B}(\mathbf{B}^2-I_B\mathbf{B}+II_B\mathbf{I})\bigg]\mathbf{B}\bigg)
=\frac{1}{det\mathbf{B}^{1/2}}\bigg[\frac{\partial \phi}{\partial I_B}\mathbf{B}-\frac{\partial \phi}{\partial II_B}\mathbf{B}^2+\frac{\partial \phi}{\partial II_B}I_B \mathbf{B}+\frac{\partial \phi}{\partial III_B}\mathbf{B}^3+\frac{\partial \phi}{\partial III_B}II_B \mathbf{B}-\frac{\partial \phi}{\partial III_B}I_B \mathbf{B}^2
+\frac{\partial \phi}{\partial I_B}\mathbf{B}-\frac{\partial \phi}{\partial II_B}\mathbf{B}^2+\frac{\partial \phi}{\partial II_B}I_B\mathbf{B}+\frac{\partial \phi}{\partial III_B}\mathbf{B}^3+\frac{\partial \phi}{\partial III_B}II_B \mathbf{B}-\frac{\partial \phi}{\partial III_B}I_B \mathbf{B}^2\bigg]
=\frac{2}{det\mathbf{B}^{1/2}}\bigg[\frac{\partial \phi}{\partial III_B}[\mathbf{B}^3+II_B\mathbf{B}-I_B\mathbf{B}^2]+(\frac{\partial \phi}{\partial I_B}+I_B\frac{\partial \phi}{\partial II_B})\mathbf{B}-\frac{\partial \phi}{\partial II_B}\mathbf{B}^2\bigg]

We know from eq. 1 in Section 1: Trace, Scalar Product, Eigenvalues: III_B\mathbf{I}=\mathbf{B}^3-I_B\mathbf{B}^2+II_B\mathbf{B}

So, \mathbf{B}^2=I_B\mathbf{B}-II_B\mathbf{I}+III_B\mathbf{B}^{-1}


\longrightarrow \boldsymbol{\sigma}=\frac{2}{det\mathbf{B}^{1/2}}\bigg[III_B\frac{\partial \phi}{\partial III_B}\mathbf{I}+\frac{\partial \phi}{\partial I_B}\mathbf{B}+\cancel{I_B\frac{\partial \phi}{\partial II_B}\mathbf{B}}-\cancel{\frac{\partial \phi}{\partial II_B}I_B\mathbf{B}}
+\frac{\partial \phi}{\partial II_B}II_B\mathbf{I} -\frac{\partial \phi}{\partial II_B}III_B\mathbf{B}^{-1}\bigg]

(9)   \begin{equation*} =\frac{2}{det\mathbf{B}^{1/2}}[(III_B\frac{\partial \phi}{\partial III_B}+\frac{\partial \phi}{\partial II_B}II_B)\mathbf{I}+(\frac{\partial \phi}{\partial I_B})\mathbf{B}-(III_B\frac{\partial \phi}{\partial II_B})\mathbf{B}^{-1}] \end{equation*}

Sometimes in other literature (see [Asaro] [Malvern] [Lubarda]) II_B is taken as \frac{1}{2}[tr(\mathbf{B^2})-(tr\mathbf{B})^2] instead of \frac{1}{2}[(tr\mathbf{B})^2-tr(\mathbf{B^2})]. If this is the case, then the last term in eq. 9 would be added instead of subtracted, the Cayley-Hamilton Theorem would become \mathbf{B}^3-I_B\mathbf{B}^2-II_B\mathbf{B}-III_B\mathbf{I}=0 instead of \mathbf{B}^3-I_B\mathbf{B}^2+II_B\mathbf{B}-III_B\mathbf{I}=0, and any strain energy functions that are a function of II_B (see examples to follow) would also need to be modified accordingly.


det\mathbf{F}=1 (III_B=1) ; \phi=\phi(I_B,II_B)

Modifying eq. 8, we get: \frac{d\phi}{d\mathbf{B}}=(\frac{\partial \phi}{\partial I_B}+I_B\frac{\partial \phi}{\partial II_B})\mathbf{I}-\frac{\partial \phi}{\partial II_B}\mathbf{B}

Substituting this simplified expression into eq. 3, we get

(10)   \begin{equation*} \boldsymbol{\sigma}=-\rho_0 \mathbf{I}+\underbrace{2\left[\left(\frac{\partial \phi}{\partial I_B}+I_B\frac{\partial \phi}{\partial II_B}\right)\mathbf{B}-\frac{\partial \phi}{\partial II_B}\mathbf{B}^2\right]}_{\text{stress for given deformation}} \end{equation*}

The underbraced expression is what you would set as the stress for a given deformation. But we can always superimpose an arbitrary pressure (e.x. hydrostatic) to the body without causing deformation. Thus, the basic constitutive expression doesn’t uniquely specify stress \therefore we introduce \rho_0.  Sometimes this hydrostatic pressure can be determined by boundary conditions.  Otherwise, an additional equation-of-state (EOS) may be needed in order to solve for the pressure.  FEA will often assume a “nearly incompressible” Poisson’s Ratio of 0.495 for rubber so that the additional EOS is not required.

We could’ve also modified eq. 9 with III_B=1 and \frac{\partial \phi}{\partial III_B}=0 and we would get

(11)   \begin{equation*} \boldsymbol{\sigma}=-\rho_0\mathbf{I}+2\left[\frac{\partial \phi}{\partial II_B}II_B\mathbf{I}+\frac{\partial \phi}{\partial I_B}\mathbf{B}-\frac{\partial \phi}{\partial II_B}\mathbf{B}^{-1}\right] \end{equation*}

Then, we can use the Cayley-Hamilton Theorem to substitute for \mathbf{B}^{-1} and arrive at the desired result (eq. 10).

Note, however, that in eq. 11, the first term in the square brackets influences all \sigma_{ii} (diagonal) terms equally just as \rho_0 does. Thus, this first term can be thought of as another pressure, and can, accordingly, be lumped together with \rho_0.


(12)   \begin{equation*} \boldsymbol{\sigma}=-\rho_0\mathbf{I}+2\left[\left(\frac{\partial \phi}{\partial I_B}\right)\mathbf{B}-\left(\frac{\partial \phi}{\partial II_B}\right)\mathbf{B}^{-1}\right] \end{equation*}

To convince yourself that eq. 12 is correct, go ahead and use either eq. 11 or eq. 12 in the following Mooney-Rivlin examples. You will see that the results do not change.

Next: 6. Phenomenological and Micromechanical Models