Ogden Materials

LAW42 uses the following strain energy density representation of the Ogden material model.(3) $$W\left({\lambda }_1,{\lambda }_2,{\lambda }_3\right)={\displaystyle \sum _{p=1}^5\frac{{\mu }_p}{{\alpha }_p}\left({\overline{\lambda }}_1{}^{{\alpha }_p}+{\overline{\lambda }}_2{}^{{\alpha }_p}+{\overline{\lambda }}_3{}^{{\alpha }_p}-3\right)}+\frac{K}{2}{\left(J-1\right)}^2$$

Where,

$$W$$

Strain energy density

${\lambda }_i$

i^th principal engineering stretch

$$J$$

Relative volume defined as: $$J={\lambda }_1\cdot {\lambda }_2\cdot {\lambda }_3=\frac{{\rho }_0}{\rho }$$

${\overline{\lambda }}_i=J^{-\frac{1}{3}}{\lambda }_i$

Deviatoric stretch

${\alpha }_p$ and ${\mu }_p$

Material constants coefficient pairs.

Up to 5 material constant pairs can be defined.

The initial shear modulus $\mu$ and bulk modulus ($$K$$) are given by:(4) $\mu =\frac{{\displaystyle \sum _{p=1}^5{\mu }_p\cdot {\alpha }_p}}{2}$

and(5) $K=\mu \cdot \frac{2\left(1+\nu \right)}{3\left(1-2\nu \right)}$

Parameters ${\alpha }_p$ and ${\mu }_p$ must be chosen so that initial shear modulus is:(6) $\mu =\frac{{\displaystyle \sum _{p=1}^5{\mu }_p\cdot {\alpha }_p}}{2}>0$

For material stability, it is required that each material constant pair (7) $${\mu }_p\cdot {\alpha }_p>0$$

A simple case of the Ogden material model is the Neo-Hooken model represented using the following equation for the strain energy density function:(8) $W=C_{10}\left(I_1-3\right)$

Where,

$I_1$

The first invariants of the right Cauchy-Green Tensor

$C_{10}$

Material constant

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:(9) $W=C_{10}\left(I_1-3\right)+C_{01}\left(I_2-3\right)$

Where,

$I_1$ and $I_2$

The first and second invariants of the right Cauchy-Green Tensor

$C_{10}$ and $C_{01}$

Material constants

If a material is truly incompressible, then $\nu =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. (10) $K=\mu \cdot \frac{2\left(1+\nu \right)}{3\left(1-2\nu \right)}=\mu \cdot \frac{2\left(1+\nu \right)}{3\left(1-2*0.5\right)}=\infty$

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

Figure 2.

In LAW42, material incompressibility is provided by using a penalty approach, which calculates the pressure proportional to a change in density.(11) $$P=K\cdot Fscal{e}_{blk}\cdot {\mathrm{f}}_{blk}\left(J\right)\cdot \left(J-1\right)$$

Where,

$$K$$

Bulk modulus

$$J=\frac{V}{V_0}=\frac{m{\rho }_0}{m\rho }=\frac{{\rho }_0}{\rho }$$

Relative volume which simplifies to relative density if mass is constant

$${\mathrm{f}}_{blk}\left(J\right)$$

Bulk coefficient scale factor versus relative volume function

$$Fscal{e}_{blk}$$

Abscissa scale factor for function $${\mathrm{f}}_{blk}\left(J\right)$$

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) then 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, i.e. 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 are modeled in LAW42 using a Maxwell model, which can be described in a simplified manner as a system of $\eta$ springs with stiffness' $$G_i$$ and dampers ${\eta }_i$:

Figure 4. Maxwell Model

The Maxwell model is represented using Prony series inputs ($G_i,{\tau }_i$). The hyperelastic initial shear modulus $\mu$ is the same as the long-term shear modulus $$G_{\infty }$$ in the Maxwell model, and ${\tau }_i$ is the relaxation time:(12) $${\tau }_i=\frac{{\eta }_i}{G_i}$$

The hyperelastic behavior in this material law is defined using the following strain energy density function:(13) $W({\lambda }_1,{\lambda }_2,{\lambda }_3)={\displaystyle \sum _{i=1}^N\frac{2{\mu }_i}{{\alpha }_i{}^2}\left({\lambda }_1{}^{{\alpha }_i}+{\lambda }_2{}^{{\alpha }_i}+{\lambda }_3{}^{{\alpha }_i}-3+\frac{1}{\beta }(J^{-{\alpha }_i\beta }-1)\right)}$

Where,

$W$

Strain energy density

${\lambda }_i$

$$i_{th}$$ principal stretch

$J$

Relative volume defined in Equation 13

$$\beta =\frac{\nu }{\left(1-2\nu \right)}$$

${\alpha }_i$ and ${\mu }_i$

Material constants coefficient pairs.

Up to 5 material constant pairs can be defined.

Poisson’s ratio must be $0<\nu <0.5$. This law can be used to model compressible or sometimes called hyperfoam materials by defining a low Poisson’s ratio value.

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 ${\eta }_i$ .

Figure 5. Maxwell Model

The Maxwell model is represented using a Prony series with inputs. The initial shear modulus is:(15) $$G_0={\displaystyle \sum _{i=1}^N{\mu }_i}$$

The sum of ${\mu }_i$ should be greater than 0.(16) $G_0=G_{\infty }+{\displaystyle \sum _i{G}_i}$

The stiffness ratio is:(17) ${\gamma }_{\infty }=\frac{G_{\infty }}{G_0}=1-{\displaystyle \sum _i{\gamma }_i}$ (18) ${\gamma }_i=\frac{G_i}{G_0}$

With, (19) ${\gamma }_i\in \left[0,1\right],{\displaystyle \sum _i{\gamma }_i}<1$

and the ground shear modulus(20) $G_0=G_{\infty }+{\displaystyle \sum _i{G}_i}$

The relative time, ${\tau }_i$ must be positive:(21) ${\tau }_i=\frac{{\eta }_i}{G_i}$

The strain energy density formulation used depends on the law_ID:

• law_ID = 1 (Ogden law):(22) $$W\left({\lambda }_1,{\lambda }_2,{\lambda }_3\right)={\displaystyle \sum _{p=1}^5\frac{{\mu }_p}{{\alpha }_p}\left({\overline{\lambda }}_1{}^{{\alpha }_p}+{\overline{\lambda }}_2{}^{{\alpha }_p}+{\overline{\lambda }}_3{}^{{\alpha }_p}-3\right)}+\frac{K}{2}{\left(J-1\right)}^2$$
• law_ID = 2 (Mooney-Rivlin law):(23) $W=C_{10}\left(I_1-3\right)+C_{01}\left(I_2-3\right)$

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.

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.(24) $${\epsilon }_c=\frac{1}{{\left({\epsilon }_b+1\right)}^2}-1$$ (25) ${\sigma }_c={\sigma }_b{\left(1+{\epsilon }_b\right)}^3$

Where,

${\epsilon }_c$

Uniaxial engineering compressive stain

${\epsilon }_b$

Equal biaxial engineering tension strain

${\sigma }_c$

Uniaxial engineering compressive stress

${\sigma }_b$

Equal biaxial engineering tension stress

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.(26) $W={\displaystyle \sum _{i=1}^N\frac{2{\mu }_i}{{\alpha }_i{}^2}}\left({\overline{\lambda }}_1{}^{{\alpha }_i}+{\overline{\lambda }}_2{}^{{\alpha }_i}+{\overline{\lambda }}_3{}^{{\alpha }_i}-3\right)+{\displaystyle \sum _{i=1}^N\frac{1}{D_i}{\left(J-1\right)}^{2i}}$

Where,

$W$

Strain energy density

$N$

Number of material constants ${\alpha }_i$, ${\mu }_i$ and $D_i$

${\overline{\lambda }}_i=J^{-\frac{1}{3}}{\lambda }_i$

Deviatoric stretch

$J$

Relative volume as defined in Equation 2

The initial shear modulus:(27) $\mu ={\displaystyle \sum _{i=1}^N{\mu }_i}$

The Bulk Modulus is calculated as $K=\frac{2}{D_1}$ based on these rules:

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

The Drücker 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.(29) ${\displaystyle \sum _{i=1}^{3}d{\tau }_{i}d{\epsilon }_i>0}$

With the change in logarithmic strain(30) $d{\epsilon }_i=\frac{d{\lambda }_i}{{\lambda }_i}$

$d{\tau }_i=J\cdot d{\sigma }_i$

The change of Kirchhoff stress

$d\tau =D:d\epsilon$

Relationship between Kirchhoff stress and logarithmic strain

The Drücker stability condition will be:(31) ${\displaystyle \sum d\epsilon :D:d\epsilon >0}$

Here $D$ is tangential material stiffness matrix and it is also the slope of stress-strain curve:(32) $D=\left[\begin{array}{ccc}{D}_{11}& D_{12}& D_{13}\\ D_{21}& D_{22}& D_{23}\\ D_{31}& D_{32}& D_{33}\end{array}\right]$

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:(33) $I_1=tr\left(D\right)=D_{11}+D_{22}+D_{33}>0$ (34) $I_2=D_{11}{D}_{22}+D_{22}{D}_{33}+D_{33}{D}_{11}-D_{23}{}^2-D_{13}{}^2-D_{12}{}^2>0$ (35) $I_3=\det \left(D\right)>0$

The Kirchhoff stress for Ogden model is:(36) ${\tau }_i={\displaystyle \sum _p{\mu }_p\left[{\overline{\lambda }}_i{}^{{\alpha }_p}-\frac{1}{3}\left({\overline{\lambda }}_1{}^{{\alpha }_p}+{\overline{\lambda }}_2{}^{{\alpha }_p}+{\overline{\lambda }}_3{}^{{\alpha }_p}\right)\right]+K\left(J^2-J\right)}$

Since $D_{ij}=\frac{\partial {\tau }_i}{\partial {\lambda }_j}$, then for a given Ogden parameter ${\alpha }_p,{\mu }_p$ with conditions $I_1>0,I_2>0\text{ and }{I}_3>0$, the strain range of material in Drücker stability could then be calculated.

Figure 6.

Then the Drücker 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

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. ⁴(37) $W=\underset{deviatoricpart}{\underbrace{{\displaystyle \sum _{i=1}^3{\displaystyle \sum _{j=1}^m\frac{{\mu }_j}{{\alpha }_j}}}\left({\overline{\lambda }}_i{}^{{\alpha }_j}-1\right)}}+\underset{sphericalpart}{\underbrace{K\left(J-1-\ln J\right)}}$

If using the damage model, the loading curves are used for both loading and unloading and the unloading stress tensor is reduced by:(38) $$\sigma =\left(1-D\right)\sigma$$

with(39) $$D=\left(1-Hys\right)\left(1-{\left(\frac{W_{cur}}{W_{\max }}\right)}^{Shape}\right)$$

Where,

$$W_{cur}$$

Current energy

$$W_{\max }$$

Maximum energy corresponding to the quasi-static behavior

$$Hys$$ and $$Shape$$

Input by user

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

Tension

Loading and Unloading

= 0

Figure 7.

Loading use loading function fct_ID_Li and ${\dot{\epsilon }}_{Li}$

Unloading use unloading function fct_ID_unL

= 1

Figure 8.

Loading and unloading all use loading function fct_ID_Li and ${\dot{\epsilon }}_{Li}$

= -1

Figure 9.

Loading use loading function fct_ID_Li and ${\dot{\epsilon }}_{Li}$

• Tension

  Unloading use unloading function fct_ID_unL

• Compression

  Unloading use loading fct_ID_Li and ${\dot{\epsilon }}_{Li}$