Ogden Materials

Hyperelastic materials can be used to model the isotropic, nonlinear elastic behavior of rubber, polymers, and similar materials. These materials are nearly incompressible in their behavior and can be stretched to very large strains.

In Radioss, material laws LAW42, LAW62, LAW69, LAW82, and LAW88 utilize different strain energy density functions of the Ogden material model ^¹ to model hyperelastic materials. ^²

Material Definition

Stretch (also called stretch ratio)

λ

is the ratio of final length and initial length. It is used for materials with large deformations. For a cube in tension:

$ε_{1} = \frac{Δ l}{l_{01}}$: Engineering strain (also called nominal strain) in direction 1
$λ_{1} = \frac{l_{1}}{l_{01}}$: Stretch in direction 1

Thus, strain and stretch are related as:

λ = 1 + ε

Principal stretch $λ_{i}$ can be used to describe the volumetric deformation by calculating the relative volume, $J$ , computed as:

J = \frac{V}{V_{0}} = \frac{l_{1} \cdot l_{2} \cdot l_{3}}{l_{01} \cdot l_{02} \cdot l_{03}} = λ_{1} \cdot λ_{2} \cdot λ_{3}

For an incompressible material the volume should not change and thus

J

=1 and thus, the stretch can be calculated for the following material tests.

Uniaxial test:
$λ_{1} = λ$ and $λ_{2}^{2} = λ_{3}^{2} = \frac{1}{λ}$
Biaxial test:
$λ_{1} = λ_{2} = λ$ and $λ_{3} = λ^{- 2}$
Planar (shear) test:
$λ_{1} = λ; λ_{3} = 1$ and $λ_{2} = λ^{- 1} $

/MAT/LAW42 (Ogden)

This material model defines a hyperelastic, viscous, and incompressible material specified using the Ogden, Neo-Hookean, or Mooney-Rivlin material models. This law is generally used to model incompressible rubbers, polymers, foams, and elastomers. This material can be used with shell and solid elements.

LAW42 uses the following strain energy density representation of the Ogden material model.

W (λ_{1}, λ_{2}, λ_{3}) = \sum_{p = 1}^{5} \frac{μ_{p}}{α_{p}} ({\bar{λ}}_{1}^{α_{p}} + {\bar{λ}}_{2}^{α_{p}} + {\bar{λ}}_{3}^{α_{p}} - 3) + \frac{K}{2} {(J - 1)}^{2}

Where,

$W$: Strain energy density
$λ_{i}$: i^th principal engineering stretch
$J$: Relative volume defined as: $J = λ_{1} \cdot λ_{2} \cdot λ_{3} = \frac{ρ_{0}}{ρ}$
${\bar{λ}}_{i} = J^{- \frac{1}{3}} λ_{i}$: Deviatoric stretch
$α_{p}$ and $μ_{p}$: Material constants coefficient pairs.; Up to 5 material constant pairs can be defined.

The initial shear modulus $μ$ and bulk modulus ( $K$ ) are given by:

μ = \frac{\sum_{p = 1}^{5} μ_{p} \cdot α_{p}}{2}

and

K = μ \cdot \frac{2 (1 + ν)}{3 (1 - 2 ν)}

Where, $ν$ is the Poisson's ratio and is only used for computing the bulk modulus.

Material Parameters

Parameters $α_{p}$ and $μ_{p}$ must be chosen so that initial shear modulus is:

μ = \frac{\sum_{p = 1}^{5} μ_{p} \cdot α_{p}}{2} > 0

For material stability, it is required that each material constant pair

μ_{p} \cdot α_{p} > 0

In general, the Ogden model can be used for strains up to 700%. The number of terms material pairs, $α_{p}$ and $μ_{p}$ , needed depends on the range of experimental data that is fit and curve fitting accuracy desired. In practice, 3 material pairs fit most data. If the material pairs are not known for a particular material, then a curve fit of uniaxial test data can be performed in Radioss using LAW69 or via separate fitting software.

Neo-Hookean Model

A simple case of the Ogden material model is the Neo-Hooken model represented using the following equation for the strain energy density function:

W = C_{10} (I_{1} - 3)

Where,

$I_{1}$: The first invariants of the right Cauchy-Green tensor
$C_{10}$: Material constant

This representation can be derived from the LAW42 Ogden strain energy density function when:

$μ_{1} = 2 \cdot C_{10}$ ; $α_{1} = 2$ , and $μ_{2} = α_{2} = 0$

The Neo-Hookean model is a simple model that is typically only accurate for strains less than 20%.

Mooney-Rivlin Model

A slightly more complex case of the LAW42 Ogden material model is the Mooney-Rivlin model, which can be represented using the following equation for the strain energy density function:

W = C_{10} (I_{1} - 3) + C_{01} (I_{2} - 3)

Where,

$I_{1}$ and $I_{2}$: The first and second invariants of the right Cauchy-Green tensor
$C_{10}$ and $C_{01}$: Material constants

This representation can be derived from the LAW42 Ogden strain energy density function when:

$μ_{1} = 2 \cdot C_{10}$ , $μ_{2} = - 2 \cdot C_{01}$ , $α_{1} = 2$ , and $α_{2} = -2$

Mooney-Rivlin constants are available from a material supplier or testing company. If they are not available, then a curve fit of uniaxial test data can be performed in Radioss using LAW69 or via separate fitting software. The Mooney-Rivlin material law is accurate for strains up to 100%.

Poisson's Ratio and Material Incompressibility

If a material is truly incompressible, then $ν = 0.5$ . However, in practice is not possible to use because that would result in an infinite bulk modulus, an infinite speed of sound, and thus an infinitely small solid element Time Step.

K = μ \cdot \frac{2 (1 + ν)}{3 (1 - 2 ν)} = μ \cdot \frac{2 (1 + ν)}{3 (1 - 2 * 0.5)} = \infty

The effect of different Poisson’s ratio input can be seen in Figure 2. The largest difference in the results is at higher amounts of strain. The results will match the test data better when $ν = 0.4997$ but this results in a time step that is 4 times lower than $ν = 0.495$ . Thus, to balance the computation time and accuracy it is recommended to use $ν = 0.495$ for incompressible rubber material.

The effect of Poisson’s ratio and Bulk modulus are similar in other Ogden material law.

Higher values of the Poisson’s ratio may lead to a very small time step or divergence for explicit simulations.

In LAW42, material incompressibility is provided by using a penalty approach, which calculates the pressure proportional to a change in density.

P = K \cdot F s c a l e_{b l k} \cdot f_{b l k} (J) \cdot (J - 1)

Where,

$K$: Bulk modulus
$J = \frac{V}{V_{0}} = \frac{m ρ_{0}}{m ρ} = \frac{ρ_{0}}{ρ}$: Relative volume which simplifies to relative density, if mass is constant
$f_{b l k} (J)$: Bulk coefficient scale factor versus relative volume function
$F s c a l e_{b l k}$: Abscissa scale factor for function $f_{b l k} (J)$

The bulk modulus (

K

) of hyperelastic materials is generally a very high value which provides the needed pressure-resistance to maintain the incompressibility condition (

J

=1). But if a material starts to compress (

J

< 1), the bulk modulus can be increased by including the fct_ID_blk input function which allows the scaling of the bulk coefficient value as a function of

J

. By default, there is no scaling and; thus, if the function identifier is zero and the value of the bulk scaling function is equal to 1. It is advisable to output (/ANIM/BRICK/DENS) and review the material density of LAW42 components to make sure that the density variation is small, that is the value of

J

is close to 1 and the material is incompressible.

Figure 3. Bulk Modulus Scale Factor Function `fct_ID`_blk

Viscous (Rate) Effects

Viscous (rate) effects are modeled in LAW42 using a Maxwell model, which can be described in a simplified manner as a system of

η

springs with stiffness'

G_{i}

and dampers

η_{i}

law82_maxwell_model — Figure 4. Maxwell Model

The Maxwell model is represented using Prony series inputs ( $G_{i}, τ_{i}$ ). The hyperelastic initial shear modulus $μ$ is the same as the long-term shear modulus $G_{\infty}$ in the Maxwell model, and $τ_{i}$ is the relaxation time:

τ_{i} = \frac{η_{i}}{G_{i}}

The $G_{i}$ and $τ_{i}$ values must be positive.

/MAT/LAW62 (VISC_HYP)

A hyper visco-elastic material law in Radioss that can be used to model polymers and elastomers.

The hyperelastic behavior in this material law is defined using the following strain energy density function:

W (λ_{1}, λ_{2}, λ_{3}) = \sum_{i = 1}^{N} \frac{2 μ_{i}}{α_{i}^{2}} (λ_{1}^{α_{i}} + λ_{2}^{α_{i}} + λ_{3}^{α_{i}} - 3 + \frac{1}{β} (J^{- α_{i} β} - 1))

Where,

$W$: Strain energy density
$λ_{i}$: $i_{t h}$ principal stretch
$J$: Relative volume defined in Equation 13
$β = \frac{ν}{(1 - 2 ν)}$
$α_{i}$ and $μ_{i}$: Material constants coefficient pairs.; Up to 5 material constant pairs can be defined.

Poisson’s ratio must be

0 < ν < 0.5

. This law can be used to model compressible or sometimes called hyperfoam materials by defining a low Poisson’s ratio value.

Note: The

μ_{i}

material coefficients are different, but can be converted using:

μ_{i}^{L A W 62} = \frac{μ_{i}^{L A W 42} \cdot α_{i}^{L A W 42}}{2}

Viscous (Rate) Effects

Viscous (rate) effects are modeled in LAW62 using a Maxwell model which can be described in a simplified manner as a system of n springs with stiffness’

G_{i}

and dampers

η_{i}

The Maxwell model is represented using a Prony series with inputs. The initial shear modulus is:

G_{0} = \sum_{i = 1}^{N} μ_{i}

The sum of $μ_{i}$ should be greater than 0.

G_{0} = G_{\infty} + \sum_{i} G_{i}

The stiffness ratio is:

γ_{\infty} = \frac{G_{\infty}}{G_{0}} = 1 - \sum_{i} γ_{i}

γ_{i} = \frac{G_{i}}{G_{0}}

With,

γ_{i} \in [0, 1], \sum_{i} γ_{i} < 1

and the ground shear modulus

G_{0} = G_{\infty} + \sum_{i} G_{i}

The relative time,

τ_{i}

must be positive:

τ_{i} = \frac{η_{i}}{G_{i}}

Note: When viscosity is included, the shear modulus in LAW62 is the initial shear modulus

G_{0} = \sum_{i = 1}^{N} μ_{i}

which includes viscosity, but in LAW42 the shear modulus is the long-term shear modulus, which does not include viscosity

G_{\infty} = \frac{\sum_{p = 1}^{5} μ_{p} \cdot α_{p}}{2}

/MAT/LAW69

This law, like /MAT/LAW42 (OGDEN) defines a hyperelastic and incompressible material specified using the Ogden or Mooney-Rivlin material models. Unlike LAW42, where the material parameters are input, this law computes the material parameters using test data from a uniaxial engineering stress-strain curve.

This material can be used with shell and solid elements.

The strain energy density formulation used depends on the law_ID:

law_ID = 1 (Ogden law):
$W (λ_{1}, λ_{2}, λ_{3}) = \sum_{p = 1}^{5} \frac{μ_{p}}{α_{p}} ({\bar{λ}}_{1}^{α_{p}} + {\bar{λ}}_{2}^{α_{p}} + {\bar{λ}}_{3}^{α_{p}} - 3) + \frac{K}{2} {(J - 1)}^{2}$
law_ID = 2 (Mooney-Rivlin law):
$W = C_{10} (I_{1} - 3) + C_{01} (I_{2} - 3)$

Material Parameters

After reading the stress-strain curve (fct_ID₁), Radioss calculates the corresponding material parameter pairs using a nonlinear least-square fitting algorithm. For the classic Ogden law, (law_ID=1), the calculated material parameter pairs are $μ_{p}$ and $α_{p}$ where the value of p is defined via the N_pair input. The maximum value is N_pair=5 with a default value of 2. Usually no more than N_pair=3 is needed for a good fit.

For the Mooney-Rivlin law (law_ID =2), the material parameter $C_{10}$ and $C_{01}$ are calculated remembering that $μ_{p}$ and $α_{p}$ for the LAW42 Ogden law can be calculated using this conversion:

$μ_{1} = 2 \cdot C_{10}$ , $μ_{2} = - 2 \cdot C_{01}$ , $α_{1} = 2$ , and $α_{2} = -2$

The minimum test data input should be a uniaxial tension engineering stress strain curve. If uniaxial compression data is available, the engineering strain should increase monotonically from a negative value in compression to a positive value in tension. In compression, the engineering strain should not be less than -1.0 since -100% strain is physically not possible.

To improve the quality of the nonlinear least square fit, it is recommended that:

The experimental data curve represents a smooth monotonically increasing function with uniform distribution of abscissa points. The number of data points in the experimental data curve should be greater than the number of parameter pairs (N_pair).
The engineering strain is negative in compression and positive in tension. For compression test data, the engineering strain should be greater than -1.0 (100% compression maximum) but tension only stress strain data can also be used.
If N_pair ≥ 3, then the test data should cover at least 100% of the tensile strain and/or 50% of the compressive strain.
N_pair should not be set to a very large value to avoid instabilities in the fitting procedure.

This material law is stable when $μ_{p} α_{p} > 0$ (with $p$ =1,…5) is satisfied for parameter pairs for all loading conditions. By default, Radioss tries to fit the curve by accounting for these conditions (I_check=2). If a proper fit cannot be found, then Radioss uses a weaker condition (I_check=1:), which ensures that the initial shear hyperelastic modulus ( $μ$ ) is positive.

To determine how well the calculated material parameters represent the input test data, the Radioss Starter outputs an “averaged error of fitting” value which is recommended to not exceed 10%. For visual comparison, the stress-strain curve calculated from the strain energy density and calculated material parameters is also output by the Radioss Starter.

Due to the friction involved in a uniaxial compression test, it is usually more accurate to take equal biaxial tension test data and convert it to uniaxial compressive data using these formulas ^³ which are valid for incompressible materials.

ε_{c} = \frac{1}{{(ε_{b} + 1)}^{2}} - 1

σ_{c} = σ_{b} {(1 + ε_{b})}^{3}

Where,

$ε_{c}$: Uniaxial engineering compressive stain
$ε_{b}$: Equal biaxial engineering tension strain
$σ_{c}$: Uniaxial engineering compressive stress
$σ_{b}$: Equal biaxial engineering tension stress

Material Incompressibility

Material LAW69 uses the same method to maintain incompressibility as LAW42. For additional information, refer to Poisson's Ratio and Material Incompressibility in LAW42.

Viscous (Rate) Effects

/VISC/PRONY must be used with LAW69 to include viscous effects. Alternatively, LAW69 could be used to extract the Ogden or Mooney-Rivlin parameters and then those parameters can be used in LAW42 with viscosity added.

/MAT/LAW82

This material model defines a hyperelastic, and incompressible material specified using the Ogden, Neo-Hookean, or Mooney-Rivlin material models. This law is generally used to model incompressible rubbers, polymers, foams, and elastomers.

This material can be used with shell and solid elements. As compared to LAW42 or LAW62, this law uses a different Ogden strain energy density formulation given in Equation 26. The LAW82 strain energy density formulation matches what is used in some other finite elements solver’s hyperelastic model and; thus, the material parameters for this form of the Ogden strain energy density are sometimes available from material suppliers or other sources.

W = \sum_{i = 1}^{N} \frac{2 μ_{i}}{α_{i}^{2}} ({\bar{λ}}_{1}^{α_{i}} + {\bar{λ}}_{2}^{α_{i}} + {\bar{λ}}_{3}^{α_{i}} - 3) + \sum_{i = 1}^{N} \frac{1}{D_{i}} {(J - 1)}^{2 i}

Where,

$W$: Strain energy density
$N$: Number of material constants $α_{i}$ , $μ_{i}$ and $D_{i}$
${\bar{λ}}_{i} = J^{- \frac{1}{3}} λ_{i}$: Deviatoric stretch
$J$: Relative volume as defined in Equation 2

The initial shear modulus:

μ = \sum_{i = 1}^{N} μ_{i}

The Bulk Modulus is calculated as

K = \frac{2}{D_{1}}

based on these rules:

If $ν = 0$ , $D_{1}$ should be entered.
If $ν \neq 0$ , $D_{1}$ input is ignored and will be recalculated and output in the Starter output using:
$D_{1} = \frac{3 (1 - 2 v)}{μ (1 + v)}$
If $ν = 0$ and $D_{1}$ =0, a default value of $ν = 0.475$ is used and $D_{1}$ is calculated using Equation 28

Neo-Hookean Model

Like LAW42, LAW82 can also be simplified to a Neo-Hookean model by using:

$μ_{1} = 2 \cdot C_{10}$ , $α_{1} = 2$ and $μ_{2} = α_{2} = 0$

Mooney-Rivlin Model

Like LAW42, LAW82 can also be simplified to a Mooney-Rivlin model by using:

$μ_{1} = 2 \cdot C_{10}$ , $μ_{2} = 2 \cdot C_{01}$ , $α_{1} = 2$ and $α_{2} = -2$

Viscous (Rate) Effects

/VISC/PRONY must be used with LAW82 to include viscous effects.

Drucker Condition Stability Check

In LAW42 and LAW69, the Drucker stability is automatically calculated by the Radioss Starter.

The Drucker stability condition checks if the change in the Kirchhoff stress corresponding to the infinitesimal change in the logarithmic strain (true strain) satisfies the following inequality.

\sum_{i = 1}^{3} d τ_{i} d ε_{i} > 0

Where, $i$ =1,2,3 principal direction.

With the change in logarithmic strain

d ε_{i} = \frac{d λ_{i}}{λ_{i}}

$d τ_{i} = J \cdot d σ_{i}$: The change of Kirchhoff stress
$d τ = D : d ε$: Relationship between Kirchhoff stress and logarithmic strain

The Drucker stability condition will be:

\sum d ε : D : d ε > 0

Here $D$ is tangential material stiffness matrix and it is also the slope of stress-strain curve:

D = [\begin{matrix} D_{11} & D_{12} & D_{13} \\ D_{21} & D_{22} & D_{23} \\ D_{31} & D_{32} & D_{33} \end{matrix}]

For a stable material, it requests tangential material stiffness $D$ be positive (slope of stress-strain curve is positive). The tangential material matrix $D$ is positive if following conditions satisfied:

I_{1} = t r (D) = D_{11} + D_{22} + D_{33} > 0

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

I_{3} = \det (D) > 0

The Kirchhoff stress for Ogden model is:

τ_{i} = \sum_{p} μ_{p} [{\bar{λ}}_{i}^{α_{p}} - \frac{1}{3} ({\bar{λ}}_{1}^{α_{p}} + {\bar{λ}}_{2}^{α_{p}} + {\bar{λ}}_{3}^{α_{p}})] + K (J^{2} - J)

Since

D_{i j} = \frac{\partial τ_{i}}{\partial λ_{j}}

, then for a given Ogden parameter

α_{p}, μ_{p}

with conditions

I_{1} > 0, I_{2} > 0 and I_{3} > 0

, the strain range of material in Drucker stability could then be calculated.

The Drucker stability criterion calculates the strain range where the material model will remain stable given a set of material parameters. This stability check cannot be made for every deformation but instead is commonly used to check material stability under uniaxial, biaxial and planar strain loading.

For example, using the following Ogden parameters:

$\begin{matrix} μ_{1} = 13.99077258830 & α_{1} = 3.788192935039 \\ μ_{2} = - 9.13454532223 & α_{2} = - 7.17617341059 \\ μ_{3} = 8.904655103235 & α_{3} = - 7.27028137148 \end{matrix}$

Then the Drucker stability will be automatically checked in Radioss Starter, and results printed in Starter output file *0.out. This shows the strain at which instability can occur for the given Ogden parameters:

CHECK THE DRUCKER PRAGER STABILITY CONDITIONS   
      -----------------------------------------------
     MATERIAL LAW = OGDEN (LAW42) 
     MATERIAL NUMBER =         1
       TEST TYPE = UNIXIAL  
        COMPRESSION:    UNSTABLE AT A NOMINAL STRAIN LESS THAN -0.3880000000000    
        TENSION:        UNSTABLE AT A NOMINAL STRAIN LARGER THAN  0.9709999999999    
       TEST TYPE = BIAXIAL  
        COMPRESSION:    UNSTABLE AT A NOMINAL STRAIN LESS THAN -0.2880000000000    
        TENSION:        UNSTABLE AT A NOMINAL STRAIN LARGER THAN  0.2780000000000    
       TEST TYPE = PLANAR (SHEAR)
        COMPRESSION:    UNSTABLE AT A NOMINAL STRAIN LESS THAN -0.3680000000000    
        TENSION:        UNSTABLE AT A NOMINAL STRAIN LARGER THAN  0.5829999999999

Note: For a Neo-Hookean material with

C_{10} > 0

(or

μ_{1} > 0

), the material is always stable and thus no critical value is found by the Drucker stability check.

For a Mooney-Rivlin material, the Drucker stability should be checked since $C_{01} or μ_{2}$ could be negative, which leads material instability.

/MAT/LAW88

This law utilizes a tabulated uniaxial tension and compression engineering stress and strain test data at different strain rates to model incompressible materials. It is only compatible with solid elements.

The material is based on the following Ogden’s strain energy density function but does not require curve fitting to extract material constants like most other hyperelastic material models. ^⁴

W = \underset{d e v i a t o r i c p a r t}{\underset{︸}{\sum_{i = 1}^{3} \sum_{j = 1}^{m} \frac{μ_{j}}{α_{j}} ({\bar{λ}}_{i}^{α_{j}} - 1)}} + \underset{s p h e r i c a l p a r t}{\underset{︸}{K (J - 1 - \ln J)}}

Instead, this law determines the Ogden function directly from the uniaxial engineering stress strain curve tabulated data.

Unlike other Ogden material laws, the Bulk Modulus must be input from either test data or extracted from Starter output of the LAW69 Ogden curve fit. When comparing results between LAW42 or LAW69 to LAW88, the same bulk modulus must be used.

Unloading Behavior

Unloading can be represented using an unloading function or by providing hysteresis and shape factor inputs to a damage model based on energy.

If using the damage model, the loading curves are used for both loading and unloading and the unloading stress tensor is reduced by:

σ = (1 - D) σ

with

D = (1 - H y s) (1 - {(\frac{W_{c u r}}{W_{\max}})}^{S h a p e})

Where,

$W_{c u r}$: Current energy
$W_{\max}$: Maximum energy corresponding to the quasi-static behavior
$H y s$ and $S h a p e$: Input by user

If an unloading curve is provided, these options are available:

Tension

Loading and Unloading

= 0

Loading use loading function fct_ID_Li and

{\dot{ε}}_{L i}

Unloading use unloading function fct_ID_unL

= 1

Loading and unloading all use loading function fct_ID_Li and

{\dot{ε}}_{L i}

= -1

Loading use loading function fct_ID_Li and

{\dot{ε}}_{L i}

Tension
Unloading use unloading function fct_ID_unL
Compression
Unloading use loading fct_ID_Li and ${\dot{ε}}_{L i}$

Viscous (Rate) Effects

Strain rate effects can be modeled by including engineering stress strain test data at different strain rates. This can be easier than calculating viscous parameters for traditional hyperelastic material models.

Conclusion

Make sure to use the material law that best fits the test data available.

For example, if minimum test data is available and the strains are not too large, the LAW42 Neo-Hookean model could be used. If the loading state is known, then it is important to have test data that represents that stress state and make sure the material model fits that test data.

References

¹ Ogden, R. W., and Non-linear Elastic Deformations. "Ellis Horwood." New York (1984)

² Miller, Kurt. "Testing Elastomers for Hyperelastic Material Models in Finite Element Analysis" Axel Products, Inc., Ann Arbor, MI (2017). Last modified April 5, 2017

http://www.axelproducts.com/downloads/TestingForHyperelastic.pdf

³ Axel Products, Inc. "Compression or Biaxial Extension", Ann Arbor, MI (2017). Last modified November 12, 2008

http://www.axelproducts.com/downloads/CompressionOrBiax.pdf

⁴ Kolling, S., P. A. Du Bois, D. J. Benson, and W. W. Feng. "A tabulated formulation of hyperelasticity with rate effects and damage." Computational Mechanics 40, no. 5 (2007): 885-899