Crash II

4. LS-DYNA Anwenderforum, Bamberg 2005

A Simplified Rubber Model with Damage

a

Stefan Kolling , Paul A. Du Bois & Dave J. Benson

DaimlerChrysler AG, EP/SPB, HPC X411, 71059 Sindelfingen, Germany

E-Mail: [email protected]
b

Consulting Engineer, Freiligrathstr. 6, 63071 Offenbach, Germany

E-Mail: [email protected]

University of California, Dept. of Mechanical and Aerospace Engineering, San Diego, USA
E-Mail: [email protected]

Abstract:
Simulation of rubber-like materials is usually based on hyperelasticity. If strain-rate dependency has to
be considered viscous dampers also have to be taken into account in the rheological model. A
disadvantage of such a material model is time-consuming parameter identification associated with the
damping constants. With MAT_SIMPLIFIED_RUBBER (Material no. 181), a material law is
implemented in LS-DYNA which allows fast generating of input data based on uniaxial static and
dynamic tensile tests at different strain rates. However, unloading, i.e. forming of a hysteresis, cannot
be modeled easily using Mat181. Thus, an extension of Mat181 based on a damage formulation was
desirable.
In this paper, we show the theoretical background and algorithmic setup of our model which has been
implemented as MAT_SIMPLIFIED_RUBBER_WITH_DAMAGE (Material no. 183) in LS-DYNA. As an
application, the validation of a soft and a hard rubber under loading and subsequent unloading at
different strain rates is shown.

Keywords:
Material Modelling, Rubber, Hyperelasticity, Strain-rate Dependency, Elastic Damage, Explicit Finite
Element Method

2005 Copyright by DYNAmore GmbH

B - II - 1

Crash II

4. LS-DYNA Anwenderforum, Bamberg 2005

Material law formulation

1.1

Basic equations

The most straight forward generalization of Hookes law to large displacements and large
deformations is hyperelasticity, see [9], [10]. A hyperelastic material is path independent and allows to
calculate the second Piola-Kirchhoff stress tensor S = 2W / C from a derivative of the energy
functional

W = W ( C) with respect to the components of the right Cauchy-Green strain tensor

C = FT F , where F = Grad x is the deformation gradient. In LS-DYNA we can distinguish two

families of hyperelastic materials. The first one is based on an energy functional expressed in the
invariants of the right Cauchy-Green tensor: W

I C = 1 : C = t r C , II C =
S=2

= W (I C , II C , III C ) . The invariants of C are given by

1 2
(I C C : C ) and III C = det C . Then, the derivative yields
2

W
W
W
III C C1 .
1+2
(I C 1 C ) + 2
I C
II C
III C

(1)
1

The Cauchy stress can now be obtained by forming = J FSF , where J = det F is the relative
volume.
The second family of hyperelastic materials is formulated in terms of principle stretch ratios. Therefore,
it is instructive to rewrite the former expressions in terms of principal stretches i . After a
T

F = RU , where RT R = 1 and U = UT , the invariants are given by

I C = + + 32 , II C = 12 22 + 22 32 + 32 12 , III C = 12 22 32 = J 2 and the Cauchy stress
and the principal engineering stress can be obtained as
decomposition
2
1

i =

2
2

1 W
j k i

j k i =: i =

W
.
i

(2)

An overview of hyperelastic laws which are implemented in LS-DYNA is given in Table 1, see also
references [1]-[6] and [8].

Law
7
2
33
27
38
83
127
181
183

Keyword
MAT_BLATZ-KO_RUBBER
MAT_ORTHOTROPIC_ELASTIC
MAT_FRAZER-NASH_RUBBER
MAT_MOONEY-RIVLIN_RUBBER
MAT_BLATZ-KO_FOAM
MAT_FU_CHANG_FOAM
MAT_ARRUDA_BOYCE_RUBBER
MAT_SIMPLIFIED_RUBBER
MAT_SIMPLIFIED_RUBBER_WITH_DAMAGE
Table 1: Overview of hyperelastic materials in LS-DYNA.

A hyperelastic law formulated in terms of principle stretches is given by Ogden:

W =

i =1

B - II - 2

j *
i 1 + K (J 1 ln J )
j =1 j
n

i =

n
p =1

p
J

*i

3
k =1

2005 Copyright by DYNAmore GmbH

* p

+K

J 1
.
J

(3)

Crash II

4. LS-DYNA Anwenderforum, Bamberg 2005

Here,

J = 12 3 and i* = iJ 1/ 3 . Note that the penalty term corresponds to

are non-integer,

a purely hydrostatic stress and, therefore, the deviatoric part of the Ogden functional generates a zero
stress in the undeformed configuration. In LS-DYNA, a tabulated version of the Ogden material law
has been implemented with the material law MAT_SIMPLIFIED_RUBBER. Here, the Ogden
functional (3) is internally determined from the uniaxial engineering stress-strain curve by defining a
tabulated function of the principal stretch ratio as follows [7]:

f ( ) =
1.2

* p

p =1

i =

1
1
f (i )
J
3

3
j =1

f ( j ) + K

J 1
.
J

(4)

Damage formulation

In rubbers the measured quasistatic loading and unloading paths are not necessarily identical.
Consequently the material cannot be considered to be hyperelastic during unloading and subsequent
reloading: the rubber has a path-dependent behaviour and a one-to-one correspondence between
stress and strain no longer exists. The current development uses a damage formulation to simulate
the rubber behaviour under cyclic loading. The implementation is limited at first to the incompressible
Ogden functional for the simulation of natural rubbers. The Ogden functional is written as

W = 1 d

W0
W0,max

i =1

j *
i 1 + U (J ) .
j =1 j
n

(5)

As can be seen in the formula above, the damage parameter d is considered to be a function of the
elastic deviatoric energy and the maximum of the elastic deviatoric energy over the deformation path :

W0 =

3
i =1

j *
i 1
j =1 j
n

W0,max = max(W0 ,W0,max )

W0
0
1
W0,max

(6)

0 d 1

The principal true stresses can no longer be computed directly from the energy functional since due to
the damage the material has become path-dependent and a one-to-one relationship between stress
and strain no longer exists. However by using the second law of thermodynamics an expression for
the principal true stresses can still be obtained:

W = 1 d

W0
W0,max

3
i =1

j *
i 1 + U ( J )
j =1 j
n

1 W
j k i
1 W0
1 U
i = (1 d )
+
j k i j k i

(7)

It is thus sufficient to compute the undamaged stress values and multiply the deviatoric part by (1-d).
The function d is tabulated and only the abscissa needs to be computed in order to perform the table
lookup. For more details concerning the damage formulation we refer to [13].

Numerical treatment

2.1

Implementation

Our implementation is a generalisation of the incompressible module in MAT_181 [7] including a

damage formulation to simulate the unloading of rubber materials. Hysteresis could only be simulated
in MAT_181 by using rate effects to fit the unloading curve. This can be cumbersome. Specific to our
formulation is that the function d of the energy function will be tabulated in such a way that a
measured unloading path from testing is exactly reproduced in the simulation. It is therefore necessary
to provide a closed-loop measurement in terms of engineering stress/strain resulting from a uniaxial
tensile (or compressive) test in input (see Figure 1).

2005 Copyright by DYNAmore GmbH

B - II - 3

Crash II

4. LS-DYNA Anwenderforum, Bamberg 2005

0
max
LCD

ULCD
UL

max

Figure 1: Definition of load curve

In Figure 1, LCD is a load curve giving the loading branch of the closed loop and ULCD is a load curve
giving the unloading branch. It is important that both load curves form a closed loop: maximum stress
and maximum strain must be identical for both load curves. In our current implementation, the
unloading curve is defined in tension as well as in compression resulting in 2 closed loops. This allows
defining a different damage evolution in tension and compression if corresponding test data are
available. If no damage is considered (ULCD=0) the original incompressible version of MAT_181 is
recovered. The function d is now internally computed from the input curves LCD and ULCD, here we
make the assumption of incompressibility: J=1 and U=0:
max

max

W0,max = V0 0,lcd d 0 =

0,lcd d 0

W0 ( 0 ) = 0,lcd d 0

(8)

W0 ( 0 )
( )
= 1 ulcd 0
W0,max
lcd ( 0 )

This allows to tabulate the damage parameter d as a function of the energy ratio. This is done during
the initialisation of the problem. For every explicit timestep, the energy ratio must be computed. For
that purpose we determine the uniaxial energy function from the quasistatic load curve (use the stressstrain curve from the table with lowest strain rate):

0 ( 0 ) is tabulated in input
= 0 +1
0

(9)

Wu ( ) = 0 d 0 = 0 d = Wu ( 0 )
This gives the energy function value for any longitudinal strain under uniaxial loading (tension and
compression must be covered) Now we specifically evaluate the deviatoric part of the Ogden function

W =

3
i =1

j *
i 1 + U ( J )
j =1 j
n

(10)

for uniaxial load conditions. With

j = k i 1 / 2 , J 1, *j = *k *i 1/ 2
we obtain

B - II - 4

2005 Copyright by DYNAmore GmbH

(11)

Crash II

4. LS-DYNA Anwenderforum, Bamberg 2005

W0,u (i ) =

n
j *
*
j
i 1 + 2
i

j =1
j =1
j
j
n

/2

(12)

and

W0,u i1/ 2 =

j *
i
j =1 j
n

/2

1 + 2

j *
i
j =1 j
n

/4

1 .

(13)

We can now define and compute an Ogden function for the energy analogously to the procedure
followed in [7] for the stress:

j *
i 1
j =1 j

( )

g *i =

( )

(14)

( )

g *i = Wu (i ) 2Wu i 1/ 2 + 4Wu 1i / 4 ...

2.2

Algorithmic setup

Using the equations derived in the last section, the damage algorithm then becomes as follows:
Step 1: calculate principal stretch ratios and engineering strains:

0i = i 1

(15)

Step 2: evaluate the Ogden functions for every stretch ratio:

( )

f = i 0 ( 0i ) +
*
i

( )

g *i = Wu ( 0i ) +

n =1

n =1

(i 1/ 2 ) 0 i(1/ 2 ) 1
n

( 2)n Wu ((i 1/ 2 )

(16)

for _ : _ (i 1/ 2 ) 1 0.01
n

Step 3: calculate undamaged true stresses and deviatoric energies:

i =

1
1
f (*i )
J
3

W0 =

3
i =1

3
j =1

f (*j ) + K

J 1
J

g (*i )

(17)

Step 4: compute the damage from the tabulated function and scale the stresses:

W0, max = max(W0, max ,W0 )

d =d

W0
W0,max

(18)

i = (1 d )( i + p ) p
3

Validation

Two rubber materials have been tested experimentally at the Ernst Mach Institute (EMI) in
Freiburg/Germany for validation of the numerical model. The experimental setup depends on the load

2005 Copyright by DYNAmore GmbH

B - II - 5

Crash II

4. LS-DYNA Anwenderforum, Bamberg 2005

direction. For compression tests, the setup consists of two pressure plates with the cube-like rubber
specimen (6x6x6mm) in between. The lower pressure plate is supported by a load cell. For tensile
tests, the specimen is fixed additionally by gluing to the pressure plates, see left hand side of Figure 2.
Compression tests with unloading on rubber cubes were performed for a hard (shore 70) and a weak
(shore 55) rubber, see right hand side of Figure 2. Moreover, dynamic tensile tests at different strain
rates (0.01/s, 1/s and 100/s) without unloading were also realised for validation.
a) Experimental setup

b) Compression test with unloading

load

pressure plate

test specimen

adhesive
bonding

stress 10-3 [GPa]

load

A hard rubber
B soft rubber

load cell
support

strain [-]

Figure 2: Experimental setup and data of compression test

In Figure 3, the results of the simulation are depicted. Since test data of the tension test with unloading
was not available, damage is assumed similar to compressive side. As can be seen, the experimental
data can be fitted exactly. As input data, the closed loop of the hysteresis in Figure 2 has been used
directly in Mat_183.

stress 10-3 [GPa]

a) Compression test

b) Compression and tensile test

A simulation
B test loading
C test unloading

A simulation tension
B simulation compression
C test tension loading
D tension unloading (hyp)
E test compression loading
F test compression unloading

strain [-]

strain [-]

Figure 3: Simulation of loading and unloading (hard rubber)

In Figure 4, the simulation of the compressive test with unloading for the hard rubber is shown.
Additionally, we chose different maximum compressions to demonstrate how the material law works if
we leave the range of defined data points. As input data, we defined the closed loop (LCD and ULCD
in Figure 1) given for maximum compression of 47.25%. In Figure 4a), we stay below this maximum
compression and go up to 40%. The loading path can be simulated, of course, exactly. However, the
unloading is an assumption: it is affine to the unloading path for 47.25% and after a certain time it is
even congruent. Although it is purely an assumption, from an engineering point of view, we detect a
healthy numerical behaviour of the material response. For the maximum compression of 47.25% in
Figure 4b), we obtain the exact curve as given in the input data. What happens for higher
compressions can be seen in Figure 4c) and d). The loading path is given as defined in the input
according to extrapolation of the curves. The unloading path is, again, given numerically by the
definition of the damage in the closed loop.

B - II - 6

2005 Copyright by DYNAmore GmbH

Crash II

4. LS-DYNA Anwenderforum, Bamberg 2005

stress 10-3 [GPa]

a) max compression 40%

A simulation
B test loading
C test unloading

b) max compression 47.25%

A simulation
B test loading
C test unloading

strain [-]

stress [GPa]

c) max compression 60%

A simulation
B test loading
C test unloading

d) max compression 80%

A simulation
B test loading
C test unloading

strain [-]

Figure 4: Compressive test (hard rubber) with unloading for different maximum compressions
Next, we compare the material laws MAT_SIMPLIFIED_RUBBER and its extension with damage
MAT_SIMPLIFIED_RUBBER_WITH_DAMAGE. First some remarks concerning LS-DYNA 9.71 and
higher:
MAT_SIMPLIFIED_RUBBER_WITH_DAMAGE is available also for strain rate dependent
materials, i.e. the load curve can be replaced by a table definition referring to the load curves
at different strain rates.
MAT_SIMPLIFIED_RUBBER has an extension with a damage model. This is, however,
intended for the simulation of failure rather than for the simulation of unloading!
a) Strain rate dependency

b) Simulation of hysteresis

stress 10-3 [GPa]

A MAT183
B MAT181

A MAT183 loading
B MAT183 unloading
C MAT181 0.01 /s
D MAT181 1 /s
E MAT181 100 /s

strain [-]

Figure 5: Comparison MAT_183 and MAT_181

For the simulation of a hysteresis loop, MAT_SIMPLIFIED_RUBBER can be used if the unloading is
determined by the definition of the strain rate. In Figure 5a), the exact curves simulating the
experiment are depicted. In Figure 5b), the strain rate sensitivity is artificially increased to obtain a

2005 Copyright by DYNAmore GmbH

B - II - 7

Crash II

4. LS-DYNA Anwenderforum, Bamberg 2005

hysteresis loop in the quasi static case. Note that in this case, the material response for higher strain
rates may be not exact anymore.
In our last validation example, we simulate the strain rate dependency of the soft rubber again in
single element tensile tests with unloading. In Figure 6, we compare the material response of MAT181
and MAT183. For each material two tensile loads (with and without rate effects) are simulated
resulting in the same maximum deformation. We then perform the analyses twice to illustrate the effect
of different settings of the unloading flag TENSION. For TENSION=-1 (Figure 6a), rate effects are
considered in the loading phase only. Consequently, the unloading path follows the quasistatic curve
in MAT_181 and dynamic and quasistatic unloading paths are identical in MAT_183. For TENSION =1
(Figure 6b) rate effects are considered in loading and unloading phase. Loading and unloading path
are then identical for MAT_181 resulting in a potentially unstable material model. With MAT_183 the
unloading paths show rate dependency and (due to the damage formulation) are always below the
loading path. The advantages of MAT_183 upon unloading are clear:

rate effects can be considered

damage formulation ensures energy dissipation and hence numerical stability

stress 10-3 [GPa]

A MAT181 qs
B MAT181 dyn
C MAT183 qs
DMAT183 dyn

b) Unloading with strain rate effect

A MAT183 qs
B MAT183 dyn
C MAT181 dyn
DMAT181 qs

stress 10-3 [GPa]

a) Unloading without strain rate effect

time [ms]

time [ms]

Figure 6: Comparison MAT_183 and MAT_181

In Figure 7, tensile tests with and without rate effects are again considered for MAT_183 only. This
time different maximum deflections (corresponding to different strain rates upon loading) are applied. It
is shown again that for TENSION=-1, quasistatic and dynamic unloading path coincide whereas for
TENSION=1, rate effects are clearly present also upon unloading.

stress 10-3 [GPa]

A dyn max strain 20%

B qs max strain 20%
C qs max strain 40%
D dyn max strain 40%

b) Unloading with strain rate effect

A dyn max strain 20%
B dyn max strain 40%
C qs max strain 20%
D qs max strain 40%

stress 10-3 [GPa]

a) Unloading without strain rate effect

time [ms]

time [ms]

Figure 7: Strain-rate dependency using MAT_183

B - II - 8

2005 Copyright by DYNAmore GmbH

4. LS-DYNA Anwenderforum, Bamberg 2005

Crash II

Summary

With the presented material formulation, exact simulation of test data for different strain rates in
tension and compression can be achieved without any parameter fitting. Furthermore, stable and
realistic unloading behaviour with energy dissipation is obtained based on a solid theoretical basis.
The implementation has been done for solid and shell elements and is currently limited to the nearly
incompressible case (Ogden functional). In the near future the damage formulation will be made
available for the general hyperelastic formulation (Hill functional) also.

References

[1]

J.O. Hallquist: LS-DYNA, Theoretical Manual, Livermore Software Technology Corporation,

Report 1018, 1991.
P.A. Du Bois: Crashworthiness Engineering Course Notes, Livermore Software Technology
Corporation, 2004.
P.J. Blatz & W.L. Ko: Application of finite elastic theory to the deformation of rubbery materials.
Trans. Soc. Rheol., 6: 223-251, 1962.
M. Mooney: A theory of large elastic deformations. J. Appl. Physics, 11:582-592, 1940.
R.S. Rivlin: Large elastic deformations of isotropic materials. Proc. Roy. Soc. London, 241: 379397, 1948.
R.W. Ogden: Large deformation isotropic elasticity: on the correlation of theory and experiment
for incompressible rubberlike solids. Proc. Roy. Soc. London, 326: 565-584, 1972.
P.A. Du Bois: A simplified approach for the simulation of rubber-like materials under dynamic
th
loading. 4 European LS-DYNA Users Conference, pp. D-I-31/46, 2003.
P.A. Du Bois, W. Fassnacht & S. Kolling: General aspects of material models in LS-DYNA. LSDYNA Forum, Bad Mergentheim, Germany 2002, V2:1-55.
R. Hill: Aspects of invariance in solid mechanics. Adv. Appl. Mech. 18: 175, 1978.
B. Storkers: On material representation and constitutive branching in finite compressible
elasticity. J. Mech. Phys. Solids, 34(2): 125-145, 1986.
P.A. Du Bois, S. Kolling, M. Koesters & T. Frank: Material modeling of polymeric materials in
rd
crashworthiness analysis. 3 Workshop for Material and Structural Behaviour at Crash
Processes (crashMAT), Freiburg, Germany 2004.
P. A. Du Bois, S. Kolling & T. Frank: Material behaviour of polymers under impact loading.
International Symposium on crashworthiness of light-weight automotive structures, Trondheim,
Norway, 2004. To appear at the International Journal of Impact Mechanics.
Timmel, M.; Kaliske, M.; Kolling, S.: Modellierung gummiartiger Materialien bei dynamischer
Beanspruchung. LS-DYNA Forum, Bamberg, Germany 2004, C-I-1/11.

[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]

2005 Copyright by DYNAmore GmbH

B - II - 9

Crash II

B - II - 10

4. LS-DYNA Anwenderforum, Bamberg 2005

2005 Copyright by DYNAmore GmbH