Fitting of hyperelastic and hyperfoam constants

Stress-strain relations for the polynomial strain energy potential

The hyperelastic polynomial form can be fitted by Abaqus up to order $N = 2$ . Since the Mooney-Rivlin potential corresponds to the case $N = 1$ , these remarks also apply by setting the higher-order coefficients to zero. The energy potential is as follows:

\begin{matrix} U = ​ & C_{10} ({\bar{I}}_{1} - 3) + C_{01} ({\bar{I}}_{2} - 3) + C_{20} {({\bar{I}}_{1} - 3)}^{2} + C_{11} ({\bar{I}}_{1} - 3) ({\bar{I}}_{2} - 3) + C_{02} {({\bar{I}}_{2} - 3)}^{2} \\ + \overset{2}{\sum_{i = 1}} \frac{1}{D_{i}} {(J_{e ℓ} - 1)}^{2 i} . \end{matrix}

The deformation modes are characterized in terms of the principal stretches. The nominal stress-strain relations are now derived for the polynomial form with $N = 2$ .

Uniaxial mode

λ_{1} = λ_{U}, λ_{2} = λ_{3} = λ_{U}^{- \frac{1}{2}}, λ_{U} = 1 + ϵ_{U}

The deviatoric strain invariants are

{\bar{I}}_{1} = λ_{U}^{2} + 2 λ_{U}^{- 1}, {\bar{I}}_{2} = λ_{U}^{- 2} + 2 λ_{U} .

We invoke the principle of virtual work to derive the nominal stress-strain relationship,

δ U = T_{U} δ λ_{U} = \frac{\partial U}{\partial λ_{U}} δ λ_{U},

and it follows that

\begin{matrix} T_{U} = ​ & \frac{\partial U}{\partial λ_{U}} = \frac{\partial U}{\partial {\bar{I}}_{1}} \frac{\partial {\bar{I}}_{1}}{\partial λ_{U}} + \frac{\partial U}{\partial {\bar{I}}_{2}} \frac{\partial {\bar{I}}_{2}}{\partial λ_{U}} \\ = & 2 (1 - λ_{U}^{- 3}) (λ_{U} \frac{\partial U}{\partial {\bar{I}}_{1}} + \frac{\partial U}{\partial {\bar{I}}_{2}}) \\ = & 2 (1 - λ_{U}^{- 3}) [C_{10} λ_{U} + C_{01} + 2 C_{20} λ_{U} ({\bar{I}}_{1} - 3) + C_{11} ({\bar{I}}_{1} - 3 + λ_{U} ({\bar{I}}_{2} - 3)) + 2 C_{02} ({\bar{I}}_{2} - 3)] . \end{matrix}

Equibiaxial mode

λ_{1} = λ_{2} = λ_{B}, λ_{3} = λ_{B}^{- 2}, λ_{B} = 1 + ϵ_{B}

The deviatoric strain invariants are

{\bar{I}}_{1} = 2 λ_{B}^{2} + λ_{B}^{- 4}, {\bar{I}}_{2} = 2 λ_{B}^{- 2} + λ_{B}^{4} .

From virtual work

δ U = 2 T_{B} δ λ_{B} = \frac{\partial U}{\partial λ_{B}} δ λ_{B},

and it follows that,

\begin{matrix} T_{B} = ​ & \frac{1}{2} \frac{\partial U}{\partial λ_{B}} = 2 (λ_{B} - λ_{B}^{- 5}) (\frac{\partial U}{\partial {\bar{I}}_{1}} + λ_{B}^{2} \frac{\partial U}{\partial {\bar{I}}_{2}}) \\ = & 2 (λ_{B} - λ_{B}^{- 5}) [C_{10} + C_{01} λ_{B}^{2} + 2 C_{20} ({\bar{I}}_{1} - 3) + C_{11} (λ_{B}^{2} ({\bar{I}}_{1} - 3) + ({\bar{I}}_{2} - 3)) \\ + 2 C_{02} λ_{B}^{2} ({\bar{I}}_{2} - 3)] . \end{matrix}

Planar (pure shear) mode

λ_{1} = λ_{S}, λ_{2} = 1, λ_{3} = λ_{S}^{- 1}, λ_{S} = 1 + ϵ_{S}

The deviatoric strain invariants are

{\bar{I}}_{1} = {\bar{I}}_{2} = λ_{S}^{2} + λ_{S}^{- 2} + 1.

From virtual work

δ U = T_{S} δ λ_{S} = \frac{\partial U}{\partial λ_{S}} δ λ_{S},

and it follows that,

\begin{matrix} T_{S} = ​ & \frac{\partial U}{\partial λ_{S}} = 2 (λ_{S} - λ_{S}^{- 3}) (\frac{\partial U}{\partial {\bar{I}}_{1}} + \frac{\partial U}{\partial {\bar{I}}_{2}}) \\ = & 2 (λ_{S} - λ_{S}^{- 3}) [C_{10} + C_{01} + 2 (C_{20} + C_{11} + C_{02}) ({\bar{I}}_{1} - 3)] . \end{matrix}

Volumetric mode

λ_{1} = λ_{2} = λ_{3} = λ_{V}, J = λ_{V}^{3}

From virtual work

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

and it follows that,

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

p = - \overset{N}{\sum_{i = 1}} 2 i \frac{1}{D_{i}} {(J - 1)}^{2 i - 1} .

Stress-strain relations for the reduced polynomial strain energy potential

The hyperelastic reduced polynomial form can be fitted by Abaqus up to order $N = 6$ . For $N = 3$ the reduced polynomial is identical to the Yeoh model, and for $N = 1$ the neo-Hookean model is retained; hence, the following also applies to these forms. The reduced polynomial energy potential is as follows:

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 the arguments in the previous section, we derive the nominal stress-strain relations for the reduced polynomial.

Uniaxial mode

λ_{1} = λ_{U}, λ_{2} = λ_{3} = λ_{U}^{- \frac{1}{2}}, λ_{U} = 1 + ϵ_{U}

T_{U} = 2 (λ_{U} - λ_{U}^{- 2}) \underset{i = 1}{\sum^{N}} i C_{i 0} {({\bar{I}}_{1} - 3)}^{i - 1} .

Equibiaxial mode

λ_{1} = λ_{2} = λ_{B}, λ_{3} = λ_{B}^{- 2}, λ_{B} = 1 + ϵ_{B}

T_{B} = 2 (λ_{B} - λ_{B}^{- 5}) \underset{i = 1}{\sum^{N}} i C_{i 0} {({\bar{I}}_{1} - 3)}^{i - 1} .

Planar (pure shear) mode

λ_{1} = λ_{S}, λ_{2} = 1, λ_{3} = λ_{S}^{- 1}, λ_{S} = 1 + ϵ_{S}

T_{S} = 2 (λ_{S} - λ_{S}^{- 3}) \underset{i = 1}{\sum^{N}} i C_{i 0} {({\bar{I}}_{1} - 3)}^{i - 1} .

Volumetric mode

λ_{1} = λ_{2} = λ_{3} = λ_{V}, J = λ_{V}^{3}

p = - \overset{N}{\sum_{i = 1}} 2 i \frac{1}{D_{i}} {(J - 1)}^{2 i - 1} .

Stress-strain relations for the hyperelastic Ogden strain energy potential

The hyperelastic Ogden form can be fitted up to order $N = 6$ :

U = \underset{i = 1}{\sum^{N}} \frac{2 μ_{i}}{α_{i}^{2}} ({\bar{λ}}_{1}^{α_{i}} + {\bar{λ}}_{2}^{α_{i}} + {\bar{λ}}_{3}^{α_{i}} - 3) + \overset{N}{\sum_{i = 1}} \frac{1}{D_{i}} {(J_{e ℓ} - 1)}^{2 i} .

Following the same approach as for the polynomial form, we can derive the nominal stress-strain equations for the Ogden form.

Uniaxial mode

λ_{1} = λ_{U}, λ_{2} = λ_{3} = λ_{U}^{- \frac{1}{2}}, λ_{U} = 1 + ϵ_{U}

T_{U} = \underset{i = 1}{\sum^{N}} \frac{2 μ_{i}}{α_{i}} (λ_{U}^{α_{i} - 1} - λ_{U}^{- \frac{1}{2} α_{i} - 1}) .

Equibiaxial mode

λ_{1} = λ_{2} = λ_{B}, λ_{3} = λ_{B}^{- 2}, λ_{B} = 1 + ϵ_{B}

T_{B} = \underset{i = 1}{\sum^{N}} \frac{2 μ_{i}}{α_{i}} (λ_{B}^{α_{i} - 1} - λ_{B}^{- 2 α_{i} - 1}) .

Planar (pure shear) mode

λ_{1} = λ_{S}, λ_{2} = 1, λ_{3} = λ_{S}^{- 1}, λ_{S} = 1 + ϵ_{S}

T_{S} = \underset{i = 1}{\sum^{N}} \frac{2 μ_{i}}{α_{i}} (λ_{S}^{α_{i} - 1} - λ_{S}^{- α_{i} - 1}) .

Volumetric mode

λ_{1} = λ_{2} = λ_{3} = λ_{V}, J = λ_{V}^{3}

p = - \overset{N}{\sum_{i = 1}} 2 i \frac{1}{D_{i}} {(J - 1)}^{2 i - 1} .

Stress-strain relations for the hyperelastic Arruda-Boyce strain energy potential

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} .

Following the same approach as for the polynomial form, we can derive the nominal stress-strain equations for the Arruda-Boyce potential.

Uniaxial mode

λ_{1} = λ_{U}, λ_{2} = λ_{3} = λ_{U}^{- \frac{1}{2}}, λ_{U} = 1 + ϵ_{U}

T_{U} = 2 μ (λ_{U} - λ_{U}^{- 2}) \sum_{i = 1}^{5} \frac{i C_{i}}{λ_{m}^{2 i - 2}} {\bar{I}}_{1}^{i - 1} .

Equibiaxial mode

λ_{1} = λ_{2} = λ_{B}, λ_{3} = λ_{B}^{- 2}, λ_{B} = 1 + ϵ_{B}

T_{B} = 2 μ (λ_{B} - λ_{B}^{- 5}) \underset{i = 1}{\sum^{5}} \frac{i C_{i}}{λ_{m}^{2 i - 2}} {\bar{I}}_{1}^{i - 1} .

Planar (pure shear) mode

λ_{1} = λ_{S}, λ_{2} = 1, λ_{3} = λ_{S}^{- 1}, λ_{S} = 1 + ϵ_{S}

T_{S} = 2 μ (λ_{S} - λ_{S}^{- 3}) \underset{i = 1}{\sum^{5}} \frac{i C_{i}}{λ_{m}^{2 i - 2}} {\bar{I}}_{1}^{i - 1} .

Volumetric mode

λ_{1} = λ_{2} = λ_{3} = λ_{V}, J = λ_{V}^{3}

p = - \frac{1}{D} (J - \frac{1}{J}) .

Stress-strain relations for the hyperelastic Van der Waals energy potential

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}} .

Following the same approach as for the polynomial form, we can derive the nominal stress-strain relations for the Van der Waals form.

Uniaxial mode

λ_{1} = λ_{U}, λ_{2} = λ_{3} = λ_{U}^{- \frac{1}{2}}, λ_{U} = 1 + ϵ_{U}

T_{U} = μ (1 - λ_{U}^{- 3}) (\frac{1}{1 - η} - a \sqrt{\frac{\tilde{I} - 3}{2}}) [λ_{U} (1 - β) + β] .

Equibiaxial mode

λ_{1} = λ_{2} = λ_{B}, λ_{3} = λ_{B}^{- 2}, λ_{B} = 1 + ϵ_{B}

T_{B} = μ (λ_{B} - λ_{B}^{- 5}) (\frac{1}{1 - η} - a \sqrt{\frac{\tilde{I} - 3}{2}}) (1 - β + β λ_{B}^{2}) .

Planar (pure shear) mode

λ_{1} = λ_{S}, λ_{2} = 1, λ_{3} = λ_{S}^{- 1}, λ_{S} = 1 + ϵ_{S}

T_{S} = μ (λ_{S} - λ_{S}^{- 3}) (\frac{1}{1 - η} - a \sqrt{\frac{\tilde{I} - 3}{2}}) .

Volumetric mode

λ_{1} = λ_{2} = λ_{3} = λ_{V}, J = λ_{V}^{3}

p = - \frac{1}{D} (J - \frac{1}{J}) .

Stress-strain relations for the hyperfoam strain energy potential

The hyperfoam potential is a modified form of the Hill strain energy potential and can be fitted up to order $N = 6$ :

U = \underset{i = 1}{\sum^{N}} \frac{2 μ_{i}}{α_{i}^{2}} [{\hat{λ}}_{1}^{α_{i}} + {\hat{λ}}_{2}^{α_{i}} + {\hat{λ}}_{3}^{α_{i}} - 3 + \frac{1}{β_{i}} (J_{e ℓ}^{- α_{i} β_{i}} - 1)] .

The deformation modes are characterized in terms of the principal stretches and the volume ratio J. The elastomeric foams are not incompressible: $J = λ_{1} λ_{2} λ_{3} \neq 1$ . The transverse stretches $λ_{2}$ and/or $λ_{3}$ are independently specified in the test data either as individual values depending on the lateral deformations or through the definition of an effective Poisson's ratio.

Uniaxial mode

λ_{1} = λ_{U}, λ_{2} = λ_{3}, J = λ_{U} λ_{2}^{2}, λ_{U} = 1 + ϵ_{U}

Equibiaxial mode

λ_{1} = λ_{2} = λ_{B}, J = λ_{B}^{2} λ_{3}^{2}, λ_{B} = 1 + ϵ_{B}

Planar mode

λ_{1} = λ_{P}, λ_{2} = 1, J = λ_{P} λ_{3}, λ_{P} = 1 + ϵ_{P}

The common nominal stress-strain relation for the three deformation modes above is

T_{j} = \frac{\partial U}{\partial λ_{j}} = \frac{2}{λ_{j}} \overset{N}{\sum_{i = 1}} \frac{μ_{i}}{α_{i}} (λ_{j}^{α_{i}} - J^{- α_{i} β_{i}}),

where $T_{j}$ is the nominal stress and $λ_{j}$ is the stretch in the direction of loading.

Simple shear mode

The simple shear deformation is described in terms of the deformation gradient,

F = [\begin{matrix} 1 & γ & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}],

where $γ$ is the shear strain. Note also that for this deformation, $J = \det (F) = 1$ . The nominal shear stress $T_{S}$ is

T_{S} = \frac{\partial U}{\partial γ} = \sum_{j = 1}^{2} {\frac{2 γ}{2 (λ_{j}^{2} - 1) - γ^{2}} \sum_{i = 1}^{N} \frac{μ_{i}}{α_{i}} (λ_{j}^{α_{i}} - 1)},

where $λ_{j}$ are the principal stretches in the plane of shearing, related to the shear strain $γ$ as follows:

λ_{1, 2} = \sqrt{1 + \frac{γ^{2}}{2} \pm γ \sqrt{1 + \frac{γ^{2}}{4}}} .

The stretch in the direction perpendicular to the plane of shearing is $λ_{3} = 1$ .

The transverse stress $T_{T}$ developed during simple shear deformation (as a result of the Poynting effect) is

T_{T} = \frac{\partial U}{\partial ε} = \overset{2}{\sum_{j = 1}} {\frac{2 (λ_{j}^{2} - 1)}{2 λ_{j}^{4} - λ_{j}^{2} (γ^{2} + 2)} \overset{N}{\sum_{i = 1}} \frac{μ_{i}}{α_{i}} (λ_{j}^{α_{i}} - 1)} .

Volumetric mode

The volumetric deformation is described as

λ_{1} = λ_{2} = λ_{3} = λ_{V}, J = λ_{V}^{3} .

The pressure p is related to volume ratio J through

p = - \frac{\partial U}{\partial J} = - \frac{2}{J} \overset{N}{\sum_{i = 1}} \frac{μ_{i}}{α_{i}} (J^{\frac{1}{3} α_{i}} - J^{- α_{i} β_{i}}) .

Least squares fit

Given experimental data, the material constants are determined through a least-squares-fit procedure, which minimizes the relative error in stress. For the n nominal stress–nominal strain data pairs, the relative error measure E is minimized,

E = \overset{n}{\sum_{i = 1}} {(1 - T_{i}^{th} / T_{i}^{test})}^{2},

where $T_{i}^{test}$ is a stress value from the test data and $T_{i}^{th}$ comes from one of the nominal stress expressions derived above.

The polynomial potential is linear in terms of the constants $C_{i j}$ ; therefore, a linear least-squares procedure can be used. The Ogden, the Arruda-Boyce, and the Van der Waals potential are nonlinear in some of their coefficients, thus necessitating the use of a nonlinear least-squares procedure.

Linear least squares fit for the polynomial model

For the full polynomial model we can rewrite the expressions for the $T_{k}^{th}$ derived above as

T_{k}^{th} = \overset{N}{\sum_{i + j = 1}} C_{i j} X_{i j} (λ_{k}), k = 1 \dots n,

where the $X_{i j} (λ_{k})$ are functions that depend on the stress state (uniaxial, biaxial, or planar), as explained above. $N = 1$ for the first-order polynomial (or Mooney-Rivlin form), and $N = 2$ for the second-order polynomial. To minimize the relative error, we need to set

\frac{\partial E}{\partial C_{i j}} = 0,

which leads to the following set of $M = \frac{1}{2} N (N + 3)$ equations:

\overset{n}{\sum_{k = 1}} \overset{N}{\sum_{i + j = 1}} \frac{X_{i j} (λ_{k}) X_{l m} (λ_{k})}{{(T_{k}^{test})}^{2}} C_{i j} = \overset{n}{\sum_{k = 1}} \frac{X_{l m} (λ_{k})}{T_{k}^{test}}, l + m = 1 \dots N .

This linear set of M equations can be solved readily to define the coefficients $C_{i j}$ .

To fit the volumetric coefficients, one needs to solve the system of N equations

\overset{n}{\sum_{k = 1}} \overset{N}{\sum_{i = 1}} \frac{X_{i} (J_{k}) X_{j} (J_{k})}{{(p_{k}^{test})}^{2}} \frac{1}{D_{i}} = \overset{n}{\sum_{k = 1}} \frac{X_{j} (J_{k})}{p_{k}^{test}}, j = 1 \dots N,

where

X_{i} (J_{k}) = - 2 i {(J_{k} - 1)}^{2 i - 1},

and

J_{k} = \frac{V_{k}}{V_{0}}

is given by the user. This system of equations can be solved readily for $D_{i}$ .

Linear least squares fit for the reduced polynomial model

For the reduced polynomial model we can rewrite the expressions for $T_{k}^{th}$ derived above as follows:

T_{k}^{th} = \overset{N}{\sum_{i = 1}} C_{i 0} X_{i} (λ_{k}),

where again the functions $X_{i} (λ_{k})$ depend on the stress state and the stretch, as outlined above, and N is the order of the reduced polynomial, which can take values up to $N = 6$ . The following also applies to the Yeoh and neo-Hookean forms since these models are special cases of the reduced polynomial, with $N = 3$ and $N = 1$ , respectively.

Following the same arguments as for the full polynomial, we arrive at the system of N equations

\overset{n}{\sum_{k = 1}} \overset{N}{\sum_{i = 1}} \frac{X_{i} (λ_{k}) X_{j} (λ_{k})}{{(T_{k}^{test})}^{2}} C_{i 0} = \overset{n}{\sum_{k = 1}} \frac{X_{j} (λ_{k})}{T_{k}^{test}}, j = 1 \dots N .

This system of equations can be solved readily for the coefficients $C_{i 0}$ . The volumetric coefficients are fitted with the same procedure as used for the general polynomial models.

Nonlinear least squares fit

The Ogden, Arruda-Boyce, and Van der Waals potentials are nonlinear in some of their coefficients; hence, a nonlinear least-squares-fit procedure is required. We use the Marquardt-Levenberg algorithm in the formulation by Twizell and Ogden (1983). Let $a_{i}$ , $i = 1 \dots m$ be the coefficients of these hyperelastic models, where m is the number of coefficients contributing to the deviatoric behavior. Specifically, $m = 2 N$ for the Ogden model, $m = 2$ for the Arruda-Boyce model, and $m = 4$ for the Van der Waals model. The coefficients are found by iterating the equation

a_{i}^{(r + 1)} = a_{i}^{(r)} - \overset{m}{\sum_{j = 1}} \overset{n}{\sum_{k = 1}} {[P_{i k}^{(r)} P_{j k}^{(r)} + γ δ_{i j}]}^{- 1} P_{j k}^{(r)} E_{k}^{(r)},

where r is the iteration count, n is the number of data points,

E_{k} = \frac{T_{k}^{test} - T_{k}^{th}}{T_{k}^{test}}

is the vector of relative errors, and

P_{i k} = \frac{\partial E_{k}}{\partial a_{i}} = - \frac{1}{T_{k}^{test}} \frac{\partial T_{k}^{th}}{\partial a_{i}}

is the derivative of the vector of relative errors with respect to the coefficients $a_{i}$ .

For $γ = 0$ , the Newton algorithm is obtained; for very large values of $γ$ , the steepest descent method is obtained. Thus, the Marquardt-Levenberg algorithm represents a compromise between these two approaches: the value of $γ$ is increased if the error grows and is reduced otherwise.

Nonlinear least squares fit for the Ogden model

After initializing the $α_{i}^{(0)}$ , the parameters $μ_{i}^{(0)}$ are found with a linear least squares fit. In the iterative procedure outlined above, the following derivatives are used:

\begin{matrix} \frac{\partial T_{k}^{th}}{\partial μ_{i}} = ​ & \frac{2}{α_{i}} (λ^{α_{i} - 1} - λ^{c α_{i} - 1}), \\ \frac{\partial T_{k}^{th}}{\partial α_{i}} = ​ & - \frac{2 μ_{i}}{α_{i}^{2}} (λ^{α_{i} - 1} - λ^{c α_{i} - 1}) - \frac{2 μ_{i}}{α_{i}} (λ^{α_{i} - 1} - c λ^{c α_{i} - 1}) \ln λ, \end{matrix}

where

c = {\begin{matrix} - \frac{1}{2}, & if
  uniaxial; \\ - 2, & if
  biaxial; \\ - 1, & if
  planar. \end{matrix}

Nonlinear least squares fit for the Arruda-Boyce model

The Arruda-Boyce model is linear in the shear modulus $μ$ but nonlinear in the locking stretch $λ_{m}$ . The locking stretch is initialized as $λ_{m}^{(0)} = max (7.0, 3.0 \times λ_{\max})$ , where $λ_{\max}$ is the maximum stretch in the user-specified test data. Given this locking stretch, the initial shear modulus, $μ^{(0)}$ , is obtained with a linear least squares fit.

In the iterative procedure outlined above, the following derivatives are used:

\frac{\partial T_{k}^{th}}{\partial μ} = {\begin{matrix} 2 (λ_{U} - λ_{U}^{- 2}) \underset{i = 1}{\sum^{5}} \frac{i C_{i}}{λ_{m}^{2 i - 2} {\bar{I}}_{1}^{i - 1}}, & if
  uniaxial; \\ 2 (λ_{B} - λ_{B}^{- 5}) \underset{i = 1}{\sum^{5}} \frac{i C_{i}}{λ_{m}^{2 i - 2} {\bar{I}}_{1}^{i - 1}}, & if
  biaxial; \\ 2 (λ_{S} - λ_{S}^{- 3}) \underset{i = 1}{\sum^{5}} \frac{i C_{i}}{λ_{m}^{2 i - 2} {\bar{I}}_{1}^{i - 1}}, & if
  planar. \end{matrix}

\frac{\partial T_{k}^{th}}{\partial λ_{m}} = {\begin{matrix} 2 μ (λ_{U} - λ_{U}^{- 2}) \underset{i = 2}{\sum^{5}} (2 - 2 i) \frac{i C_{i}}{λ_{m}^{2 i - 1} {\bar{I}}_{1}^{i - 1}}, & if
  uniaxial; \\ 2 μ (λ_{B} - λ_{B}^{- 5}) \underset{i = 2}{\sum^{5}} (2 - 2 i) \frac{i C_{i}}{λ_{m}^{2 i - 1} {\bar{I}}_{1}^{i - 1}}, & if
  biaxial; \\ 2 μ (λ_{S} - λ_{S}^{- 3}) \underset{i = 2}{\sum^{5}} (2 - 2 i) \frac{i C_{i}}{λ_{m}^{2 i - 1} {\bar{I}}_{1}^{i - 1}}, & if
  planar. \end{matrix}

Nonlinear least squares fit for the Van der Waals model

The Van der Waals model is linear in the shear modulus $μ$ but nonlinear in the locking stretch $λ_{m}$ , the global interaction parameter a, and the mixture parameter $β$ . The locking stretch is initialized as $λ_{m}^{(0)} = max (10.0, 3.0 \times λ_{\max})$ , where $λ_{\max}$ is the maximum stretch in the user-specified test data. Given this guess for the locking stretch, we make use of an expression proposed by Kilian et al. (1986) to initialize the global interaction parameter

a^{(0)} = \frac{{(λ_{m}^{(0)})}^{2}}{{(λ_{m}^{(0)})}^{3} - 1} .

The invariant mixture parameter is initialized to $β^{(0)} = 0$ . Given these initial values, the shear modulus, $μ^{(0)}$ , is initialized using a linear least-squares-fit procedure.

In the iterative procedure outlined above, the following derivatives are used:

\frac{\partial T_{k}^{th}}{\partial μ} = {\begin{matrix} (1 - λ_{U}^{- 3}) (\frac{1}{1 - η} - a \sqrt{\frac{\tilde{I} - 3}{2}}) [λ_{U} (1 - β) + β], & if
  uniaxial; \\ (λ_{B} - λ_{B}^{- 5}) (\frac{1}{1 - η} - a \sqrt{\frac{\tilde{I} - 3}{2}}) (1 - β + β λ_{B}^{2}), & if
  biaxial; \\ (λ_{S} - λ_{S}^{- 3}) (\frac{1}{1 - η} - a \sqrt{\frac{\tilde{I} - 3}{2}}), & if
  planar. \end{matrix}

\frac{\partial T_{k}^{th}}{\partial λ_{m}} = {\begin{matrix} - μ (1 - λ_{U}^{- 3}) \frac{η λ_{m}}{{(1 - η)}^{2} (λ_{m}^{2} - 3)} [λ_{U} (1 - β) + β], & if
  uniaxial; \\ - μ (λ_{B} - λ_{B}^{- 5}) \frac{η λ_{m}}{{(1 - η)}^{2} (λ_{m}^{2} - 3)} (1 - β + β λ_{B}^{2}), & if
  biaxial; \\ - μ (λ_{S} - λ_{S}^{- 3}) \frac{η λ_{m}}{{(1 - η)}^{2} (λ_{m}^{2} - 3)}, & if
  planar. \end{matrix}

\frac{\partial T_{k}^{th}}{\partial a} = {\begin{matrix} - μ (1 - λ_{U}^{- 3}) \sqrt{\frac{\tilde{I} - 3}{2}} [λ_{U} (1 - β) + β], & if
  uniaxial; \\ - μ (λ_{B} - λ_{B}^{- 5}) \sqrt{\frac{\tilde{I} - 3}{2}} (1 - β + β λ_{B}^{2}), & if
  biaxial; \\ - μ (λ_{S} - λ_{S}^{- 3}) \sqrt{\frac{\tilde{I} - 3}{2}}, & if
  planar. \end{matrix}

\frac{\partial T_{k}^{th}}{\partial β} = {\begin{matrix} 2 (1 - λ_{U}^{- 3}) [\frac{\partial U^{'}}{\partial β} (λ_{U} (1 - β) + β) + U^{'} (1 - λ_{U})], & if
  uniaxial; \\ 2 (λ_{B} - λ_{B}^{- 5}) [\frac{\partial U^{'}}{\partial β} (1 - β + β λ_{B}^{2}) + U^{'} (λ_{B}^{2} - 1)], & if
  biaxial; \\ 0, & if
  planar. \end{matrix}

In the derivatives of $\partial T_{k}^{th} / \partial β$

U^{'} = \frac{\partial U}{\partial \tilde{I}} = \frac{μ}{2} (\frac{1}{1 - η} - a \sqrt{\frac{\tilde{I} - 3}{2}}),

and

\frac{\partial U^{'}}{\partial β} = μ [\frac{1}{4} \frac{1}{y - 3} \frac{1}{η (1 - η^{2})} - \frac{a}{8} \sqrt{\frac{2}{\tilde{I} - 3}}] ({\bar{I}}_{1} - {\bar{I}}_{2}) .

In the planar case ${\bar{I}}_{1} = {\bar{I}}_{2}$ ; hence, $\partial T_{k}^{th} / \partial β$ vanishes.

Drucker stability check

Abaqus checks the Drucker stability of the material for the first three modes of deformation described above. The Drucker stability condition requires that the change in the Kirchhoff stress $d τ$ following from an infinitesimal change in the logarithmic strain $d ε$ satisfies the inequality

d τ : d ε > 0.

Using $d τ = D \cdot d ε$ , the inequality becomes

d ε \cdot D : d ε > 0,

thus requiring the tangential material stiffness $D$ to be positive definite for material stability to be satisfied.

For the isotropic elastic formulation considered here, the inequality can be represented in terms of the principal stresses and strains:

d τ_{1} d ε_{1} + d τ_{2} d ε_{2} + d τ_{3} d ε_{3} > 0.

Polynomial form

With the incompressibility assumption for the two hyperelastic models, the Kirchhoff stress is equal to the Cauchy stress: $τ = σ$ and, thus,

d σ : d ε > 0.

In addition, we can choose any value for the hydrostatic pressure without affecting the strains. For the stability calculation a convenient choice is $σ_{3} = d σ_{3} = 0$ , which gives us

d σ_{1} d ε_{1} + d σ_{2} d ε_{2} > 0.

The infinitesimal strain changes are related to the changes in stretch ratios by the equations

d ε_{1} = \frac{d λ_{1}}{λ_{1}}, d ε_{2} = \frac{d λ_{2}}{λ_{2}} .

The stresses follow from the strain energy, which in turn follow from the changes in the strain invariants or in the stretches.

The relation between changes in the stress and changes in strain are described by the matrix equation

{\begin{matrix} d σ_{1} \\ d σ_{2} \end{matrix}} = [\begin{matrix} D_{11} & D_{12} \\ D_{21} & D_{22} \end{matrix}] {\begin{matrix} d ε_{1} \\ d ε_{2} \end{matrix}},

where

D_{11} = 4 (λ_{1}^{2} + λ_{3}^{2}) (\frac{\partial U}{\partial {\bar{I}}_{1}} + λ_{2}^{2} \frac{\partial U}{\partial {\bar{I}}_{2}}) + 4 {(λ_{1}^{2} - λ_{3}^{2})}^{2} (\frac{\partial^{2} U}{\partial {\bar{I}}_{1}^{2}} + 2 λ_{2}^{2} \frac{\partial^{2} U}{\partial {\bar{I}}_{1} \partial {\bar{I}}_{2}} + λ_{2}^{4} \frac{\partial^{2} U}{\partial {\bar{I}}_{2}^{2}}),

D_{22} = 4 (λ_{2}^{2} + λ_{3}^{2}) (\frac{\partial U}{\partial {\bar{I}}_{1}} + λ_{1}^{2} \frac{\partial U}{\partial {\bar{I}}_{2}}) + 4 {(λ_{2}^{2} - λ_{3}^{2})}^{2} (\frac{\partial^{2} U}{\partial {\bar{I}}_{1}^{2}} + 2 λ_{1}^{2} \frac{\partial^{2} U}{\partial {\bar{I}}_{1} \partial {\bar{I}}_{2}} + λ_{1}^{4} \frac{\partial^{2} U}{\partial {\bar{I}}_{2}^{2}}),

D_{12} = D_{21} = 4 λ_{3}^{2} \frac{\partial U}{\partial {\bar{I}}_{1}} + 4 λ_{3}^{- 2} \frac{\partial U}{\partial {\bar{I}}_{2}} + 4 (λ_{1}^{2} - λ_{3}^{2}) (λ_{2}^{2} - λ_{3}^{2}) (\frac{\partial^{2} U}{\partial {\bar{I}}_{1}^{2}} + (λ_{1}^{2} + λ_{2}^{2}) \frac{\partial^{2} U}{\partial {\bar{I}}_{1} \partial {\bar{I}}_{2}} + λ_{1}^{2} λ_{2}^{2} \frac{\partial^{2} U}{\partial {\bar{I}}_{2}^{2}}) .

For material stability $D$ must be positive definite; thus, it is necessary that

D_{11} + D_{22} > 0,

D_{11} D_{22} - D_{12} D_{21} > 0,

for all relevant values of $λ_{1}$ , $λ_{2}$ , and $λ_{3}$ .

Ogden form

For the Ogden form we follow the same approach as the polynomial form. Using $λ_{3} = λ_{1}^{- 1} λ_{2}^{- 1}$ , we have

U = \underset{i = 1}{\sum^{N}} \frac{2 μ_{i}}{α_{i}^{2}} (λ_{1}^{α_{i}} + λ_{2}^{α_{i}} + λ_{1}^{- α_{i}} λ_{2}^{- α_{i}} - 3),

σ_{1} = λ_{1} \frac{\partial U}{\partial λ_{1}} = \underset{i = 1}{\sum^{N}} \frac{2 μ_{i}}{α_{i}} (λ_{1}^{α_{i}} - λ_{1}^{- α_{i}} λ_{2}^{- α_{i}}),

σ_{2} = λ_{2} \frac{\partial U}{\partial λ_{2}} = \underset{i = 1}{\sum^{N}} \frac{2 μ_{i}}{α_{i}} (λ_{2}^{α_{i}} - λ_{1}^{- α_{i}} λ_{2}^{- α_{i}}),

and the material stiffness $D$ that we check for positive definiteness is

[\begin{matrix} D_{11} & D_{12} \\ D_{21} & D_{22} \end{matrix}] = \underset{i = 1}{\sum^{N}} 2 μ_{i} λ_{1}^{- α_{i}} λ_{2}^{- α_{i}} [\begin{matrix} λ_{1}^{2 α_{i}} λ_{2}^{α_{i}} + 1 & 1 \\ 1 & λ_{1}^{α_{i}} λ_{2}^{2 α_{i}} + 1 \end{matrix}] .

Arruda-Boyce form

For positive values of the shear modulus, $μ$ , and the locking stretch, $λ_{m}$ , the Arruda-Boyce form is always stable. Hence, it suffices to check the coefficients to determine whether the material satisfies Drucker stability.

Van der Waals form

When the Van der Waals model is employed in its admissible stretch range given by $\tilde{I} < λ_{m}^{2}$ , its stability depends on the global interaction parameter, a, for positive values of the initial shear modulus, $μ$ , and the locking stretch, $λ_{m}$ . To verify the Drucker stability of the Van der Waals model, we can employ the equations derived for the polynomial models by making use of the fact that

\frac{\partial U}{\partial {\bar{I}}_{1}} = (1 - β) \frac{\partial U}{\partial \tilde{I}} and \frac{\partial U}{\partial {\bar{I}}_{2}} = β \frac{\partial U}{\partial \tilde{I}} .

To determine the admissible stretch range, we need to find the two positive real-valued roots neighboring $λ = 1$ of the equation

\tilde{I} (λ) - λ_{m}^{2} = 0

for each of the three stress states—uniaxial, biaxial, and planar—by using a simple bisection method.

Hyperfoam

The Kirchhoff stress-strain relation for the uniaxial, biaxial, planar, and volumetric deformation modes is

τ_{j} = λ_{j} \frac{\partial U}{\partial λ_{j}} = 2 \overset{N}{\sum_{i = 1}} \frac{μ_{i}}{α_{i}} [λ_{j}^{α_{i}} - J^{- α_{i} β_{i}}] .

Taking the total differential of $τ_{j}$ and using $d ε_{j} = d λ_{j} / λ_{j}$ ,

d τ_{j} = 2 \overset{N}{\sum_{i = 1}} μ_{i} [λ_{j}^{α_{i}} d ε_{j} + β_{i} J^{- α_{i} β_{i}} (d ε_{1} + d ε_{2} + d ε_{3})], j = 1, 2, 3.

Since we cannot use the incompressibility assumption, we have to use all three principal stress and strain components and a $3 \times 3$ $D$ matrix,

{\begin{matrix} d τ_{1} \\ d τ_{2} \\ d τ_{3} \end{matrix}} = [\begin{matrix} D_{11} & D_{12} & D_{13} \\ D_{21} & D_{22} & D_{23} \\ D_{31} & D_{32} & D_{33} \end{matrix}] {\begin{matrix} d ε_{1} \\ d ε_{2} \\ d ε_{3} \end{matrix}} .

Specifically,

{\begin{matrix} d τ_{1} \\ d τ_{2} \\ d τ_{3} \end{matrix}} = 2 \underset{i = 1}{\sum^{N}} μ_{i} [\begin{matrix} λ_{1}^{α_{i}} + A_{i} & A_{i} & A_{i} \\ A_{i} & λ_{2}^{α_{i}} + A_{i} & A_{i} \\ A_{i} & A_{i} & λ_{3}^{α_{i}} + A_{i} \end{matrix}] {\begin{matrix} d ε_{1} \\ d ε_{2} \\ d ε_{3} \end{matrix}},

where $A_{i} = β_{i} J^{- α_{i} β_{i}}$ ,

For the simple shear case the principal stretches $λ_{1}$ and $λ_{2}$ are computed from the shear strain $γ$ (as given in an earlier expression). Thus, the same form of equations is used in checking material stability during simple shear deformation.

For material stability (i.e., for $D$ to be positive definite) the following conditions must be satisfied:

D_{11} + D_{22} + D_{33} > 0,

D_{11} D_{22} + D_{22} D_{33} + D_{33} D_{11} - D_{23}^{2} - D_{13}^{2} - D_{12}^{2} > 0,

\det (D) > 0.