A material is said to be hyperelastic if there exists an elastic potential function W (or strain-energy density function) which is a scalar function of one of the strain or deformation tensors, whose derivative with respect to a strain component determines the corresponding stress component. This can be expressed by:

(4–166)

where:

The Lagrangian strain can be expressed as follows:

(4–167)

where:

The nominal strain in the current configuration is defined by:

(4–168)

The deformation tensor C_ij is composed of the products of the deformation gradients F_ij:

(4–169)

where:

and the right stretch tensor is .

The Kirchhoff stress is defined:

(4–170)

and the Cauchy stress is obtained by:

(4–171)

The eigenvalues (principal stretch ratios) of C_ij are , , and , and exist only if:

(4–172)

which can be re-expressed as:

(4–173)

where:

(4–174)

and

(4–175)

J is also the ratio of the deformed elastic volume over the reference (undeformed) volume of materials (Ogden([296]) and Crisfield([295])).

When there is thermal volume strain, the volume ratio J is replaced by the elastic volume ratio J_el which is defined as the total volume ratio J over thermal volume ratio J_th, as:

(4–176)

and the thermal volume ratio J_th is:

(4–177)

where:

Under the assumption that material response is isotropic, it is convenient to express the strain-energy function in terms of strain invariants or principal stretches (Simo and Hughes([253])).

(4–178)

or

(4–179)

Define the volume-preserving part of the deformation gradient, , as:

(4–180)

and thus

(4–181)

The modified principal stretch ratios and invariants are then:

(4–182)

(4–183)

The strain-energy potential can then be defined as:

(4–184)

Following are several forms of strain-energy potential (W) provided (as options TBOPT in TB,HYPER) for the simulation of incompressible or nearly incompressible hyperelastic materials.

4.6.3.1.  Arruda-Boyce Model

The form of the strain-energy potential for Arruda-Boyce model is:

(4–185)

where:

μ = initial shear modulus of material (input on TBDATA commands with TB,HYPER)

λL = limiting network stretch (input on TBDATA commands with TB,HYPER)

d = material incompressibility parameter (input on TBDATA commands with TB,HYPER)

The initial bulk modulus is:

(4–186)

As the parameter λ_L goes to infinity, the model is converted to Neo-Hookean form.

4.6.3.2.  Blatz-Ko Model

The form of strain-energy potential for the Blatz-Ko model is:

(4–187)

where:

μ = initial shear modulus of material (input on TBDATA commands with TB,HYPER)

The initial bulk modulus is defined as:

(4–188)

4.6.3.3.  Extended Tube Model

The extended-tube model is a physics-based polymer model which introduces the physical consideration on the molecular scale into the formulation of the strain-energy potential. The model considers the topological constraints as well as the limited chain extensibility of network chains in the filled rubbers.

The elastic strain-energy potential consists of the elastic energy from the crosslinked network and the constraint network as well as the volumetric strain energy:

(4–189)

where the initial shear modulus is given by G = G_c + G_e, and:

Ge = constraint contribution to modulus

Gc = crosslinked contribution to modulus

δ = extensibility parameter

β = empirical parameter (0 ≤ β ≤1)

d1 = material incompressibility parameter

The model is equivalent to a two-term Ogden model with the following parameters:

4.6.3.4.  Gent Model

The form of the strain-energy potential for the Gent model is:

(4–190)

where:

μ = initial shear modulus of material (input on TBDATA commands with TB,HYPER)

Jm = limiting value of (input on TBDATA commands with TB,HYPER)

d = material incompressibility parameter (input on TBDATA commands with TB,HYPER)

The initial bulk modulus is:

(4–191)

As the parameter J_m goes to infinity, the model is converted to Neo-Hookean form.

4.6.3.5.  Mooney-Rivlin

This option includes two-, three-, five-, and nine-term Mooney-Rivlin models. The form of the strain-energy potential for a two-parameter Mooney-Rivlin model is:

(4–192)

where:

c10, c01, d = material constants (input on TBDATA commands with TB,HYPER)

The form of the strain-energy potential for a three-parameter Mooney-Rivlin model is

(4–193)

where:

c10, c01, c11, d = material constants (input on TBDATA commands with TB,HYPER)

The form of the strain-energy potential for five-parameter Mooney-Rivlin model is:

(4–194)

where:

c10, c01, c20, c11, c02, d = material constants (input on TBDATA commands with TB,HYPER)

The form of the strain-energy potential for nine-parameter Mooney-Rivlin model is:

(4–195)

where:

c10, c01, c20, c11, c02, c30, c21, c12, c03, d = material constants (input on TBDATA commands with TB,HYPER)

The initial shear modulus is given by:

(4–196)

The initial bulk modulus is:

(4–197)

4.6.3.6.  Neo-Hookean

The form Neo-Hookean strain-energy potential is:

(4–198)

where:

μ = initial shear modulus of materials (input on TBDATA commands with TB,HYPER)

d = material incompressibility parameter (input on TBDATA commands with TB,HYPER)

The initial bulk modulus is related to the material incompressibility parameter by:

(4–199)

where:

K = initial bulk modulus

4.6.3.7.  Ogden Compressible Foam Model

The strain-energy potential of the Ogden compressible foam model is based on the principal stretches of left-Cauchy strain tensor, which has the form:

(4–200)

where:

N = material constant (input as NPTS on TB,HYPER)

μi, α i, βi = material constants (input on TBDATA commands with TB,HYPER)

The initial shear modulus, μ, is given as:

(4–201)

The initial bulk modulus K is defined by:

(4–202)

For N = 1, α ₁ = -2, μ₁= –μ, and β = 0.5, the Ogden option is equivalent to the Blatz-Ko option.

4.6.3.8.  Ogden Potential

The Ogden form of strain-energy potential is based on the principal stretches of left-Cauchy strain tensor, which has the form:

(4–203)

where:

N = material constant (input as NPTS on TB,HYPER)

μi, α i, dk = material constants (input on TBDATA commands with TB,HYPER)

Similar to the Polynomial form, there is no limitation on N. A higher N can provide better fit the exact solution, however, it may, on the other hand, cause numerical difficulty in fitting the material constants and also it requests to have enough data to cover the entire range of interest of the deformation. Therefore a value of N > 3 is not usually recommended.

The initial shear modulus, μ, is given as:

(4–204)

The initial bulk modulus is:

(4–205)

For N = 1 and α ₁ = 2, the Ogden potential is equivalent to the Neo-Hookean potential. For N = 2, α ₁ = 2 and α ₂ = -2, the Ogden potential can be converted to the 2 parameter Mooney-Rivlin model.

4.6.3.9.  Polynomial Form

The polynomial form of strain-energy potential is

(4–206)

where:

N = material constant (input as NPTS on TB,HYPER)

cij, dk = material constants (input on TBDATA commands with TB,HYPER)

In general, there is no limitation on N. A higher N may provide better fit the exact solution, however, it may, on the other hand, cause numerical difficulty in fitting the material constants and requires enough data to cover the entire range of interest of deformation. Therefore a very higher N value is not usually recommended.

The Neo-Hookean model can be obtained by setting N = 1 and c₀₁ = 0. Also for N = 1, the two parameters Mooney-Rivlin model is obtained, for N = 2, the five parameters Mooney-Rivlin model is obtained and for N = 3, the nine parameters Mooney-Rivlin model is obtained.

The initial shear modulus is defined:

(4–207)

The initial bulk modulus is:

(4–208)

4.6.3.10.  Yeoh Model

The Yeoh model is also called the reduced polynomial form. The strain-energy potential is:

(4–209)

where:

N = material constant (input as NPTS on TB,HYPER)

Ci0 = material constants (input on TBDATA commands with TB,HYPER)

dk = material constants (input on TBDATA commands with TB,HYPER)

The Neo-Hookean model can be obtained by setting N = 1.

The initial shear modulus is defined:

(4–210)

The initial bulk modulus is:

(4–211)

The anisotropic constitutive strain-energy density function W is:

(4–212)

where:

It is assumed that the material is nearly incompressible or purely incompressible. The volumetric part W_v is absolutely independent of the isochoric part W_d.

The volumetric part, W_v, is assumed to be only function of J as:

(4–213)

The isochoric part W_d is a function of the invariants of the isochoric part of the right Cauchy Green tensor and the two constitutive material directions A, B in the undeformed configuration. The material directions yield so-called structural tensors of the microstructure of the material. Thus, the strain-energy density yields:

(4–214)

where:

Two strain energy potentials, as forms of polynomial or exponential functions, are available for characterizing the isochoric part of the strain-energy potential.

The third invariant is ignored due to the incompressible assumption. The parameter ς is defined as:

(4–215)

In Equation 4–214 the irreducible basis of invariants:

(4–216)

The exponential-function-based strain energy potential function is given as:

(4–217)

For information about the splitting of the volumetric and deviatoric parts of the strain-density function, see Finite Strain Elasticity. For further information about this material model, see Anisotropic Hyperelasticity in the Mechanical APDL Material Reference.

The option of user subroutine allows users to define their own strain-energy potential. Use of subroutine userhyper is necessary to provide the derivatives of the strain-energy potential with respect to the strain invariants. See the Guide to ANSYS User Programmable Features for more information about writing a user hyperelasticity subroutine.

Stresses (output quantities S) are true (Cauchy) stresses in the global coordinate system. They are computed from the second Piola-Kirchhoff stresses using:

(4–218)

where:

Strains (output as EPEL) are the Hencky (logarithmic) strains (see Equation 3–6). They are in the global coordinate system. Thermal strain (output as EPTH) is reported as:

(4–219)

The hyperelastic constants in the strain-energy density function of a material determine its mechanical response. Therefore, in order to obtain successful results during a hyperelastic analysis, it is necessary to accurately assess the material constants of the materials being examined. Material constants are generally derived for a material using experimental stress-strain data. It is recommended that this test data be taken from several modes of deformation over a wide range of strain values. In fact, it has been observed that to achieve stability, the material constants should be fit using test data in at least as many deformation states as will be experienced in the analysis.

For hyperelastic materials, simple deformation tests (consisting of six deformation modes) can be used to accurately characterize the material constants. (See Material Curve Fitting in the Mechanical APDL Material Reference.) All available laboratory test data will be used to determine the hyperelastic material constants. The six different deformation modes are graphically illustrated in Figure 4.25: Illustration of Deformation Modes. Combinations of data from multiple tests will enhance the characterization of the hyperelastic behavior of a material.

Figure 4.25:  Illustration of Deformation Modes

Although the algorithm accepts up to six different deformation states, it can be shown that apparently different loading conditions have identical deformations, and are thus equivalent. Superposition of tensile or compressive hydrostatic stresses on a loaded incompressible body results in different stresses, but does not alter deformation of a material. As depicted in Figure 4.26: Equivalent Deformation Modes, we find that upon the addition of hydrostatic stresses, the following modes of deformation are identical:

With several equivalent modes of testing, we are left with only three independent deformation states for which one can obtain experimental data.

Figure 4.26:  Equivalent Deformation Modes

The following sections outline the development of hyperelastic stress relationships for each independent testing mode. In the analyses, the coordinate system is chosen to coincide with the principal directions of deformation. Thus, the right Cauchy-Green strain tensor can be written in matrix form by:

(4–220)

where:

λi = 1 + εi principal stretch ratio in the ith direction

εi = principal value of the engineering strain tensor in the ith direction

The principal invariants of C_ij are:

(4–221)

(4–222)

(4–223)

For each mode of deformation, fully incompressible material behavior is also assumed so that third principal invariant, I₃, is identically one:

(4–224)

Finally, the hyperelastic Piola-Kirchhoff stress tensor, Equation 4–166 can be algebraically manipulated to determine components of the Cauchy (true) stress tensor. In terms of the left Cauchy-Green strain tensor, the Cauchy stress components for a volumetrically constrained material can be shown to be:

(4–225)

where:

p = pressure

4.6.7.1. Uniaxial Tension (Equivalently, Equibiaxial Compression)

As shown in Figure 4.25: Illustration of Deformation Modes, a hyperelastic specimen is loaded along one of its axis during a uniaxial tension test. For this deformation state, the principal stretch ratios in the directions orthogonal to the 'pulling' axis will be identical. Therefore, during uniaxial tension, the principal stretches, λ_i, are given by:

(4–226)

(4–227)

Due to incompressibility Equation 4–224:

(4–228)

and with Equation 4–227,

(4–229)

For uniaxial tension, the first and second strain invariants then become:

(4–230)

and

(4–231)

Substituting the uniaxial tension principal stretch ratio values into the Equation 4–225, we obtain the following stresses in the 1 and 2 directions:

(4–232)

and

(4–233)

Subtracting Equation 4–233 from Equation 4–232, we obtain the principal true stress for uniaxial tension:

(4–234)

The corresponding engineering stress is:

(4–235)

4.6.7.2. Equibiaxial Tension (Equivalently, Uniaxial Compression)

During an equibiaxial tension test, a hyperelastic specimen is equally loaded along two of its axes, as shown in Figure 4.25: Illustration of Deformation Modes. For this case, the principal stretch ratios in the directions being loaded are identical. Hence, for equibiaxial tension, the principal stretches, λ_i, are given by:

(4–236)

(4–237)

Utilizing incompressibility Equation 4–224, we find:

(4–238)

For equibiaxial tension, the first and second strain invariants then become:

(4–239)

and

(4–240)

Substituting the principal stretch ratio values for equibiaxial tension into the Cauchy stress Equation 4–225, we obtain the stresses in the 1 and 3 directions:

(4–241)

and

(4–242)

Subtracting Equation 4–242 from Equation 4–241, we obtain the principal true stress for equibiaxial tension:

(4–243)

The corresponding engineering stress is:

(4–244)

4.6.7.3. Pure Shear

(Uniaxial Tension and Uniaxial Compression in Orthogonal Directions)

Pure shear deformation experiments on hyperelastic materials are generally performed by loading thin, short and wide rectangular specimens, as shown in Figure 4.27: Pure Shear from Direct Components. For pure shear, plane strain is generally assumed so that there is no deformation in the 'wide' direction of the specimen: λ₂ = 1.

Figure 4.27:  Pure Shear from Direct Components

Due to incompressibility Equation 4–224, it is found that:

(4–245)

For pure shear, the first and second strain invariants are:

(4–246)

and

(4–247)

Substituting the principal stretch ratio values for pure shear into the Cauchy stress Equation 4–225, we obtain the following stresses in the 1 and 3 directions:

(4–248)

and

(4–249)

Subtracting Equation 4–249 from Equation 4–248, we obtain the principal pure shear true stress equation:

(4–250)

The corresponding engineering stress is:

(4–251)

4.6.7.4. Volumetric Deformation

The volumetric deformation is described as:

(4–252)

As nearly incompressible is assumed, we have:

(4–253)

The pressure, P, is directly related to the volume ratio J through:

(4–254)

4.6.7.5. Least Squares Fit Analysis

By performing a least squares fit analysis the Mooney-Rivlin constants can be determined from experimental stress-strain data and Equation 4–233, Equation 4–243, and Equation 4–250. Briefly, the least squares fit minimizes the sum of squared error between experimental and Cauchy predicted stress values. The sum of the squared error is defined by:

(4–255)

where:

E = least squares residual error

Ti(Cj) = engineering stress values (function of hyperelastic material constants

n = number of experimental data points

Equation 4–255 is minimized by setting the variation of the squared error to zero: δ E² = 0. This yields a set of simultaneous equations which can be used to solve for the hyperelastic constants:

(4–256)

It should be noted that for the pure shear case, the hyperelastic constants cannot be uniquely determined from Equation 4–250. In this case, the shear data must by supplemented by either or both of the other two types of test data to determine the constants.

From Equation 4–225, the deviatoric stress is determined solely by the deformation and the derivatives of the elastic potential function , , and . These derivatives are called the response functions and are determined analytically for the hyperelastic potentials in Isotropic Hyperelasticity and Anisotropic Hyperelasticity. These response functions can also be determined directly from experimental data, bypassing the need to fit the potential function's parameters to the experimental data.

The constitutive Equation 4–234, Equation 4–243, and Equation 4–250 give an explicit relationship between the stress, deformation, and response functions. The experimental data consists of the measured deformation and stress so that the only unknowns in the constitutive equations are the response functions.

The volumetric response function ( ) is determined either analytically (from the polynomial volumetric potential function given as the second term on the right side of Equation 4–206) or as experimental data of volume ratio versus pressure. Given the experimental data, the volume ratio for a general deformation is used to determine the pressure and hence the volumetric response function.

The response functions for the first and second deformation invariant are determined from the experimental data from uniaxial tension, equibiaxial tension, pure shear or combined uniaxial tension and compression experiments. Additionally, for incompressible materials, uniaxial compression experiments are equivalent to equibiaxial tension and can be used in place of equibiaxial data to determine the response functions.

Combined uniaxial tension plus compression data cannot be combined with other data sets (except pressure-volume), and gives only a material behavior that depends on the first invariant response function. Combinations of the other experimental data sets can be used to determine the first and second invariant response functions. Any combination of uniaxial tension, equibiaxial tension or pure shear data can be used.

For a single set of experimental data, the respective constitutive equation has two unknown response functions. Setting the second invariant response function ( ) to zero, the constitutive equation can be solved for the first invariant response function ( ).

Given two sets of experimental data for independent deformations (uniaxial-equibiaxial, uniaxial-pure shear, or equibiaxial-pure shear), the resulting two equations can be solved for the two unknown response functions. With experimental data for uniaxial tension, equibiaxial tension, and pure shear deformations, the resulting three constitutive equations can be solved for a best fit of the two response functions.

The response functions can therefore be determined from experimental data over the range of experimental deformation; however, the constitutive response is needed for any general deformation. The constitutive response is obtained by determining an experimental deformation that is nearest to the general deformation, where the definition of nearest is that the general deformation and the experimental deformation have the same first invariant I₁.

For example, given any arbitrary deformation state, Equation 4–221 gives the equation for the first invariant. Equating this equation to the uniaxial deformation first invariant of Equation 4–230 gives

(4–257)

where:

Solving the resulting equation for the tensile value of λ gives the uniaxial deformation that is nearest to the general deformation. Then, λ is used to determine the response function from the uniaxial experimental data, which corresponds to the required response function for the general deformation.

In general, Equation 4–257 has two valid solutions, one in tension and one in compression. For use with the uniaxial tension experimental data, the tensile solution is used. For use with the combined uniaxial tension and compression data, the tension or compression λ value is chosen so that the difference between the uniaxial deformation gradient's eigenvalues are closest to the eigenvalues of the actual deformation gradient.

The nearest equibiaxial and pure shear deformations are determined in a similar manner by equating the first invariant for the general deformation with the first invariant for the equibiaxial or pure shear deformation. Then the experimental stretch is determined by choosing the tensile value from the resulting valid solutions for λ.

Stability checks are provided for the Mooney-Rivlin hyperelastic materials. A nonlinear material is stable if the secondary work required for an arbitrary change in the deformation is always positive. Mathematically, this is equivalent to:

(4–258)

where:

Since the change in stress is related to the change in strain through the material stiffness tensor, checking for stability of a material can be more conveniently accomplished by checking for the positive definiteness of the material stiffness.

The material stability checks are done at the end of preprocessing but before an analysis actually begins. At that time, the program checks for the loss of stability for six typical stress paths (uniaxial tension and compression, equibiaxial tension and compression, and planar tension and compression). The range of the stretch ratio over which the stability is checked is chosen from 0.1 to 10. If the material is stable over the range then no message will appear. Otherwise, a warning message appears that lists the Mooney-Rivlin constants and the critical values of the nominal strains where the material first becomes unstable.