Hyperelastic material behavior

Definitions and basic kinematic results

We first introduce some definitions and basic kinematic results that will be used in this section. Some of these items have already been discussed in Introduction and Basic Equations: they are repeated here for convenience.

Writing the current position of a material point as $x$ and the reference position of the same point as $X$ , the deformation gradient is

F \overset{def}{=} \frac{\partial x}{\partial X} .

Then J, the total volume change at the point, is

J \overset{def}{=} det (F) .

For simplicity, we define

\bar{F} \overset{def}{=} J^{- \frac{1}{3}} F

as the deformation gradient with the volume change eliminated.

We then introduce the deviatoric stretch matrix (the left Cauchy-Green strain tensor) of $\bar{F}$ as

\bar{B} \overset{def}{=} \bar{F} \cdot {\bar{F}}^{T}

so that we can define the first strain invariant as

\begin{matrix} (1) & {\bar{I}}_{1} \overset{def}{=} trace \bar{B} = I : \bar{B}, \end{matrix}

where $I$ is a unit matrix, and the second strain invariant as

\begin{matrix} (2) & {\bar{I}}_{2} \overset{def}{=} \frac{1}{2} ({\bar{I}}_{1}^{2} - trace (\bar{B} \cdot \bar{B})) = \frac{1}{2} ({\bar{I}}_{1}^{2} - I : \bar{B} \cdot \bar{B}) . \end{matrix}

The variations of $\bar{B}$ , $\bar{B} \cdot \bar{B}$ , ${\bar{I}}_{1}$ , ${\bar{I}}_{2}$ , and J will be required during the remainder of the development. We first define some variations of basic kinematic quantities that will be needed to write these results.

The gradient of the displacement variation with respect to current position is written as

δ L \overset{def}{=} \frac{\partial δ u}{\partial x} .

The virtual rate of deformation is the symmetric part of $δ L$ :

δ D \overset{def}{=} sym (δ L) = \frac{1}{2} (δ L + δ L^{T}),

which we decompose into the virtual rate of change of volume per current volume (the “virtual volumetric strain rate”),

δ ε^{vol} \overset{def}{=} I : δ D,

and the virtual deviatoric strain rate,

δ e \overset{def}{=} δ D - \frac{1}{3} δ ε^{vol} I .

The virtual rate of spin of the material is the antisymmetric part of $δ L$ :

δ W \overset{def}{=} asym (δ L) = \frac{1}{2} (δ L - δ L^{T}) .

The variations of $\bar{B}$ , $\bar{B} \cdot \bar{B}$ , ${\bar{I}}_{1}$ , ${\bar{I}}_{2}$ , and J are obtained directly from their definitions above in terms of these quantities as

δ \bar{B} = δ e \cdot \bar{B} + \bar{B} \cdot δ e + δ W \cdot \bar{B} - \bar{B} \cdot δ W = H_{1} : δ e + δ W \cdot \bar{B} - \bar{B} \cdot δ W,

where

{(H_{1})}_{i j k l} \overset{def}{=} \frac{1}{2} (δ_{i k} {\bar{B}}_{j l} + {\bar{B}}_{i k} δ_{j l} + δ_{i l} {\bar{B}}_{j k} + {\bar{B}}_{i l} δ_{j k});

\begin{matrix} δ (\bar{B} \cdot \bar{B}) & ​ = δ e \cdot \bar{B} \cdot \bar{B} + \bar{B} \cdot \bar{B} \cdot δ e + 2 \bar{B} \cdot δ e \cdot \bar{B} + δ W \cdot \bar{B} \cdot \bar{B} - \bar{B} \cdot \bar{B} \cdot δ W, \\ = H_{2} : δ e + δ W \cdot \bar{B} \cdot \bar{B} - \bar{B} \cdot \bar{B} \cdot δ W, \end{matrix}

where

{(H_{2})}_{i j k l} \overset{def}{=} \frac{1}{2} (δ_{i k} {\bar{B}}_{j p} {\bar{B}}_{p l} + {\bar{B}}_{i p} {\bar{B}}_{p k} δ_{j l} + δ_{i l} {\bar{B}}_{j p} {\bar{B}}_{p k} + {\bar{B}}_{i p} {\bar{B}}_{p l} δ_{j k}) + {\bar{B}}_{i k} {\bar{B}}_{j l} + {\bar{B}}_{i l} {\bar{B}}_{j k};

\begin{matrix} (3) & δ {\bar{I}}_{1} = 2 \bar{B} : δ e; \end{matrix}

\begin{matrix} (4) & δ {\bar{I}}_{2} = 2 ({\bar{I}}_{1} \bar{B} - \bar{B} \cdot \bar{B}) : δ e; \end{matrix}

and

\begin{matrix} (5) & δ J = J δ ε^{vol} . \end{matrix}

The Cauchy (“true”) stress components are defined from the strain energy potential as follows. From the virtual work principal the internal energy variation is

δ W_{I} = \int_{V} σ : δ D d V = \int_{V^{0}} J σ : δ D d V^{0},

where $σ$ are the components of the Cauchy (“true”) stress, V is the current volume, and $V^{0}$ is the reference volume.

We decompose the stress into the equivalent pressure stress,

p \overset{def}{=} - \frac{1}{3} I : σ,

and the deviatoric stress,

S \overset{def}{=} σ + p I,

so that the internal energy variation can be written

δ W_{I} = \int_{V^{0}} J (S : δ e - p δ ε^{vol}) d V^{0} .

For isotropic, compressible materials the strain energy, U, is a function of ${\bar{I}}_{1}$ , ${\bar{I}}_{2}$ , and J:

U = U ({\bar{I}}_{1}, {\bar{I}}_{2}, J),

so that

δ U = \frac{\partial U}{\partial {\bar{I}}_{1}} δ {\bar{I}}_{1} + \frac{\partial U}{\partial {\bar{I}}_{2}} δ {\bar{I}}_{2} + \frac{\partial U}{\partial J} δ J .

Hence, using Equation 3, Equation 4, and Equation 5,

\begin{matrix} (6) & δ U = 2 [(\frac{\partial U}{\partial {\bar{I}}_{1}} + {\bar{I}}_{1} \frac{\partial U}{\partial {\bar{I}}_{2}}) \bar{B} - \frac{\partial U}{\partial {\bar{I}}_{2}} \bar{B} \cdot \bar{B}] : δ e + J \frac{\partial U}{\partial J} δ ε^{vol} . \end{matrix}

Since the variation of the strain energy potential is, by definition, the internal virtual work per reference volume, $δ W_{I}$ , we have

\begin{matrix} (7) & δ W_{I} = \int_{V^{0}} J (S : δ e - p δ ε^{vol}) d V^{0} = \int_{V^{0}} δ U d V^{0} . \end{matrix}

For a compressible material the strain variations are arbitrary, so this equation defines the stress components for such a material as

\begin{matrix} (8) & S = \frac{2}{J} DEV [(\frac{\partial U}{\partial {\bar{I}}_{1}} + {\bar{I}}_{1} \frac{\partial U}{\partial {\bar{I}}_{2}}) \bar{B} - \frac{\partial U}{\partial {\bar{I}}_{2}} \bar{B} \cdot \bar{B}] \end{matrix}

and

\begin{matrix} (9) & p = - \frac{\partial U}{\partial J} . \end{matrix}

When the material response is almost incompressible, the pure displacement formulation, in which the strain invariants are computed from the kinematic variables of the finite element model, can behave poorly. One difficulty is that from a numerical point of view the stiffness matrix is almost singular because the effective bulk modulus of the material is so large compared to its effective shear modulus, thus causing difficulties with the solution of the discretized equilibrium equations. Similarly, in Abaqus/Explicit the high bulk modulus increases the dilatational wave speed, thus reducing the stable time increment substantially. Another problem is that, unless reduced-integration techniques are used, the stresses calculated at the numerical integration points show large oscillations in the pressure stress values, because—in general—the elements cannot respond accurately and still have small volume changes at all numerical integration points. To avoid such problems, Abaqus/Standard offers a “mixed” formulation for such cases. The concept is to introduce a variable, $\hat{J}$ , that is used in place of the volume change, J, in the definition of the strain energy potential. The internal energy integral, $W_{I}$ , is augmented with the constraint that $J - \hat{J} = 0$ , imposed by the use of a Lagrange multiplier, $\hat{p}$ , and integrated over the volume:

{(W_{I})}^{A} \overset{def}{=} \int_{V^{0}} [U ({\bar{I}}_{1}, {\bar{I}}_{2}, \hat{J}) - \hat{p} (J - \hat{J})] d V^{0} .

Taking the variation of this definition,

δ {(W_{I})}^{A} = \int_{V^{0}} [J S : δ e + (\frac{\partial U}{\partial \hat{J}} + \hat{p}) δ \hat{J} - J \hat{p} δ ε^{vol} - (J - \hat{J}) δ \hat{p}] d V^{0} .

Since $δ \hat{J}$ is an independent variation in this expression, the Lagrange multiplier is

\hat{p} = - \frac{\partial U}{\partial \hat{J}},

and its variation is

δ \hat{p} = - \hat{Q} : δ e + \frac{\partial \hat{p}}{\partial \hat{J}} δ \hat{J},

where

\hat{Q} \overset{def}{=} J \frac{\partial S}{\partial \hat{J}} .

These results allow us to write the augmented internal energy variation as

\begin{matrix} (10) & δ {(W_{I})}^{A} = \int_{V^{0}} [(J S + (J - \hat{J}) \hat{Q}) : δ e - J \hat{p} δ ε^{vol} - (J - \hat{J}) \frac{\partial \hat{p}}{\partial \hat{J}} δ \hat{J}] d V^{0}, \end{matrix}

which implies that continuity of the $\hat{J}$ interpolation across elements is not required.

This augmented formulation can be used for any value of compressibility except fully incompressible behavior. For most element types $\hat{J}$ is interpolated independently in each element: Abaqus uses constant $\hat{J}$ in most first-order elements and linear variation of $\hat{J}$ with respect to position in second-order elements. The only element type where continuity of the $\hat{J}$ interpolation across elements is enforced is C3D4H: a first-order tetrahedron with a linear interpolation of $\hat{J}$ continuous across elements.

When the material is fully incompressible, U is a function of the first and second strain invariants— ${\bar{I}}_{1}$ and ${\bar{I}}_{2}$ —only, and we write the internal energy in the augmented form,

{(W_{I})}^{A} \overset{def}{=} \int_{V^{0}} [U - \hat{p} (J - 1)] d V^{0},

where $\hat{p}$ is again a Lagrange multiplier introduced to impose the constraint $J - 1 = 0$ in such a way that the variation of ${(W_{I})}^{A}$ can be taken with respect to all kinematic variables, thus giving

\begin{matrix} (11) & δ {(W_{I})}^{A} = \int_{V^{0}} [J S : δ e - J \hat{p} δ ε^{vol} - (J - 1) δ \hat{p}] d V^{0} . \end{matrix}

The Lagrange multiplier $\hat{p}$ is interpolated in the same way as $\hat{J}$ is interpolated in the augmented formulation for almost incompressible behavior; that is, $\hat{p}$ is assumed to be constant in most first-order elements and to vary linearly with respect to position in second-order elements. The only exception is element type C3D4H, which is a first-order tetrahedron with a linear interpolation of $\hat{p}$ continuous across elements.

Rate of change of the internal virtual work

The rate of change of the internal virtual work is required for use in the Newton method, which is generally used in Abaqus/Standard to solve the nonlinear equilibrium equations (after discretization by finite elements). It will also be used when we extend the elasticity model to viscoelastic behavior for small (linearized) vibrations about a predeformed state.

When the pure displacement formulation is used for the compressible case, the deviatoric stress components, $S$ , are defined by Equation 8, from which we can show that

d (J S) = J (C^{S} : d e + Q d ε^{vol} + d W \cdot S - S \cdot d W),

where the “effective deviatoric elasticity” of the material, $C^{S}$ , is defined as

\begin{matrix} C^{S} \overset{def}{=} & \frac{2}{J} (\frac{\partial U}{\partial {\bar{I}}_{1}} + {\bar{I}}_{1} \frac{\partial U}{\partial {\bar{I}}_{2}}) H_{1} - \frac{2}{J} \frac{\partial U}{\partial {\bar{I}}_{2}} H_{2} + \frac{4}{J} (\frac{\partial^{2} U}{\partial {\bar{I}}_{1}^{2}} + \frac{\partial U}{\partial {\bar{I}}_{2}} + 2 {\bar{I}}_{1} \frac{\partial^{2} U}{\partial {\bar{I}}_{1} \partial {\bar{I}}_{2}} + {\bar{I}}_{1}^{2} \frac{\partial^{2} U}{\partial {\bar{I}}_{2}^{2}}) \bar{B} \bar{B} \\ - \frac{4}{J} (\frac{\partial^{2} U}{\partial {\bar{I}}_{1} \partial {\bar{I}}_{2}} + {\bar{I}}_{1} \frac{\partial^{2} U}{\partial {\bar{I}}_{2}^{2}}) (\bar{B} \cdot \bar{B} \bar{B} + \bar{B} \bar{B} \cdot \bar{B}) + \frac{4}{J} \frac{\partial^{2} U}{\partial {\bar{I}}_{2}^{2}} \bar{B} \cdot \bar{B} \bar{B} \cdot \bar{B} \\ - \frac{4}{3 J} [\frac{\partial U}{\partial {\bar{I}}_{1}} + 2 {\bar{I}}_{1} \frac{\partial U}{\partial {\bar{I}}_{2}} + {\bar{I}}_{1} \frac{\partial^{2} U}{\partial {\bar{I}}_{1}^{2}} + ({\bar{I}}_{1}^{2} + 2 {\bar{I}}_{2}) \frac{\partial^{2} U}{\partial {\bar{I}}_{1} \partial {\bar{I}}_{2}} + 2 {\bar{I}}_{1} {\bar{I}}_{2} \frac{\partial^{2} U}{\partial {\bar{I}}_{2}^{2}}] (I \bar{B} + \bar{B} I) \\ + \frac{4}{3 J} (2 \frac{\partial U}{\partial {\bar{I}}_{2}} + {\bar{I}}_{1} \frac{\partial^{2} U}{\partial {\bar{I}}_{1} \partial {\bar{I}}_{2}} + 2 {\bar{I}}_{2} \frac{\partial^{2} U}{\partial {\bar{I}}_{2}^{2}}) (I \bar{B} \cdot \bar{B} + \bar{B} \cdot \bar{B} I), \end{matrix}

and the deviatoric stress rate-volumetric strain rate coupling term, $Q$ , is

Q \overset{def}{=} \frac{\partial (J S)}{\partial J} = 2 (\frac{\partial^{2} U}{\partial {\bar{I}}_{1} \partial J} + {\bar{I}}_{1} \frac{\partial^{2} U}{\partial {\bar{I}}_{2} \partial J}) \bar{B} - 2 \frac{\partial^{2} U}{\partial {\bar{I}}_{2} \partial J} \bar{B} \cdot \bar{B} - \frac{2}{3} ({\bar{I}}_{1} \frac{\partial^{2} U}{\partial {\bar{I}}_{1} \partial J} + 2 {\bar{I}}_{2} \frac{\partial^{2} U}{\partial {\bar{I}}_{2} \partial J}) I .

From Equation 9 it can be shown that

d (J p) = - J (Q : d e + K d ε^{vol}),

where K is the effective bulk modulus of the material,

K \overset{def}{=} - (J \frac{\partial p}{\partial J} + p) = J \frac{\partial^{2} U}{\partial J^{2}} + \frac{\partial U}{\partial J} .

Thus,

\begin{matrix} (12) & d δ W_{I} = \int_{V} [⌊ \begin{matrix} δ e & δ ε^{vol} \end{matrix} ⌋ : [\begin{matrix} C^{S} & Q \\ Q & K \end{matrix}] : {\begin{matrix} d e \\ d ε^{vol} \end{matrix}} - σ : (2 δ ε \cdot d ε - δ L^{T} \cdot d L)] d V, \end{matrix}

since

δ e : (d W \cdot S - S \cdot d W) + S : d δ e - p d δ ε^{vol} = - σ : (2 δ ε \cdot d ε - δ L^{T} \cdot d L) .

For the mixed formulation introduced for almost incompressible materials, the rate of change of the augmented variation of internal energy, Equation 10, is

\begin{matrix} d δ {(W_{I})}^{A} = \int_{V} [ & ⌊ \begin{matrix} δ e & δ ε^{vol} & δ \hat{J} \end{matrix} ⌋ [\begin{matrix} {\tilde{C}}^{S} & \hat{Q} & \hat{A} \\ \hat{Q} & - \hat{p} & - \partial \hat{p} / \partial \hat{J} \\ \hat{A} & - \partial \hat{p} / \partial \hat{J} & - \tilde{K} \end{matrix}] {\begin{matrix} d e \\ d ε^{vol} \\ d \hat{J} \end{matrix}} \\ - (σ + \frac{1}{J} (J - \hat{J}) \hat{Q}) : (2 δ ε \cdot d ε - δ L^{T} \cdot d L)] d V, \end{matrix}

where

{\tilde{C}}^{S} \overset{def}{=} C^{S} + (J - \hat{J}) \frac{\partial C^{S}}{\partial \hat{J}},

\hat{Q} \overset{def}{=} J \frac{\partial S}{\partial \hat{J}},

\hat{A} \overset{def}{=} \frac{1}{J} (J - \hat{J}) \frac{\partial \hat{Q}}{\partial \hat{J}},

and

\tilde{K} \overset{def}{=} \frac{1}{J^{2}} (\hat{K} - (J - \hat{J}) \frac{\partial \hat{K}}{\partial \hat{J}}),

in which

\hat{K} \overset{def}{=} - J \frac{\partial \hat{p}}{\partial \hat{J}} .

For the case of incompressible materials the rate of change of the augmented variation of internal energy is similarly obtained from Equation 11 as

\begin{matrix} d δ {(W_{I})}^{A} = \int_{V} [ & ⌊ \begin{matrix} δ e & δ ε^{vol} & δ \hat{p} \end{matrix} ⌋ [\begin{matrix} C^{S} & 0 & 0 \\ 0 & - \hat{p} & - 1 \\ 0 & - 1 & 0 \end{matrix}] {\begin{matrix} d e \\ d ε^{vol} \\ d \hat{p} \end{matrix}} \\ - σ : (2 δ ε \cdot d ε - δ L^{T} \cdot d L)] d V . \end{matrix}

Particular forms of the strain energy potential

Several particular forms of the strain energy potential are available in Abaqus. The incompressible or almost incompressible models make up:

the polynomial form and its particular cases—the reduced polynomial form, the neo-Hookean form, the Mooney-Rivlin form, and the Yeoh form;
the Ogden form;
the Arruda-Boyce form; and
the Van der Waals form.

In addition, a hyperelastic model for highly compressible, elastic materials is offered.

Polynomial form and particular cases

Given isotropy and additive decomposition of the deviatoric and volumetric strain energy contributions in the presence of incompressible or almost incompressible behavior, we can write the potential, in general, as

U = f ({\bar{I}}_{1} - 3, {\bar{I}}_{2} - 3) + g (J_{e ℓ} - 1) .

Setting $g = \sum_{i = 1}^{N} \frac{1}{D_{i}} {(J_{e ℓ} - 1)}^{2 i}$ and expanding $f ({\bar{I}}_{1} - 3, {\bar{I}}_{2} - 3)$ in a Taylor series, we arrive at

\begin{matrix} (13) & U \overset{def}{=} \sum_{i + j = 1}^{N} C_{i j} {({\bar{I}}_{1} - 3)}^{i} {({\bar{I}}_{2} - 3)}^{j} + \sum_{i = 1}^{N} \frac{1}{D_{i}} {(J_{e ℓ} - 1)}^{2 i} . \end{matrix}

This form is the polynomial representation of the strain energy in Abaqus. The parameter N can take values up to six; however, values of N greater than 2 are rarely used when both the first and second invariants are taken into account. $C_{i j}$ and $D_{i}$ are temperature-dependent material parameters. The value of N and tables giving the $C_{i j}$ and $D_{i}$ values as functions of temperature are specified by the user. The elastic volume strain, $J_{e ℓ}$ , follows from the total volume strain, J, and the thermal volume strain, $J_{t h}$ , with the relation

J_{e ℓ} = \frac{J}{J_{t h}},

and $J_{t h}$ follows from the linear thermal expansion, $ε_{t h}$ , with

J_{t h} = {(1 + ε_{t h})}^{3},

where $ε_{t h}$ follows from the temperature and the isotropic thermal expansion coefficient defined by the user.

The $D_{i}$ values determine the compressibility of the material: if all the $D_{i}$ are zero, the material is taken as fully incompressible. If $D_{1} = 0$ , all $D_{i}$ must be zero.

Regardless of the value of N, the initial shear modulus, $μ_{0},$ and the bulk modulus, $k_{0},$ depend only on the polynomial coefficients of order $N = 1$ :

μ_{0} = 2 (C_{10} + C_{01}), k_{0} = \frac{2}{D_{1}} .

If $N = 1$ , so that only the linear terms in the deviatoric strain energy are retained, the Mooney-Rivlin form is recovered:

U = C_{10} ({\bar{I}}_{1} - 3) + C_{01} ({\bar{I}}_{2} - 3) + \frac{1}{D_{1}} {(J_{e ℓ} - 1)}^{2} .

The Mooney-Rivlin form can be viewed as an extension of the neo-Hookean form (discussed below) in that it adds a term that depends on the second invariant of the left Cauchy-Green tensor. In some cases this form will give a more accurate fit to the experimental data than the neo-Hookean form; in general, however, both models give similar accuracy since they use only linear functions of the invariants. These functions do not allow representation of the “upturn” at higher strain levels in the stress-strain curve.

Particular forms of the polynomial model can also be obtained by setting specific coefficients to zero. If all $C_{i j}$ with $j \neq 0$ are set to zero, the reduced polynomial form is obtained:

U = \overset{N}{\sum_{i = 1}} C_{i 0} {({\bar{I}}_{1} - 3)}^{i} + \overset{N}{\sum_{i = 1}} \frac{1}{D_{i}} {(J_{e ℓ} - 1)}^{2 i} .

Following Yeoh (1993) the justification for reducing the general polynomial series expansion by omitting the dependence on the second invariant arises from the following observations. The sensitivity of the strain energy function to changes in the second invariant is generally much smaller than the sensitivity to changes in the first invariant. In addition, the ${\bar{I}}_{2}$ -dependence is difficult to measure, so it might be preferable to neglect it rather than to calculate it based on potentially inaccurate measurements. Finally, it appears that omitting the dependence on the second invariant if data for only a particular mode of deformation are known might enhance the prediction for other deformation states. This conjecture is supported by investigating the so-called reduced stresses in the presence of almost incompressible behavior:

\begin{matrix} \frac{σ_{1} - σ_{2}}{λ_{1}^{2} - λ_{2}^{2}} = ​ & 2 (\frac{\partial U}{\partial {\bar{I}}_{1}} + λ_{3}^{2} \frac{\partial U}{\partial {\bar{I}}_{2}}), \\ \frac{σ_{1} - σ_{3}}{λ_{1}^{2} - λ_{3}^{2}} = ​ & 2 (\frac{\partial U}{\partial {\bar{I}}_{1}} + λ_{2}^{2} \frac{\partial U}{\partial {\bar{I}}_{2}}), \\ \frac{σ_{2} - σ_{3}}{λ_{2}^{2} - λ_{3}^{2}} = ​ & 2 (\frac{\partial U}{\partial {\bar{I}}_{1}} + λ_{1}^{2} \frac{\partial U}{\partial {\bar{I}}_{2}}), \end{matrix}

where the $σ_{i}$ , $i = 1 \dots 3$ represent the principal Cauchy (“true”) stresses. If the derivatives with respect to ${\bar{I}}_{2}$ are omitted and different stress states—uniaxial, biaxial, and planar—are considered, the reduced stresses have the same form regardless of the stress state.

Measurements of the ${\bar{I}}_{2}$ -dependence of carbon-black reinforced rubber vulcanizates confirming these findings can be found in Kawabata, Yamashita, et al. (1995). The paper of Kaliske and Rothert (1997) also supports the notion that often the prediction of general deformation states based on a uniaxial measurement can be enhanced only by ignoring the ${\bar{I}}_{2}$ -dependence.

In this context it is worth noting that the mathematical structure of the Arruda-Boyce model can be viewed as a fifth-order reduced polynomial, where the five coefficients $C_{10} \dots C_{50}$ are implicit nonlinear functions of the two parameters $μ$ and $λ_{m}$ in the Arruda-Boyce form. However, the Arruda-Boyce model offers a physical interpretation of the parameters, which the general fifth-order reduced polynomial fails to provide.

The Yeoh form (Yeoh, 1993) can be viewed as a special case of the reduced polynomial with $N = 3$ :

U = \overset{3}{\sum_{i = 1}} C_{i 0} {({\bar{I}}_{1} - 3)}^{i} + \overset{3}{\sum_{i = 1}} \frac{1}{D_{i}} {(J_{e ℓ} - 1)}^{2 i} .

Typically, if $C_{10} = O (1)$ , the second coefficient will be negative and one to two orders of magnitude smaller [i.e., $C_{20}$ is $- O (0.1)$ to $- O (0.01)$ ], while the third coefficient $C_{30}$ is again one to two orders of magnitude smaller but positive [i.e., $C_{30}$ is $+ O (1. E - 2)$ to $+ O (1. E - 4)$ ]. These magnitudes will create the typical S-shape of the stress-strain behavior of rubber; at low strains $C_{10}$ represents the initial shear modulus, which softens at moderate strains due to the effect of the negative second coefficient $C_{20}$ and is followed by an upturn at large strains due to the positive third coefficient $C_{30}$ . Thus, the Yeoh model often provides an accurate fit over a large strain range.

If the reduced-polynomial strain-energy function is simplified further by setting $N = 1$ , the neo-Hookean form is obtained:

U = C_{10} ({\bar{I}}_{1} - 3) + \frac{1}{D_{1}} {(J_{e ℓ} - 1)}^{2} .

This form is the simplest hyperelastic model and often serves as a prototype for elastomeric materials in the absence of accurate material data. It also has some theoretical relevance since the mathematical representation is analogous to that of an ideal gas: the neo-Hookean potential represents the Helmholtz free energy of a molecular network with Gaussian chain-length distribution (see Treloar, 1975).

The user can request that Abaqus calculate the $C_{i j}$ and $D_{i}$ values from measurements of nominal stress and strain in simple experiments. The basis of this calculation is described in Fitting of hyperelastic and hyperfoam constants.

Ogden form

The Ogden strain energy potential is expressed in terms of the principal stretches. In Abaqus the following formulation is used:

\begin{matrix} (14) & U \overset{def}{=} \sum_{i = 1}^{N} \frac{2 μ_{i}}{α_{i}^{2}} ({\bar{λ}}_{1}^{α_{i}} + {\bar{λ}}_{2}^{α_{i}} + {\bar{λ}}_{3}^{α_{i}} - 3) + \sum_{i = 1}^{N} \frac{1}{D_{i}} {(J_{e ℓ} - 1)}^{2 i}, \end{matrix}

where

{\bar{λ}}_{i} = J^{- \frac{1}{3}} λ_{i} \to {\bar{λ}}_{1} {\bar{λ}}_{2} {\bar{λ}}_{3} = 1.

Hence, the first part of Ogden's strain energy function depends only on ${\bar{I}}_{1}$ and ${\bar{I}}_{2}$ . Ogden's energy function cannot be written explicitly in terms of ${\bar{I}}_{1}$ and ${\bar{I}}_{2}$ . It is, however, possible to obtain closed-form expressions for the derivatives of U with respect to ${\bar{I}}_{1}$ and ${\bar{I}}_{2}$ .

The value of N and tables giving the $μ_{i}$ and $α_{i}$ values as functions of temperature are specified by the user. If $N = 2$ , $α_{1} = 2$ , and $α_{2} = - 2$ , the Mooney-Rivlin model is obtained. If $N = 1$ and $α_{1} = 2$ , Ogden's model degenerates to the neo-Hookean material model. In the Ogden form the initial shear modulus, $μ_{0}$ , depends on all coefficients:

μ_{0} = \underset{i = 1}{\sum^{N}} μ_{i},

and the initial bulk modulus, $k_{0}$ , depends on $D_{1}$ as before. The user can request that Abaqus calculate the $μ_{i}$ and $α_{i}$ values from measurements of nominal stress and strain.

Arruda-Boyce form

The hyperelastic Arruda-Boyce potential has the following form:

U = μ \overset{5}{\sum_{i = 1}} \frac{C_{i}}{λ_{m}^{2 i - 2}} ({\bar{I}}_{1}^{i} - 3^{i}) + \frac{1}{D} (\frac{J_{e ℓ}^{2} - 1}{2} - \ln J_{e ℓ}),

where

C_{1} = \frac{1}{2}, C_{2} = \frac{1}{20}, C_{3} = \frac{11}{1050}, C_{4} = \frac{19}{7000}, and C_{5} = \frac{519}{673750} .

The deviatoric part of the strain energy density comes from Arruda and Boyce (1993). This model is also known as the eight-chain model, since it was developed starting out from a representative volume element where eight springs emanate from the center of a cube to its corners. The values of the coefficients $C_{1} \dots C_{5}$ arise from a series expansion of the inverse Langevin function, which arises in the statistical treatment of non-Gaussian chains. The series expansion is truncated after the fifth term. The coefficient $λ_{m}$ is referred to as the locking stretch. Approximately at this stretch the slope of the stress-strain curve will rise significantly. The initial shear modulus, $μ_{0}$ , is related to $μ$ with the expression

μ_{0} = μ (1 + \frac{3}{5 λ_{m}^{2}} + \frac{99}{175 λ_{m}^{4}} + \frac{513}{875 λ_{m}^{6}} + \frac{42039}{67375 λ_{m}^{8}}) .

A typical value of $λ_{m}$ is 7, for which $μ_{0} = 1.0125 μ$ .

The initial bulk modulus is obtained as $K_{0} = 2 / D$ . To the deviatoric part of the strain energy density we add a simplified representation of the volumetric strain energy density, which requires only one material parameter, so that all material parameters can be estimated easily even with limited knowledge of the material behavior. This volumetric representation has been used successfully by Kaliske and Rothert (1997) and provides sufficient accuracy for most engineering elastomeric materials.

The Arruda-Boyce potential depends on the first invariant only. The physical interpretation is that the eight chains are stretched equally under the action of a general deformation state. The first invariant, ${\bar{I}}_{1} = λ_{1}^{2} + λ_{2}^{2} + λ_{3}^{2}$ , directly represents this elongation.

When the Arruda-Boyce form is chosen, the user can specify the coefficients as functions of temperature; alternatively, Abaqus can perform a fit of the test data specified by the user to determine the coefficients.

Van der Waals form

The hyperelastic Van der Waals potential, also known as the Kilian model, has the following form:

U = μ {- (λ_{m}^{2} - 3) [\ln (1 - η) + η] - \frac{2}{3} a {(\frac{\tilde{I} - 3}{2})}^{\frac{3}{2}}} + \frac{1}{D} (\frac{J_{e ℓ}^{2} - 1}{2} - \ln J_{e ℓ}),

where

\tilde{I} = (1 - β) {\bar{I}}_{1} + β {\bar{I}}_{2} and η = \sqrt{\frac{\tilde{I} - 3}{λ_{m}^{2} - 3}} .

The name “Van der Waals” draws on the analogy in the thermodynamic interpretation of the equations of state for rubber and gas. While the neo-Hookean model can be compared with an ideal gas in that it starts out from a Gaussian network with no mutual interaction between the “quasi-particles” (Kilian, 1981), the Van der Waals strain energy potential is analogous to the equations of state of a real gas. This introduces two additional material parameters: the locking stretch, $λ_{m}$ , and the global interaction parameter, a. (Similarly, the Van der Waals equation for a real gas introduces two parameters to account for excluded volume and modified exchange of momentum between the particles.)

The locking stretch, $λ_{m}$ , accounts for finite extendability of the non-Gaussian chain network. In contrast to the Arruda-Boyce model the mathematical structure of the Van der Waals potential is such that the strain energy tends to infinity as the locking stretch, $λ_{m}$ , is reached; more precisely, as $\tilde{I} \to λ_{m}^{2}$ . Thus, the Van der Waals potential cannot be used at stretches larger than the locking stretch.

The global interaction parameter, a, models the interaction between the chains; it is difficult to estimate. Kilian et al. (1986) point out that, given Mooney-Rivlin coefficients and a locking stretch $λ_{m}$ , a suitable value for the global interaction parameter is

a = \frac{2 C_{01}}{3 μ} + \frac{λ_{m}^{2}}{λ_{m}^{3} - 1},

where $μ$ is the initial shear modulus at low strains and $C_{01}$ is the second Mooney-Rivlin parameter. Given a positive initial shear modulus, $μ$ , and locking stretch, $λ_{m}$ , too large a positive interaction parameter, a, will lead to Drucker instability in the tensile range. Realistic values of the global interaction parameter, a, will contribute to the characteristic S-shape of tensile stress-strain curves of rubber in the middle strain range before the final upturn as the locking stretch is approached, without causing instability.

The parameter $β$ represents a linear mixture parameter combining both invariants ${\bar{I}}_{1}$ and ${\bar{I}}_{2}$ into $\tilde{I}$ ; for $β = 0.0$ , the Van der Waals potential will be dependent on the first invariant only. Admissible values for this parameter are $0.0 \leq β \leq 1.0$ .

When the Van der Waals potential is chosen, the user can specify the coefficients as functions of temperature; alternatively, the parameters can be fitted from user-defined test data. The data fitting procedure will not necessarily yield a value of $β$ between zero and one. If during the curve fitting procedure the parameter $β$ leaves the admissible range, the curve fitting procedure is aborted and restarted with a fixed value of $β = 0.0$ . Alternatively, the curve fitting procedure can be used with a user-defined value of $β$ .

Strain energy potential for highly compressible elastomers

While the previous forms are intended for incompressible or almost incompressible materials, the elastic foam energy function is designed for describing highly compressible elastomers (Storåkers, 1986). This energy function has the form

\begin{matrix} (15) & U \overset{def}{=} \sum_{i = 1}^{N} \frac{2 μ_{i}}{α_{i}^{2}} [{\hat{λ}}_{1}^{α_{i}} + {\hat{λ}}_{2}^{α_{i}} + {\hat{λ}}_{3}^{α_{i}} - 3 + \frac{1}{β_{i}} (J_{e ℓ}^{- α_{i} β_{i}} - 1)], \end{matrix}

where

{\hat{λ}}_{i} = J_{t h}^{- \frac{1}{3}} λ_{i} \to {\hat{λ}}_{1} {\hat{λ}}_{2} {\hat{λ}}_{3} = J_{e ℓ} .

The volumetric and the deviatoric contributions are coupled in this expression, which can be demonstrated clearly by writing the expression in the form

\begin{matrix} (16) & U = \sum_{i = 1}^{N} \frac{2 μ_{i}}{α_{i}^{2}} [J_{e ℓ}^{\frac{1}{3} α_{i}} ({\bar{λ}}_{1}^{α_{i}} + {\bar{λ}}_{2}^{α_{i}} + {\bar{λ}}_{3}^{α_{i}} - 3) + 3 (J_{e ℓ}^{\frac{1}{3} α_{i}} - 1) + \frac{1}{β_{i}} (J_{e ℓ}^{- α_{i} β_{i}} - 1)] . \end{matrix}

Series expansion of the last two terms in terms of $J_{e ℓ} - 1$ shows that the first-order terms vanish and that the coefficients of the second-order terms are equal to $\frac{1}{2} α_{i}^{2} (\frac{1}{3} + β_{i}) {(J_{e ℓ} - 1)}^{2}$ . Hence, a stable material is obtained if $β_{i} > - \frac{1}{3}$ . For each term in the energy function, the coefficient $β_{i}$ determines the degree of compressibility. $β_{i}$ is related to the Poisson's ratio, $ν_{i}$ , by the expressions

β_{i} = \frac{ν_{i}}{1 - 2 ν_{i}}, ν_{i} = \frac{β_{i}}{1 + 2 β_{i}} .

Thus, if $β_{i}$ is the same for all terms, we have a single effective Poisson's ratio, $ν$ . The value of N and tables giving the $μ_{i}$ , $α_{i}$ , and $ν_{i}$ values as functions of temperature are specified by the user for the hyperfoam material model. The Poisson's ratio is valid for finite values of the logarithmic principal strains $e_{1}, e_{2}, e_{3}$ ; $e_{2} = e_{3} = - ν e_{1}$ in uniaxial tension. For $β_{i} = ν_{i} = 0$ there is no Poisson's effect. The initial shear modulus, $μ_{0}$ , again follows from

μ_{0} = \underset{i = 1}{\sum^{N}} μ_{i},

and the initial bulk modulus follows from

k_{0} = \underset{i = 1}{\sum^{N}} 2 μ_{i} (\frac{1}{3} + β_{i}) .

If Poisson's ratio is constant and known, Abaqus can calculate the $μ_{i}$ and $α_{i}$ from measurements of nominal stress and stretch as before. If Poisson's ratio depends on the level of straining, Abaqus can also calculate the $β_{i}$ from the nominal lateral strains.

Subroutine UHYPER

Abaqus/Standard also allows other forms of strain energy potentials to be defined for isotropic materials via user subroutine UHYPER by programming the first and second derivatives of U with respect to ${\bar{I}}_{1}$ , ${\bar{I}}_{2}$ , and J in that subroutine.