COMSOL 6.2 - Hyperelastic Materials

A hyperelastic material is defined by its elastic strain energy density Ws, which is a function of the elastic strain state. It is often referred to as the energy density. The hyperelastic formulation normally gives a nonlinear relation between stress and strain, as opposed to Hooke’s law in linear elasticity.

Most of the time, the right Cauchy–Green deformation tensor C is used to describe the current state of strain (although one could use the left Cauchy–Green tensor B, the deformation gradient tensor F, and so forth), so the strain energy density is written as Ws(C).

For isotropic hyperelastic materials, any state of strain can be described in terms of three independent variables — common choices are the invariants of the right Cauchy–Green tensor C, the invariants of the Green–Lagrange strain tensor, or the principal stretches.

Once the strain energy density is defined, the second Piola–Kirchhoff stress is computed as

In the general case, the expression for the energy Ws is symbolically evaluated down to the components of C using the invariants definitions prior to the calculations of the components of the second Piola–Kirchhoff stress tensor. The differentiation is performed in components on the local coordinate system.

	In Equation View, the definitions of the stress components are shown as solid.Sl11 = 2*d(solid.Ws,solid.Cl11), solid.Sl12 = d(solid.Ws,solid.Cl12), and so on. The factor 2 in front of the differentiation operator for the shear stresses is omitted, since the symmetry in the Cauchy–Green tensor will cause two equal contributions.

	Modeling Geometric Nonlinearity

For hyperelastic materials, the decomposition between elastic and inelastic deformation is made using a multiplicative decomposition of the deformation gradient

Here, the inelastic deformation tensor Fin depends on the inelastic process, such as thermal expansion, hygroscopic swelling, or plasticity.

In this case, the strain energy density depends on the elastic deformation only, Ws(Cel), and the second Piola–Kirchhoff stress can be written in terms of its elastic counterpart

where Jin = det(Fin) and

If thermal expansion is present, a stress-free volume change occurs. This is a pure volumetric change, so the multiplicative decomposition of the deformation gradient tensor in Equation 3-4 implies

Here, the thermal volume ratio, Jth, depends on the thermal stretch λth, which for linear thermal expansion in isotropic materials can be written in terms of the isotropic coefficient of thermal expansion, αiso, and the absolute change in temperature

Here, the term αiso(T − Tref) is the thermal strain. The isotropic thermal gradient is therefore a diagonal tensor defined as

	The internal variables for the thermal stretch and the thermal volume ratio are named solid.stchth and solid.Jth.

When the coefficient of thermal expansion α is anisotropic, the thermal strain is computed from

and the anisotropic thermal gradient is defined as

The tangent coefficient of thermal expansion can be computed then as

On the opposite, when the tangent coefficient of thermal expansion is specified as a function of temperature, the thermal strain can be found as

where a matrix exponent is assumed, and the Hencky thermal strain tensor is defined as

When the thermoelastic damping effect needs to be taken into account, the heat source term can be computed as

where the Mandel stress tensor is introduced as

Hygroscopic swelling is an internal strain caused by changes in moisture content. This strain depends linearly on the moisture content

where βh is the coefficient of hygroscopic swelling, cmo is the moisture concentration, and cmo,ref is the strain-free reference concentration. The coefficient of hygroscopic swelling can represent isotropic or anisotropic swelling. The anisotropic hygroscopic gradient is defined as

It is possible to combine the hyperelastic material models with plasticity, viscoplasticity, and creep. Since these models are primarily used for large strain applications, only the large strain formulation is available. The decomposition between elastic and plastic deformation is made using a multiplicative decomposition of the deformation gradient tensor,

Here, the plastic deformation tensor Fpl depends on the plastic flow rule, yield function, and plastic potential.

For some classes of hyperelastic materials it is convenient to split the strain energy density into volumetric (also called dilatational) and isochoric (also called distortional or volume-preserving) contributions. The elastic deformation tensor is then multiplicatively decomposed into the volumetric and isochoric components

with Fel,vol as the volumetric elastic deformation (a diagonal tensor) and

the isochoric elastic deformation gradient. Isochoric deformation means that the volume ratio is kept constant during deformation, so the isochoric elastic deformation is computed by scaling it by the elastic volume ratio. The elastic volume ratio is defined by

and the volumetric deformation as

By using Jel it is possible to define the isochoric-elastic deformation gradient

the isochoric-elastic right Cauchy–Green tensor

and the isochoric-elastic Green–Lagrange strain tensor

This scaling changes the eigenvalues of the tensor, but not its principal directions, so the original and isochoric tensors remain coaxial to each other.

Some authors call

and

the modified tensors. Note that

	The internal variables for the isochoric-elastic Cauchy–Green deformation tensor in local coordinate system are named solid.CIel11, solid.CIel12, and so on.

The other two invariants normally used together with Jel are the first and second invariant of the isochoric-elastic right Cauchy–Green deformation tensor

Since

, the third invariant is never explicitly used.

	The internal variables for the invariants Jel, , and are named solid.Jel, solid.I1CIel, and solid.I2CIel.

The invariants of the isochoric (modified) elastic Green–Lagrange strain tensor are related to the invariants of the isochoric-elastic right Cauchy–Green deformation tensor

	The internal variables for the invariants of the isochoric elastic Green–Lagrange strain tensor are named solid.I1eIel, solid.I2eIel, and solid.I3eIel.

The isochoric and volumetric responses of most materials are coupled (Ref. 12), in this case, the strain energy density is written as a function of the elastic deformation

however, in many situations it is possible to assume that the pressure and deviatoric stress are uncoupled, and therefore for numerous hyperelastic materials it is possible to define the strain energy density as the sum of isochoric and volumetric counterparts

here,

is the isochoric elastic right Cauchy–Green deformation tensor, and Jel is the elastic volume ratio.

The volumetric strain energy density, Wvol, is defined as an expression of the elastic volumetric deformation Jel and the bulk modulus κ.

The volumetric strain energy density is defined as:

The volumetric strain energy density (Ref. 12) is defined as:

The volumetric strain energy density (Ref. 13) is defined as:

The volumetric strain energy density (Ref. 15) is defined as:

and the volumetric strain energy density (Ref. 16) is defined as:

For a given volumetric strain energy density, the volumetric stress (pressure) is calculated as

When the quadratic expression in Equation 3-40 is used, the pressure becomes linearly related to the volume change:

but for other volumetric representations the pressure becomes a nonlinear function of the volume change. For instance, if the expression in Equation 3-42 is used instead of the quadratic function, the pressure reads

A numerical scheme is said to exhibit locking if the accuracy of the approximation deteriorates as a parameter tends to a limiting value (Ref. 17). Finite elements in solid mechanics are said to “lock” when exhibiting an unphysical response to deformation (Ref. 18). Locking can occur for many different reasons. For linear elastic materials, this typically happens as Poisson’s ratio tends to 0.5, or the bulk modulus is much larger than the shear modulus. Numerical errors arise because the shape functions are unable to properly describe the volume preserving deformation.

To avoid the locking problem in computations, the mixed formulation replaces pm in Equation 3-45 with a corresponding interpolated pressure help variable pw, which adds an extra degree of freedom to the ones defined by the displacement vector u.

The general procedure is the same as described in Mixed Formulation for Linear Elastic materials.

For incompressible hyperelastic materials, the volumetric strain energy density Wvol is not defined at all, and the strain energy density Ws consists only of the isochoric contribution and the incompressibility constraint

Instead of deriving the pressure from the volumetric strain energy, a weak constraint is added to account for the incompressibility condition

The auxiliary pressure variable, pw, acts as Lagrange multiplier to enforce the constraint Jel = 1. This variable, positive in compression, is then used as the pressure when computing stresses.

The contribution to the virtual work is

The second Piola–Kirchhoff stress is then given by

and the Cauchy stress tensor by

For nearly incompressible hyperelastic materials, the strain energy density Ws is decoupled into isochoric and volumetric counterparts:

here, Wiso is the isochoric strain energy density and Wvol is the volumetric strain energy density.

The isochoric strain energy density depends on the material model, and it can be an expression involving the following:

•

For User Defined hyperelastic materials, it is also possible to define the strain energy density in terms of components of the isochoric elastic deformation gradient

The volumetric strain energy density can be any of the expression described in Volumetric Response. The default value for the bulk modulus is set to 100 times the equivalent shear modulus at infinitesimal deformation, which gives an initial Poisson’s ratio of approximately ν = 0.495.

An auxiliary variable, pw, is added to map the pressure pm derived from the volumetric strain energy density as described in Equation 3-45.

Then the variational problem is computed by the so-called perturbed Lagrangian method (Ref. 19), so the contribution to the virtual work reads

The second Piola–Kirchhoff stress is then defined as

where the extra pressure variable satisfies the weak constraint

the Cauchy stress tensor then reads

	Using the Quadratic volumetric strain energy density results in the only mixed formulation that returns a symmetric coupled stiffness matrix.

Different hyperelastic material models are constructed by specifying different elastic strain energy expressions. There are several predefined material models, as well as the option to enter user defined expressions for the isochoric and volumetric strain energy densities.

The strain energy density for the version of the Neo-Hookean material is written in terms of the elastic volume ratio Jel and the first invariant of the elastic right Cauchy–Green deformation tensor I1(Cel) (Ref. 11, 14)

Here, λ and μ are the Lamé parameters (SI unit: Pa).

In Ref. 8 the coupled strain energy density is defined from Lamé parameters

The and versions use the isochoric invariant

to define the isochoric strain energy density

where the volumetric strain energy Wvol can use any of the expressions described in Volumetric Response. See Nearly Incompressible Hyperelastic Materials for details.

The option uses the same isochoric strain energy, but an extra variable is added to enforce the incompressibility condition Jel = 1, see Incompressible Hyperelastic Materials.

	It is possible to define the isotropic hyperelastic material from other pair of elastic moduli, see the conversion from one set of elastic moduli to Lamé parameters in Table 3-1. See also the description of the Neo-Hookean material model in the Solid Mechanics interface documentation.

One of the simplest hyperelastic material models is the St Venant–Kirchhoff material, which is an extension of a linear elastic material into the hyperelastic regime.

For the response, the elastic strain energy density is written with two parameters and two invariants of the elastic Green–Lagrange strain tensor, I1(εel) and I2(εel)

Here, λ and μ are Lamé parameters (SI unit: Pa). The bulk modulus κ is calculated from κ = λ + 2μ/3.

The and the versions use the isochoric invariants

and

to define the isochoric strain energy density

The elastic volume ratio Jel and the bulk modulus κ are used to define the volumetric strain energy density Wvol, see Volumetric Response and Nearly Incompressible Hyperelastic Materials.

	It is possible to define the isotropic hyperelastic material from other pair of elastic moduli, see the conversion from one set of elastic moduli to Lamé parameters in Table 3-1. See also the description of the St Venant–Kirchhoff material model in the Solid Mechanics interface documentation.

For the and versions, the isochoric strain energy density is written in terms of the two isochoric invariants of the elastic right Cauchy–Green deformation tensors

and

The material parameters C10 and C01 (SI unit: Pa) are related to the Lamé parameter (shear modulus) μ = 2(C10 + C01).

The elastic volume ratio Jel and the bulk modulus κ are used to define the volumetric strain energy density Wvol, see Volumetric Response and Nearly Incompressible Hyperelastic Materials.

The option uses the same isochoric strain energy, but an extra variable is added to enforce the incompressibility condition Jel = 1, see Incompressible Hyperelastic Materials.

	See also the description of the Mooney–Rivlin, Two Parameters material model in the Solid Mechanics interface documentation.

Rivlin and Saunders (Ref. 2) proposed a phenomenological model for small deformations in rubber-based materials on a polynomial expansion of the first two invariants of the elastic right Cauchy–Green deformation, so the strain energy density is written as an infinite series

with C00 = 0. This material model is sometimes also called polynomial hyperelastic material.

In the first-order approximation, the material model recovers the Mooney–Rivlin strain energy density

while the second-order approximation incorporates second-order terms

The and the versions use the isochoric invariants of the elastic right Cauchy–Green deformation tensors

and

to define the isochoric strain energy density

The elastic volume ratio Jel and the bulk modulus κ are used to define the volumetric strain energy density Wvol, see Volumetric Response and Nearly Incompressible Hyperelastic Materials.

The incompressible option uses the same isochoric strain energy, but an extra variable is added to enforce the incompressibility condition Jel = 1, see Incompressible Hyperelastic Materials.

	See also the description of the Mooney–Rivlin, Five Parameters material model in the Solid Mechanics interface documentation.

The Mooney–Rivlin, nine parameters material model is an extension of the polynomial expression to third order terms and the strain energy density is written as

where

and

are the isochoric invariants of the elastic right Cauchy–Green deformation tensors.

The elastic volume ratio Jel and the bulk modulus κ are used to define the volumetric strain energy density Wvol, see Volumetric Response and Nearly Incompressible Hyperelastic Materials.

The option uses the same isochoric strain energy, but an extra variable is added to enforce the incompressibility condition Jel = 1, see Volumetric Response and Incompressible Hyperelastic Materials.

	See also the description of the Mooney–Rivlin, Nine Parameters material model in the Solid Mechanics interface documentation.

Yeoh proposed (Ref. 1) a phenomenological model in order to fit experimental data of filled rubbers, where Mooney–Rivlin and Neo-Hookean models were to simple to describe the stiffening effect in the large strain regime. The strain energy was fitted to experimental data by means of three parameters (SI unit: Pa), and the first invariant of the elastic right Cauchy–Green deformation tensors I1(Cel)

The shear modulus depends on the deformation, and it is calculated as

This imposes a restriction on the coefficients c1, c2, c3, since μ > 0.

The and the versions use the isochoric invariant of the elastic right Cauchy–Green deformation tensor

to define the isochoric strain energy density

The elastic volume ratio Jel and the bulk modulus κ are used to define the volumetric strain energy density Wvol, see Volumetric Response and Nearly Incompressible Hyperelastic Materials.

The option uses the same isochoric strain energy, but an extra variable is added to enforce the incompressibility condition Jel = 1, see Incompressible Hyperelastic Materials.

	See also the description of the Yeoh material model in the Solid Mechanics interface documentation.

The Neo-Hookean material model usually fits well to experimental data at moderate strains but fails to model hyperelastic deformations at high strains. In order to model rubber-like materials at high strains, Ogden adapted (Ref. 1) the energy of a Neo-Hookean material to

Here μp (SI unit: Pa) and αp (dimensionless) are material parameters, and λel1, λel2, and λel3 are the principal elastic stretches such as Jel = λel1λel2λel3.

The Ogden model is empirical, in the sense that it does not relate the material parameters αp and μp to physical phenomena. The parameters αp and μp are obtained by curve-fitting measured data, which can be difficult for N > 2. The most common implementation of Ogden material is with N = 2, so four parameters are needed.

The and the versions use the isochoric elastic stretches

so the isochoric strain energy density is defined as

The isochoric elastic stretches define a volume preserving deformation, since

The shear modulus at infinitesimal deformation is then defined from

The elastic volume ratio Jel and the bulk modulus κ are used to define the volumetric strain energy density Wvol, see Volumetric Response and Nearly Incompressible Hyperelastic Materials.

The option uses the same isochoric strain energy, but an extra variable is added to enforce the incompressibility condition Jel = 1, see Incompressible Hyperelastic Materials.

	Sometimes a slightly different definition of the strain energy function for the Ogden material is used. If you have given material data for an Ogden material, be careful that the definitions used are the same. If not, you will need to rescale the μp coefficients.

	See also the description of the Ogden material model in the Solid Mechanics interface documentation.

The Storakers material (Ref. 18 and Ref. 25) is used to model highly compressible foams. The strain energy density is written in a similar fashion as in Ogden material:

The initial shear and bulk moduli are computed from the parameters μk (SI unit: Pa) and βk (dimensionless) as

for constant parameters βk = β, the initial bulk modulus becomes κ = 2μ(β + 1/3), so a stable material requires μ > 0 and β > −1/3. In this case, the Poisson’s ratio is given by ν = β/(2β + 2/3), which means that for a Poisson’s ratio larger than −1, β > −2/9 is needed.

	See also the description of the Storakers material model in the Solid Mechanics interface documentation.

The Varga material model (Ref. 1) describes the strain energy in terms of the elastic stretches as

The and versions use the isochoric elastic stretches defined as

to define the isochoric strain energy density:

The simplest Varga model is obtained by setting c1 = μ and c2 = 0:

The elastic volume ratio Jel and the bulk modulus κ are used to define the volumetric strain energy density Wvol, see Volumetric Response and Nearly Incompressible Hyperelastic Materials.

The option uses the same isochoric strain energy, but an extra variable is added to enforce the incompressibility condition Jel = 1, see Incompressible Hyperelastic Materials.

	See also the description of the Varga material model in the Solid Mechanics interface documentation.

The other hyperelastic materials described are phenomenological models in the sense that they do not relate the different material parameters (normally obtained by curve-fitting experimental data) to physical phenomena.

Arruda and Boyce (Ref. 3) derived a material model based on Langevin statistics of polymer chains.

The strain energy density is defined by

here, n is the chain density, kB is Boltzmann constant, T is the absolute temperature, and β is computed from the inverse Langevin function β = L−1(λchain/λlock).

The number of segments in the polymeric chain, N, is used to derive the locking stretch of the chain,

, and the equivalent chain stretch λchain is computed from the first invariant of the elastic right Cauchy–Green deformation tensor I1(Cel),

The macroscopic shear modulus is defined as μ0 = nkBT. This material model is sometimes also called the eight-chain model.

In Ref. 3, a series expansion of the inverse Langevin function is used to simplify the expression for the strain energy density. In this case, the strain energy density reads

Here, μ0 is the initial macroscopic shear modulus, I1(Cel) is the first invariant of the elastic right Cauchy–Green deformation tensor, and the coefficients cp are obtained by series expansion of the inverse Langevin function.

Arruda and Boyce truncated the series and used only the first five terms. The coefficients cp (dimensionless) are listed in Table 3-2:

Table 3-2: First Five Coefficients of Arruda–Boyce Material Model.
c1	c2	c3	c4	c5
1/2	1/20	11/1050	19/7000	519/673750

Other authors (Ref. 1) use only the first three terms of the series.

The and versions use the isochoric invariant

to define the isochoric strain energy density, so when using the inverse Langevin function the chain stretch is defined as

, and in the five terms approximation the isochoric strain energy is defined as

Then, the elastic volume ratio Jel and the bulk modulus κ are used to define the volumetric strain energy density Wvol, see Volumetric Response and Nearly Incompressible Hyperelastic Materials.

The option uses the same isochoric strain energy, but an extra variable is added to enforce the incompressibility condition Jel = 1, see Incompressible Hyperelastic Materials.

	See also the description of the Arruda–Boyce material model in the Solid Mechanics interface documentation.

Many hyperelastic material models are difficult to fit to experimental data. The Gent material (Ref. 20 and Ref. 21) is a simple phenomenological constitutive model based on only two parameters, μ and jm, which defines the strain energy density as

Here, μ (SI unit: Pa) is the shear modulus and jm (dimensionless) is a limiting value for I1 − 3, which takes care of the limiting polymeric chain extensibility of the material.

Since the strain energy density does not depend on the second invariant I2, the Gent model is often classified as a generalized Neo-Hookean material. The strain energy density tends to be the one of incompressible Neo-Hookean material as jm → ∞. The Gent material is the simplest model of the limiting chain extensibility family.

The and formulations use the isochoric invariants

to define the isochoric strain energy density

The elastic volume ratio Jel and the bulk modulus κ are used to define the volumetric strain energy density Wvol, see Volumetric Response and Nearly Incompressible Hyperelastic Materials.

The option uses the same isochoric strain energy, but an extra variable is added to enforce the incompressibility condition Jel = 1, see Incompressible Hyperelastic Materials.

	See also the description of the Gent material model in the Solid Mechanics interface documentation.

Kilian and co-workers formulated a van der Waals equation of state for real networks of polymer chains from the analogy of an ideal gas (Ref. 22–24). The constitutive model is based on four parameters which define the averaged invariant

in order to define the isochoric strain energy density

Here, μ (SI unit: Pa) is the shear modulus. The maximum chain stretch λm (dimensionless) represents a limiting value for the averaged invariant

, which accounts for the maximum chain extensibility in networks with finite chain lengths. The dimensionless parameter β phenomenologically averages the two isochoric invariants

and

, and the parameter α (dimensionless) represents the global interaction between polymer chains.

The and formulations use the elastic volume ratio Jel and the bulk modulus κ to define the volumetric strain energy density Wvol, see Volumetric Response and Nearly Incompressible Hyperelastic Materials.

The option uses the same isochoric strain energy, but an extra variable is added to enforce the incompressibility condition Jel = 1; see Incompressible Hyperelastic Materials

	See also the description of the van der Waals material model in the Solid Mechanics interface documentation.

The Blatz–Ko material model was developed for foamed elastomers and polyurethane rubbers, and it is used to model compressible isotropic hyperelastic materials (Ref. 1, 26).

The elastic strain energy density is written with three parameters and the three invariants of the elastic right Cauchy–Green deformation tensor, I1(Cel), I2(Cel), and I3(Cel) = Jel2

Here, ϕ is an interpolation parameter bounded to 0 < ϕ < 1, the parameter μ (SI unit: Pa) is the shear modulus, and β is an expression of Poisson’s ratio, β = ν/(1 − 2ν), or Lamé parameters, β = λ/2μ.

When the parameter β → ∞, or equivalently, Poisson’s ratio tends to 0.5; the strain energy simplifies to a similar form of the Mooney–Rivlin, Two Parameters model

In the special case of ϕ = 1, the strain energy reduces to a compressible Neo-Hookean model (use 2β = λ/μ = 2ν/(1 − 2ν)),

This expression is also equivalent to the compressible Storakers model consisting in one term, and the material parameter defined as α1 = 2.

	See also the description of the Blatz–Ko material model in the Solid Mechanics interface documentation.

Gao proposed (Ref. 27) a simple compressible hyperelastic material where the strain energy density is defined by two parameters, a (SI unit: Pa) and n, and two invariants of the elastic right Cauchy–Green deformation tensors Cel:

Here, the invariant I−1(Cel) is calculated as in the Blatz–Ko model:

Gao proposed that the material is unconditionally stable when the parameters are bounded to 1 < n < 3 and 0 < a, and related these parameters under small strain to the Young’s modulus and Poisson’s ratio by

Since n = (1 + ν)/(1 − 2ν) and it is bounded to 1 < n < 3, this material is stable for an initial Poisson’s ratio in the range of 0 < ν < 2/7.

	See also the description of the Gao material model in the Solid Mechanics interface documentation.

The Murnaghan strain energy density is used in nonlinear acoustoelasticity. Most conveniently, it is expressed in terms of the three invariants of the elastic Green–Lagrange strain tensor, I1(εel), I2(εel), and I3(εel):

Here, l, m, and n (SI unit: Pa) are the Murnaghan third-order elastic moduli, which can be found experimentally for many commonly encountered materials such as steel and aluminum, and λ and μ are the Lamé parameters. Setting the Murnaghan moduli to zero, l = m = n = 0, recovers a compressible St Venant–Kirchhoff material.

	See also the description of the Murnaghan material model in the Solid Mechanics interface documentation.

Delfino proposed a simple hyperelastic material for modeling carotid arteries (Ref. 29). The strain energy density is defined by two parameters, a (SI unit: Pa) and b (dimensionless), and the first invariant of the isochoric elastic right Cauchy–Green deformation tensors

The parameter a plays the role of shear modulus at the small strain limit.

The option uses the isochoric strain energy, and an extra variable is added to enforce the incompressibility condition Jel = 1, see Incompressible Hyperelastic Materials.

Fung (Ref. 29) proposed one of the most popular strain energy functions to describe the deformation in soft biological tissues and arteries. For this compressible hyperelastic material, the strain energy density reads

the parameter c (SI unit: Pa) scales the global stiffness, and the quadratic form Q (dimensionless) depends on the Green–Lagrange strain

here, A is a 6-by-6 symmetric matrix (dimensionless), and E is a vector which contains the elements of the elastic Green–Lagrange strain tensor sorted in either standard or Voigt order, see Orthotropic and Anisotropic Materials.

For modeling the deformation in arteries, Fung proposed that the principal directions of the strain tensor coincide with the radial, circumferential, and axial directions of the artery, so the coefficient matrix A is described by nine independent parameters (Ref. 29), see Orthotropic Material to see the structure of such matrix.

The and formulations use the components of the isochoric elastic Green–Lagrange strain

to define the isochoric strain energy density

also, these formulations use the elastic volume ratio Jel and the bulk modulus κ to define the volumetric strain energy density Wvol, see Volumetric Response and Nearly Incompressible Hyperelastic Materials.

The option uses the same isochoric strain energy, but an extra variable is added to enforce the incompressibility condition Jel = 1, see Incompressible Hyperelastic Materials.

The extended tube model is a micromechanics-inspired material model that considers the network constraints from molecular chains and the limited chain extensibility (Ref. 33–35). The isochoric strain energy density of the model consists of two terms

The term Wc represents the energy from the cross-linking of the network, and We represents the energy from confining tube constrains. These are written as

Here, Gc and Ge are material parameters (SI unit: Pa),

is the first invariant of the elastic right Cauchy–Green deformation tensor,

are the isochoric principal stretches, and α and β are dimensionless coefficients.

The parameters Gc and Ge plays the role of shear modulus at the small strain limit, G0 = Gc + Ge. The parameter α corresponds to the maximum finite chain extensibility, as the strain energy Wc is singular for stretch values such as

The option uses the isochoric strain energy, and it adds an extra variable to enforce the incompressibility condition Jel = 1, see Incompressible Hyperelastic Materials.

Note that the material might be unstable for certain combinations of material parameters, see Ref. 33-34 for a valid range of parameter values.

It is possible to define the strain energy density for compressible, nearly incompressible or incompressible hyperelastic materials when selecting the user defined option.

For hyperelastic materials, enter an expression for the elastic strain energy Ws, which can include any expressions involving the following:

•

Components of Cel, the elastic right Cauchy–Green deformation tensor in the local material coordinate system.

	The internal variables for these invariants are named solid.I1Cel, solid.I2Cel, and solid.I3Cel.

•

Principal elastic stretches λel1, λel2, and λel3, which are the square-root of the eigenvalues of the elastic right Cauchy–Green deformation tensor Cel.

	The internal variables for the principal elastic stretches are named solid.stchelp1, solid.stchelp2, and solid.stchelp3.

	The internal variables for these invariants are named solid.I1eel, solid.I2eel, and solid.I3eel.

When the option is selected, the elastic strain energy is decoupled into the volumetric and isochoric components:

•

The isochoric strain energy, Wiso, as an expression involving the invariants of the isochoric elastic right Cauchy–Green tensor

and

; the invariants of the isochoric elastic Green–Lagrange strain

, and

; or the principal isochoric elastic stretches defined as

When the option is selected, an extra variable is added to map the pressure, see The Locking Problem and Nearly Incompressible Hyperelastic Materials for details.

	The internal variables for Jel, , and are named solid.Jel, solid.I1CIel, and solid.I2CIel. The internal variables for , , and are named solid.I1eIel, solid.I2eIel, and solid.I3eIel.

When the option is selected, enter an expression for the isochoric elastic strain energy Wiso, as is done for the option. An extra variable is added to enforce the incompressibility condition Jel = 1, see Incompressible Hyperelastic Materials and Mixed Formulation.

	The strain energy density must not contain any other expressions involving displacement or their derivatives. Examples of such expressions are components of the displacement gradient ∇u and deformation gradient F = ∇u + I tensors, their transpose, inversions, as well as the global material system components of C and ε. If they occur, such variables are treated as constants during symbolic differentiations.

	Select the Use elastic deformation gradient check box to define the elastic strain energy density Ws, the isochoric strain energy density, Wiso, and the volumetric strain energy, Wvol, in terms of the components the elastic deformation gradient Fel. Then, the second Piola–Kirchhoff stress is defined from .

	See also the description of the User Defined material model in the Solid Mechanics interface documentation.

Some nonlinear effects observed in rubbers, such as hysteresis in stress-stretch curves, residual strains, and stress softening effects, are not accounted in the formulation of common hyperelastic materials. The Mullins effect (Ref. 43-45) describes the stress-softening phenomenon observed under cyclic loading in elastomers and biological materials.

Ogden and Roxburgh (Ref. 46) used an additional state variable to model the Mullins effect. The state variable η is introduced to memorize the microstructural damage on reinforced rubber after repeated loading-unloading cycles. The modified isochoric strain energy density reads

here, Wiso is the isochoric strain energy of the undamaged material, and ϕ(η) is referred as the damage function. The choice of the damage function ϕ(η) is completely arbitrary as long as some constrains are fulfilled. The authors (Ref. 46) proposed a state function based on the error function, which defined how the state variable η varies as a function of the isochoric strain energy

here, erf(.) is the error function, r and m are positive parameters, and Wmax is the maximum attained value of the isochoric strain energy density on the loading path.

The associated microstructural damage is computed from

and the isochoric strain energy modified by the Mullins effect reads

Over the years others authors (Ref. 47-48) have proposed different flavors of Ogden–Roxburgh model, the version implemented in COMSOL Multiphysics uses by default the hyperbolic tangent function, and a parameter to define the maximum allowed damage d∞ = 1/r. The microstructural damage is then computed from

Here, d∞, α, and Wsat are positive parameters. It is also possible to use the error function as in the original formulation (Ref. 46), in which case the microstructural damage is computed from