Fitting of hyperelastic and hyperfoam constants

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:

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

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:

Following the arguments in the previous section, we derive the nominal stress-strain relations for the reduced polynomial.

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

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

The hyperelastic Arruda-Boyce potential has the following form:

where

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

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

where

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

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

The deformation modes are characterized in terms of the principal stretches and the volume ratio J. The elastomeric foams are not incompressible: $J={\lambda }_1{\lambda }_2{\lambda }_3\ne 1$. The transverse stretches ${\lambda }_2$ and/or ${\lambda }_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.

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,

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

The polynomial potential is linear in terms of the constants $C_{ij}$; 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.

For the full polynomial model we can rewrite the expressions for the $T_k^{\mathrm{th}}$ derived above as

where the $X_{ij}\left({\lambda }_k\right)$ 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

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

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

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

where

and

is given by the user. This system of equations can be solved readily for $D_i$.

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

where again the functions $X_i\left({\lambda }_k\right)$ 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

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

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\mathrm{\dots }m$ be the coefficients of these hyperelastic models, where m is the number of coefficients contributing to the deviatoric behavior. Specifically, $m=2N$ 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

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

is the vector of relative errors, and

is the derivative of the vector of relative errors with respect to the coefficients $a_i$.

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

After initializing the ${\alpha }_i^{\left(0\right)}$, the parameters ${\mu }_i^{\left(0\right)}$ are found with a linear least squares fit. In the iterative procedure outlined above, the following derivatives are used:

where

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

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

The Van der Waals model is linear in the shear modulus $\mu$ but nonlinear in the locking stretch ${\lambda }_m$, the global interaction parameter a, and the mixture parameter $\beta$. The locking stretch is initialized as ${\lambda }_m^{\left(0\right)}=max\left(10.0,3.0\times {\lambda }_{\max }\right)$, where ${\lambda }_{\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

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

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

In the derivatives of $\partial T_k^{\mathrm{th}}/\partial \beta$

and

In the planar case ${\overline{I}}_1={\overline{I}}_2$; hence, $\partial T_k^{\mathrm{th}}/\partial \beta$ vanishes.

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\mathit{\tau }$ following from an infinitesimal change in the logarithmic strain $d\mathit{\epsilon }$ satisfies the inequality

Using $d\mathit{\tau }=\mathbf{D}\cdot d\mathit{\epsilon }$, the inequality becomes

thus requiring the tangential material stiffness $\mathbf{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: