Hyperelasticity refers to a constitutive response that is derivable from an elastic free energy potential and is typically used for materials which experience large elastic deformation. Applications for elastomers such as vulcanized rubber and synthetic polymers, along with some biological materials, often fall into this category.
The microstructure of polymer solids consists of chain-like molecules. The flexibility of these molecules allows for an irregular molecular arrangement and, as a result, the behavior is very complex. Polymers are usually isotropic at small deformation and anisotropic at larger deformation as the molecule chains realign to the loading direction. Under an essentially monotonic loading condition, however, many polymer materials can be approximated as isotropic, which has historically been popular in the modeling of polymers.
Some classes of hyperelastic materials cannot be modeled as isotropic. An example is fiber reinforced polymer composites. Typical fiber patterns include unidirectional and bidirectional, and the fibers can have a stiffness that is 50-1000 times that of the polymer matrix, resulting in a strongly anisotropic material behavior. Another class of anisotropic materials that can experience large deformation is biomaterials, such as muscles and arteries, in which the anisotropic behavior is due to their fibrous structure.
The typical volumetric behavior of hyperelastic materials can be grouped into two classes. Materials such as polymers typically have small volumetric changes during deformation and these are incompressible or nearly-incompressible materials. An example of the second class of materials is foams, which can experience large volumetric changes during deformation, and these are compressible materials.
The available hyperelastic material constitutive models are derived from strain-energy potentials that are functions of the deformation invariants. An exception is the response function model which obtains the constitutive response functions directly from experimental data. The hyperelastic material models are defined through data tables (TB,HYPER or TB,AHYPER) and include:
Incompressible or nearly-incompressible isotropic models: Neo-Hookean, Mooney-Rivlin, Polynomial Form, Ogden Potential, Arruda-Boyce, Gent, Yeoh, and Extended Tube. These models work with following elements following elements: SHELL181, PLANE182, PLANE183, SOLID185, SOLID186 , SOLID187, SOLID272, SOLID273, SOLID285, SOLSH190, SHELL208, SHELL209, SHELL281, PIPE288, PIPE289, and ELBOW290.
Compressible isotropic models: Blatz-Ko and Ogden Compressible Foam options are applicable to compressible foam or foam-type materials. These models work with the following elements: SHELL181, PLANE182, PLANE183, SOLID185, SOLID186 , SOLID187, SOLID272, SOLID273, SOLID285, SOLSH190, SHELL208, SHELL209, SHELL281, PIPE288, PIPE289, and ELBOW290.
Incompressible or nearly-incompressible isotropic response function hyperelastic model. This model (TB,EXPE) uses experimental data and works with the following elements: SHELL181, PLANE182, PLANE183, SOLID185, SOLID186 , SOLID187, SOLID272, SOLID273, SOLID285, SOLSH190, SHELL208, SHELL209, SHELL281, PIPE288, PIPE289, and ELBOW290.
Invariant-based anisotropic strain-energy potential. This option (TB,AHYPER) works with the following elements: PLANE182 and PLANE183 with plane strain and axisymmetric options, and SOLID185, SOLID186, SOLID187, SOLID272, SOLID273, SOLID285, and SOLSH190.
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:
S_{ij} = components of the second Piola-Kirchhoff stress tensor |
W = strain-energy function per unit undeformed volume |
E_{ij} = components of the Lagrangian strain tensor |
C_{ij} = components of the right Cauchy-Green deformation tensor |
The Lagrangian strain can be expressed as follows:
(4–167) |
where:
δ_{ij} = Kronecker delta (δ_{ij} = 1, i = j; δ_{ij} = 0, i ≠ j) |
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:
F_{ij} = components of the deformation gradient tensor |
X_{i} = undeformed position of a point in direction i |
x_{i} = X_{i} + u_{i} = deformed position of a point in direction i |
u_{i} = displacement of a point in direction i |
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:
I_{1}, I_{2}, and I_{3} = invariants of C_{ij}, |
(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:
α = coefficient of the thermal expansion |
ΔT = temperature difference about the reference temperature |
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.
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.
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) |
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:
G_{e} = constraint contribution to modulus |
G_{c} = crosslinked contribution to modulus |
δ = extensibility parameter |
β = empirical parameter (0 ≤ β ≤1) |
d_{1} = material incompressibility parameter |
The model is equivalent to a two-term Ogden model with the following parameters:
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) |
J_{m} = 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.
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:
c_{10}, c_{01}, 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:
c_{10}, c_{01}, c_{11}, 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:
c_{10}, c_{01}, c_{20}, c_{11}, c_{02}, 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:
c_{10}, c_{01}, c_{20}, c_{11}, c_{02}, c_{30}, c_{21}, c_{12}, c_{03}, 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) |
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 |
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, α _{1} = -2, μ_{1}= –μ, and β = 0.5, the Ogden option is equivalent to the Blatz-Ko option.
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}, d_{k} = 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 α _{1} = 2, the Ogden potential is equivalent to the Neo-Hookean potential. For N = 2, α _{1} = 2 and α _{2} = -2, the Ogden potential can be converted to the 2 parameter Mooney-Rivlin model.
The polynomial form of strain-energy potential is
(4–206) |
where:
N = material constant (input as NPTS on TB,HYPER) |
c_{ij}, d_{k} = 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_{01} = 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) |
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) |
C_{i0} = material constants (input on TBDATA commands with TB,HYPER) |
d_{k} = 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:
W_{v} = volumetric part of the strain energy |
W_{d} = deviatoric part of strain energy (often called isochoric part of the strain energy) |
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:
ρ, ρ_{o} = mass densities in the current and initial configurations |
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.
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:
Uniaxial Tension and Equibiaxial Compression.
Uniaxial Compression and Equibiaxial Tension.
Planar Tension and Planar Compression.
With several equivalent modes of testing, we are left with only three independent deformation states for which one can obtain experimental data.
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_{3}, 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 |
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) |
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) |
(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: λ_{2} = 1.
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) |
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) |
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 |
T_{i}(C_{j}) = 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^{2} = 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_{1}.
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:
λ_{1},λ_{2}, and λ_{3} = stretches for the general deformation |
λ = unknown uniaxial stretch that yields an equivalent value of the first invariant I_{1}. |
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:
dσ = change in the Cauchy stress tensor corresponding to a change in the logarithmic strain |
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.