Figure 1 The energy dissipation is a primary feature to be predicted in many engineering applications such
as the estimate of the rolling resistance of tires and hysteretic losses in biological tissues and the
design of vibration absorbers.
Outline
In this work, the behavior of carbon blacklled rubber is thoroughly analyzed with the inten
tion of developing a constitutive model able to reproduce both static and dynamic material
responses. However, due to the similarities between rubber and biological tissues, most of the
results presented can be applied in a wider context.
The objectives of this thesis can be summarized as follows:
to dene adhoc experimental procedures for the evaluation of viscoelastic constitutive
parameters;
to formulate a constitutive model able to predict hysteresis losses at nite strain and
low frequencies;
to solve the numerical problems related to the identication of viscoelastic constitutive
parameters.
Several nonlinear viscoelastic models have been examined thoroughly and for each of them
advantages and disadvantages are highlighted.
A series of experiments concerning both static and dynamic tests were performed aimed
at measuring all the relevant nonlinear eects. Temperature and strain rate dependencies
were investigated and discussed. The standard methodology was applied to perform both
tensile and compressive quasistatic tests. Some shortcomings of this procedure, resulting
in a unreliable stressstrain constitutive curve around the undeformed conguration, were
identied. This lead to the design of a nonstandard cylindrical specimen able to bear both
tensile and compressive loadings. Consequently, the inuence of the shape factor was removed
and the same boundary conditions, in tension and compression, was applied. This allowed
the stiness around the undeformed conguration to be evaluated in detail.
The quasistatic experimental results also allowed the inuence of the Mullins eect on
the quasistatic response to be investigated: during the loading cycles, there is a signicant
reduction in the stress at a given level of strain, which is a consequence of the internal material
rearrangement, i.e., the Mullins eect. This damage phenomenon is sometimes reported to
induce transverse isotropy in the material, which is usually assumed to be isotropic. The
results of the experiments have claried this issue.
vii
To measure the dynamic properties of the rubber compound, creep, relaxation and cyclic
tests at dierent strain rates were performed. While relaxation and creep experiments incor
porate the longterm material response, the stress arising from loading/unloading cycles at
dierent strainrates involves the shortest characteristic times, thus the highest characteristic
frequencies. This methodology allows the accurate reproduction of the actual material opera
tive range. In addition, the incorporation of creep tests into loading/unloading cycles proved
to be the most eective methodology to evaluate nonlinear viscoelastic parameters.
Finally, to allow a comparison with the literature results, the usual harmonic testing
procedure was also applied.
Thereafter, all the constitutive equations under consideration have been analyzed thor
oughly in terms of their capabilities of describing the collated experimental data.
The material coecients were initially identied by a procedure which relies on the sep
arate identication of the instantaneous (elastic) term and of the dissipative (inelastic)
part. By means of high strain rate loading path, the elastic moduli of the material were
identied. Thereafter, relaxation and oscillatory data was used to obtain the characteristic
times and the dissipative moduli. It has been proved that, to guarantee the wellconditioning
of the resulting leastsquares problem, i.e., to avoid relevant numerical error, relaxation tests
should be applied for the viscoelastic kernel identication rather than the standard oscillatory
tests. However, the results of this procedure were unsatisfactory for some of the models.
To overcome these limitations, a joint identication of the elastic and dissipative terms
was introduced. A common choice to deal with the numerical diculties, related to the tran
scendental dependence of the constitutive function upon the characteristic times, is to x
them apriori, e.g., one or two per decades of the experimental time range. This choice, while
often used by many authors, laed to unsatisfactory results when dealing with lled rubber. As
a matter of fact, it does not account for the wellknown property of carbon blacklled elas
tomers of having characteristic times very close to each other. As a consequence, an enhanced
iterative scheme, which actually allows a more accurate estimate of the characteristic times
clusters, has been introduced. Moreover, this iterative scheme deals with the requirement of
keeping the number of constitutive parameters to a minimum in order to avoid nonuniqueness
of their determination and to provide a clear physical interpretation for each of them.
The resulting identication problem requires the minimization of a nonlinear functional,
which was solved numerically. The penchant of local optimization algorithms to become
trapped in local minima on such landscapes required the adoption of a nonlocal optimization
strategy based on a Pattern search algorithm.
Although all the models under consideration were able to reproduce the relaxation or the
cyclic test, they failed to extrapolate the material behavior under dierent kinds of defor
mation. Furthermore, they encounter diculties in the prediction of the dynamic response
especially at low frequencies and low strain rates. This should be considered an important
drawback as the frequencies up to 1015 Hz are often the most signicant in many applications,
e.g., rolling tire at 100 km/h or human heart rate.
The extension of the denition of dynamic moduli to the nonlinear case provided new
insights into the dynamic response of a broad class of constitutive equations. In view of
the novel denition, it has been possible to prove that every nonlinear viscoelastic equation
cannot account for the linear frequency dependence of the storage modulus observed at low
frequencies. Therefore, a onedimensional constitutive model, based on hysteretic damping
and formulated in the frequency domain, has been proposed. The corresponding timedomain
representation of the model required the introduction of the Hilbert transform.
viii Chapter 0. Summary
Original Contributions
The main original contributions to the existing literature could be summarized as follows.
New experimental evidence concerning the deformation of lled rubber at nite strain
and nite strain rate.
New experimental evidence concerning the transverse isotropy induced by the Mullins
eect.
A critical analysis of the literature proposals for nonlinear viscoelastic models.
A study of the numerical conditioning of the minimization problem resulting from the
identication of nonlinear viscoelastic constitutive equations.
A novel iterative identication technique, which has allowed the characteristic times of
the material to be evaluated more accurately.
An experimental/numerical comparison of the dierent nonlinear viscoelastic models
through the proposed identication technique.
The extension of the denition of dynamic moduli to the nonlinear case.
The proof of the incapability of most literature proposals to describe the dynamic moduli
at low frequencies.
A onedimensional constitutive model based on hysteretic damping able to t the dy
namic behavior of lled rubber at low frequencies.
Further Developments
Having considered all the outlined results, the following points are still under development.
The extension of the onedimensional model proposed to the threedimensional case.
The introduction of a constitutive equation in the time domain, able to reproduce the
low frequency behavior of the dynamic moduli.
The introduction of transversely isotropic and anisotropic nonlinear viscoelastic consti
tutive models.
Experimental procedures to identify the anisotropic viscoelastic models.
Fiber reinforced polymer, collagen ber bundles, and human ligaments are materials which,
nowadays, play a signicant role in numerous engineering applications. All of them display
strong anisotropic response due to ber orientation. While the constitutive theory of isotropic
materials has reached a certain level of completion, substantial eort must be made to de
velop viscoelastic models for anisotropic materials at nite strain and to place them within a
consistent thermodynamic framework.
Further complications arise from the identication of these anisotropic nonlinear viscoelas
tic models, since a number of independent quasistatic and dynamic tests are necessary to
identify separately the material coecients. Very often, these experiments are dicult to
perform, particularly for in vivo softtissues. Therefore, the problem of identifying correctly
the material parameters is still a challenging engineering task.
Finally, nite element multiscale codes are increasingly used to simulate complex biological
tissues and rubber structures. In view of the signicant computational eort required by
ix
the simulation of an entire human organ or a vehicle tire, the computational eciency of
these numerical tools must a primary requirement. Parallel computing techniques should be
considered and new parallel algorithms based on GPU computing investigated.
Structure of the Thesis
Chapter 1. The behavior of carbon black lledrubber in relation to quasistatic and
dynamic responses is examined in detail. New insights regarding the material behavior
are provided by examining the results from the performed experiments.
Chapter 2. The main aspects of the nonlinear theory of elasticity are discussed.
Chapter 3. The main approaches followed to model nonlinear viscoelastic solids during
isothermal deformation are thoroughly described. Some of the most common integral
models are reviewed and advantages and disadvantages of each are highlighted. The
concept of dynamic moduli, introduced in linear viscoelasticity and referred to as storage
and loss moduli, is extended, in a consistent manner, to the nonlinear case. Finally, a
onedimensional model based on hysteretic damping is introduced.
Chapter 4. The standard identication procedure of the material parameters for a
nonlinear viscoelastic (NLV) constitutive equation is analyzed in view of the collated
experimental results. The main feature of this approach are evaluated by considering
Fungs constitutive model. Thereafter, a joint identication of the elastic and of the
dissipative term is introduced.
Chapter 5. The behavior of isotropic, almostincompressible, nonlinear elastic and
viscoelastic materials is simulated by means of the ABAQUS FEA code. Simple de
formations are considered and the numerical results are compared with the analytical
solutions. Finally, shortcomings of the ABAQUS nite viscoelasticity model are high
lighted and discussed.
x Chapter 0. Summary
Chapter 1
Rubber Phenomenology
Chapter Outline. In this chapter the behavior of carbon blacklled rubber in relation to quasistatic and
dynamic responses is examined in detail. In particular, the main features of the microstructure of the material
and their inuence on the macromechanical response are highlighted. The eects of strain, strainrate and
temperature on the constitutive response are discussed. New insights regarding the material phenomenological
behavior are provided by examining the results of the experiments carried out in this work. Mullins and Payne
eects, which are peculiar in the behavior of lled elastomers, are reviewed and new results are shown.
2 Chapter 1. Rubber Phenomenology
1.1 Material Description 3
1.1 Material Description
Owing to its unique physical properties, rubber plays a keyrole in countless industrial appli
cations. Tires, vibration absorbers and shoe soles are only but a few of the myriad uses of
rubber in an industry which in 2009 had an estimated market value of two billion euro.
The term rubber is actually misleading: it is used both to indicate the material, technically
referred to as natural rubber, and the broad class of synthetic elastomers which share with
natural rubber some fundamental chemical properties. Indeed, the majority of rubber used
for industrial applications are synthetically produced and derived from petroleum
1
.
Rubber, or elastomer, has an internal structure which consists of exible, long chain
molecules that intertwine with each other and continually change contour due to thermal
agitation. Elastomers are polymers with long chains (Ferry, 1980). The morphology of an
elastomer can be described in terms of convolution, curls and kinks. Convolutions represent
the longrange contour of an entire molecular chain, which forms entanglements (knots). Curls
are shorter range molecular contours that develop between entanglements and crosslinks,
and kinks are molecular bonds within a curl. Each molecular bond has rotational freedom
that allows the direction of the chain molecule to change at every bond. Thus the entire
molecular chain can twist, spiral and tangle with itself or with adjacent chains. This basic
morphology is shared among all the fty thousand compounds used in the market today
and generically referred to by the term rubber. Despite this intricate internal structure, the
random orientation of the molecular chains results in a material which is externally isotropic
and homogeneous.
Prior of using, the neat elastomer is subjected to physical/chemical treatments to enhance
its mechanical properties. One of these treatments consists of the addition, through heating, of
sulfurbased curatives which create crosslinks among the macromolecules chains; this process
is commonly called vulcanization (see, e.g., Callister, 2007).
Figure 1.1 highlights the dierent behavior of a vulcanized and a nonvulcanized rub
ber specimen subjected to a tensile loading. Initially, both of the elastomers have a similar
intertwined internal structure. When stretched, the macromolecules of the nonvulcanized
compound disentangle themselves according to the direction of the applied force. This mi
crostructural change results in a more ordered internal state with a subsequent reduction in
entropy. Thereafter, the macrobrownian motions of the macromolecules cause the chains to
slide back, one onto each other, to the disordered state. Finally, once the external load is
removed (step d in the gure), each macromolecule maintains its state of maximum entropy.
Therefore, the initial overall shape is not recovered: all the energy externally supplied to
stretch the specimen is dissipated by the viscous friction among the macromolecules.
A dierent microstructural response occur during the deformation of the vulcanized spec
imen. Indeed, when subjected to an external traction, the molecular chains dispose parallel
to the macrodisplacement and, because of the crosslinks introduced
by the vulcanization, they cannot slide back to the initial disordered state. By removing
the external loading, the system tends towards the initial state of maximum entropy and
the specimen recover the initial length. In this case, the external supplied energy is totally
recovered.
The behavior of a real elastomer slightly diers from this simplied description. Indeed,
even if the elastomer is vulcanized, the macromolecules can partially slide one onto each other
with a dissipation of the mechanical energy.
At the end of the vulcanization process for some specic applications, such as in tires,
reinforcing ller, usually carbon black, is added to the compound. This carbon based curative
1
The uncontrolled growth of the petroleum price in 2008 has forced many producers to substitute partially
the synthetically produced compounds with natural rubber.
4 Chapter 1. Rubber Phenomenology
initial
state
a b c d
e f g h
stretched stretched
final
state
Figure 1.1 Eects of stretching on a nonvulcanized (above) and a vulcanized (below) elastomer.
lends to the material the black color typical of tires.
The tensile strength of rubber increases with increasing ller content up to a certain level.
Beyond this level, the tensile strength decreases with higher ller concentrations. Goldberg
et al. (1989) suggested that this is because high amounts of carbon black llers cause the
carbon black to agglomerate into large clusters and these clusters impart aws that can easily
create cracks and lead to a catastrophic failure. The quantity of ller present in the elastomer
is measured in phr, parts per hundred by weight of elastomer; the concentration at which
maximum tensile strength is obtained, varies with the type of carbon black. For carbon
black llers with smaller particle size, the maximum tensile strength is attained at lower
concentrations than those for large particle sized carbon black llers.
The resulting mechanical characteristics such as strength, tear and abrasion resistance,
along with stiness, considerably increase with respect to the neat elastomer. The addition
of ller contributes also to alter greatly the viscous behavior and temperature dependence.
For example, unlled elastomers exhibit a linear viscoelastic behavior for shear strains up to
20 % or more, while a carbon blacklled elastomer shows a pronounced nonlinear behavior
at shear strains as low as 0.5 % (Chazeau et al., 2000).
In the next section the standard phenomenology of carbon blacklled rubber will be pre
sented and the inuence on the constitutive response of temperature and ller concentration
will be discussed. Although the focus is on traditional vulcanized rubber, other thermoplas
tic elastomers show similar mechanical properties even if their chemical composition is quite
dierent. Moreover, from a macroscopic point of view, the behavior of such materials is very
close to the behavior of some biological soft tissues, such as ligaments and tendons, for what
concerns both their static and dynamic responses.
1.2 Standard Phenomenology 5
1.0 1.2 1.4 1.6 1.8 2.0
0.5
0.0
0.5
1.0
1.5
2.0
2.5
Nominal Strain
N
o
m
i
n
a
l
S
t
r
e
s
s
M
P
a
Shear
Tension
Figure 1.2 Experimental data on carbon blackreinforced styrene butadiene rubber for tensile (circle) and pure
shear (square) tests (Drozdov, 2007). The ratio of the tangent stiness around the undeformed
conguration, i.e., nominal strain equal to 1, is approximately equal to 3.
1.2 Standard Phenomenology
1.2.1 Quasistatic
The behavior of carbon blacklled elastomers can be primarily described as hyperelastic:
under static or quasistatic loading dissipative eects are negligible.
There have been numerous experimental studies addressing the response of rubber un
der quasistatic loading conditions, including uniaxial tension/compression, shear, equibiaxial
tension (Drozdov & Dorfmann, 2003; LarabaAbbes et al., 2003; Przybylo & Arruda, 1998;
Treloar, 2005). In all these experimental conditions, the resulting constitutive curves are
strongly nonlinear. However, constitutive nonlinearities coupled with heterogeneous strain
eld could lead to experimental results which are very dicult to analyze
1
. Thus, displace
ment elds leading to homogeneous deformation should be opted for. A typical example is
the equibiaxial (twodimensional) extension test which is preferred to the equivalent uniaxial
compression, because the diculties related to the bulging of the specimen under compressive
loading are avoided (Treloar, 2005).
The typical stressstrain constitutive curves of a carbon blacklled elastomer are shown
in Fig. 1.2 (Drozdov, 2007). The material is subjected to uniaxial tension/compression, and
pure shear. In the typical working range (0.8 2.0) the constitutive nonlinearities are
evident; indeed, as the breaking point is approached, the material stiness rapidly increases so
that the slope of the experimental curves begins to rise. As a consequence of the intertwining
internal structure, during compression, high levels of loading force are suddenly reached, i.e.,
the material is much stier with a nonsymmetric behavior between tensile and compressive
stresses.
From Fig. 1.2, it is evident that the shear modulus G around the undeformed conguration,
1
In recent years the use of digital image correlation techniques to evaluate heterogeneous strain elds
has spread rapidly (see, e.g., Chevalier et al., 2001; Sasso et al., 2008). However, for technical reasons, these
techniques are limited to experiments involving very low strain rate.
6 Chapter 1. Rubber Phenomenology
V
0
Figure 1.3 Volume dilatation for a rubber specimen undergoing a uniaxial tensile experiment (results provided
by Reichert et al., 1987). The volume change remains limited over a wide strain range.
i.e., nominal strain equal to 1, has a lower value compared to the Young modulus E in tensile
experiments. The ratio E/G is approximately equal to 3, which corresponds to a Poisson
function in the undeformed conguration equal to = 0.5, meaning that the material is
incompressible.
The incompressibility of carbon blacklled rubber has been conrmed by a number of
dierent researchers over the years (Bischo et al., 2001; MacKnight, 1966; Ogden, 1976; Penn,
1970; Reichert et al., 1987). Experiments by Reichert et al. (1987) in Fig. 1.3 show a limited
volume variation (V/V
0
0.01) at large strain ( 4) corroborating the incompressibility
constraint introduced in many constitutive equations (see also Mott & Roland, 2010; Mott
et al., 2008; Voinovich, 2010).
The eects upon the quasistatic response of an increasing quantity of reinforcing ller have
been studied and results have been provided in (Yeoh & Fleming, 1997) for pure shear tests
(see Fig. 1.4). The addition of carbon black produces higher value of the initial stiness (i.e.,
tangent modulus around the undeformed conguration) with respect to the neat elastomer,
while it makes the compound more sensitive to temperature variations. Indeed, the same
qualitative behavior has been reported whatever the loading conditions.
The inuence of the temperature on the stressstrain curve is shown in Fig. 1.5. At very
low temperatures, the polymer will behave like glass and exhibit a high modulus. As the
temperature is increased, the polymer will undergo a transition from a hard glassy state to
a soft rubbery state in which the modulus can be several orders of magnitude lower than
it was in the glassy state. The transition from glassy to rubbery behavior is continuous and
the transition zone is often referred to as the leathery zone. The onset temperature of the
transition zone, moving from glassy to rubbery, is known as the glass transition temperature,
or T
g
.
The stiness reduction produced by the temperature is strongly aected by the amount of
ller. Results in Fig. 1.6 by ChanliauBlanot et al. (1989) show that a compound with ller
content of 45 phr has a percentage variation of the stiness with the temperature higher than
a compound with 0 phr of ller.
The majority of rubber compounds currently used in industry have a glass transition
temperature T
g
much less than 0 . Hence, within the common operative range the material
1.2 Standard Phenomenology 7
0 1 2 3 4
0.0
0.5
1.0
1.5
M
P
a
2.0 phr
1.5 phr
1.0 phr
0.5 phr
Figure 1.4 Results of shear tests on rubber specimens with an increasing ller concentration (Yeoh & Fleming,
1997). The initial material stiness shows a monotonic growth for higher value of ller content in
the range {0.5, 1.0, 1.5, 2.0} phr.
behaves as elastic and the eects due to a glassy state are avoided.
1.2.2 Dynamic
The material behavior abovedescribed refers to the quasistatic response. However, elas
tomers subjected to real world loading conditions possess uidlike characteristics typical of
a viscoelastic material. When loaded by means of a stepwise strain, they stressrelax, i.e., the
reaction force resulting from the application of an initial peak falls to an asymptotic value,
which is theoretically reached after an innite time (see the experiments by Khan et al., 2006,
shown in Fig. 4.8). Moreover, if an external force is suddenly applied, creep is observed and
the strain begins to change slowly towards a limiting value.
Both these phenomena are caused by the complex geometrical entanglements between
chains, which produce a local enhancement of the residual (Van der Walls) force. Under
prolonged loading, such entanglementcohesion will slowly breakdown, giving rise to the
phenomena of stressrelaxation and creep described above (Treloar, 2005). For shorter times
of stressing, these eects are limited and the elastic contribution is predominant.
This behavior provides evidences of the fading memory property of the material. There
fore, the entire strain (and temperature) history must aect the constitutive behavior of lled
rubber elastomers. While the strainrate sensitivity and the failure time dependency are rec
ognized and welldocumented in the case of other materials such metals, the incorporation of
historydependent properties of elastomers requires further clarication.
A frequently employed characterization of elastomers is achieved through sinusoidal strain
histories of frequency . This type of material characterization is frequently referred to as dy
namic meaning that it implicates moving parts, diering from methods leading to quasistatic
response. Therefore, in this context, the adjective dynamic is not reserved to phenomena
involving inertia (e.g., wave propagation) which can be neglected in most of the experimental
conditions.
Under the action of dynamic loading, the deformation of rubber, like other viscoelas
tic solids, occurs with a certain delay owing to viscous friction inside the material. Under
harmonic deformation, this delay manifests itself by a phase shift between the applied dis
8 Chapter 1. Rubber Phenomenology
100 50 0 50 100
0
1000
2000
3000
4000
T C
E
M
P
a
Figure 1.5 Young modulus as function of temperature for polyamide6. The thick dashed lines indicate the
transition leathery zone (Shan et al., 2007).
0 20 40 60 80 100 120
0
200
400
600
800
1000
T C
E
M
P
a
45 phr
29 phr
0 phr
Figure 1.6 Young modulus temperature dependence of a rubberpolyethylene blend for a ller content in the
range {0, 29, 45} phr (ChanliauBlanot et al., 1989).
1.2 Standard Phenomenology 9
0 1000 2000 3000 4000 5000 6000
2.6
2.8
3.0
3.2
3.4
t s
N
o
m
i
n
a
l
S
t
r
e
s
s
M
P
a
Figure 1.7 Nominal stress as function of time for a relaxation experiment on AdipreneL100 (Khan et al.,
2006).
placement and the load (Boiko et al., 2010). This shift is proportional to the viscous losses.
In order to explain thoroughly elastomers behavior under oscillatory deformation, let
u(t) = u
0
+ u sin(t) (1.1)
be the longitudinal displacement in an uniaxial deformation from which the nonlinear La
grangian strain follows as
(t) =
0
+
1
sin(t) (1.2)
obtained by dividing u by the length l
0
of the undeformed specimen. The imposed strain
function (1.2) implies, in the nonlinear case, the timedependent nominal stress response (t),
i.e., force applied to the specimen divided by the initial area, whose steady state response is
assumed to have the Fourier series
(t) =
b
0
2
+
k=1
[a
k
sin(kt) + b
k
cos(kt)] . (1.3)
Here,
S(
0
, ,
1
) :=
1
1
a
1
(
0
, ,
1
), (1.4)
L(
0
, ,
1
) :=
1
1
b
1
(
0
, ,
1
) (1.5)
are the storage and loss moduli, also generically referred to as complex moduli. In general,
neither S nor L depend on
1
, if [
1
[ is small (small strain). On the contrary, the
aforementioned moduli for carbon blackreinforced rubber show a rather strong dependence
on
1
in the case [
1
[ is large. This nonlinear amplitude dependence is called the Payne
eect (see Sec. 1.3).
The storage and loss moduli frequency dependence bears no special name, but it is of
fundamental importance to understand the dynamic behavior of elastomers.
Figure 1.8 outlines the dynamic moduli as function of the frequency for dierent values of
static prestrain
0
(Lee & Kim, 2001). At lower frequencies ( 0) the storage modulus tends
to a nite nonzero value with a nonzero derivative. As it will be shown in a following chapter,
10 Chapter 1. Rubber Phenomenology
0 50 100 150
3.4
3.6
3.8
4.0
4.2
4.4
4.6
4.8
f Hz
S
M
P
a
0
0.95
0
0.75
0
0.65
0 50 100 150
0.20
0.25
0.30
0.35
0.40
0.45
0.50
0.55
f Hz
L
M
P
a
0
0.95
0
0.75
0
0.65
Figure 1.8 Storage S and loss L moduli as functions of the frequency in the range [0, 1200] Hz for
dierent values of static prestrain
0
{0.65, 0.75, 0.95}. The amplitude value was
1
= 0.63 for
all the experiments (Lee & Kim, 2001).
this behavior cannot be described by (linear or nonlinear) standard viscoelastic constitutive
equations. The data collated by Lee & Kim (2001) suggest a nonmonotonic dependence of
the storage modulus upon the static prestrain
0
: from
0
= 0.65 to
0
= 0.75, the storage
modulus S considerably decreases, but it increases again at
0
= 0.95. A similar, but less
accentuated, trend is shown by the loss modulus. Experiments collated in (Gottenberg &
Christensen, 1972; Osanaiye, 1996) and more recently in (Luo et al., 2010) are in agreement
with Lee & Kims results.
As in the static case, the dynamic behavior of elastomers also exhibits a very strong
temperature dependence. This eect is much more pronounced than in the comparable types
of tests conducted upon metals, where the mechanical properties could reasonably be taken
as temperature independent within the common working range.
A standard assumption made in the modeling of lled elastomers, which can be corrobo
rated by experimental data, is the socalled thermorehologically simple behavior. Within this
context, the basic postulate is that a viscoelastic mechanical property  relaxation function,
creep function or complex moduli  at a series of dierent temperatures, when plotted against
the logarithm of time or frequency can be superimposed to form a single curve (Pipkin, 1986;
1.2 Standard Phenomenology 11
3 2 1 0 1 2 3
5
10
15
20
Log
10
f Hz
S
M
P
a
49 C
38 C
26 C
16 C
5 C
7 C
3 2 1 0 1 2
0
1
2
3
4
Log
10
f Hz
L
M
P
a
49 C
38 C
26 C
16 C
5 C
7 C
Figure 1.9 Storage S and loss L moduli plotted against logarithmic frequency (log
10
f) for six dierent
temperatures T {7, 5, 16, 26, 38, 49} (Gottenberg & Christensen, 1972). The same results
are shown in a linear frequency scale in Fig. 1.10.
12 Chapter 1. Rubber Phenomenology
0 2 4 6 8 10
4
6
8
10
12
14
16
f Hz
S
M
P
a
49 C
38 C
26 C
16 C
5 C
7 C
0 2 4 6 8 10
1
2
3
4
f Hz
L
M
P
a
49 C
38 C
26 C
16 C
5 C
7 C
Figure 1.10 Storage S and loss L moduli plotted against frequency in the range f [0, 15] Hz for six dierent
temperatures T {7, 5, 16, 26, 38, 49} (Gottenberg & Christensen, 1972).
4 3 2 1 0 1 2
5
10
15
20
Log
10
t s
K
M
P
a
49 C
26 C
5 C
Figure 1.11 Relaxation function plotted against logarithmic time (log
10
t) for three dierent temperatures
T {5, 26, 49} (Christensen, 2003).
1.3 Other Nonlinear Eects 13
Figure 1.12 Preconditioning cycles of a particlereinforced dumbbell specimen with 20 phr (left) and 60 phr
(right) of carbon black with maximum stretch := /
0
= 3 (Dorfmann & Ogden, 2004).
Williams et al., 1955), shifting the various curves at dierent temperatures along the time
or frequency axis. Such a temperature dependence is schematically shown in Fig. 1.9 for
the storage S and loss L moduli. Similar temperature dependence is shown for the relaxation
function in Fig. 1.11 (Christensen, 2003). Materials obeying this empirical principle are called
thermorehological simple. Even for thermorehological simple materials such a procedure can
be expected to be valid only over a limited time and temperature range, primarily above the
glass transition temperature (see, e.g., Singh et al., 2006).
1.3 Other Nonlinear Eects
Apart from the standard phenomenology described in the previous section, carbon black
lled elastomers present some eects peculiar of this class of materials. These eects are the
Mullins eect, which concerns the quasistatic behavior, and the Payne eect, dealing with
the dynamic response.
1.3.1 Mullins eect
The Mullins eect (Mullins, 1947) is a strain induced softening phenomenon, which is as
sociated mainly with a signicant reduction in the stress at a given level of strain during
the unloading path as compared with the stress on initial loading in stressstrain cyclic tests
(Dorfmann & Ogden, 2003).
In lled rubber this phenomenon is due to the mechanical hysteresis from ller particles
debonding from each other or from the polymer chains caused by the stretching. Owing to this,
highly reinforced elastomers suer a more pronounced stiness reduction than those with low
ller content. After the rst few loading/unloading cycles the internal microstructure reaches
a permanent state and changes in stiness become no more signicant. Figure 1.12 represents
typical loading/unloading curves for a rubber specimen subjected to multiple cycles of uniaxial
stretching (Dorfmann & Ogden, 2004). Although this anelastic eect is irreversible for a xed
temperature, an increase in the temperature of the specimen could result in a partial recovery
of the previously broken bonds and, consequently, on a recovery of the material stiness.
14 Chapter 1. Rubber Phenomenology
For the sake of completeness, in the following the main approaches used to describe the
Mullins eect will be reviewed (see, e.g. Ouyang, 2006; Vakada, 2005, and references therein).
The rst attempt to develop a quantitative theory to account for the softening resulting
from rubber stretching was developed by Blanchard & Parkinson (1952). They considered
that value of the shear modulus G is a measure of the total number of crosslinks within
rubber and a reection of the chemical crosslinks produced within vulcanization as well as
linkages between rubber and ller. They suggested that the decrease in G was due to the
breakdown of linkages between ller and rubber. Their interpretation has provided a useful
starting point for the work of other researchers.
One of the other early investigations was done by Mullins & Tobin (1957) who considered
the lled rubber as a heterogeneous system comprising hard and soft phases. The hard phase
was considered to be inextensible and the soft phase to have the characteristics of gum rubber.
During deformation, hard regions are broken down and transformed into soft regions. Then
the fraction of the soft region becomes greater with the increasing tension which in turn is
responsible for the reduced material stiness. However, Mullins & Tobin did not provide a
direct physical interpretation for their model.
More recently, new insights into Mullins eect have been obtained and many researchers
proposed their own constitutive model (Dorfmann & Ogden, 2004; Govindjee & Simo, 1992;
Horgan et al., 2004; Ogden & Roxburgh, 1999; Qi & Boyce, 2004).
In Govindjee & Simo (1992) a micromechanically based continuum damage model for
carbonblack lled elastomers was introduced. The keypoint of the paper was to incorporate
both a damage induced phenomenon such as Mullins eect and the viscous behavior of a
theory of viscoelasticity. Within the framework of damage elasticity, relaxation processes in
the material are described via stresslike convected internal variables, governed by dissipative
evolution equations (see Chap. 3); they are interpreted as the nonequilibrium interaction
stresses between the polymer chains in the network.
Ogden & Roxburgh (1999) proposed to account for the Mullins eect with a phenomeno
logical model based on the theory of incompressible isotropic elasticity amended by the incor
poration of a single continuous damage parameter. The dissipation is measured by a damage
function which depends only on the damage parameter and on the point of the primary loading
path from which unloading begins. A specic form of this function with two adjustable ma
terial constants, coupled with standard forms of the (incompressible, isotropic) strainenergy
function, was used to illustrate the qualitative features of the Mullins eect in both simple
tension and pure shear. However any eects of residual strain were not incorporated.
Dorfmann & Ogden (2004) introduced a constitutive model for the Mullins eect with
permanent set in particlereinforced rubber. The theory of pseudoelasticity has been used
for this model, the basis of which is the inclusion of two variables in the energy function in
order to capture separately the stress softening and residual strain eects. The dissipation of
energy i.e. the dierence between the energy input during loading and the energy returned
on unloading is also accounted for in the model by the use of a dissipation function, which
evolves with deformation history.
A phenomenological model based on the limiting chain extensibility associated with the
Gent model of rubber elasticity has been proposed by Horgan et al. (2004). The Gent strain
energy function (Gent, 1996) was modied to incorporate stress softening characteristics typ
ical of the Mullins eect. Although the Gent model is phenomenological in nature, a micro
scopically based interpretation was given to all of its constitutive parameters. In this way,
it has been possible to develop a model for the Mullins eect based on the alteration of the
polymeric network. Indeed they showed that their approach is a particular case of the more
general framework of pseudoelasticity developed in (Ogden & Roxburgh, 1999).
1.3 Other Nonlinear Eects 15
1
S
M
P
a
0 phr
10 phr
20 phr
30 phr
40 phr
50 phr
60 phr
70 phr
1
L
M
P
a
0 phr
10 phr
20 phr
30 phr
40 phr
50 phr
60 phr
70 phr
Figure 1.13 Strain dependence of the storage and loss moduli (Payne eect) at 70 and 10 Hz for a rubber
compound with dierent concentration of carbon black ller (Wang, 1999). The graphs suggest a
monotonic dependence of the dynamic moduli on the ller content in the range [0, 70] phr.
The Payne eect becomes unnoticeable for low reinforced elastomers ( {0, 10} phr).
1.3.2 Payne eect
Another softening phenomena which manifests the dependence of the stress upon the entire
history of deformation is the socalled Payne eect. Like the Mullins eect, this is a softening
phenomena but it concerns the behavior of carbon blacklled rubber subjected to oscillatory
displacement. Indeed, the dynamic part of the stress response presents a rather strong non
linear amplitude dependence, which is actually the Payne eect (Chazeau et al., 2000; Huber
et al., 1996; Payne, 1962).
For a dynamic strain arising from a harmonic displacement (1.1), the storage and loss
moduli depends nonlinearly upon the strain amplitude
1
as shown in Fig. 1.13 for a strain
amplitude in the range [0.1, 0.6] and a frequency f = 2/ = 10 Hz.
There have been several attempts to explain the Payne eect by macroscale mechanism
based models. Chazeau et al. (2000) classify them as (i) llerstructure models, (ii) matrix
ller bonding and debonding models and (iii) phenomenological or nonlinear network models.
They also state: Payne himself suggested qualitatively that the amplitude dependence of the
storage and loss moduli were due to a ller network in which the ller contacts depended on
the strain amplitude. At lower amplitudes, he argued that the ller contacts are largely intact
16 Chapter 1. Rubber Phenomenology
and contribute to the high value of the modulus [moduli, the author]. Conversely, at higher
amplitudes the ller structure has broken down and does not have time to reform. Therefore,
Paynes explanation is of class (ii).
Following the work of Payne, Kraus (1984) proposed an empirical model based on the
agglomeration/deagglomeration kinetics of ller aggregates, assuming a Van der Waals type
interaction between the particles. In a paper addressing universal properties in the dynamic
deformation of lled rubbers, Huber et al. (1996) introduced the rheological model of Zener
with a nonlinear and linear spring and a dashpot to corroborate the phenomenologically based
formula
G
0
G
=
1
1 + (
1
/a
c
)
2m
, (1.6)
where G
0
the corresponding
value at very small strain. Moreover, a
c
is a constant and m 0.6 is nearly universal, i.e.
to a large extent independent of temperature, frequency, ller content and type of carbon.
Whilst Huber et al. (1996) call (1.6) a theoretical result, it is still based on a rheological model.
Chazeau et al. (2000) stress this eect in their paper, and so it qualies no, or no much, better
than the phenomenological approach of continuum mechanicians (see Lion & Kardelky, 2004,
for references) who postulate nonlinear stress strain behavior. In those approaches the matrix
ller bonding and debonding is formulated considering the dependence upon the entire stress
history with the debonding modeled by the appropriate irreversibility properties.
In 1999 Wang (Wang, 1999) investigated the impact of the ller network, both its strength
and architecture on the dynamic modulus and hysteresis during dynamic strain. It was found
that the ller network can substantially increase the eective volume of the ller due to rubber
trapped in the agglomerates, leading to high elastic modulus. During the cyclic strain, while
the stable ller network can reduce the hysteresis of the lled rubber, the breakdown and
reformation of the ller network would cause an additional energy dissipation resulting in
the higher hysteresis. The experiments, shown in Fig. 1.13, were done at strain amplitudes
1
[0.1, 60] % with a frequency of 10 Hz under a constant temperature of 70. The results
show that the dependence of the storage and loss moduli upon
1
is strongly inuenced by the
quantity of carbon black ller used in the compound, vanishing for low reinforced elastomers.
Therefore, higher hysteresis at low temperature and low hysteresis at high temperature could
be achieved by depressing ller network formation.
Even though the Payne eect has been known for more than 40 years, a model able to
describe such a phenomenon in the relevant frequency and amplitude range is still missing
1
.
1.4 Experimental Techniques
1.4.1 Testing Procedures
The constitutive nonlinearities of carbon blacklled elastomers must be treated warily by any
experimenter in his exploration of the material properties. Indeed, nonlinearities coupled with
nonuniform strain elds could lead to experimental results very dicult to analyze because
the strain nonuniformities can easily mask the actual nonlinear behavior of the material
(Beatty, 1996). Moreover, the ability of lled rubber to undergo nite strains is a compelling
reason to characterize the material through displacement elds for which the relation between
stress components and the position vector is known at any point of the body.
1
Recently Hfer & Lion (2009) have proposed a new model which seems able to describe the Payne eect;
their constitutive relationship is based on Volterratype fading memory scalar internal variables, but they need
a lot of them to obtain a convincing match of the storage and loss moduli with the experiments.
1.4 Experimental Techniques 17
Among the solutions of the balance equation, every deformation, in equilibrium with zero
body force and supported by suitable surface traction alone, is called a controllable solution.
A universal solution, or universal relation, is a controllable solution valid for all the materials
in a given class (Ericksen, 1954; Pucci & Saccomandi, 1997; Saccomandi, 2001).
Among universal solutions, homogeneous strain are the preferred way to test the behavior
of nonlinear homogeneous isotropic materials (Beatty & Hayes, 1992a). Since the strain eld
is uniform, conducting and measuring displacement and forces are equivalent to control and
measure strains and stresses (Haupt, 2002).
In solid mechanics the testing procedures usually rely upon tension/compression, torsion
and shear experiments.
Thin specimens with a constant crosssection are the preferred way to test the material
behavior under tensile loadings. Indeed, a homogeneous uniaxial state of stress prevails in
the central thin shaft. Thick cylinders are commonly employed either for compressive loading
or torsion testing. In the rst case, barrel deformation of the lateral mantle is avoided by
a proper lubrication of the platelets, while torsion test can be interpreted in terms of a
universal relation. A torsion test carried out on a cylindrical tube produces a homogeneous
shear stress distribution if the wall is thin enough. Shear of a short cylindrical specimen
results in a homogeneous strain eld, as far as the heighttodiameter ratio remains limited
and the bulging of the lateral surface can be ignored.
Nonuniform strain eld could also be used (Beatty & Hayes, 1992b). However, two
dierent situations could occur: nonhomogeneities of the deformation eld are ignored and
the stress eld is interpreted as homogeneous (Przybylo & Arruda, 1998), while combined
mechanicaloptical techniques are used to decouple constitutive nonlinearities and strain in
homogeneities (Chevalier et al., 2001; Sasso et al., 2008). The rst solution could lead to
an approximation error which could become relevant depending upon the testing conditions.
The latter requires hard data processing and its applicability is limited to low strain rate
processes.
1.4.2 Specimen Geometry
The material properties inferred from an experiment can be strongly inuenced by the speci
men geometry. Therefore, a lot of care should be taken to ensure that strain elds within the
specimen reect the ideal homogeneous deformation state. Owing to this, the study of new
specimens and clamps is an active research area (Castellucci et al., 2008; Rittel et al., 2002;
Zhao et al., 2009).
In Tables 1.1 and 1.2 several specimens commonly employed to characterize the stress
strain properties of lled rubber are reviewed. These include: cylindrical specimen for com
pression / tension (e.g., Bergstrm & Boyce, 1998; Lion, 1998), dumbbell (e.g., Drozdov &
Dorfmann, 2003; Kar & Bhowmick, 1997; Yoshida et al., 2004), rubber strip (e.g., Przybylo
& Arruda, 1998; Ramorino et al., 2003), cylindrical (or rectangular) double shear specimen
(e.g., Chazeau et al., 2000; Dorfmann et al., 2002) and compression tension hourglass (e.g.,
Haupt & Sedlan, 2001). In particular, Tab. 1.1 shows those specimen actually employed
for the material characterization carried out in this work, while Tab. 1.2 presents some new
literature proposals.
A Dumbbell specimen is represented in Tab. 1.1a. Its typical dogbone shape is stan
dardized together with molding techniques and dimensions in the ASTM D412 norm (ASTM,
1998, 2003). The proper size of the specimen depends upon the load cell equipped in the
testing machine. The enlarged boundaries are intended to increase the contact area between
the holders and the material, allowing a wider range of strains to be imposed. Moreover, this
particular shape contributes to the prevention of the onset of a fracture near the clamps where
18 Chapter 1. Rubber Phenomenology
the maximum stress is located. Even if a heterogeneous strain eld is present, the deformation
could reasonably be assumed as homogeneous in the middle of the shaft. Therefore, while
the displacement is exerted controlling the distance between the clamps, the deformation is
monitored in the central region. External measurement equipments such as straingauges,
extensometers or fast cameras are commonly employed.
All these devices can accurately measure displacement elds arising from low strain rate
processes. At higher frequencies, however, the additional inertia due to strain gauges and
extensometers can overcome the inertia of the specimen resulting in incorrect measurements.
Moreover, fast cameras produce a large amount of data, which is dicult to handle and
analyze.
Because of these limitations, dumbbell specimens are rarely used for high frequency testing,
while thin rubber strips (tension) or thick cylindrical specimens (compression) are preferred
1
.
In both cases, the deformation is assumed as proportional to the distance between the platelets
by neglecting nonhomogeneities of the deformation eld.
The length of the rubber strips used in dynamic tests is chosen according to the physical
dimensions of the DMA load cell (Dynamic Mechanical Analyzer). Based on this, very small
specimens are also used (e.g., L = 2 mm and H = 2 cm), requiring very careful handling of the
interfaces between the specimen and the load device. Pressure controlled holders can avoid
slipping, but for technical reasons, they are rarely used in testing machines. An alternative
solution could be gluing the specimen into the clamps.
Cylindrical shaped specimens are commonly employed for compression experiments. Dur
ing the test, the displacement is exerted controlling the distance between the platelets. If the
plates are accurately lubricated, the bulging of the lateral surface is avoided and the nonho
mogeneities of the strain eld can be considered negligible. To this end temperature inert
silicon or graphite based lubricants are commonly employed.
Dimensions and shapes of cylindrical specimens are standardized in the ASTM D575
norm (ASTM, 2003). Indeed, the diametertoheight ratio is limited by the occurrence of
buckling instabilities. In particular, by considering a hingedhinged isotropic homogeneous
beam subjected to a compressive load, the critical value of the stretch
c
can be calculated in
terms of the beam shape factor (Fig. 1.14). Based on this, a tall cylindrical specimen, with
H = 54 mm and diameter D = 18 mm, was engineered. Our intention was to use a unique
specimen for both compressive and tensile loads avoiding any eects due to the dierent
shape factors and clamping conditions. Although the chosen heighttodiameter ratio results
in
c
= 0.7, to safely avoid any buckling instabilities, the deformation was limited to = 0.8.
A knife extensometer was used to monitor the strain in the central shaft, while the edges were
glued to the plates for the transmission of the tensile force. An epoxy adhesive, capable of
resisting high temperatures, allowed us to reach tensile stresses of the order of 10 MPa.
The results of these tension/compression tests will be discussed in the next section.
1.5 Experimental Evidences
The experiments presented in this section were conducted by the author at the Bridgestone
Technical Center Europe s.p.a. in Rome and at Dipartimento di Ingegneria Chimica e Mate
riali di Sapienza Universit di Roma. The material tested was a carbon blacklled rubber.
Three dierent compounds were used: a weakly, a medium and a heavily reinforced com
pound, which are indicated, respectively, as A, B and C in the following. Further details on
1
At higher frequencies, specimen geometry and clamping conditions become critical and dierent setup
are used (ASTM, 2003; Lakes, 2004).
1.5 Experimental Evidences 19
Table 1.1 Schematic representation of the rubber specimens used in this work.
Shape Remarks
(a)
L
W
T
clamp
Test: Tension. Shape, size and
molding techniques are standardized
in the ASTM D412 norm (ASTM,
2003). Dimensions: W = 7 mm,
L = 30 mm, T = 2 mm.
(b)
D
H lubricant
Test: Compression. Dimensions:
H = 25 mm, D = 19 mm (ASTM
D575 norm, ASTM, 2003). Bound
ary Conditions: platelets lubricated
with graphite.
(c)
D
H
glue
Test: TensionCompression. Di
mensions: H = 54 mm,
D = 19 mm. Boundary Con
ditions: specimen glued to the
platelets.
(d)
L
W
T
clamp
Test: Tension. Dimensions:
L = 50 mm, W = 5 mm.
Boundary Conditions: in absence
of pressure controlled clamps, the
specimen could be glued to the grips
to avoid slipping at the interface.
20 Chapter 1. Rubber Phenomenology
Table 1.2 Schematic representation of specimens used in the literature for elastomers testing.
Shape Remarks
(a)
D
H
specimen fixtures
glue
Test: Tension  Compression 
Torsion. Typical dimensions:
H = 20 mm, D = 30 mm.
Boundary Conditions: the specimen
is molded directly into the grips
(see, e.g., Haupt & Sedlan, 2001).
(b)
H
H
D
fixed
fixed
Test: (double) Shear. Typical
dimensions: H = 2.5 mm,
D = 19 mm. For low strain,
the state of deformation can be rea
sonably interpreted as a pure shear
state (see, e.g., Castellucci et al.,
2008).
(c)
W
T
H
Test: Shear. Typical dimensions:
H = 2 mm, W = 12 mm,
T = 2 mm. The specimen can sus
tain larger strain than the standard
double shear specimen. Moreover,
this test does not necessitate of any
special holders (see, e.g., Zhao et al.,
2009).
1.5 Experimental Evidences 21
0.3 0.4 0.5 0.6 0.7 0.8 0.9
0
1
2
3
4
5
6
c
H
D
Figure 1.14 Heighttodiameter ratio, H/D, as function of the critical strain
c
.
0 500 1000 1500
1.0
1.1
1.2
1.3
1.4
1.5
t s
N
o
m
i
n
a
l
S
t
r
e
s
s
M
P
a
1.45 1.5
1.5
1.65
Figure 1.16 Nominal stress plotted against stretch () for the preconditioning cycles performed on a dumbbell
specimen. The inset shows the leftdrift of the rst three cycles of the constitutive curve which
is actually a manifestation of the Mullins eect. The sixth an seventh repetitions overlap.
1.0 1.1 1.2 1.3 1.4
0.0
0.2
0.4
0.6
0.8
1.0
N
o
m
i
n
a
l
S
t
r
e
s
s
M
P
a
Figure 1.17 Nominal stress plotted against stretch () for the preconditioning cycles performed on a rubber
strip (only the loading paths are shown). The stiness reduction caused by the Mullins eect is
evident.
1.5 Experimental Evidences 23
90 mm
clamp
3 mm
60 mm
Figure 1.18 Rubber pad used to investigate the transverse isotropy induced by the Mullins eect. 6 dierent
specimens, 3 in the direction of the displacement, 3 in the orthogonal direction, were cut after
the preconditioning procedure.
to a velocity of the crossing bar of 50 mm/min (Figs. 1.15 and 1.16). The number of repeti
tions was established by observing that after seven deformation cycles, the material reached
a permanent state and the Mullins eect is no longer signicant. This behavior was observed
for all the compounds.
The softening phenomenon manifestation of the Mullins eect is due to the change in the
microstructure caused by deformation. This internal damage could induce a preferred direc
tion resulting in a dierent material symmetry with respect to the neat elastomer (Dorfmann
& Ogden, 2004; Horgan et al., 2004).
In order to provide new insights into this phenomenon, a thin pad of rubber with dimen
sions W = 60 mm, L = 90 mm and T = 3 mm was subjected to an uniaxial displacement
leading to the strain shown in Fig. 1.15. Thereafter, six specimens were cut, three in the
direction of the displacement (Specimen V in Fig. 1.18) and three in the orthogonal direction
(Specimen H). If the transverse isotropy induced by the Mullins eect were relevant, the two
classes of specimens would have dierent preferred directions, thus dierent uniaxial behavior.
However, experiments in Fig. 1.19 show that the stressstrain curves of the vertical and
horizontal preconditioned specimens are suciently close to one another. Moreover, after the
rst few loading/unloading cycles the dierences between the experimental curves become
negligible.
In view of these results, the microstructural changes caused by the preconditioning can be
ignored and the material can be assumed reasonably as isotropic.
In the following, all the results have been obtained with preconditioned specimens.
1.5.2 Quasistatic
The quasistatic tests were conducted with a Zwick/Roell z010 for both tensile and compres
sive loads.
A dumbbell specimen (Tab. 1.1a) with size L = 30 mm, W = 7 mm and T = 2 mm,
clamped at the upper and lower ends for the transmission of the tensile force was used. The
experiment was monitored by exerting a displacement with a static load corresponding to ten
sion up to the stretch = 1.5. The corresponding strain was obtained from the displacement
measured in the thin shaft of the specimen, through a contact extensometer equipped with
knife edges.
24 Chapter 1. Rubber Phenomenology
1.0 1.1 1.2 1.3 1.4 1.5
0.0
0.5
1.0
1.5
N
o
m
i
n
a
l
S
t
r
e
s
s
M
P
a
Cycle 1
Specimen V.
Specimen H.
N
o
m
i
n
a
l
S
t
r
e
s
s
M
P
a
Cycle 3
Specimen V.
Specimen H.
N
o
m
i
n
a
l
S
t
r
e
s
s
M
P
a
Cycle 7
Specimen V.
Specimen H.
Figure 1.19 Nominal stress plotted against stretch () in the range [1.0, 1.5] for the horizontally H and
vertically V prestretched specimens. After few loadingunloading cycles the two curves overlap.
1.5 Experimental Evidences 25
1.0 1.1 1.2 1.3 1.4 1.5
0.0
0.5
1.0
1.5
N
o
m
i
n
a
l
S
t
r
e
s
s
M
P
a
65C
45C
25C
Figure 1.20 Nominal stress plotted against stretch in the range [1.0, 1.5]. The experiment was repeated
for three dierent temperatures T {25, 45, 65} . The velocity of the crossbar for the force
transmission was 10 mm/min.
For the compression tests, cylindrical specimens with diametertoheight ratio D/H = 0.76
(Tab. 1.1b) were used. The top and bottom surfaces of the cylinder were lubricated with
graphite to guarantee uniform lateral displacement over the height and, consequently, avoid
bulging of the mantle surface. Central transfer of the load was very accurate so that bending
and torsional deformations, if present, were negligible. Both tension and compression tests
were repeated for three temperatures in the range T 25, 45, 65 and three velocities of
the crossbar v 10, 30, 50 mm/min. To ensure the homogeneity of the temperature eld
inside the sample, each specimen was kept at a constant temperature for one hour. Moreover,
after this heating process, the preconditioning procedure was repeated to take account of the
rebonded physical crosslinks, responsible for the Mullins eect, caused by the temperature
increase.
The tensile and compressive stressstrain constitutive curves are shown in Figs. 1.20 
1.23. Strain rate eects on the material stiness are very limited and the loading curves
almost overlap (Figs. 1.22 and 1.23). Moreover, as seen in Figs. 1.20 and 1.21, the changes in
the compound stiness for a temperature range of 25  65 are negligible.
All the results up to this point, have dealt with quasistatic compression and tension tests.
They have revealed a very dierent material response to tensile and compressive loads. In
particular, at lower strains, the stressstrain curve shows an inection point and, consequently,
a change of the tangent stiness around the undeformed conguration (Figs. 1.22  1.23 and
1.20  1.21). This behavior could be interpreted in terms of the dierent micromechanical
phenomena undergoing compression, e.g., macromolecules entanglement grows, and tension,
e.g., macromolecular chains disentangle themselves. However, for very low strain, a much
smoother change of the stiness has been reported (Mott & Roland, 1995, 1996; Roland et al.,
1999). Therefore, the reliability of the measurements in the initial region of the constitutive
curves can be questioned, e.g., a non accurate transmission of the load may result in an
overestimate of the initial stiness.
To investigate thoroughly this point, a tall cylindrical specimen with sizes and shape
reported in Tab. 1.1c was molded. This nonstandard specimen was used for both tensile and
26 Chapter 1. Rubber Phenomenology
0.5 0.6 0.7 0.8 0.9 1.0
1.5
1.0
0.5
0.0
N
o
m
i
n
a
l
S
t
r
e
s
s
M
P
a
65C
45C
25C
Figure 1.21 Nominal stress plotted against stretch in the range [0.4, 1.0] for a cylindrical specimen. The
experiment was repeated for three dierent temperatures T {25, 45, 65} . The velocity of the
crossbar for the force transmission was 10 mm/min.
1.0 1.1 1.2 1.3 1.4 1.5
0.0
0.5
1.0
1.5
N
o
m
i
n
a
l
S
t
r
e
s
s
M
P
a
50 mmmin
30 mmmin
10 mmmin
Figure 1.22 Nominal stress plotted against stretch in the range [1.0, 1.5]. The experiment was repeated
for three velocities of the crossbar v {10, 30, 50} mm/min. The temperature was held constant
at 25 .
1.5 Experimental Evidences 27
0.6 0.7 0.8 0.9 1.0
3
2
1
0
N
o
m
i
n
a
l
S
t
r
e
s
s
M
P
a
50 mmmin
30 mmmin
10 mmmin
Figure 1.23 Nominal stress plotted against stretch in the range [0.6, 1.0] for cylindrical specimen. The
experiment was repeated for three velocities of the crossing bar v {10, 30, 50} mm/min. The
temperature was held constant at 25 .
compressive loading, avoiding the inuence on the measurements of the clamping conditions
and of the shape factor.
The results of these nonstandard tests are shown in Fig. 1.24. The change in the tangent
stiness between compression (thick cylinder) and tension (dumbbell) is noticeable, while a
much smoother change is observed with the tall cylindrical specimen. Moreover, the material
stiness measured with dumbbell and tall cylindrical specimens is comparable, while it diers
considerably from that measured with the thick cylindrical shape. This suggests a low relia
bility at lower strains of the measurements, in particular during compression, owing to both
an inaccurate transfer of the load and an improper lubrication. By using a unique specimen,
these experimental diculties have been overcome and the material behavior in the range
[0.8, 1.2], i.e., the actual working range, have been properly described.
1.5.3 Dynamic
The carbon blacklled rubber behavior cannot be inferred by means of static experiments
only. Indeed, real world loading conditions imply loading rates which are outside the range of
quasistatic tests (v 10 mm/min). Dynamic tests for elastomers are usually conducted by
statically stretching the specimen to a large value of strain and then making it oscillate with
a small amplitude sinusoidal time law (Cho & Youn, 2006; Darvish & Crandall, 2001; Knauss
et al., 2008; White et al., 2000). However, this procedure does not allow the triggering of all
the nonlinearities of the dynamic response.
In this work, with the intention of reproducing loading conditions in agreement with the
actual operative range, oscillatory tests at nite strain, relaxation, creep and cyclic experi
ments were conducted. While the oscillatory tests require a very precise Dynamic Mechanical
Analyzer (DMA) to reach higher frequencies, the other experiments do not require expensive
testing machinery. Therefore, apart from the material characterization itself, the intention
was to investigate the possibility of inferring dynamic material properties through relaxation
and creep rather than loading and unloading cycles at nite strain or standard harmonic tests.
28 Chapter 1. Rubber Phenomenology
"
0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4
3
2
1
0
1
N
o
m
i
n
a
l
S
t
r
e
s
s
M
P
a
Figure 1.24 Nominal stress as functions of strain in the range [0.7, 1.4] for three dierent specimen shapes:
dumbbell, thick cylinder and tall cylinder. The results show a dierent value of the tangent
stiness around the undeformed conguration ( = 1) in tension with respect to compression,
probably caused by the lubrication of the platelets. These dierences disappear by using the
same specimen for both tensile and compressive loads.
The stress relaxation tests were performed on both cylindrical and dumbbell specimens
(Tab. 1.1 a and b).
The cylindrical specimen was compressed, starting from the initial undeformed congu
ration, = 1, up to the nal strain = 0.83 (17 %) in t = 0.7 s with a constant strain
rate. Thereafter, the deformation was held xed for 100 s. In particular, the strain rate of the
initial ramp was
0.24 s
1
, corresponding to a velocity of the crossbar for the transmission
of the load v 370 mm/min. The resulting stress relaxation curves are shown in Fig. 1.25
for three temperatures in the range T 0, 25, 65 . It can be seen from the graphs that
an increase in temperature results in a reduction of the material stiness and, consequently,
in a lower (absolute) value of the maximum force.
Theoretically, the same deformation history with an innite strain rate, i.e., t 0,
would have allowed the direct measurement of the viscoelastic properties. However, this is
not possible when dealing with laboratory equipment. Indeed, not accounting for the nite
strain rate of the initial ramp would result in an underestimate of the material characteristic
times (Antonakakis et al., 2006).
The laboratory environments normally imposes a range for the observable time scales.
The highest sampling rate the acquisition channel can reach determines the shortest achiev
able time; besides the duration of the experiment is an upper bound for time scales. The
Zwick/Roell z010 equipment is able to acquire data up to the frequency 10 kHz; 1 kHz was
used, whereas the experiments lasts 30 s. This choice was a compromise between the minimum
observable time scale and number of data samples recorded.
Relaxation tests permit the capturing of the material behavior involving larger charac
teristic times. Since in many engineering applications (e.g., tires, shock absorbers, etc.), the
shortest intrinsic times are also signicant, loading/unloading cycles at high strain rate were
performed. As shown in Figs. 1.26 and 1.27, the loading/unloading path was repeated for
four dierent speeds of the rising ramp in the range
min
= 0.14 s
1
to
max
= 1.09 s
1
. All
1.5 Experimental Evidences 29
0 5 10 15 20 25 30
1.5
1.0
0.5
0.0
t s
N
o
m
i
n
a
l
S
t
r
e
s
s
M
P
a
65 C
25 C
0 C
Figure 1.25 Nominal stress plotted against time in the range t [0, 30] for the relaxation experiment re
peated for three dierent temperatures T {0, 25, 65} . The strain rate during the initial ramp
was
0.24 s
1
, corresponding to a velocity of the crossbar for the transmission of the load
v 370 mm/min.
the loading paths from the undeformed conguration, = 1, to the maximum strain = 0.83
were displacement controlled; all the unloading paths were force controlled up to the zero
force. This has allowed us to perform after each cycle a three seconds creep test to recover the
undeformed, stressfree initial conguration. The timerate of the force controlled unloading
paths were proportional to those of the loading ramps.
Figure 1.27 outlines the stressstrain constitutive curves related to the strain history rep
resented in Fig. 1.26. The stiness growth for increasing values of the strain rate
is evident.
Moreover, the dissipated energy over a cycle, which is proportional to the area of the cycle,
shows a monotonic dependence upon the strain rate. Indeed, this behavior is shown by all
the viscoelastic materials either solids or uids.
Standard oscillatory test, both in tension and in compression, were also performed. They
were obtained through a sinusoidal displacement of amplitude
1
overimposed on a static
stretch
0
as outlined in Fig. 1.28a in the case of a compression test with
0
= 0.17,
1
= 0.1 and f
1
= 5 Hz. The experiment was repeated for
1
ranging in 0.01, 0.05, 0.1
and frequency f
1
[0, 70] Hz.
The time history of the stress is shown in Fig. 1.28b. It is evident from the graph that
the material response involves both short and longterm contributions. In particular, the
relaxation phenomenon associated with larger characteristic times is evident for t < 4 s, while
the stress settles for t > 5 s. The steady state response, corresponding to the last few cycles
in Fig. 1.28, was used to extract the dynamic moduli through Eqs. (1.4) and (1.5). Figures
1.30, 1.31 and 1.32 outline the variation of the dynamic moduli with the frequency f
1
and
their dependence upon the temperature, the compound type and the dynamic amplitude. In
particular, Fig. 1.32 displays a much stronger frequency dependence of a heavily reinforced
compound (Compound C) with respect to a weakly reinforced one (Compound A).
Dynamic tension tests were conducted with a GABO Eplexor 500N testing machine on
the rubber strip represented in Tab. 1.1d. Two prestrain were used,
0
= 0.20 and
0
= 0.40
respectively and the strain amplitude
1
ranges from 0.01 to 0.13. It is seen from Fig. 1.29
30 Chapter 1. Rubber Phenomenology
0 5 10 15 20 25 30 35
0.8
0.85
0.9
0.95
1
t (s)
0 5 10 15 20 25 30 35
1.5
1
0.5
0
t (s)
N
o
m
i
n
a
l
S
t
r
e
s
s
(
M
P
a
)
(a)
(b)
Figure 1.26 (a) Strain, , and (b) nominal stress plotted against time in the range t [0, 35] for the
loadingunloading experiment. The loading path was repeated for four strain rate
in the range
[0.03, 0.3].
0.8 0.82 0.84 0.86 0.88 0.9 0.92 0.94 0.96 0.98 1
1.5
1
0.5
0
N
o
m
i
n
a
l
S
t
r
e
s
s
(
M
P
a
)
Figure 1.27 Nominal stress plotted against strain, , for four dierent strain rate
in the range
[0.03, 0.3].
The arrow highlights the stiness increase due to the increasing strain rate.
1.5 Experimental Evidences 31
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
0.75
0.8
0.85
0.9
0.95
1
t (s)
1
0.002. Below this value the graphs indicate a peak value at positive
1
and a small
drop for even smaller
1
. Such behavior has also been measured by Hfer & Lion (2009).
Figure 1.33 displays storage, S, and loss, L, moduli plotted against
1
for dierent values
of static prestrain
0
as indicated in the insets. It is evident that both moduli react to the
static prestress and that, in the considered range (
0
0.2, 0.4), S and L increase with
growing
0
. This results partially conrmed the measurements of (Lee & Kim, 2001), even if
in those results a clear monotonic dependence upon
0
was absent at lower strain.
Figure 1.30 shows the analogous behavior of the storage and loss moduli as functions
of frequency for dierent values of
1
. These graphs show a monotonic growth of the
storage and loss moduli and give no hint that the values for S and L would saturate at larger
frequencies.
To see whether the dependencies of storage and loss moduli of Figs. 1.30 and 1.34 are
reliably secured experimental results, the experiments of Fig. 1.34 were repeated several times
under the same conditions. We would be able to show the results by reproducing the graph
of Fig. 1.34, but in order to amplify the dierence between several curves we produced Fig.
1.35, in which the storage and loss moduli are shown in doubly logarithmic representation
for 0.0005
1
0.2. Figure 1.35 shows results for the storage and loss moduli of ve
fold repetitions of the experiments as shown in Figs. 1.30 and 1.34. It is clearly seen that
32 Chapter 1. Rubber Phenomenology
1.0 1.1 1.2 1.3 1.4
0.2
0.4
0.6
0.8
1.0
1.2
N
o
m
i
n
a
l
S
t
r
e
s
s
M
P
a
Figure 1.29 Nominal stress as function of the stretch for dierent values of the dynamic strain amplitudes
1
[0.01, 0.13] in the case of a tensile test (
0
= 0.2, f
1
= 15 Hz) on a rubber strip. For the
largest
1
, the stressstrain curve is no more elliptic, meaning that the nonlinearities become
relevant. The slope variation of the ellipse major axis, which indicates a reduction of the Storage
modulus, is a manifestation of the Payne eect.
0 10 20 30 40 50 60 70
0.8
1
1.2
1.4
1.6
1.8
2
2.2
x 10
7
f
1
(Hz)
S
(
P
a
)
1
= 0.01
1
= 0.05
1
= 0.1
0 10 20 30 40 50 60 70
1
2
3
4
5
6
7
8
9
10
x 10
6
f
1
(Hz)
L
(
P
a
)
1
= 0.01
1
= 0.05
1
= 0.1
Figure 1.30 Storage, S, and loss, L, moduli, plotted against the frequency f
1
in the range f
1
[0, 70] Hz for
three dierent strain amplitude
1
{0.01, 0.05, 0.1}. The temperature was T = 25 and
prestrain was
0
= 0.15 (compression) for all the frequencies. The reduction of the dynamic
moduli with increasing dynamic amplitudes is a manifestation of the Payne eect.
1.5 Experimental Evidences 33
0 10 20 30 40 50 60 70
0.8
1
1.2
1.4
1.6
1.8
2
2.2
2.4
x 10
7
f
1
(Hz)
S
(
P
a
)
T = 0 C
T = 25 C
T = 65 C
0 10 20 30 40 50 60 70
0
2
4
6
8
10
12
14
16
18
x 10
6
f
1
(Hz)
L
(
P
a
)
T = 0 C
T = 25 C
T = 65 C
Figure 1.31 Storage, S, and loss, L, moduli, plotted against the frequency f
1
in the range f
1
[0, 70] Hz
for three dierent temperature T {0, 25, 65} . The dynamic amplitude was
1
= 0.05 and
prestrain was
0
= 0.15 (compression) for all the frequencies.
0 10 20 30 40 50 60 70
0
1
2
3
4
5
6
x 10
7
f
1
(Hz)
S
(
P
a
)
Compound A
Compound B
Compound C
0 10 20 30 40 50 60 70
0
0.5
1
1.5
2
2.5
x 10
7
f
1
(Hz)
L
(
P
a
)
Compound A
Compound B
Compound C
Figure 1.32 Storage, S, and loss, L, moduli, as functions of the frequency f
1
in the range f
1
[0, 70] Hz for the
three dierent compounds tested. The content of ller monotonically increases from compound
A to C. The dynamic amplitude was
1
= 0.05, the prestrain
0
= 0.15 (compression) and
the temperature T = 25 for all the frequencies.
34 Chapter 1. Rubber Phenomenology
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
x 10
7
1
S
(
P
a
)
0
= 0.2
0
= 0.4
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14
1.5
2
2.5
3
3.5
4
x 10
6
1
L
(
P
a
)
0
= 0.2
0
= 0.4
(b) (a)
Figure 1.33 (a) Storage, S, and (b) loss, L, moduli, plotted against the strain amplitude
1
in the range
1
[0, 0.13] for the static prestrains listed in the insets. The graphs suggest a dependence on
the static prestrain. The value of the frequency was f
1
= 10 Hz.
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
x 10
7
1
S
(
P
a
)
10 Hz
15 Hz
50 Hz
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14
1.5
2
2.5
3
3.5
4
4.5
5
x 10
6
1
L
(
P
a
)
10 Hz
15 Hz
50 Hz
(a) (b)
Figure 1.34 (a) Storage, S, and (b) loss, L, moduli, plotted against the strain amplitude
1
in the range
1
[0, 0.13] for three dierent driving frequencies f = 10, 15 and 50 Hz. The variation of S
and L with
1
are manifestations of the Payne eect. The prestrain was
0
= 0.2 (tension) for
all the frequencies.
reproducibility for
1
0.005 is excellent. Below this value, the measured points show
a somewhat enlarged spreading, more for the loss than for the storage modulus, but still
suciently close to one another that the robustness of the results can be trusted.
1.5 Experimental Evidences 35
10
3
10
2
10
1
10
7
1
S
(
P
a
)
10 Hz
15 Hz
50 Hz
10
3
10
2
10
1
10
6.3
10
6.4
10
6.5
10
6.6
1
L
(
P
a
)
10 Hz
15 Hz
50 Hz
(a)
(b)
Figure 1.35 Doubly logarithmic representation of (a) storage, S, and (b) loss, L, moduli as function of
1
for three dierent frequencies f = 10, 15 and 50 Hz. For each of these, ve repetitions of the
experiment are shown. The value of the prestrains was
0
= 0.2 (tension).
36 Chapter 1. Rubber Phenomenology
Chapter 2
Nonlinear Elasticity
Chapter Outline. In this chapter the main aspects of the nonlinear theory of elasticity are presented. As
nonlinear elasticity, and, in particular, hyperelasticity, is such a useful tool in the description of the behavior of
carbon blacklled rubber undergoing quasistatic loadings, the main methodologies for describing the behavior
of materials subjected to large strains are introduced. Some of the results herein presented will be applied to
nonlinear viscoelastic constitutive models and discussed in subsequent chapters.
38 Chapter 2. Nonlinear Elasticity
2.1 Kinematics 39
2.1 Kinematics
In the following section, the basic concepts used to describe the (nite) deformation of a
simple material are briey presented. A comprehensive introduction of nite elasticity can
be found, for instance, in (Beatty, 1996; Holzapfel, 2000; Ogden, 1997).
We consider a continuous body which occupies a connected open subset of a three
dimensional Euclidean point space, and we refer to such a subset as a conguration of the
body. We identify an arbitrary conguration as a reference conguration and denote this by
B
r
.
Let points in B
r
be labelled by their position vectors X relative to an arbitrarily chosen
origin and let B
r
denote the boundary of B
r
. Now suppose that the body is deformed
quasistatically from B
r
so that it occupies a new conguration, B say, with boundary B.
We refer to B as the current or deformed conguration of the body. The deformation is
represented by the mapping : B
r
B which takes points X in B
r
to points x in B. Thus,
x = (X), X B
r
, (2.1)
where x is the position vector of the point X in B. The mapping is called the deformation
from B
r
to B and is required to be onetoone. Its inverse
1
satises
X =
1
(x), x B. (2.2)
Both and its inverse are assumed to satisfy proper regularity conditions, e.g., C
2
(B
r
) C
0
(B
r
).
For simplicity we consider only Cartesian coordinate systems and let X and x respectively
have coordinates X
and x
i
, where , i 1, 2, 3, so that x
i
=
i
(X
x
i
, (2.4)
where grad is the gradient operator in B. With the use of the notation dened by
J = det F. (2.5)
The equation
dx = FdX (2.6)
(in components dx
i
= F
i
dX
G
(2.10)
=J
1
F
T
= J
1
F
G
F
T
(2.11)
where J = det F, T is the so called nominal stress tensor and is the Cauchy stress tensor.
The mechanical interpretation of these stress measures are:
the second PiolaKirchho stress tensor represents a contact force density measured in
the reference conguration per unit of reference area;
the Cauchy stress tensor represents a contact force density measured in the current
conguration per unit of current area;
T
is called rst PiolaKirchho ; it expresses the contact force density in the reference
frame per unit of current area.
We remark that the only assumption used to introduce denitions (2.9)  (2.11) is that
a strain energy density function can be dened in the reference conguration. Indeed, this
is the most general way of describing nonlinear elastic simple materials. Moreover, since
the function depends only on the CauchyGreen strain tensor G, and it is dened in the
reference conguration, it is not aected by any change of observer. The previous requirement
is mechanically known as Principle of Frame Invariance and for the constitutive relations (2.9)
 (2.11) is automatically fullled (Liu, 2004; Murdoch, 2005; Rivlin, 2002, 2005; Truesdell &
Noll, 1965).
In order to clarify this assertion, let us suppose that a rigidbody motion, i.e.,
x
= Qx +c
is superimposed on the deformation x = (X), where Q and c are constant with respect to
the position X (c is the translation vector). Q belongs to the class of the orthogonal tensor,
which we will call Orth
3
. The resulting deformation gradient, say F
, is given by
F
= QF
and
G
:=
1
2
_
F
T
F
II
_
= G
2.3 Restrictions on the Strain Energy Function 41
Therefore, using equation (2.11) the following relation for the Cauchy stress tensor holds for
each deformation gradient F and for all Q Orth
3
:
= QQ
T
. (2.12)
Relation (2.12) expresses the fact that the constitutive law (2.11) is objective. In essence, it
means that the material properties are independent on the superimposed rigidbody motions.
2.3 Restrictions on the Strain Energy Function
The form of the constitutive law can be simplied if the material is characterized by some
symmetry properties. From the physical point of view this means that there exist a change
in the reference placement such that after this change the material is indistinguishable.
The set of all material symmetry transformations at a material point X depends on the
selected reference conguration and for hyperelastic materials can be dened as
g
R,X
=
_
H Lin
+
[ det H = 1
_
H
T
GH
_
= (G)
_
(2.13)
where Lin
+
is the space of the positive denite tensor. The set g
R,X
is a group (see, e.g.,
Ogden, 1997, for the completed proof) and it is called material symmetry group.
Materials which are indistinguishable after every rotation of the reference frame are called
isotropic materials. In such a case, it results
Orth
+
3
g
R,X
(2.14)
i.e., the symmetry group contains the class of all the rotations of the reference frame. Con
stitutive equations for isotropic materials are, actually, the simplest ones.
According to (2.13), for isotropic materials the energy function must fulll the condition:
Q Orth
+
3
,
_
Q
T
GQ
_
= (G) (2.15)
Orth
+
3
being the class of all the rotations.
Every scalar function of a symmetric tensor Gwhich satises (2.15) is called an Isotropic
Tensor Function of G. An isotropic scalarvalued function of Gis also called a scalar invariant
of G. It may easily be checked that the principal invariants of G, dened by
I
1
(G) = trG, (2.16)
I
2
(G) =
1
2
_
I
2
1
(G) trG
2
, (2.17)
I
3
(G) = det G (2.18)
are scalar invariants in accordance with denition (2.15).
Rivlin & Ericksen (1955) showed that a scalarvalued function of a symmetric tensor G is
isotropic if and only if it is expressible as a function of I
1
(G), I
2
(G) and I
3
(G).
Hence, for isotropic materials the strain energy density function takes the form:
(G) =
(
I
1
(G),
I
2
(G),
I
3
(G)) (2.19)
or, since
G =
1
2
(CI) , (2.20)
it is natural to express the strain energy density function in terms of the invariants of the
strain tensor C, i.e.,
(C) =
(I
1
, I
2
, I
3
), (2.21)
42 Chapter 2. Nonlinear Elasticity
where
I
1
=
I
1
(C), I
2
=
I
2
(C), I
3
=
I
3
(C). (2.22)
According to equation (2.19) and denitions (2.16)  (2.18), the most general form of the
second PiolaKirchho stress tensor for an isotropic and hyperelastic material is:
T =2
C
=2
_
I
1
+ I
1
I
2
_
I 2
I
2
C+ 2 I
3
I
3
C
1
, (2.23)
where the following equalities have been used:
I
1
C
= I,
I
2
C
= I
1
I C,
I
3
C
= I
3
C
1
. (2.24)
From equation (2.11), the relation between the Cauchy stress and the strain invariants
follows
=2 I
1/2
3
I
3
I + 2I
1/2
3
_
I
1
+ I
1
I
2
_
B2I
1/2
3
I
2
B
2
=
0
(I
1
, I
2
, I
3
) I +
1
(I
1
, I
2
, I
3
) B+
2
(I
1
, I
2
, I
3
) B
2
. (2.25)
By applying the CayleyHamilton theorem, the previous equation can rewritten as
=
0
(I
1
, I
2
, I
3
) I +
1
(I
1
, I
2
, I
3
) B+
1
(I
1
, I
2
, I
3
) B
1
, (2.26)
being
_
0
=
0
I
2
2
,
1
=
1
+ I
1
2
,
1
= I
3
2
(2.27)
Assuming that the stress vanishes in the reference conguration (T(I) = 0), one gets the
following restriction on the strain energy :
I
1
C=I
+ 2
I
2
C=I
+
I
3
C=I
= 0 (2.28)
A stressfree reference conguration is commonly called a natural state.
2.4 Compressibility
A typical choice to model compressible materials is to decompose the left CauchyGreen strain
tensor into a pure isochoric and a pure volumetric part (Flory, 1961; Sansour, 2008):
C = C(I
1/3
3
I) (2.29)
so that det C = 1.
Furthermore, the rst and the second modied invariants are introduced as the invariants
of C in the same manner of those of C in (2.16)(2.17):
I
1
= I
1
(C) = I
1/3
3
I
1
, I
2
= I
2
(C) = I
2/3
3
I
2
. (2.30)
According to (2.29), the relation I
3
= 1 holds for all deformations.
2.4 Compressibility 43
In the eld of nonlinear mechanics, an ansatz assumed by several researchers is that the
strain energy function is additively decomposed as
(I
1
, I
2
, I
3
) =
I
(I
1
, I
2
) +
V
(I
3
), (2.31)
where
I
depends only upon the isochoric part of the deformation and
V
depends on changes
in volume (Ogden, 1997; Sansour, 2008; Simo et al., 1985). This choice could eventually leads
to a nonphysical behavior at large strains (Eihlers & Eppers, 1998). From equation (2.25),
the Cauchy stress tensor becomes
=2 I
1/2
3
_
I
I
3
+
V
I
3
_
I + 2I
1/2
3
_
I
I
1
+ I
1
I
I
2
_
B
2I
1/2
3
I
I
2
B
2
(2.32)
and from denitions (2.30) one gets the following derivatives of the the strain energy function
with respect to the modied invariants I
1
, I
2
and I
3
,
I
I
1
=I
1/3
3
I
I
1
,
I
I
2
=I
2/3
3
I
I
2
,
I
I
3
=
1
3
I
1
3
_
I
I
1
I
1
+ 2
I
I
2
I
2
_
.
(2.33)
By substituting equations (2.33) into (2.32), one gets:
D
(B
D
) = 2 I
1/2
3
_
I
I
1
+ I
1
I
I
2
_
B
D
2 I
1/2
3
I
I
2
_
B
2
_
D
(2.34)
V
(I
3
I) = 2 I
1/2
3
V
I
3
I. (2.35)
in which the subscripts D and V represent the deviatoric and the volumetric part of the
tensor.
Equations (2.34)(2.35) state that:
1. a superimposed deviatoric stress doesnt produce any volume changes, but only shape
changes (tr
D
= 0);
2. a superimposed pressure uniquely produces a volume change.
Points 1 and 2 are the extension to the nonlinear case of the deviatoric/volumetric stress
decomposition introduced in the framework of linear elasticity (and are consistent with that).
Equation (2.35) implies that the reference conguration is stressfree if the following re
striction on the function
V
is valid:
V
I
3
I
3
=1
= 0 . (2.36)
Therefore, assuming a sucient regularity for the function
V
, expression (2.31) can be
expanded by Taylor series around the undeformed conguration (I
3
= 1) as:
V
(I
3
) =
i=2
1
D
i
(I
3
1)
i
(2.37)
44 Chapter 2. Nonlinear Elasticity
where
(D
i
)
1
=
1
i!
V
I
i
3
I
3
=1
. (2.38)
With similar assumptions, one obtain the following expansion around the reference con
guration (I
1
= 3, I
2
= 3) for the function
I
:
I
(I
1
, I
2
) =
i,j=0
c
ij
_
I
1
3
_
i
_
I
2
3
_
j
(2.39)
where
c
ij
=
_
_
0, if (i, j) = (0, 0),
1
i!j!
i+j
i
I
1
j
I
2
I
1
=3, I
2
=3
, otherwise.
(2.40)
2.5 Incompressibility
Experimental evidences shown in Chapter 1 have revealed a negligible change in volume occur
ing during the deformation. This behavior allows the modeling of rubber as an incompressible
material. From one hand, this assumption simplies the determination of equilibrium solu
tions, but, on the other hand, it makes the constitutive relation hard to implement in a
numerical code. Therefore, both nearincompressible and incompressible materials will be
considered in the following.
Every deformation allowed in a hyperelastic incompressible material must satisfy:
I
3
= det C = 1 (2.41)
The constraint (2.41) denes an hypersurface in the space of the deformation gradients.
Any stress normal to this surface, i.e., in the direction det C/F, does not expend work
on any (virtual) incremental deformation x compatible with the constraint. The stress is,
hence, determined by the constitutive law unless a vector parallel to det C/F. From an
energetic point of view this is tantamount to assume the strain energy function as
(I
1
, I
2
, I
3
) = (I
1
, I
2
, 1) p (I
3
1) , (2.42)
where p(I
1
, I
2
) is the Lagrange multiplier associated to the constraint (2.41) which depends
upon boundary conditions.
Once more, assuming a sucient regularity of the function
I
, one obtains the following
Taylor expansion around the reference conguration
I
(I
1
, I
2
) =
N
i=1
c
ij
(I
1
3)
i
(I
2
3)
j
. (2.43)
Equation (2.43) was rstly introduced by Rivlin & Saunders (1952), and for this reason it is
sometimes called RivlinSaunderss expansion. Some of the most used material models and
the respective parameters c
ij
are reported in table 2.1.
The following expression of the Cauchy stress for an incompressible material follows from
equations (2.25) and (2.42):
= p I + 2
_
I
1
+ I
1
I
2
_
B2
I
2
B
2
. (2.44)
In order to verify the stressfree condition ((I) = 0) in the reference conguration, the
unknown pressure eld p must satisfy:
p(I
1
, I
2
)[
C=I
= p(3, 3) = 2
I
1
C=I
+ 4
I
2
C=I
(2.45)
2.6 Homogeneous deformations 45
Table 2.1 Material models based on Rivlins expansion (Hartmann, 2001a).
Author/Model
Mooney Rivlin c
10
c
01
Isihara et al. c
10
c
01
c
20
NeoHooke c
10
Yeoh c
10
c
20
c
30
James et al. c
10
c
01
c
11
c
20
c
30
Biderman c
10
c
01
c
20
c
30
Tschoegl c
10
c
01
c
11
Tschoegl c
10
c
01
c
22
Lion c
10
c
01
c
50
Haupt/Sedlan c
10
c
01
c
11
c
02
c
30
2.6 Homogeneous deformations
In the following we analyze some elementary problems in which the deformation is homoge
neous, i.e. the deformation gradient F is constant in whole body. Homogeneous deformations
are equilibrium solution for all the class of hyperelastic materials; for this reason they are
called universal solutions (see, e.g., Pucci & Saccomandi, 1997).
2.6.1 Simple Tension
In the case of simple tension or simple compression the deformation is given by
x
1
=
1
X
1
, x
2
=
2
X
2
, x
3
=
2
X
3
, (2.46)
hence the deformation gradient is a diagonal matrix
F = Diag
1
,
2
,
2
, (2.47)
where
1
,
2
are called principal stretches, which are constant because the deformation is
homogeneous. Here the
1
direction is the direction of the external load. Since in a simple
tension test, the lateral surfaces of the specimen are supposed to be unloaded, the principal
stresses corresponding to the directions 2 and 3 vanish.
For deformation (2.47) the strain invariants become:
I
1
=
2
1
+ 2
2
2
,
I
2
=
2
1
+ 2
2
2
. (2.48)
Incompressible Materials
In the case of incompressibility, the constraint (2.41) must be satised, and the following
relation between the stretches must hold:
2
2
= 1
2
=
1/2
1
, (2.49)
46 Chapter 2. Nonlinear Elasticity
Table 2.2 Uniaxial Cauchy stress arising from a stretch for some of the most used incompressible material
models.
Model
11
NeoHooke
11
= 2c
10
_
1
_
MooneyRivlin
11
= 2 (c
10
+ c
01
)
_
2
_
Yeoh (c
30
= 0)
11
= 2
_
2
_ _
c
10
+ 2 c
20
(
2
3 + 2)
I
1
=
2
1
+ 2
1
1
,
I
2
=
2
1
+ 2
1
(2.50)
The unknown pressure eld can be determined from the condition that the stresses
22
and
33
vanish
p = 2
_
I
1
+ I
1
I
2
_
I
1
=
I
1
,I
2
=
I
2
1/2
1
2
_
I
2
_
I
1
=
I
1
,I
2
=
I
2
1
1
. (2.51)
In the following the solution of the simple tension problem will be presented for three of
the most used nonlinear elastic models, viz. NeoHooke, MooneyRivlin and Yeoh model.
The strain energy function in the NeoHookean model is
I
= c
10
(I
1
3) (2.52)
thus, from equation (2.44):
= p I +
0
B, (2.53)
where
0
= 2 c
10
is the so called initial shear modulus (shear modulus in the reference cong
uration).
From equation (2.51), the expression of the pressure eld follows:
p =
0
1
(2.54)
and the stress in the direction 1 becomes:
11
=
0
_
_
. (2.55)
In the same manner by means of equation (2.44) the unidimensional stressstrain relation
can be obtained for MooneyRivlin and Yeoh model. These results are shown in table 2.2.
Compressible Materials
In the compressible case the relation (2.49) is not anymore valid, thus the orthogonal stretch
2
must be derived from the implicit relation:
22
(
1
,
2
) = 0 (2.56)
or equivalently
33
(
1
,
2
) = 0, which follows from the condition that the components of the
stress
22
and
33
vanish. Equation (2.56) has to be solved numerically for each
1
.
Numerical results in the cases of NeoHooke, MooneyRivlin and Yeoh models will be
discussed in Chapter 5.
2.6 Homogeneous deformations 47
Model
12
NeoHooke
12
= 2c
10
MooneyRivlin
12
= 2 (c
10
c
01
)
Yeoh (c
30
= 0)
12
= 2 [c
10
+ 2c
20
(I
1
3)]
Table 2.3 Simple shear test results for some of the most used material models.
2.6.2 Simple Shear
Another example of a homogeneous deformation state is a simple shear dened by
x
1
= X
1
+ X
2
, x
2
= X
2
, x
3
= X
3
, (2.57)
where is the amount of shear strain. Thus,
F =
_
_
1 0
0 1 0
0 0 1
_
_
, B =
_
_
1 +
2
0
1 0
0 0 1
_
_
, C =
_
_
1 0
1 +
2
0
0 0 1
_
_
(2.58)
and
I
1
(C) = 3 +
2
, I
2
(C) = 3 +
2
, I
3
(C) = 1. (2.59)
that is the strain invariants are even function of the shear strain.
The simple shear is an isochoric deformation that is possible in every compressible, homo
geneous, and isotropic hyperelastic material. The constitutive relation (2.26) shows that the
shear stress related to the shear strain is given by:
12
= (I
1
, I
2
, 1), (2.60)
wherein the generalized shear response function is dened by
(I
1
, I
2
, 1) =
1
(I
1
, I
2
, 1)
1
(I
1
, I
2
, 1) (2.61)
From (2.26) one gets
_
0
(I
1
, I
2
, 1) = 2
I
3
2 I
2
I
2
1
(I
1
, I
2
, 1) = 2
I
1
1
(I
1
, I
2
, 1) = 2
I
2
(2.62)
thus,
12
= 2
_
I
1
I
2
_
(2.63)
It is seen that the shear stress is an odd function of the amount of shear.
Furthermore notice that the shear stress is in the direction of the shear if and only if
(I
1
, I
2
, 1) > 0. However, the only presence of shear stress cannot produce a simple shear
state.
It also follows from (2.32) that the scalar function
1
and
1
are determined by the
normal stress dierences; we have:
11
33
=
1
2
, (2.64)
22
33
=
1
2
, (2.65)
48 Chapter 2. Nonlinear Elasticity
where
33
=
0
+
1
+
1
= (
2
)
2
(2.66)
where since
i
are even function of , this dependence has been explicitly shown in the last
term. Moreover the last relation allows to determine (
2
), hence
0
. We highlight that the
normal stresses are unchanged if the shear is reverted. If these are not furnished, the block
will tend to contract or to expand. Such normal stress eects are typical of problems in nite
elasticity.
The most striking feature of the simple shear problem is that the results (2.64)(2.66) do
not involve the shear stress. On the contrary, the shear stress is determined by the dierence
of the normal stresses:
12
=
11
22
, (2.67)
and it is determined in the same way for every homogeneous, isotropic hyperelastic material,
regardless of the form of the response functions. The formula (2.67) is an example of an
universal relation in the nite elasticity theory.
Chapter 3
Nonlinear Viscoelasticity
Chapter Outline. In the rst section of this chapter the main approaches followed in modeling nonlinear
viscoelastic solids during isothermal deformation are thoroughly described. Since lled rubber can reasonably
be assumed isotropic, only isotropic constitutive relations are considered. Thereafter, the attention is focused
on integral viscoelastic models. In particular, some of the most common integral models are reviewed and
advantages and disadvantages of each of them are highlighted. Uniaxial stretch histories are investigated by
reducing the general threedimensional model to the onedimensional case. In this context, the concept of
dynamic moduli, introduced in linear viscoelasticity and referred to as storage and loss moduli, is extended,
in a consistent manner, to the nonlinear case. Moreover, for some of the constitutive equations examined,
the analytical value of the storage and loss moduli is computed and their frequency behavior discussed with
reference to the collated experiments on carbon blacklled rubber. Finally the predicting capabilities of the
classical linear integral viscoelastic model and of a new model based on hysteretic damping are analyzed.
50 Chapter 3. Nonlinear Viscoelasticity
3.1 Nonlinear Theory of Viscoelasticity 51
3.1 Nonlinear Theory of Viscoelasticity
The theory of viscoelasticity is crucial when describing materials, such as rubber, which exhibit
time dependent stressstrain behavior.
Indeed, carbon blacklled rubber, when loaded with timedependent external forces, suf
fers a state of stress which is the superposition of two dierent aspects: a time independent,
longterm, behavior (sometimes improperly called hyperelastic) opposed to a time depen
dent, shortterm, response. Stepstrain relaxation tests suggest that short term stresses are
larger than the long term or quasistatic ones (Johnson et al., 1995). Moreover, oscillatory
(sinusoidal) tests indicate that dissipative anelastic eects are signicant, which leads to the
consideration of a constitutive relation which depends not only on the current value of the
strain but on the entire strain history. This assumption must be in accordance with some
principles which restrict the class of reliable constitutive equations. These restrictions can be
classied as physical and constitutive. The former are restrictions to which every rational
physical theory must be subjected to, e.g., frame indierence. The latter, on the other hand,
depends upon the material under consideration, e.g., internal symmetries.
The principle of determinism (Truesdell & Noll, 1965) belongs to the rst class. It states
that: the stress at a given material point is determined by the entire past history of the motion
in a neighborhood of the considered point.
Another basic assumption in every rational constitutive theory is the principle of frame in
dierence: the response of the system must be the same for all observers (Liu, 2004; Murdoch,
2005; Rivlin, 2002, 2005; Truesdell & Noll, 1965).
Then one is lead to consider constitutive restrictions. One important assumption, proposed
by Noll (1958), postulates that the material is simple, which means that the stress at a given
material point depends only on the history of the rst order spatial gradient of the deformation,
in a small neighborhood of the material point. Therefore, the inuence of the higher order
spatial gradients is ignored. Although this assumption is usually considered nonconstitutive,
it legitimately belongs to the class of constitutive restrictions because there exhist materials,
e.g., porous media and functionally graded, leading to an unacceptable approximation if higher
order spatial gradients are ignored
1
.
Other simplications of the constitutive laws are obtained by assuming that the material
is nonaging, which means that the microscopic changes at the timescale of experimental
test can be ignored, which indeed, complies with the experimental observations discussed in
Sec. 1.2 (see, e.g., Wineman, 2009). An additional assumption, which is also corroborated by
the experimental data, is the isotropy of the material, i.e., the material at a given point in one
reference conguration is indistinguishable in its response from the same material after it has
been statically rotated into another reference conguration. Then, the constitutive equation
is simplied accordingly. As expected, the constitutive laws for isotropic materials are, by
far, the simplest ones. Finally, the internal material constraints of incompressibility is yet an
other way to restrict and simplify the constitutive laws.
All the results presented are consistent with these principles.
In the following, we are motivated to develop a general nonlinear theory of viscoelastic
ity because, in the practical application of tire industry, rubber materials are used under
conditions which do not comply with the innitesimal deformation assumptions of the linear
theory. For these materials, the range of deformation beyond which superposition and thereby
linearity holds is extremely limited. Anyway, one of the rst requirements for a nonlinear con
stitutive law is that, for a very small deformation, the model reduces to the corresponding
1
The dependence upon higher order spatial gradient may have signicant applications in the modeling
of solids subjected to concentrated strains and fractured. As this, the applications on rubber modeling are
forthcoming (see, e.g., dellIsola et al., 2009).
52 Chapter 3. Nonlinear Viscoelasticity
linear model (see, e.g., Quintanilla & Saccomandi, 2007).
In the following section a review of the constitutive equations used to model nonlinear
viscoelastic solids undergoing isothermal deformation is provided.
3.2 State of the Art
Pioneering Works
Polymeric materials, such as rubber, exhibit a mechanical response which can not be properly
described by means of elastic and viscous eects. In particular, elastic eects account for
materials which are able to store mechanical energy with no dissipation. On the other hand,
a viscous uid in a hydrostatic stress state dissipates energy, but is unable to store it. As
the experimental results reported in Chapter 1 have shown, lled rubber present both the
characteristics of a viscous uid and of an elastic solid. Viscoelastic constitutive relations
have been introduced with the intent of describing the behavior of materials able to both
store and dissipate mechanical energy.
The origin of the theory of viscoelasticity may be traced to various isolated researchers in
the last decades of the 19
th
Century. This early stage of development is essentially due to the
work of Maxwell, Kelvin and Voigt who independently studied the one dimensional response
of such materials. The linear constitutive relationships introduced therein are the base of
rheological models which are still used in many applications (see, e.g., Malkin, 1995). Their
works led to Boltzmanns (Boltzmann, 1874) rst formulation of three dimensional theory for
the isotropic medium, which was generalized later to the anisotropic case by Volterra (1912).
However, in lled rubber, constitutive nonlinearities turn out to be signicant even at
small strains because of their internal entangled structure. In addition, when the material
is able to bear large strains, geometrical nonlinearities also become relevant. Therefore an
approach considering both constitutive and geometrical nonlinearities must be pursued.
The early stage of development, following the works of Maxwell, Kelvin and Voight, has
constituted the starting point for many other researchers. In particular in the mid 40s dierent
modeling strategies have actually been used to describe nonlinear viscoelastic solids.
Internalvariables Formulations
A general approach, introduced by Coleman & Noll (1961), is to formulate the constitutive
equation in terms of thermodynamic statevariables: the internal energy is expressed as func
tion of both the current values of strain (stress) and the socalled internal state variables
(Coleman, 1964; Coleman & Gurtin, 1967; Coleman & Noll, 1963)
1
. The latter may be iden
tied with local microstructural quantities, e.g., ller content (Schapery, 1997), or suitably
dened averages (Lion, 1997). Rate eects are introduced through evolution equations, which
usually relate time ratesofchange of internal variables to thermodynamic forces, which are
the derivatives of the internal energy with respect to each internal variable.
Simo (1987) proposed a constitutive equation based on an internal variables formulation
which has provided a starting point for many successive works (Govindjee & Simo, 1992;
Holzapfel, 1996; Holzapfel & Gasser, 2001; Yoshida et al., 2004). In Simos approach, the
internal energy is split according to the multiplicative decomposition of the deformation gra
1
In the isothermal and adiabatic context considered here, the internal energy coincides with the Helmotz
free energy; furthermore, the ClausiusDuhem inequality reduces to the ClausiusPlanck inequality.
3.2 State of the Art 53
dient into dilatational and volumepreserving parts
2
Even though this choice might lead to nonphysical results at nite strains (Eihlers &
Eppers, 1998), Simos model is able to reproduce the hysteretic behaviors of carbonlled
rubber, incorporating also the Mullins eect. In this case the internal energy is split as
the sum of three dierent parts: (i) a volumetric term depending upon the volume change
J = det F, (ii) a term depending upon the isochoric deformation F = J
1
F and (iii) a term
relying on the internal variable q representing the nonequilibrium part of the stress. The
evolution equation of q is postulated assuming the generalized force is proportional to the
derivative of the internal energy with respect to the isochoric strain. This approach is found
to be computationally very ecient, and thus adopted in many commercial nite element
codes. However, Simos model has not been conclusively proven to satisfy the second law
of thermodynamics for all the admissible processes. Hence, it is essentially restricted to
viscoelastic response for strain states near the elastic equilibrium.
Govindjee & Simo (1992) developed a similar model on the basis of the micromechanical
structures of the carbon black particles and rubber matrices: the relaxation processes in the
material were described through stresslike internal variables. These variables are governed
by dissipative evolution equations, and interpreted as the nonequilibrium stresses due to the
interaction between the polymer chains. Holzapfel (1996) proposed a model in which the
internal energy is additively partitioned in the standard volumetric term plus an isochoric
part. The latter depends both on the isochoric strain and a set of internal variables
that can be regarded as an internal strain tensor. This approach generalizes the additive
decomposition introduced by Simo.
An advantage of the statevariable formulations is that, in contrast to the other approaches,
it is not restricted to isotropic responses. Anisotropic eects could be easily taken into account,
e.g., by introducing state variables depending upon ber orientations as in Holzapfel & Gasser
(2001). Moreover physical theories, such as dislocation models, may be introduced directly in
the formulation of the evolution equations. However, it has been reported that the viscosity
alone is not enough to reproduce the large hysteretic energy behavior of highdamping rub
ber, used in vibration absorbers. Thus, Yoshida et al. (2004) proposed a constitutive model
consisting of two parts: an elastoplastic term with a straindependent isotropic hardening
law, representing the energy dissipation of the material, and a second part consisting of a hy
perelastic body with a damage model, which expresses the evolutional direction of the stress
tensor.
Additive Decomposition of
A decomposition of the deformation gradient into elastic and inelastic terms leads to alter
native formulations of the strain energy function (Meggyes, 2001). This decomposition was
rst proposed by (Sidoro, 1974) and later by Lubliner (1985) who extended the pioneering
work of Green & Tobolsky (1946). Although in the framework of elastoplasticity the decom
position of the deformation gradient, into elastic and plastic terms, relies on clear physical
assumptions, there is a lack of evidence in the context of viscoelasticity. However, it has been
successfully applied in many nonlinear constitutive equations (e.g., Bonet, 2001; Hasanpour
et al., 2009; Haupt & Sedlan, 2001; Hoo Fatt & AlQuraishi, 2008; Huber & Tsakmakis, 2000;
Lion, 1997, and many others).
In this context, it is assumed that the deformation gradient can be decomposed as
F = F
e
F
i
. (3.1)
2
The multiplicative decomposition of the deformation gradient into volumetric and isochoric parts was
rst introduced by Flory (1961) (see, e.g., Sansour, 2008).
54 Chapter 3. Nonlinear Viscoelasticity
The inelastic term F
i
, sometimes called viscous term F
v
, introduces an intermediate con
guration. However, the decomposition (3.1) is a conceptual one, and cannot generally be
determined experimentally since neither F
e
nor F
i
are observable quantities (Vidoli & Sciarra,
2002). The inelastic term in (3.1) was also extended to three, four or more deformation parts
F
i
= F
(1)
i
...F
(N)
i
(3.2)
and was adopted and studied, for example, in elastoplasticity and viscoelasticity (see, e.g.,
Bonet, 2001; Haupt, 1985; Haupt & Sedlan, 2001; Meggyes, 2001, and references therein).
The decomposition (3.1) is generally followed by the ansatz on the internal energy for
which is split as the sum of an equilibrium part and an overstress term, i.e.,
=
e
(C) +
o
(C
e
) . (3.3)
Here C is the rightCauchGreen strain tensor and C
e
= F
T
e
F
e
is the elastic strain in the
intermediate conguration. This form of the internal energy has been postulated by several
researchers (Hasanpour et al., 2009; Haupt & Sedlan, 2001; Hoo Fatt & AlQuraishi, 2008;
Lion, 1997).
Bonet (2001) proposed a more general equation for , e.g.,
=
e
(C) +
o
(C, C
i
) , (3.4)
in which the overstress term depends both on the whole strain tensor C and on the inelastic
(viscous) strain C
i
= F
T
i
F
i
; assuming that the viscous components is proportional to the long
term expression, e.g.,
o
(C, C
i
) =
e
(C
e
) , (3.5)
equation (3.4) reduces to (3.3).
All these constitutive choices for the free energy lead to dierent expressions of the
stress in terms of the deformation gradient. By applying the Coleman and Noll procedure
(Coleman & Noll, 1963), i.e., by restricting the form of the stress tensor in such a way that
the ClausiusPlank inequality is veried for every admissible process , the Piola symmetric
stress tensor is shown below
T = T
e
+T
i
, (3.6)
where T
e
is the equilibrium stress and T
i
the overstress. In particular, for an internal energy
of the form (3.3), the following relations between
e
,
o
and T
e
, T
i
are valid:
T
e
=
e
C
, T
i
=
o
C
i
. (3.7)
Equation (3.7) is not sucient to determine the behavior of the material. In order to
complete the description, the evolution equations (or ow rules) of the internal variables F
e
and/or F
i
, which determine the way viscoelastic processes evolve, must be dened. Often
the evolution equations are suitably dened to be ecient with respect to time integration
algorithms (Holzapfel & Gasser, 2001). A common choice for the ow rule is to apply a gener
alization of the onedimensional linear Maxwellmodel to the threedimensional and nonlinear
regime. In this case the evolution equations are assumed to be linear, and the overstress term
arising from them is the generalization of the extrastress arising in Maxwell element (see,
e.g, Holzapfel, 2000; Huber et al., 1996). Lion (1997) proposed nonlinear evolution equations
based on strain, time and temperature. Bonet (2001) also used nonlinear evolution equations
of rate type for the internal variables. These are based on a particular linear relaxation form
of the Maxwell model which leads to a viscoelastic formulation that can be seen as a particular
case of a large strain viscoplastic model. A variational formulation of Bonets model has been
developed in (Fancello et al., 2008).
3.2 State of the Art 55
Integral Formulation
The internalvariables formulation is not the only way to dene the internal energy. Following
the seminal work of Boltzmann (1874), Green & Rivlin (1957) and successively Coleman & Noll
(1961) proposed constitutive relations for which the stress a time t depends upon the entire
history of deformation up to the current time instant. However, the denition of the internal
energy was developed accordingly and in agreement with the fading memory properties, i.e.,
strains which occurred in the distant past have less inuence on the present value of than
those which occurred in the more recent past.
To mathematically express the fading memory property, Coleman & Noll (1961) introduced
the following inner product in the space of deformation histories
C
t
1
(s) : C
t
2
(s) :=
_
0
tr
_
C
t
1
(s)C
t
2
(s)
_
h
2
(s)ds , (3.8)
which induces the norm
_
_
C
t
(s)
_
_
H
t
:=
_
C
t
(s) : C
t
(s)
_
1/2
. (3.9)
Here, C
t
represents the history of the rightGreen strain tensor up to time t, i.e.,
C
t
(s) = C(t s) , s [0, ) . (3.10)
Moreover h(t) is called obliviator of order r and it satises the following conditions (Trues
dell & Noll, 1965):
1. h(s) is dened for 0 s < and has a positive real value: h(s) > 0;
2. h(s) is normalized by the condition h(0) = 1;
3. h(s) decays to zero monotonically for large s in such a way that
lim
s
s
r
h(s) = 0. (3.11)
The norm (3.9) equipped with a function h satisfying properties 1, 2 and 3 is called fading
memory norm; in this topology, two deformation histories are distant if they are distant in
the recent past, i.e., deformations which occur in the recent past have more weight than those
which occurred in the distant past.
Fichera (1979) was the rst to point out that in order to assure the existence of equilib
rium solutions the inuence function has to satisfy the decreasing condition (3.11) (see the
contribution of Fichera in Carroll & Hayes, 1996).
The mathematical assumptions behind the theory of fading memory have been recently
reviewed by Drapaca et al. (2007). Denition (3.9) leads to the socalled Strong Principle
of Fading Memory (Truesdell & Noll, 1965) (SPFM)
1
, which denes the class of admissible
internal energy functionals:
There exists an obliviator h(t) of order greater than n+1/2 such that the constitutive function
is dened and ntimes Frchetdierentiable in a neighborhood of the zero strain history.
The SPFM alone is not sucient to dene properly an internal energy, but rather restrictive
conditions must be satised to assure the existence of a stationary point of (see, e.g.,
Fabrizio et al., 1995).
Over the years, many researchers have dealt with a proper denition of internal energy
accounting for deformation histories (Del Piero & Deseri, 1997; Fabrizio & Morro, 1992;
1
Less restrictive conditions on are stated in the socalled Weak Principle of Fading Memory (WPFM)
(Drapaca et al., 2007).
56 Chapter 3. Nonlinear Viscoelasticity
Golden, 2005, and references therein). It is well known that the internal energy and entropy
of a material with memory is generally not uniquely dened (Golden, 2001). A fundamental
result in this area is due to Gurtin & Hrusa (1988) who obtained a necessary and sucient
condition for the existence of the internal energy arising from a stressstrain constitutive
relation of single integral type. Moreover they were able to develop the following explicit
formula for ,
(C
t
) =
1
(C(t)) +
_
0
(, C(), C(t ))d . (3.12)
The majority of models obtained from an a priori internal energy ts within this single hered
itary framework (Haupt & Lion, 2002; Hfer & Lion, 2009; Lion & Kardelky, 2004). The
stress arising from the constitutive assumption (3.12) involves a single hereditary integral of
a nonlinear function of the strain. In this wide sense, single integral constitutive relations
encompass a class of viscoelastic models equivalent to dierential and fractional dierential
models (Adolfsson et al., 2005; Hanyga, 2007; Hanyga & Seredynska, 2007).
The theory of single integral constitutive equations developed by Gurtin & Hrusa was
extended to multiple integral functionals by Hanyga & Seredynska (2007). In this context,
the internal energy reads as
(C
t
) =
1
(C(t))+
+
N
n=1
_
0
...
_
0
(, C(), C(t
1
), ..., C(t
n
))d
1
... d
n
,
(3.13)
where N is a positive integer. The models introduced by Green & Rivlin (1957), Pipkin
& Rogers (1968), and Hassani et al. (1998) t into this enlarged framework
1
. However, in
general, their applicability to describe the behavior of real materials is questioned since many
parameters are necessary to t the experimental data (see Chap. 5 of Lockett, 1972).
Single Integral Formulation
In applied viscoelasticity not all the constitutive equations are formulated by an apriori
dened internal energy , but the constitutive model is expressed directly by the functional
relation between the stress and the strain through an hereditary integral. In rheology this class
of constitutive models is called RivlinSawyers models; Fungs (Fung, 1972), Fosdick and Yus
(Fosdick & Yu, 1998) and many other models currently used belong to this constitutive class.
Because RivlinSawyers models are not obtained with the Coleman and Nolls procedure, their
thermodynamic consistency must be veried a posteriori.
Single hereditary formulation has proven to reproduce all the crucial aspects of rubber
behavior (hysteresis, relaxation and creep). In the simplest situation, the current value of the
stress is the sum of two dierent contributions: a purely elastic term depending on the current
value of the strain and a hereditary integral depending on the whole strain history. In the
linear model of viscoelasticity, introduced by Bernstein et al. (1963), the stress dependence
on the strain history is assumed to be linear, i.e.,
(t) =2 G(t) + tr G(t) I
+
_
t
0
n
k
t
n
0 , (3.14)
The condition of innite dierentiability can be dropped obtaining a much weaker require
ments for k (see, e.g., Lion & Kardelky, 2004),
k(t) 0, k
(t) 0, k
(t) 0 . (3.15)
In order to express the reduced relaxation function k(t), a discrete relaxation spectrum,
whose form was derived by various molecular models (Pipkin, 1986), is generally used. For
mally,
k(t) =
N
i=1
k
i
+
N
i=1
(1 k
i
) e
i
,
N
i=1
k
i
< 1 . (3.16)
Equation (3.16) is commonly referred to as Pronys series.
Recently fractional calculus has started to play an increasing role in polymer rheology
(Adolfsson & Enelund, 2003; Adolfsson et al., 2005; GilNegrete et al., 2009; Hanyga, 2007;
Haupt & Lion, 2002; Haupt et al., 2000; Lion & Kardelky, 2004). This is due to the fact
that the frequency dependence of the dynamic moduli of llerreinforced rubber is fairly weak
and essentially of the powerlaw type (Lion & Kardelky, 2004). As shown in the literature
such behavior can be represented with a minimum of material constants using the fractional
calculus.
In terms of fractional derivatives, the relaxation function can be dened as
k(t) = 1 +
1
(1 )
_
t
0
(t s)
(s)ds , (3.17)
where 0 < < 1, (x) =
_
0
z
x1
e
z
dz is the Eulerian Gamma function and (s) is a
suitable function such that lim
t
k(t) < . The time derivative of k becomes the leftsided
RiemannLiouville fractional derivative, i.e.,
0
D
t
(t) =
1
(1 )
d
dt
_
t
0
(t s)
(s)ds . (3.18)
The behavior of the dynamic moduli arising from a fractional order viscoelastic kernel,
like (3.17), has been studied in (Rogers, 1983).
Another way to introduce fractional derivatives is through rheological models of fractional
order. In particular, the fractional Maxwell element corresponds to a spring in series with a
fractional damper. The onedimensional linear stress, , versus strain, , relation of a spring
in parallel with the fractional Maxwell element can expressed in terms of fractional derivatives
(Metzler & Nonnenmacher, 2003), e.g.,
(t) =
_
t
0
_
eq
+
0v
E
(t s)
__
(s)ds (3.19)
where the kernel
E
(t) =
i=0
t
i
(1 + i)
, (3.20)
58 Chapter 3. Nonlinear Viscoelasticity
is the MittagLeer function (see Haupt et al., 2000; Metzler & Nonnenmacher, 2003, and
references therein). Since E
1
(t) equals the exponential function, the Mittag Leer function
is also known as the fractional exponential function.
Fractional order models present a relevant drawback due to the diculty in handling nu
merically constitutive equations of fractional order, in particular of dierential type (Adolfsson
et al., 2005). Moreover, the identication of the constitutive parameters relies on a strongly
illconditioned minimization problem (see Chapter 4). For this reason, they are rarely imple
mented in commercial codes and therefore their use is very limited.
QuasiLinear Viscoelasticity
Because of the inherent nonlinear behavior exhibited by most of carbon blacklled rubber,
the linear formulation is not applicable in general.
In the context of nonlinear viscoelasticity, one of the simplest model is the QuasiLinear
viscoelastic model proposed by Fung (Fung, 1972). In one dimension, he suggested the fol
lowing relationship for the second PiolaKirchho stress T in terms of the stretch := /
0
,
T(t) =
_
t
k(t s)
T
e
[(s)]
s
ds , (3.21)
that is, the tensile stress at time t is the sum of contributions of all the past changes, each
governed by the same relaxation function. T
e
(), a function of alone, is the nonlinearly
elastic response.
Rewriting Eq. (3.21) in the form
T(t) =
_
t
k(t s)
T
e
(s)ds , (3.22)
we see that the stress response depends linearly upon the nonlinear function of the strain
T
e
(), from which the name QuasiLinear derives. If the material is in the natural state for
t < 0, Eq. (3.22) reduces to
T(t) = T
e
[(t)] +
_
t
0
k(t s)
t s
T
e
[(s)] ds , (3.23)
since k(0) = 1 and all the functions are smooth in 0 t < . Equation (3.23) states that
the tensile stress at time t is equal to the instantaneous elastic response T
e
decreased by an
amount depending on the past history, since
k(t) is negative.
Recently many investigators have proposed their own nonlinear viscoelastic constitutive
relationship. Among them, the predictive capabilities of the models introduced in (Fosdick
& Yu, 1998; Hallquist, 1998; Hibbit et al., 2007; Shim et al., 2004; Yang et al., 2000) will be
analyzed in Chap. 4 with respect to the experimental data shown in Chap. 1.
There are several applications of the viscoelastic theory concerning the behavior of carbon
blacklled elastomers at high strain rates (10
2
10
3
s
1
) (Hoo Fatt & Ouyang, 2007; Shim
et al., 2004; Yang et al., 2000). In all these models the time derivative of the strain explicitly
appears in the hereditary term. Hoo Fatt & Ouyangs model is developed from the BKZ
constitutive equation (Bernstein et al., 1963) and is reported to be able to capture the high
modulus due to high strain rates. However, it shows some shortcomings owing to a zero
Youngs modulus in the undeformed conguration and, hence, it will not be considered in the
comparison addressed in Chap. 4.
In particular, for this model, the Cauchy stress arising from a constant strain rate test,
say
0
, such that = 1 +
0
t, is
= 2
1
(I
1
3)
2
_
__
1
_
2
1
k(
0
)
_
2
+
2
3
__
1
2
__
d
3.3 QuasiLinear Viscoelasticity 59
and, hence, the Youngs modulus around the undeformed conguration is zero ([/]
=1
= 0),
which contrasts with the experimental evidences.
Advanced nite element codes are often called upon to simulate tires and biological soft
tissues because of the complex behavior of these NLV materials. Among many, two of the
most used FEA codes, which includes a nite viscoelasticity model, are the Abauqs FEA and
the LSDyna code. Both of these numerical tools are used in dierent branches of engineering
(e.g, aeronautical, automotive, structural).
In particular, the LSDyna nite viscoelastic relationship (Hallquist, 1998) takes into ac
count rate eects through linear viscoelasticity by a convolution integral. The model corre
sponds to a Maxwell uid consisting of dampers and springs in series. The Abaqus FEA model
is reminiscent of, and similar to, a wellestablished model of nite viscoelasticity, namely the
PipkinRogers model (Pipkin & Rogers, 1968). This model, with an appropriate choice of
the constitutive parameters, reduces to the Fung (QLV) model (Ciambella et al., 2009; Hibbit
et al., 2007).
Dierential Viscoelasticity
Finally, it is worth mentioning another approach used to describe nonlinear viscoelastic solids:
nonlinear dierential viscoelasticity (Biot, 1954; Schapery, 1997; Tvedt, 2008). This theory
has been successfully applied to model nite amplitude waves propagation (Destrade & Sac
comandi, 2004; Destrade et al., 2009; Hayes & Saccomandi, 2000). It is the generalization
to the threedimensional nonlinear case of the rheological element composed by a dashpot in
series with a spring. Thus in the simplest case, the stress depends upon the current values of
strain and strain rate only. In this sense, it can account for the nonlinear shortterm response
and the creep behavior, but it fails to reproduce the longterm material response (e.g., relax
ation tests). The socalled MooneyRivlin viscoelastic material (Beatty & Zhou, 1991) and
the incompressible version of the model proposed by Landau & Lifshitz (1986) belong to this
class.
The substantive grade 1 is generally used referring to such models for remarking the de
pendence of stress on the strain rate only (see Truesdell & Noll, 1965). Constitutive equations
with higher order time derivatives are also used (see, e.g., Fosdick & Yu, 1996).
3.3 QuasiLinear Viscoelasticity
3.3.1 Fungs Model
A quite general integral series representation of the internal energy was proposed by Pipkin
& Rogers (1968). Dai et al. (1992) used the rst term of such an integral series to describe
the nonhomogeneous deformation of a nonlinearly viscoelastic slab. The constitutive relation
they obtained is
T(t) = R[C(t), 0] +
_
t
0
(t s)
(R[C(s), t s])ds. (3.24)
Here R[C(t), 0] represents the stress due to an instantaneous deformation occurring at time
t = 0, while R[C(s), ] is a strain dependent tensorial relaxation function which in the case
of isotropy has the form
R[C(), ] =
0
(, I
i
())I +
1
(, I
i
())C() +
1
(, I
i
())C
1
(), (3.25)
60 Chapter 3. Nonlinear Viscoelasticity
where
0
,
1
,
1
are scalar functions of time and of the principal strain invariants I
1
, I
2
and I
3
(2.16)(2.18) at time . The expression given be eq. (3.24) incorporates the assumption
that there has been no deformation prior to time t = 0.
Johnson et al. (1996a) have shown that, if the relaxation property can be described by a
scalar function k(t), the single integral representation (3.24) is equivalent to the QuasiLinear
Viscoelastic (QLV) model rst introduced by Fung (1972), i.e.,
T(t) =
0
(t)I +
1
(t)C(t) +
1
(t)C
1
(t)
+
_
t
0
k(t s)
_
0
(s)I +
1
(s)C(s) +
1
(s)C
1
(s)
_
ds (3.26)
which follows from (3.24)(3.25) with the identication
i
(, I
i
()) = k()
i
(I
1
(), I
2
(), I
3
()), i 1, 1, 0 (3.27)
where k(t) is the reduced viscoelastic kernel. In the next, the dependence of
i
upon the strain
invariants I
1
, I
2
, I
3
will be specied only when necessary.
The QuasiLinear Viscoelastic (QLV) model has proven to be a successful phenomenologi
cal model for describing the nonlinear viscoelastic behavior of solids (see references in Johnson
et al., 1996b; Rajagopal & Wineman, 2008).
If the reference conguration is stress free, the coecients
0
,
1
,
1
cannot be arbitrary
assigned, but the following restriction
0
(3, 3, 1) +
1
(3, 3, 1) +
1
(3, 3, 1) = 0 (3.28)
must hold.
In the case of incompressibility Fungs model reads as
T(t) =p(t)C
1
(t) +
0
(t)I +
1
(t)C(t)
+
_
t
0
k(t s)
0
(s)I +
1
(s)C(s) ds, (3.29)
where p(t) is the Lagrange multiplier associated to the incompressibility constraint, which in
the dynamic case reads as
t, det C(t) = 1. (3.30)
As in the static case, p(t) must be determined from equilibrium equations and boundary
conditions.
Equation (3.29) states that the stress at time t is equal to the instantaneous response
decreased by an amount depending on the past history, since
k(t) is generally negative valued.
In other words, the QLV model reects strain history dependent stress and fading memory.
In order to express relaxation properties, Pronys series (3.16) might be used, i.e.,
k(t) =
N
i=1
0
+
N
i=1
_
1
i
0
_
e
i
,
N
i=1
i
=
, (3.31)
where
0
is a real constant representing the shear modulus in the reference conguration,
is the ultimate value to which the shear modulus settles after an innite time and
i
are the
characteristic time constants.
To emphasize some limits of Fungs model, we henceforth focus on an incompressible
viscoelastic solid for which the instantaneous response is modeled by a NeoHookean stress
strain relationship, i.e.,
T
e
= pC
1
+
0
I (3.32)
3.3 QuasiLinear Viscoelasticity 61
Then from (3.29) we have the identication
0
=
0
,
1
=
1
= 0, that yields:
T(t) = p(t)C
1
(t) +
0
k(t)I, (3.33)
since k(0) = 1, or for the Cauchy stress tensor:
(t) = p(t)I +
0
k(t)B(t). (3.34)
Let us consider an uniaxial deformation described by
x
1
= (t)X
1
, x
2
= (t)
1/2
X
2
, x
3
= (t)
1/2
X
3
, (3.35)
where (t) is the stretch ratio in the direction of the extension. The resulting deformation
gradient has the diagonal form
F(t) = Diag
_
(t), (t)
1/2
, (t)
1/2
_
. (3.36)
Assuming that the uniaxial deformation arises from a uniaxial tension with
11
,= 0,
22
=
33
= 0 enables us to compute the Lagrange multiplier p(t). The resulting nonzero
component of the Cauchy stress is
11
(t) =
_
+
N
i=1
(
0
i
) exp
t/
i
_
[
2
(t)
1
(t)], (3.37)
hence, system response in the steady state (t max
1
, ...,
N
) is
SS
11
(t) =
[
2
(t)
1
(t)], (3.38)
which is the response of a purely elastic nondissipative material.
3.3.2 Fosdick and Yus model
In the framework of nonlinear viscoelasticity, Fosdick & Yu (1998) proposed their own con
stitutive equation.
They assumed that the second PiolaKirchho stress tensor is given by
T(t) =
0
(t)I +
1
(t)C(t) +
1
(t)C
1
(t)
+ J(t)F
1
(t)
__
0
k(s) [C
t
(t s) I] ds
_
F
T
(t), (3.39)
where C
t
(s) is the relative right CauchyGreen strain tensor:
C
t
(s) = F
T
t
(s)F
t
(s) = F
T
(s)F
T
(s)F(s)F
1
(t), (3.40)
being F
t
(s) = F(s)F
1
(t).
If a strain C(t) is suddenly applied at time t = 0, e.g.,
C(t) =
_
I, if t 0
C
+
(t) ,= I, if t > 0
(3.41)
equation (3.39) takes the following form
T(t) = T
e
(t) + J(t)F
1
(t)F
T
(t)
__
t
0
k(t s)C(s)ds
_
F
1
(t)F
T
(t), (3.44)
or in terms of Cauchy stress,
(t) = q(t)I +
0
(t)B(t) +
1
(t)B
2
(t) +F
T
(t)
__
t
0
k(t s)C(s)ds
_
F
1
(t). (3.45)
Fosdick and Yus model has been successfully applied to describe nite amplitude wave
propagation (Destrade & Saccomandi, 2006; Hayes & Saccomandi, 2000; Salvatori & Sanchini,
2005).
By comparing (3.44) with Fungs constitutive equation (3.29), the dierences between the
two models appear:
the instantaneous part of the stress is the same in the two models;
the dissipative term of Fosdicks model diers from that of the QLV model (3.29)
since it represents the history of the symmetric PiolaKirchho stress transformed by
pushforward and pullback deformations.
To investigate thoroughly the behavior of Fosdick and Yus model, let us consider a simple
shear deformation of amount (t) in the plane 12, e.g,
F(t) =
_
_
1 (t) 0
0 1 0
0 0 1
_
_
. (3.46)
Hence, the rightCauchyGreen strain tensor reads as
C(t) =
_
_
1 (t) 0
(t) 1 +
2
(t) 0
0 0 1
_
_
. (3.47)
If there is no traction acting on the lateral surfaces, it results
33
(t) = 0, then the Lagrangian
multiplier p(t) into (3.45) can be computed. The Cauchy stress resulting from Eq. (3.39) has
the following components:
11
(t) =
2
(t)[
0
(t) + 3
1
(t) +
1
(t)
2
(t)],
12
(t) = (t) [
0
(t) + 2
1
(t) +
1
(t)
2
(t)] +
_
0
22
(t) =
1
(t)
2
(t) +
_
0
k(t s)F
1
(s)
e
(s)F(s)ds
_
F
1
(t)
_
, (3.48)
where
e
is the instantaneous elastic Cauchy stress response (elastic response at very short
times), k is the socalled viscoelastic kernel, which characterizes the stress relaxation and
satises k(0) = 1.
Also, SYM represents the symmetric part of the bracketed term. The constitutive re
lation (3.48) is valid for compressible as well as incompressible solids. In the latter case the
hydrostatic term pI in
e
(where p is a Lagrange multiplier) is a workless constraint stress
in both the instantaneous response and in the history term, as expected. Indeed, for an
incompressible solid, it results J(t) = 1, for each t, and
e
has the general form:
e
= pI +
1
B+
2
B
2
, (3.49)
where
1
,
2
are scalar functions of time and of the rst and second principal invariants,
I
1
, I
2
, of C.
Then (3.48) reduces to
(t) = p(t)I +
1
(t)B(t) +
2
(t)B(t)
2
+
2
i=1
SYM
_
F(t)
__
t
0
k(t s)
i
(s)C(s)
i
ds
_
F
1
(t)
_
,
(3.50)
where p(t) = p(t) +
_
t
0
k(t s)
i
(s)C(s)
i
ds
_
F
T
(t),
(3.51)
with p(t) the arbitrary pressure. This is then a special case within the model (3.29).
A numerical comparison between the QLV and the Abaqus FEA model for the simple shear
and uniaxial extension case has been performed to further highlight the dierences between
the two models. The results of the simulations are reported in Section 5.4.
3.4 Linear Viscoelasticity
3.4.1 Reduction from nonlinear theory
In order to prove the compatibility of Fungs QLV model with the linear theory of innitesimal
viscoelasticity (see, e.g., Coleman & Noll, 1961; Gurtin & Sternberg, 1962) we will reduce
the QLV constitutive equation to (3.14) by ignoring higher order terms in the deformation
gradient.
The base assumption of linear elasticity, thus viscoelasticity, is that the second order term
in a formal series expansion (with respect to a small scalar evolution parameter ) of the
displacement gradient can be ignored.
To make the meaning of small precise, let assume that F = F() is the deformation
gradient at time , and let put
H = F I. (3.52)
The tensor H is the gradient of the vector u = u() of the displacement from the reference
conguration.
We put
H
= H (3.53)
where = [[H
= I + H+ o(
2
)
C
= F
T
= I + 2 E+ o(
2
)
C
1
= I 2 E+ o(
2
)
(3.54)
where E = (H+H
T
)/2 is the symmetric part of the displacement gradient, while the linearized
strain invariants are
I
1
= 3 + 2 trE + o(
2
), (3.55)
I
2
= 3 + 4 trE + o(
2
), (3.56)
I
3
= 1 + 2 trE + o(
2
). (3.57)
Let recall the QLV model introduced by Fung
T(t) = T
e
(t) +
_
t
0
k(t s)T
e
(s)ds,
T
e
(t) =
0
(t)I +
1
(t)C(t) +
1
(t)C
1
(t),
(3.58)
3.4 Linear Viscoelasticity 65
where
(I
1
, I
2
, I
3
), = 1, 0, 1, is a function of the strain invariants I
1
, I
2
and I
3
at time t.
For small strains, the coecients
(I
1
, I
2
, I
3
) =
0
+ 2 trE
_
1
+ 2
2
+
3
_
=0
,
0
e
(t) =
_
0
1
+
0
0
+
0
1
I + 2
_
0
1
0
1
E+
+ 2 trE I
=1,0,1
_
1
+ 2
2
+
3
_
=0
.
(3.60)
If the reference conguration is a natural state, then
0
1
+
0
0
+
0
1
= 0 and the QLV
constitutive equation reduces to the linear viscoelastic model, i.e.,
(t) =2 E(t) + tr E(t) I+
+
_
t
0
follows, e.g.,
=
0
1
0
1
= 2
=1,0,1
_
1
+ 2
2
+
3
_
=0
(3.62)
An improved model in linear viscoelasticity can be obtained from equation (3.61) in the
hypothesis of dierent relaxation properties for the deviatoric and spheric parts of the strain.
It means that the previous equation can be rewritten as:
D
(t) = G(0) E
D
(t) +
_
t
0
E
D
(t s)
G(s) ds (3.63)
V
(t) = H(0) trE(t)I +
_
t
0
trE(t s)I
H(s) ds (3.64)
where a subscript ()
D
represents the deviatoric part of a tensor, G(t) and H(t) are the
viscoelastic kernels of the deviatoric and volumetric part of the stress, respectively. They are
usually called shear and bulk modulus.
3.4.2 Dynamic Moduli
There are practical situations in which viscoelastic bodies may be subjected to steady state
oscillatory conditions. The fading memory property guarantees that, for such displacement
law, the stress gets into a steady state, hence, it is reasonable to expect that a special form
of the stress might arise (see, e.g., Christensen, 2003). This situation will here be analyzed.
Let us consider the case of isotropic materials subjected to an uniaxial deformation, and
let the stressstrain relation
(t) = k
(0) (t) +
_
t
0
(t s)
k
(s) ds (3.65)
66 Chapter 3. Nonlinear Viscoelasticity
designate either the deviatoric or the volumetric part of the stress in equations (3.63) and
(3.64), depending upon whether = 1 or = 2.
With integration by parts, equation (3.65) becomes
(t) = k
(t) (0) +
_
t
0
k
(t s) (s) ds (3.66)
We let the strain history be specied as being a harmonic function of time according to:
(t) =
0
+
1
sin (t) (3.67)
hence
(t) = k
(t)
0
+
1
_
t
0
k
(t s) cos ( s) ds
= k
(t)
0
+
1
_
t
0
k
(t)
0
+
1
cos ( t)
_
t
0
k
(s) cos ( s) ds
+
1
sin ( t)
_
t
0
k
(s) sin ( s) ds
(3.68)
It is useful to decompose the viscoelastic kernel k
(t) as:
k
(t) =
k
(t) (3.69)
where in accordance to fading memory properties, it results:
(t) 0 as t (3.70)
By means of (3.69) and (3.70), the steady state response of the system (t ) is:
s
(t)
0
+
1
sin( t)
_
+
_
0
(t) sin( s) ds
_
+ cos( t)
_
0
(t) cos( s) ds
=
k
0
+
1
[S
() sin( t) + L
() =
+
_
0
sin( s) ds
L
() =
_
0
cos( s) ds
(3.72)
are generally referred to as the storage and the loss moduli. Moreover the complex function:
k
() = S
() + j L
(), j =
1 (3.73)
is called dynamic modulus.
From equation (3.71) we can obtain the inverse relation between the stress and the complex
moduli
S
() :=
2
2
_ 2
s
(t) sin( t) dt = (
2
1
)
1
_ 2
s
(t)
(t) dt (3.74)
L
() :=
2
2
_ 2
SS
(t) cos( t) dt = (
2
1
)
1
_ 2
SS
(t) (t) dt (3.75)
3.4 Linear Viscoelasticity 67
which allows us to calculate these moduli from experimental data.
Pronys series form is generally chosen for the viscoelastic kernels, i.e.
k
(t)
k
0
)
1
= k
(t) = 1 +
N
i=1
k
i
_
exp
i
1
_
(3.76)
where k
0
is the modulus (shear or bulk) in the reference conguration. The normalized storage
and loss moduli become:
S
()
k
0
=1
N
i=1
k
i
+
N
i=1
k
i
2
i
2
1 +
2
i
2
(3.77)
L
()
k
0
=
N
i=1
k
i
1 +
2
i
2
(3.78)
3.4.3 Some Remarks on Energy Dissipation
Dissipation in materials with memory is strictly related to the structure of the viscoelastic
kernel.
To investigate more thoroughly this point, let us consider an incompressible isotropic
material following the constitutive relation (3.61) subjected to the dynamic oscillation (3.67).
For such a material the volumetric part of the stress is not constitutively assigned and the
average rate of working per unit volume in the steady state oscillatory condition is:
E
d
(,
1
) =
2
_ 2
s
(t) (t) dt
=
1
2
2
1
L() (3.79)
Equation (3.79) furnishes an alternative denition of loss modulus and contributes to explain
better the relation with the dissipation.
Assuming, once more, for the viscoelastic kernel the Pronys series expansion, one gets:
E
d
(,
1
) =
1
2
2
1
k
0
N
i=1
k
i
2
1 +
2
i
2
(3.80)
hence E
d
(,
1
) attends its maximum E
M
for , that is:
E
M
= lim
E
d
(,
1
) =
1
2
2
1
k
0
N
i=1
k
i
i
(3.81)
The value of the maximum amount of energy is strictly related to the time behavior of
the viscoelastic kernel k(t); indeed, since
k(t) = k
0
_
1 +
N
i=1
k
i
_
exp
i
1
_
_
, (3.82)
the time derivative of (3.82) around the initial time instant is:
k(t)
t
t=0
= G
0
N
i=1
k
i
i
(3.83)
68 Chapter 3. Nonlinear Viscoelasticity
0 1 2 3 4 5
0.0
0.2
0.4
0.6
0.8
1.0
s
1
E
d
k
0
Differential
Integral
Figure 3.1 Average rate of working per unit volume normalized with respect to k
0
for k
1
= 1,
1
= 1, = 1
and
1
= 1.
0 1 2 3 4 5
0.0
0.2
0.4
0.6
0.8
1.0
s
1
k
0
1
Differential
Integral
Figure 3.2 Storage modulus normalized with respect to k
0
for k
1
= 1,
1
= 1 and = 1.
Comparison between equation (3.80) and (3.83) shows us that for linear viscoelastic ma
terials the average rate of working increases as fast as the memory fades.
This result can be compared with the case of dierential viscoelasticity, which can be con
sidered the limit case of a material with instantaneous memory. The stressstrain constitutive
relation is:
(t) = k
0
(t) + (t) (3.84)
where k
0
, (k
0
> 0, > 0) are material constants. In this case the average rate of working
is:
E
d
(,
1
) =
2
2
1
2
(3.85)
that diers considerably from (3.80), as reported in gure (3.1).
Dierences between (3.80) and (3.85) reect on the dierent storage and loss moduli
behaviors (Figs. (3.2) and (3.3)).
3.5 Nonlinear Dynamic Moduli 69
0 1 2 3 4 5
0.0
0.2
0.4
0.6
0.8
1.0
s
1
k
0
Differential
Integral
Figure 3.3 Loss modulus normalized with respect to k
0
for k
1
= 1,
1
= 1 and = 1.
3.5 Nonlinear Dynamic Moduli
The relation between loss modulus and dissipated energy indicates how to extend the def
initions of dynamic moduli to nonlinear viscoelastic models. For small strain and assumed
linear viscoelastic stressdeformation model, the stress response is monochromatic at the same
frequency of sinusoidal input. However, real world loading conditions do not comply with this
linear assumption, since modes at multiple frequencies are also excited. Therefore, the as
sumptions underlying equation (3.71) are no longer valid, requiring the introduction of a more
general denition.
Experimental evidence showing the variation of the dynamic moduli with respect to the
frequency for carbon blacklled rubber (Osanaiye, 1996), a polyurethane matrix (Gottenberg
& Christensen, 1964) and mozzarella cheese (Singh et al., 2006) are displayed in Fig. 3.4.
When assessed at the lowest frequencies, both the storage and loss moduli of all the
materials considered depend linearly on the frequency. For carbon blacklled rubber this
behavior has been reported by a number of dierent researchers (Lee & Kim, 2001; Luo
et al., 2010). In particular, Luo et al. use the empirical Kraus model endowed with a set of
parameters varying linearly with the frequency, in order to achieve an adeguate match with
the experimental data. However, such an empirical relation cannot be derived from any known
constitutive equation. Consequently, the possibility of introducing a stressstrain constitutive
relation able to reproduce the linear dependence of the dynamic moduli at low frequencies
has been investigated. Because, the linear viscoelastic model (3.65) is unable to reproduce
this behavior, a simple nonlinear model may be appropriate?
In particular, three dierent classes of nonlinear constitutive equation have been consid
ered: (i) elastic, (ii) dierential viscoelastic and (iii) integral viscoelastic. Although (i) cannot
account for dissipative eects, the derivation of the dynamic moduli in this case was useful for
the calculation involving models (ii) and (iii) which, instead, are generally used to describe
materials showing relevant hysteresis losses. Commonly employed constitutive equations be
longing to the class (i) are the NeoHooke and MooneyRivlin hyperelastic models (see, e.g.,
Hartmann, 2001b). The socalled MooneyRivlin viscoelastic material (Beatty & Zhou, 1991)
and the incompressible version of the model proposed by Landau & Lifshitz (1986) belong to
class (ii). Finally, most of the nonlinear integral viscoelastic models t within class (iii) (see
70 Chapter 3. Nonlinear Viscoelasticity
0 1 2 3 4 5
0
1
2
3
4
s
1
M
P
a
0 1 2 3 4 5
0.0
0.2
0.4
0.6
0.8
1.0
s
1
M
P
a
0 1 2 3 4 5
3.0
3.5
4.0
4.5
5.0
5.5
6.0
6.5
s
1
M
P
a
0 1 2 3 4 5
0.2
0.4
0.6
0.8
1.0
1.2
1.4
s
1
M
P
a
0 1 2 3 4 5
0.05
0.10
0.15
0.20
0.25
0.30
s
1
M
P
a
0 1 2 3 4 5
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.10
s
1
M
P
a
f
Figure 3.4 Frequency dependence of the storage and loss moduli for dierent materials: (a)(b) carbon
blacklled rubber (Osanaiye, 1996); (c)(d) polyurethane matrix containing salt and aluminum
powder (Gottenberg & Christensen, 1964); (e)(f ) mozzarella cheese (Singh et al., 2006).
3.5 Nonlinear Dynamic Moduli 71
Chap. 4).
3.5.1 Denitions
Most of the solid materials, under a suddenly applied deformation, exhibit a decaying stress
as function of time, while they react to a fast load with a time increasing deformation. Both
these phenomena are consequences of the fading memory properties of the material.
Fading memory mathematically translates into a constitutive relation for which the stress
at time t, say (t), depends on the deformation history up to time t, say
t
, i.e.,
(t) = F[
t
]. (3.86)
It should be noted that the function F depends upon the time only implicitly through .
Indeed, this corresponds to the assumption that the material is nonaging.
The idea of fading memory is rendered mathematically through the relaxation property
(see, for instance, Fabrizio et al., 1995). If
1
and
2
are two deformations such that
1
(t) =
2
(t) for each t t
0
, then the relaxation property states that
lim
t
_
F[
t
1
] F[
t
2
]
_
= 0 . (3.87)
In the case of dynamic motion, let
u(t) = u
0
+ u
1
sin(t) (3.88)
be the longitudinal displacement in an uniaxial deformation from which the nonlinear La
grangian strain is as follows
(t) =
0
+
1
sin(t) (3.89)
obtained by dividing u(t) by the reference length l
0
. Therefore, by assuming into (3.87)
1
(t) = (t) and
2
(t) = (t + T), with T = 2/, then for each t,
1
(t) =
2
(t) and
Eq. (3.87); thus,
lim
t
_
F[
t
] F[
t+T
]
_
= 0 , (3.90)
and the steady state stress, say
s
(t), arising from deformation (3.89) is Tperiodic.
Therefore, fading memory properties assure that the material response to a harmonic
deformation gets into a periodic steady state.
Under weak regularity assumptions, e.g.,
s
(t) L
2
_
s
(t) =
C
0
2
+
i=1
_
S
i
sin(i t) +
C
i
cos(i t)
, (3.91)
being
S
i
and
C
i
the Fourier coecients of
s
(t) dened as
S
i
=
_
/
/
s
(t) sin(i t) dt , i = 1, 2, ...,
C
i
=
_
/
/
s
(t) cos(i t) dt , i = 1, 2, ... .
(3.92)
For materials whose constitutive response is linear, only the rst coecients, i.e., i = 1,
in the series (3.91) are non zero. In the nonlinear case, however, the stress arising from the
harmonic displacement (3.88) is no longer monochromatic since components at multiple of the
72 Chapter 3. Nonlinear Viscoelasticity
input frequency are excited. Therefore, by considering the steady state stress response (3.91),
the denition of dynamic moduli can be extended to the case of a nonlinear stressstrain
relation, viz.
S =
S(
0
,
1
, ) :=
1
_
/
/
s
() sin() ds =
S
1
1
,
L =
L(
0
,
1
, ) :=
1
_
/
/
s
() cos() ds =
C
1
1
.
(3.93)
The denitions given in Eqs. (3.93) are consistent with those introduced in the linear context
and, hence, the same nomenclature of storage and loss moduli has been adopted.
By introducing the stationary average rate of working over a cycle, i.e.,
=
(
0
,
1
, ) := lim
t
_
t+2/
t
() ()d , (3.94)
the loss modulus can be expressed in terms of as L = 2/(
2
1
).
Elasticity
In the general framework of materials with memory, nonlinear elasticity can be considered a
subset of this theory aimed at describing materials for which memory eects can be ignored.
Indeed, stress depends only upon the current value of strain. Although elastic constitutive
equations are not used to describe the dynamic behavior of elastomers, since they cannot
account for dissipated energy, for the sake of completeness the dynamic moduli have also
been derived.
The constitutive equation for a Cauchy elastic material reads as
(t) = f[(t)] (3.95)
where f is a nonlinear smooth function of (t).
If (t) is Tperiodic, i.e., periodic of period T, then (t) is also Tperiodic:
(t + T) = f[(t + T)] = f[(t)] = (t) , (3.96)
and, therefore, (t) =
s
(t) being the Fourier series (3.91) convergent.
Through the change of coordinates s = /, the denitions (3.93) become
S =
S(
0
,
1
) =
1
1
_
1
1
f[
0
+
1
sin(s)] sin(s) ds ,
L =
L(
0
,
1
) =
1
1
_
1
1
f[
0
+
1
sin(s)] cos(s) ds ,
(3.97)
which explicate that, for elastic constitutive equations, both the storage and loss moduli are
frequency independent.
Under stronger regularity assumptions on f, e.g., f is an analytic function of , the Taylors
series around the prestrained conguration =
0
(
1
= 0) converges to f, i.e.,
s
(t) = f[
0
+
1
sin(t)] =
i=0
i
1
i!
f
i
(
0
) sin
i
(t) , (3.98)
3.5 Nonlinear Dynamic Moduli 73
with
f
i
(
0
) =
_
i
f/
i
1
=0
. Therefore, Eqs. (3.97) become
S =
i=0
i1
1
i!
f
i
(
0
)
_
1
1
sin
i+1
(s) ds =
i=0
f
i
(
0
)
i1
1
i!
1 + (1)
i+1
(1 +
i
2
)
(
3+i
2
)
L =
i=0
i1
1
i!
f
i
(
0
)
_
1
1
sin
i
(s) cos(s) ds = 0
(3.99)
where () is the Euler gamma function. It should be noted that to write Eqs. (3.99), the
innite summation has been taken out of the integral since the series (3.98) satises the
hypothesis of the Dominated Convergence Theorem (see, e.g., Kolmogorov, 1999).
Equations (3.99) state that in the purely elastic case, i.e., for a material described by the
constitutive relation (3.95), the loss modulus and the dissipated energy are zero. Indeed,
this result is in accordance with the denition of hyperelastic materials for which the work
done by body and surface forces is converted into kinetic energy and stored elastic energy,
without dissipation.
Equations (3.99) allow the storage modulus to be obtained simply by computing the
derivatives of f whatever constitutive equation in the form (3.95) is chosen.
Dierential Viscoelasticity
The constitutive relation of a viscoelastic material is:
(t) = g[, ] . (3.100)
Hereafter, a superimposed dot will represent the rstorder time derivative.
The constitutive equation (3.100) is generally referred to as dierential viscoelasticity to
distinguish it from the integral formulation introduced in the following. The substantive
grade 1 is also used referring to (3.100) for emphasizing the dependence of on the strain
rate only (and not on the higher order time derivatives).
As in the case of a purely elastic material, if (t) is Tperiodic, also (t) is Tperiodic,
e.g.,
(t + T) = g[(t + T), (t + T)] = g[(t), (t)] = (t) ; (3.101)
therefore
s
(t) = (t).
If the function g is analytic with respect to
1
, Taylors series around the prestrained
conguration =
0
(
1
= 0) coincides with g (see the Appendix); therefore, the storage and
loss modulus become
S =
S(
0
,
1
, ) =
i=0
i
n=0
in
i
1
n!(i n)!
g
i
(
0
)
_
1
1
sin(s)
n+1
cos(s)
in
ds ,
L =
L(
0
,
1
, ) =
i=0
i
n=0
in
i
1
n!(i n)!
g
i
(
0
)
_
1
1
sin(s)
n
cos(s)
in+1
ds ,
(3.102)
with g
i
(
0
) =
_
(
n
g/
n
) (
in
g/
in
)
1
=0
.
The derivatives with respect to the frequency for 0 are
S
=0
=
i=1
i
1
(i 1)!
g
i
(
0
)
_
1
1
sin(s)
i
cos(s) ds = 0 ,
L
=0
=
i=1
1
2
2i1
1
(2i 2)!
g
2i1
(
0
)
(i
1
2
)
(i + 1)
.
(3.103)
74 Chapter 3. Nonlinear Viscoelasticity
Integral Viscoelasticity
To take into account the fading memory properties of the material, a hereditary integral should
be introduced in the constitutive equation (see, e.g., Christensen, 2003; Drapaca et al., 2007;
Tschoegl, 1989). A convenient formulation to represent most of the nonlinear viscoelastic
models used in the literature is
(t) = f[, ] + l[, ]
_
t
< +, lim
t
k(t) = 0
(3.105)
and, hence, k
< k
0
, with k
0
= k(0).
Constitutive equation (3.104) assures that the stress arising from a Tperiodic strain is
Tperiodic, viz.
(t + T) =f[(t + T), (t + T)]
+ l[(t + T), (t + T)]
_
+
0
k()h[(t + T ), (t + T )] d =
=f[(t), (t)] + l[(t), (t)]
_
+
0
s
(t) =
0
2
+
+
i=1
_
S
i
sin(i t) +
C
i
cos(i t)
(3.107)
with
0
=
0
(f
0
, l
0
, H
0
) ,
S
i
=
S
i
(f
S
i
, f
C
i
, l
S
i
, l
C
i
, H
S
i
, H
C
i
; , t) ,
C
i
=
C
i
(f
S
i
, f
C
i
, l
S
i
, l
C
i
, H
S
i
, H
C
i
; , t) ,
(3.108)
and
H
0
= h
0
(k
k
0
), k
0
< ,
H
S
i
= h
S
i
_
+
0
k(s) cos(is)ds h
C
i
_
+
0
k(s) sin(is)ds,
H
C
i
= h
S
i
_
+
0
k(s) sin(is)ds + h
C
i
_
+
0
k(s) cos(is)ds .
(3.109)
The existence of the integrals in (3.109) is assured by the asymptotic properties (3.105) of the
viscoelastic kernel.
In Eqs. (3.108)(3.109) f
S
i
, f
C
i
, l
S
i
, l
C
i
, h
S
i
and h
C
i
are the Fourier coecients of the
functions f, l and h obtained by Eqs. (3.92); therefore alle these coecients depend on the
frequency, but this dependence has been omitted in (3.109) for the sake of brevity.
3.5 Nonlinear Dynamic Moduli 75
Equation (3.107) explicates the dependence of the harmonic stress upon the oscillatory
strain for the nonlinear constitutive equation (3.104). The dynamic moduli are obtained by a
standard projection of (3.107) over sin(t) and cos(t). After some manipulations, the result
is
S =
S(
0
,
1
, ) =
1
1
_
f
S
1
+
l
0
2
H
S
1
+
l
S
1
2
H
0
_
+
1
2
1
+
i=1
_
l
S
i
H
C
i+1
+ l
S
i+1
H
C
i
+
1
2
1
+
i=1
_
l
C
i
H
S
i+1
l
C
i+1
H
S
i
L =
L(
0
,
1
, ) =
1
1
_
f
C
1
+
l
0
2
H
C
1
+
l
C
1
2
H
0
_
+
1
2
1
+
i=1
_
l
S
i
H
S
i+1
+ l
S
i+1
H
C
i
+
1
2
1
+
i=1
_
l
C
i
H
S
i+1
+ l
C
i+1
H
S
i
(3.110)
If f, l and h are analytic functions of
1
, Equations (3.110) allow the derivative of the
storage modulus for 0 to be calculated, viz.
S
=0
=
k
k
0
2
1
i=1
_
l
S
i
h
C
i+1
l
S
i
h
C
i+1
+
l
S
i+1
h
C
i
+ l
S
i+1
h
C
i
_
=0
+
k
k
0
2
1
i=1
_
l
C
i
h
S
i+1
+ l
C
i
h
S
i+1
l
C
i+1
h
S
i
l
C
i+1
h
S
i
_
=0
.
(3.111)
It can be easily proven that the coecients h
S
i
(0) and l
S
i
(0) (h
C
i
(0) and l
C
i
(0)) are zero if
i is even (odd); furthermore, the derivatives
h
S
i
(0) and
l
S
i
(0) (
h
C
i
(0) and
l
C
i
(0) = 0.
3.5.2 Onedimensional model
In the previous paragraph some inconsistencies between the dynamic moduli, following from
denitions (3.93), and the experimental data, collated from the literature, were highlighted.
In particular, while the results of experiments display a nonzero derivative of the storage
modulus for which tends to zero, all the constitutive equations considered result in an
horizontal tangent at = 0. This frequency behavior produces a relevant error in the tting
at lowest frequencies and is an important drawback as in many operative conditions the
material is subjected to strain rate corresponding to frequencies lower than 200 s
1
.
In order to obtain an accurate match with the experimental data, a constitutive equation
not belonging to the classes (i), (ii) and (iii) must be considered.
As the basis for an appropriate constitutive equation, the linear viscoelastic model already
introduced in Sec. 3.4 was used, i.e.,
(t) = k
0
(t) +
_
t
k
1
1
k
N
N
k
1
~
(a) (b)
k
N
~
Figure 3.5 Panel (a) displays the Nterms generalized Maxwell element while panel (b) displays the modied
Maxwell element with hysteretic dampers.
From Eqs. (3.110), the dynamic moduli are obtained as
S() = k
0
+
_
+
0
k(s) cos(s) ds ,
L() =
_
+
0
k(s) sin(s) ds .
(3.114)
If the Pronys series is used to represent the viscoelastic kernel k(t), viz.
k(t) = k
0
_
1
N
i=1
k
i
+
N
i=1
k
i
e
t/
i
_
,
N
i=1
k
i
< 1 , (3.115)
Eqs. (3.114) can be expressed in terms of characteristic amplitudes k
i
and times
i
:
S =
S(, k
0
, k
i
,
i
) = k
0
_
1
N
i=1
k
i
1 +
2
i
2
_
,
L =
L(, k
0
, k
i
,
i
) = k
0
_
N
i=1
k
i
1 +
2
i
2
_
,
(3.116)
and Eqs. (3.77) and (3.78) are recovered. It should be noted that in the linear case neither
moduli depend on
0
or on
1
, but only on the frequency .
Figure 3.6 shows the frequency dependence of the storage and loss moduli for a carbon
blacklled elastomer (experiments carried out by the author). In the same gure the results
of the tting obtained with model (3.116) are shown.
In order to match adequately the experimental data, some modications of the linear
model must be considered. To this end it is convenient to reformulate Eq. (3.112) in terms of
its equilibrium and overstress parts. The stress can be split into two contributions, i.e.,
(t) =
el
(t) +
ov
(t) = k
(t) +
N
i=1
(i)
ov
(t) , (3.117)
with k
= k
0
(1
N
i=1
k
i
). The evolution equations of the overstress terms are:
(i)
ov
(t) +
1
(i)
ov
(t) = k
0
k
i
(t) i = 1, ..., N . (3.118)
3.5 Nonlinear Dynamic Moduli 77
0 100 200 300 400
10
10.5
11
11.5
12
12.5
13
13.5
14
14.5
(s
1
)
S
(
M
P
a
)
Experimental
Standard Kernel
Complex Kernel
0 100 200 300 400
1
2
3
4
5
6
(s
1
)
L
(
M
P
a
)
Experimental
Standard Kernel
Complex Kernel
(a) (b)
Figure 3.6 Fitting of the storage, (a), and loss, (b), moduli with the standard linear viscoelastic model
(Standard Kernel ) and with the model proposed (Complex Kernel ).
Equations (3.117)  (3.118) are standard constitutive equations in material rehology. They
correspond to a generalized Maxwell element in parallel to a spring (Fig. 3.5a). By integrating
(3.118) and substituting into (3.117), the linear viscoelastic model (3.112) is recovered.
Appropriate models to describe the linear frequency dependence of dynamic moduli can
be formulated by considering the constitutive equation of a linear hysteretic damper (Inaudi
& Makris, 1996), which in Fourier domain reads as
(j ) = k [1 j sign()] (j ) (3.119)
In Eq. (3.119), , are the Fourier transforms of stress and defomation, while j =
1 and
is a frequency independent loss factor; sign() is the signum function, i.e., sign(x) = 1 if
x > 0, sign(x) = 1 if x < 0, sign(x) = 0 if x = 0. A timedomain representation of the
element (3.119) would involve the Hilberttransform (Inaudi & Kelly, 1995).
In order to describe the steady state response of mechanical systems subjected to har
monic deformations and, hence, attention has been focused on the frequency domain. The
generalized Maxwell model, governed by the evolution equation (3.118), can be modied by
introducing N hysteretic blocks in parallel, as shown in Fig. 3.5b. The constitutive equation
of the ith block is given by:
i
ov
(j)
(j)
= j k
0
(1 +
i
i
)
k
i
1 +
2
2
i
+ k
0
(
i
+
i
)
i
k
i
1 +
2
2
i
, (3.120)
which arises from the series of a hysteretic spring and dashpot. The constitutive equation of
the modied Maxwell element shown in Fig. 3.5 is
(j) =
el
(j) +
N
i=1
i
ov
(j) , (3.121)
with k
= k
0
_
1
N
i=1
k
i
_
. Indeed, the constitutive assumptions (3.120) are equivalent to
introduce a complex stiness
k into Eq. (3.112). As a consequence, the hypotheses underlying
Eqs. (3.105) are no longer valid.
In the case of linear stressstrain constitutive model, the storage and loss moduli are
obtained as the real and imaginary parts of the frequency response function (Inaudi & Makris,
78 Chapter 3. Nonlinear Viscoelasticity
1996), i.e.,
S() = k
0
_
1
N
i=1
k
i
+
N
i=1
(
i
+
i
)
i
k
i
1 +
2
2
i
_
,
L() = k
0
N
i=1
_
(1 +
i
i
)
k
i
1 +
2
2
i
_
,
(3.122)
which, for
i
0, coincides with the standard dynamic moduli (3.116).
The results of the tting are shown in Fig. 3.6 for N = 4 hysteretic dampers in parallel
with a spring. A more accurate match with the experimental data is achieved than with the
standard model. Moreover, the initial slope of the storage modulus is reproduced even with
a small number of paramenters.
Chapter 4
Model Identication
Chapter Outline. In the rst section, the standard identication procedure of the material parameters
for a nonlinear viscoelastic (NLV) constitutive equation is discussed thoroughly. This routine relies on the
separate identication of the instantaneous (elastic) and dissipative (inelastic) parts. The advantages and
disadvantages of this approach are evaluated by considering Fungs constitutive model. Thereafter, a joint
identication of the elastic and dissipative parts is introduced. To this end, the constitutive equation is
rewritten in a form suitable for encompassing all the single hereditary models. Therefore, the identication
of the constitutive coecients is reduced to the solution of a nonlinear optimization problem. An iterative
technique allowed the rening of the relaxation times around the most signicant time constants. The results
of the compression tests with cylindrical carbon blacklled rubber specimens, already summarized in Chapter
1, are here used to compare the main features of each model.
80 Chapter 4. Model Identication
4.1 Position of the Problem 81
4.1 Position of the Problem
Owing to the importance of elastomer components in engineering applications, the accurate
prediction of their mechanical behavior under operational conditions is a relevant subject in
industry. Often the material is exposed to dierent short and longtime loads simultaneously.
Thus it is necessary to provide constitutive models able to predict the material response for
several loading conditions. Moreover, the models must be endowed with a set of parameters
which accurately represents the material behaviour over the entire working range.
Despite the Mullins eect, in most rubberlike materials no further damage can be ob
served, even at large deformations. The consideration of predamaged materials, which are
in a stable stationary state with respect to the Mullins eect, leads to the choice of a nite
viscoelastic constitutive equation.
In the following, the QLV model introduced by Fung (1972) is used to illustrate a possible
procedure for the material parameters identication. Fungs model reads as (see Section 3.3.1)
T(t) = T
e
(t) +
_
t
0
k(t s)T
e
(s) ds,
T
e
(t) =
0
(t)I +
1
(t)C(t) +
1
(t)C
1
(t),
(4.1)
where
(t) =
(I
1
(t), I
2
(t), I
3
(t)), = 1, 0, 1, is a function of the strain invariants I
1
, I
2
and I
3
at time t.
The standard identication algorithm relies on three successive steps:
1. the identication of the instantaneous response T
e
(t);
2. the solution of the linear integral equation for k(t);
3. the identication of the viscoelastic kernel parameters, e.g., the coecients g
i
,
i
in the
case of the Prony series for k(t).
In the next sections each of the steps 1, 2 and 3 are discussed thoroughly. However,
this procedure requires experimental measurements which are very dicult to perform: a
very fast loading ramp is necessary for the identication of the instantaneous stress T
e
,
while the viscoelastic kernel necessitates very long relaxation tests. As a consequence of the
restrictions imposed by the experimental equipment on the maximum exerting strain rate, the
experimental procedure actually employed in the laboratory leads to approximations which
are not satisfactory when nonlinear viscoelastic materials are involved.
In order to overcome these limitations, a combined experimentalnumerical analysis based
on an identication routine, both of the long and shorttime material response, is presented.
This procedure actually generalizes the results in Knauss & Zhao (2007) to nonlinear consti
tutive equations.
4.2 Standard Identication Procedure
4.2.1 Instantaneous Response
The instantaneous stress T
e
(t) in equation (4.1) is by denition the tensile stress instanta
neously generated when a step function is imposed on the specimen.
Strict laboratory measurements of T
e
(t) according to this denition are dicult, because
at a sudden application of the load, transient waves are induced in the specimen and the
resulting stress response measurement will be confused by these elastic waves. However, if
82 Chapter 4. Model Identication
the relaxation function k(t) is a continuous function, then T
e
(t) might be estimated by the
tensile stress response in a loading experiment with a suciently high rate of loading.
In order to justify this procedure, let us consider an uniaxial deformation (t) produced
by an uniaxial stress state (T
11
(t) = T(t), T
22
(t) = T
33
(t) = 0). Now, if by some monotonic
process, (t) is increased in a time interval from 0 to () =
, then by integrating by parts
Eq. (4.1) at time t = , the result is:
T() = T
e
_
+
_
0
k(s)T
e
[( s)] ds. (4.2)
Since the relaxation function k(t) is a continuously varying decreasing function, as s
increases from 0 to , the integrand never changes sign; by the mean theorem there exists a
[0, ] such that
T() = T
e
_
+
k() T
e
[( )] . (4.3)
Therefore, being k(t) C
1
(IR) and 0 as 0, one gets
k()
T() T
e
_
, (4.4)
which allows the integral term into (4.2) to be ignored. This result conrms that, in order
to model correctly the elastic part T
e
, very fast loading/unloading experiments should be
preferred to quasistatic tests. Indeed, in the latter case, the material response would involve
both elastic and dissipative eects.
If k(t) is represented in terms of the Prony series, e.g.,
k(t) =
N
i=1
k
i
+
N
i=1
(1 k
i
) e
i
,
N
i=1
k
i
< 1 , (4.5)
the inequality in (4.4) yields the maximum value of , i.e., the minimum strain rate /,
for which T
e
can be estimated as the tensile stress response T. In particular, since k(t) is a
completely monotonic decreasing function, the result is
k()
k(0)
=
N
i=1
k
i
i
<
N
min
, (4.6)
which provides a sucient condition for the inequality (4.4) to be valid, i.e.,
min
N
. (4.7)
Inequality (4.7) states that the duration in time taken for the ramp to reach the stretch
should be much less than the minimum characteristic time of the material. However, since
rubber is wellknown to have characteristic times of the order of a few milliseconds, or even
less, laboratory measurements of T
e
are very dicult.
For modeling purposes, a common choice is to assume for T
e
an hyperelastic constitutive
equation (Hoo Fatt & Ouyang, 2007; Pena et al., 2007; Shim et al., 2004; Yang & Shim,
2004). The relation between the stress and the strain is indirectly specied through the
relation between the internal energy and the strain (Rivlin & Saunders, 1952), i.e.,
W(I
1
, I
2
) =
M
m=0
N
n=0
c
mn
(I
1
3)
m
(I
2
3)
n
, (4.8)
therefore, specifying the constitutive parameters c
mn
into (4.8). As a consequence of the in
compressibility of rubber, the functional W depends only upon the rst and second invariants
of the rightCauchy Green strain tensor.
4.2 Standard Identication Procedure 83
In the case of uniaxial stretch of a thin sheet of material, the constitutive law for T
e
can
be expressed in terms of the stretch and of the strain energy function W as follows (see
Chap. 2),
T
e
=
1
W
; (4.9)
however, from an experimental point of view, it is more convenient to express (4.9) in terms
of the nominal stress, i.e., force over reference area, which is the actual measured quantity,
viz.
e
=
W
=
M
m=0
N
n=0
c
mn
(m, n, ) , (4.10)
where is a nonlinear function of and of the indices m, n, e.g.,
(m, n, ) =2 m
_
2
+
2
3
_
m1
_
1
2
__
2 +
1
2
3
_
n
+
+ 2 n
_
2
+
2
3
_
m
_
2 +
1
2
3
_
n1
_
1
1
3
_
.
(4.11)
Note that, whatever choice of constitutive parameters used in (4.8), the relation between the
stress and the material coecients c
mn
is linear.
In order to achieve a more useful form of (4.10) with respect to the identication process,
it is possible to introduce the vector c IR
mn
of the material coecients, i.e.,
c
T
= c
01
, ..., c
0n
, c
10
, ..., c
1n
, ..., c
m0
, ..., c
mn
, (4.12)
and the matrix M
k mn
(IR),
=
_
_
(0, 1,
1
) (0, 2,
1
) ... ... (m, 0,
1
) ... (m, n,
1
)
(0, 1,
2
) (0, 2,
2
) ... ... (m, 0,
2
) ... (m, n,
2
)
... ... ... ... ... ... ...
(0, 1,
k
) (0, 2,
k
) ... ... (m, 0,
k
) ... (m, n,
k
)
_
_
, (4.13)
where
i
= (i ) is the stretch at time t = i . In the previous equations, the coecient c
00
has been set at zero, which implies that the strain energy is zero in the reference conguration.
With the above notation, the constitutive equation (4.10) can be rewritten as the linear
system
e
= c, (4.14)
where
e
=
e
(t
1
), ...,
e
(t
k
) is the discretetime stress vector.
Since the material is assumed to be incompressible, only two stretches can be varied
separately and, hence, biaxial tests would suce to determine the form of the strain energy
function (4.14), i.e., the constitutive relation (4.10). However, in the present investigation,
only uniaxial tests were carried out, because of the limited availability of biaxial specimens.
Consequently, the resulting identication problem is expected to be illconditioned (Ogden
et al., 2004).
The minimization problem arising from the identication of (4.14) is
min
cIR
mn
_
_
_c
e
_
_
_
2
2
(4.15)
where
e
is the vector of the stress as recorded and sampled by the acquisition equipment,
while  
2
is the standard l
2
norm. A LeastSquare method (see, e.g., Golub, 1996) has been
used to solve (4.15) in Sec. 4.2.4, where results and discussions are reported.
84 Chapter 4. Model Identication
4.2.2 Viscoelastic Kernel
Once the instantaneous response T
e
(t) in equation (4.1) has been identied, the problem of
determining the (normalized) viscoelastic kernel k(t), i.e., solving the integral equation (4.1),
can be addressed.
The onedimensional constitutive equation (4.2) can be rewritten through integration by
parts as
T(t) = k(t)T
e
(0) +
_
t
0
k(s)
T
e
(t s) ds . (4.16)
If T
e
(0) ,= 0, introducing the normalized quantities
T(t) = T(t)/T
e
(0) and D(t) =
T
e
(t)/T
e
(0),
one gets
k(t) =
T(t)
_
t
0
k(s)D(t s) ds, (4.17)
which is a Volterra equation of the second kind for k(t) (Linz, 1985).
A more useful form of (4.17) can be achieved by a standard timediscretization. Indeed,
during an experiment, the displacement and the force are recorded at uniform time intervals,
e.g.,
t
i
= i , i = 0, 1, ..., N, =
t
N
. (4.18)
Since k and T
e
are smooth functions of the time variable t, for small the integral can be
estimated by a nite sum through a quadrature rule (see, e.g., Press, 2007). For the uniform
mesh (4.18) a simple scheme is the trapezoidal rule, which gives:
_
t
i
0
D(t
i
s) k(s) ds =
2
D
i
k
0
+
i1
j=1
D
ij
k
j
+
2
D
0
k
i
; (4.19)
therefore equation (4.17) has the following discrete form:
_
T
0
= k
0
,
T
i
=
_
1 +
2
D
0
_
k
i
+
2
D
i
k
0
+
i1
j=1
D
ij
k
j
, i = 1, ..., N
(4.20)
Since k(0) = k
0
= 1,
T
i
=
T
i
D
i
/2 and the result is
T
i
=
_
1 +
2
D
0
_
k
i
+
i1
j=1
D
ij
k
j
, i = 1, ..., N , (4.21)
which can be easily expressed in terms of the linear system
T = k, (4.22)
where
T =
T
1
, ...,
T
N
, k = k
1
, ..., k
N
and
=
_
_
1/ + D
0
/2 0 0 ... ... ...
D
1
1/ + D
0
/2 0 ... ... ...
D
2
D
1
1/ + D
0
/2 0 ... ...
... ... ... ... ... ...
D
N
D
N1
... D
2
D
1
1/ + D
0
/2
_
_
(4.23)
is a lower triangular matrix.
Through the (4.20), the identication of the viscoelastic kernel k(t) has been reduced to
the solution of the linear system (4.22), which requires the inversion of the matrix .
4.2 Standard Identication Procedure 85
The structural properties of strongly depend upon the experimental test performed to
measure T(t). In the case of a stepstrain relaxation with an innite speed of the initial ramp,
the strain is suddenly increased to the value , kept constant thereafter, e.g.,
(t) =
_
1, t < 0
, t 0
(4.24)
In this case, the derivative of instantaneous stress T
e
is vanishing for t 0,
T
e
(t) =
T
e
()
= 0, t 0. (4.25)
Therefore, for each l D
l
= 0, and the matrix reduces to the identity matrix. In this ideal
situation, measuring the stress response, i.e., measuring
T, is equivalent to measuring the
viscoelastic kernel k(t). However, to catch the time dependence of k for t it would be
necessary to carry out an innitely long experimental test.
In the case of lled rubber, an accurate prediction of k can be obtained with a relaxation
experiment lasting ten minutes or less (see, e.g., Antonakakis et al., 2006). It should be noted
that a premature cut o may give an erroneous limiting value of k(t) for t .
In the laboratory the prescription of (t) in the form (4.24) is necessarily supplanted by a
ramp history wherein the strain increases along a constant strain rate until the predetermined
value is reached at time t
0
, e.g.,
(t) =
_
_
1, t < 0
1 + ( 1)
t
t
0
, 0 t < t
0
, t t
0
. (4.26)
Thus, if the undeformed conguration is stressfree, then the instantaneous stress T
e
0
vanishes
and the equation (4.17) becomes
T(t) =
_
t
0
k(s)
T
e
(t s)ds , (4.27)
which is a Volterra equation of the rst kind (Linz, 1985).
In analogy with the previous case, by using the trapezoidal quadrature rule, the discrete
form of (4.27) can be introduced
_
_
_
T
0
= 0,
T
i
=
2
k
i
+
i1
j=1
k
j
D
ij
, i = 1, ..., N
, (4.28)
where
T
i
= T
i
/
T
e
0
/2
D
i
and
D
i
=
T
e
i
/
T
e
0
. In matrix notation, the last equation can be
rewritten as
T = k (4.29)
where
T =
T
1
, ...,
T
N
, k = k
1
, ..., k
N
and
=
_
_
1/2 0 0 ... ... ...
D
1
1/2 0 ... ... ...
D
2
D
1
1/2 0 ... ...
... ... ... ... ... ...
D
N1
D
N2
... D
2
D
1
1/2
_
_
. (4.30)
86 Chapter 4. Model Identication
Dierentiating Eq. (4.26), the result is
dT
e
(t)
dt
=
T
e
()
=
_
_
0, t < 0
1
t
0
T
e
()
, 0 t < t
0
0, t t
0
, (4.31)
hence,
T
e
(t) ,= 0 if 0 t < t
0
. The bandwidth of the matrix is b = t
0
/, where  is
the oor function. Therefore, the inuence of the history upon the current value of the stress
is limited to the previous b samples.
The 1norm of the matrix is

1
= max
1jN
N
i=1
[
i j
[ =
1
2
+
b
i=1
[D
i
[ (4.32)
and, hence, its condition number
1
depends upon the bandwidth b. For a fastly increasing
ramp, the resulting matrix has a lower condition number and the linear system (4.29) is well
conditioned. On the opposite, a slowly increasing ramp results in a larger condition number
of and measurement errors on T(t) may propagate to the solution k(t) and eventually
invalidate the result (see, e.g., Golub, 1996). In any case, errors due to a wrong tting of the
instantaneous stress T
e
(t) will propagate to k(t).
A frequently employed characterization of the viscoelastic kernel is achieved through si
nusoidal strain histories of frequency overimposed on a static displacement (Lee & Kim,
2001), e.g.,
(t) =
_
_
1, t < 0
1 + ( 1)
t
t
0
, 0 t < t
0
+
1
sin((t t
0
)), t t
0
, (4.33)
In analogy with the previous case, the relation between the viscoelastic kernel and the stress
is governed by the linear system (4.29); however, the derivative of T
e
with respect to the time
does not vanish for t > t
0
, i.e.,
dT
e
(t)
dt
=
T
e
()
=
_
_
0, t < 0
1
t
0
T
e
()
, 0 t < t
0
1
cos((t t
0
))
T
e
()
, t t
0
. (4.34)
In this case, hence, the matrix is a fully lower triangular matrix.
4.2.3 Fitting of the Prony Series
Once the viscoelastic kernel has been obtained by inverting the matrices or , the functional
form of k(t) must be chosen and the constitutive parameters identied accordingly.
A common assumption for k is
k(t; k, ) = 1 +
N
i=1
k
i
_
e
t/
i
1
_
, k = k
1
, ..., k
N
, =
1
, ...,
N
, (4.35)
1
The condition number of an invertible matrix A is dened as product between the norm of A and the
norm of its inverse A
1
(see, e.g., Golub, 1996).
4.2 Standard Identication Procedure 87
which is the wellknown Prony series. This functional form, derived from various molecular
models, has been successfully applied to model relaxation phenomena in lled rubber (e.g.,
Antonakakis et al., 2006; Haupt & Lion, 2002; Johnson et al., 1994; Park & Schapery, 1999;
Quigley et al., 1995).
The minimization problem arising from the identication of a Nterms Prony series is
min
(k,)IR
N
IR
N
_
_
_K(k, )
K
_
_
_
2
2
, (4.36)
where K(k, ) = k(t
1
; k, ), ..., k(t
M
; k, ),
K is the kernel obtained by inverting the matri
ces or and t
1
, ..., t
M
are the discrete time instants.
Ten, twenty or even more terms in the series (4.35) are sometimes used to increase the
accuracy of the tting model. However, this choice leads to illconditioned solutions of the
optimization problem, which as a result are strongly aected by the choice of the initial point.
Owing to the transcendental dependence of k upon the characteristic times, the mini
mization problem (4.36) is nonlinear, thus, nonconvex. A typical procedure, applied in the
literature for reducing the computational eort due to the minimization of a nonconvex func
tional, is to x apriori the relaxation times
i
(see, e.g., Knauss & Zhao, 2007). The result
is a linear optimization problem in the reduced set of variables k.
If =
1
, ...,
N
are the xed characteristic times, then (4.36) reduces to
min
kIR
N
_
_
_Hk
K
_
_
_
2
2
, (4.37)
where
K
m
= (
K
m
1) and h is the M N matrix
h =
_
_
e
t
1
/
1
1 e
t
1
/
2
1 ... e
t
1
/
N
1
e
t
2
/
1
1 e
t
2
/
2
1 ... e
t
2
/
N
1
... ... ... ...
e
t
M
/
1
1 e
t
M
/
2
1 ... e
t
M
/
N
1
_
_
. (4.38)
If
i
,=
j
for each i ,= j, the functions e
t/
1
, ..., e
t/
N
are linearly independent. Hence,
the matrix h is full column rank and h
T
h is invertible. The solution of (4.37) can be achieved
with the standard linear Least Squares algorithm, i.e., by inverting the matrix h
T
h.
The time interval, within which the characteristic times are initially selected, is chosen
according to the duration of the experiment, i.e., the upper limit for
i
, and the sampling
rate of the acquisition channel, i.e., the lower limit. Within this range, the relaxation times are
initially xed at equal intervals along the logarithmic time scale with one or two increments
per decade (Knauss & Zhao, 2007). However, this choice seems to lead to unsatisfactory
results especially when nonlinear viscoelastic models are involved. Therefore, an iterative
tting procedure for rening the estimate of the relaxation times is introduced in Section 4.3.
4.2.4 Identication Results
In the following the identication results obtained through the steps 1,2 and 3 of the above
described procedure are reported.
Compression tests on thick cylindrical specimens were used for the tting (see Section
1.5). In particular, the instantaneous response was obtained exerting a displacement with
a constant strain rate up to the stretch 0.8. The velocity was 500 mm/min which is
approximately fty times the velocity of a standard quasistatic test.
88 Chapter 4. Model Identication
Table 4.1 Identication results of instantaneous term for three dierent material models ( is the condition
number of the matrix
T
).
Model Name c
10
c
01
c
20
c
30
NeoHookean 1.1426 10
6
   1
MooneyRivlin 1.1426 10
6
0   2910
Yeoh 1.4543 10
6
 0.3987 10
7
1.3956 10
7
184880
Yeoh 2 1.3312 10
6
 1.4101 10
6
 315
Three models were considered for the instantaneous term T
e
, respectively with one (c
10

NeoHookean), two (c
10
, c
01
 MooneyRivlin) and three (c
10
, c
20
, c
30
 Yeoh) material coe
cients respectively.
The function lsqlin in the MATLAB Optimization Toolbox was employed to solve nu
merically the problem (4.15). Actually, to guarantee the physical admissibility of the tting
models, the following constraints have been introduced (see pag. 243 in Hartmann, 2001b;
Holzapfel, 2000):
c
10
> 0, c
01
0, c
20
0, c
30
> 0. (4.39)
The results of the identication of T
e
are shown in Fig. 4.1. All the curves can roughly
describe the experimental data. However, to allow a comparison among the dierent models,
the relative errors are also plotted in the gure. In order to avoid a division by small values
of the force when stretches are close to one, the relative error is computed as proposed in
(Ogden et al., 2004), i.e.,
err
i
:=
(c)
i
e
_
i
max
_
0.5,
_
e
_
i
_. (4.40)
The optimum values of the coecients c
ij
together with the condition numbers of the matrix
T
are shown in Tab. 4.1.
For NeoHooke and MooneyRivlin models, the relative error at low strains is larger than
10 %. The NeoHookean and MooneyRivlin curves coincide because from the optimization
problem the result is c
01
= 0. The 3parameters Yeoh model shows a good agreement with
the experimental data with an error less than 20 % for strain lower than 0.95.
It can be seen from Tab. 4.1 that as the number of parameters increases, the condition
number increases greatly, meaning that some of the directions in the space of the parameters
c are not signicant. In the case of the Yeoh model, by a singular value decomposition of the
matrix
T
, the results is that the direction corresponding to the parameter c
03
is associated
with the smallest eigenvalue, meaning that c
03
is not signicant in describing the experimental
data.
To conrm this result, a 2 parameters Yeoh model was also considered (c
30
= 0). While
the relative error remains limited the condition number is drastically reduced. However, this
model with the choice of the coecients reported in Table 4.1 leads to a nonphysical behavior
at larger strains. Therefore this solution must be discarded aposteriori.
After the instantaneous response of the material has been estimated, the relaxation
kernel k can be evaluated by inverting the matrix in (4.29). As a result of the discussed
wellconditioning of , the errors for computing its inverse
1
remain limited.
The matrix is assembled from the relaxation tests (see Sec. 1.5). In this case the cylinder
was compressed to the stretch = 0.15 with a constant strain rate of 100 mm/min. Thereafter
4.2 Standard Identication Procedure 89
M
P
a
data
Y c
30
0
MR
NH
e
r
r
b
Y c
30
0
MR
NH
Figure 4.1 Nominal stress, , plotted against the stretch for NeoHookean, MooneyRivlin and Yeoh
material models (Panel (a)). Since c
01
= 0 in Tab. 4.1, NeoHookean and MooneyRivlin curves
coincide. Panel (b) shows the corresponding relative errors.
90 Chapter 4. Model Identication
0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5
1
0
1
2
3
4
5
M
P
a
Y c
30
0
Y
MR
NH
Figure 4.2 Identication results of the instantaneous response for three dierent material models. A very
dierent behavior is shown outside the range of the experimental data (represented in the gure
as bigdots). The 2parameters Yeoh model exhibits a nonphysical behavior.
the strain was held xed for 100 s. The sampling rate of the external force and displacement
was 100 Hz, corresponding to a period = 0.01 s. Therefore the characteristic times in the
Prony series were chosen equally spaced in logarithmic scale in the range
i
[0.01, 100] s.
After k has been calculated, the minimization problem (4.37) can be addressed.
In this case, the identication results are reported in Figure 4.4 for an increasing num
ber of exponential terms in the Prony series. With N = 10 coecients there is a perfect
agreement between the tting model and the experimental curve, with a relative error below
2%. However, the condition number of the matrix h
T
h increases dramatically (
5
= 127.2,
10
= 2327 and
20
= 4.371 10
7
). Therefore, there might be multiple sets of parameters
k resulting in the same value of the objective function, but which produce a very dierent
behavior outside the observed time range.
The robustness of the solution k has also been investigated. Because of the high values
of the condition number, it is expected that a small error either in the experimental data,
e.g.,
K, or in the model, e.g., h, produces a large variation in the identied parameters. In
particular, by considering the perturbed optimization problem:
min
kIR
N
_
_
_(h + dh(
1
)) (k + dk)
K d
K(
2
)
_
_
_
2
2
, (4.41)
the sensitivity of the minimum can be investigated. In (4.41), dh(
1
) is the variation in the
model matrix induced by a random error of maximum amplitude
1
upon the characteristic
times
i
, i.e.,
1
= (1 +
1
r)
1
, ...,
N
= (1 +
1
r)
N
, (4.42)
with 1 r 1. Moreover, the vector d
K
i
=
2
r
K
i
. (4.43)
Table 4.2 outlines the variation of the minimum solution of (4.41) for increasing error
amplitudes
1
and
2
in the range 0.01, 0.3. It is evident that the solution is strongly
aected by the error amplitudes. Even though the dierent values of k correspond to the
4.2 Standard Identication Procedure 91
T
a
b
l
e
4
.
2
V
a
r
i
a
t
i
o
n
o
f
t
h
e
m
i
n
i
m
u
m
s
o
l
u
t
i
o
n
i
n
t
h
e
c
a
s
e
o
f
N
=
2
0
P
r
o
n
y
s
c
o
e
c
i
e
n
t
s
f
o
r
i
n
c
r
e
a
s
i
n
g
e
r
r
o
r
a
m
p
l
i
t
u
d
e
s
1
a
n
d
2
i
n
t
h
e
r
a
n
g
e
{
0
.
0
1
,
0
.
1
,
0
.
3
}
.
O
n
l
y
t
h
e
t
i
m
e
c
o
n
s
t
a
n
t
s
a
s
s
o
c
i
a
t
e
d
w
i
t
h
t
h
e
l
a
r
g
e
s
t
t
h
r
e
e
c
o
e
c
i
e
n
t
s
k
a
r
e
s
h
o
w
n
i
n
t
h
e
t
a
b
l
e
.
d
h
d
k
k
1
1
k
2
2
k
3
3
[
%
]
[
%
]
[
%
]
[
%
]
[
%
]
0
0
0
0
0
1
2
.
9
1
.
0
3
0
6
0
.
0
2
5
7
1
2
.
9
1
5
0
.
0
1
5
9
8
.
4
7
4
3
0
1
0
0
.
7
1
1
.
4
0
.
1
1
9
2
1
.
0
3
0
6
0
.
0
2
1
2
1
2
.
9
1
5
0
.
0
2
1
0
8
.
4
7
4
3
0
1
0
0
7
1
2
4
.
4
0
.
1
0
3
5
0
.
6
7
6
2
0
.
0
4
6
0
5
.
5
6
0
3
0
.
0
1
3
5
3
0
0
3
0
0
2
1
1
1
8
.
3
0
.
0
6
2
7
0
.
4
4
3
7
0
.
0
5
4
7
0
.
6
7
6
2
0
.
0
2
8
5
8
.
4
7
4
3
1
0
0
.
1
0
2
.
6
0
.
1
2
7
0
1
.
0
4
0
3
0
.
0
2
5
7
1
3
.
0
2
3
0
.
0
1
5
8
8
.
4
6
1
1
1
0
.
1
0
.
7
1
3
.
4
0
.
1
1
7
2
1
.
0
4
0
3
0
.
0
2
1
3
1
3
.
0
2
3
0
.
0
2
0
9
8
.
4
6
1
1
1
0
0
.
1
7
1
2
3
.
3
0
.
1
0
2
3
0
.
6
7
1
5
5
0
.
0
4
5
9
5
.
5
2
0
4
0
.
0
1
3
9
3
0
.
2
7
6
1
3
0
0
.
1
2
1
1
1
8
.
3
0
.
0
6
3
3
0
.
4
4
7
7
9
0
.
0
5
4
0
0
.
6
7
1
5
0
.
0
2
9
0
8
.
4
6
1
1
0
0
1
0
1
8
.
6
0
.
1
1
2
7
1
.
1
2
7
5
0
.
0
2
6
5
1
3
.
9
8
9
0
.
0
2
0
5
0
.
6
2
9
9
1
0
1
1
0
.
7
2
6
.
1
0
.
1
0
5
1
1
.
1
2
7
5
0
.
0
2
6
9
0
.
6
2
9
9
0
.
0
2
2
7
1
3
.
9
8
9
1
0
1
0
1
7
1
0
6
.
3
0
.
0
9
0
4
0
.
6
2
9
9
0
.
0
2
9
6
5
.
1
6
2
0
0
.
0
2
6
4
1
.
1
2
7
5
1
0
3
0
1
2
1
1
1
8
.
8
0
.
0
6
9
7
0
.
4
8
4
9
0
.
0
4
7
4
0
.
6
2
9
9
0
.
0
3
3
6
8
.
3
4
1
7
3
0
0
2
.
8
0
3
6
.
1
0
.
0
9
9
8
1
.
3
2
1
5
0
.
0
3
6
8
0
.
5
6
7
4
0
.
0
2
9
8
1
6
.
1
3
7
3
0
1
2
.
8
0
.
7
4
0
.
8
0
.
0
9
4
2
1
.
3
2
1
5
0
.
0
4
1
1
0
.
5
6
7
4
0
.
0
2
6
7
1
6
.
1
3
7
3
0
1
0
2
.
8
7
9
2
.
8
0
.
0
7
9
8
0
.
5
6
7
4
0
.
0
4
3
9
8
.
0
7
6
5
0
.
0
4
3
5
1
.
3
2
1
5
3
0
3
0
2
.
8
2
1
1
1
9
.
1
0
.
0
6
5
0
0
.
5
6
7
4
0
.
0
5
3
7
0
.
5
3
7
3
0
.
0
4
3
7
8
.
0
7
6
5
92 Chapter 4. Model Identication
1 100 200 300
1
100
200
300
1 100 200 300
1
100
200
300
Figure 4.3 Contour plot of matrix
mn
with m, n 300 for the performed relaxation test. The bandwidth
of the matrix is evident in the gure and corresponds to a rising time of the initial ramp t
0
approximately equal to 0.5 s.
same value of the objective function, they lead to a completely dierent frequency behavior,
as shown in Fig. 4.5.
The material behavior arising from the identied parameters c and k is shown in Fig. 4.6
compared to the experimental data. It is seen that the aforementioned procedure leads to
unsatisfactory results because the material stiness is heavily underestimated and the relax
ation behavior around the stresskink is completely missed. Therefore, in the next section,
an alternative identication method for both the instantaneous and dissipative parts of the
constitutive equation is introduced.
4.3 Joint Identication
In this section a joint identication of the instantaneous and dissipative terms is presented. To
this end a single hereditary formulation for the constitutive equation of a nonlinear viscoelastic
solid is introduced which enables us to include also dierential and fractional dierential
models (see for instance Adolfsson et al., 2005; Lion & Kardelky, 2004).
The hereditary formulation for the current value of the second PiolaKirchho stress T(t)
is introduced as follows:
T(t) = T
(e)
(C(t)) +(C(t))
_
t
0
k(t s)
(t s)
(C(s), C(t))ds, (4.44)
where T
(e)
represents the instantaneous stress as a function of the current value of deforma
tion, , and k account for the memory properties of the material. In particular, is a
tensorvalued function of the current strain value C(t), whilst is a tensorvalued function
which depends on both the current strain C(t) and the strain history C(s).
4.3 Joint Identication 93
0 5 10 15 20 25 30
0.70
0.75
0.80
0.85
0.90
0.95
1.00
t s
k
a
data
N 20
N 10
N 5
0 5 10 15 20 25 30
0.0
0.5
1.0
1.5
2.0
t s
e
r
r
b
N 20
N 10
N 5
Figure 4.4 Viscoelastic kernel k plotted against time t for an increasing number of exponential terms in the
Prony series (Panel (a)). Even with the lowest number of parameters N = 5 the relative error is
less than 2% (Panel (b)).
94 Chapter 4. Model Identication
0 1 2 3 4 5
0.85
0.90
0.95
1.00
s
1
1
0.3,
2
0.3
1
0.01,
2
0.1
1
0,
2
0
0 1 2 3 4 5
0.00
0.01
0.02
0.03
0.04
0.05
0.06
0.07
s
1
1
0.3,
2
0.3
1
0.01,
2
0.1
1
0,
2
0
Figure 4.5 Normalized storage, S, and normalized loss, L, moduli plotted against the frequency for dierent
values of the error amplitudes
1
and
2
in the case of N = 20 Pronys parameters.
4.3 Joint Identication 95
0 5 10 15 20 25 30
1.2
1.0
0.8
0.6
0.4
0.2
0.0
t s
M
P
a
data
fitting
Figure 4.6 Nominal stress plotted against time in the range t [0, 30] s for the collected experimental data
and the Fungs constitutive equation with a NeoHookean elastic part. The model coecients have
been identied through the discussed identication procedure. The material stiness is heavily
underestimated and the relaxation behavior around the stresskink is completely missed
From an experimental point of view, it is convenient to express the force recorded during
the experiment (divided by the reference crosssection area) in terms of the (nominal) strain.
To this aim, the previous constitutive relation can be equivalently expressed in terms of the
rst PiolaKirchho stress tensor:
(t) =
(e)
(t) +F(t) (t)
_
t
0
k(t s)
(t s)
(s, t)ds , (4.45)
with = FT and
(e)
= FT
(e)
.
To express the viscoelastic kernel, the discrete relaxation spectrum, leading to the Pronys
series is used:
k(t) = 1 +
N
i=1
k
i
(e
t/
i
1) ,
N
i
k
i
< 1 ; (4.46)
The number N of exponential terms in (4.46) is generally selected to increase the accuracy of
the tting model. In the literature, up to N = 20 relaxation times were used (Antonakakis
et al., 2006; Knauss & Zhao, 2007); however a large value of N could lead to illconditioned
identication problems where identied parameters show high sensitivity to errors in the
experimental data.
Beside dierential viscoelasticity, accounted for with a discontinuous kernel, the constitu
tive relation (4.44) can also include fractional dierential models. In this case the kernel k is
expressed as (see for instance Adolfsson et al., 2005; Metzler & Nonnenmacher, 2003):
k(t) = 1 +
N
i=1
k
i
(E
i
[(t/
i
)
i
] 1), (4.47)
with
N
i=1
k
i
< 1,
i
> 0, t 0 ,
96 Chapter 4. Model Identication
where E
(u) =
j=0
u
j
(1 + j )
, > 0, (4.48)
and is the Eulerian Gamma function. Thermodynamics laws implies the following restric
tions for the coecients
i
:
0 <
i
1, i = 1, ..., N (4.49)
Remark that N additional constitutive parameters
i
are introduced with respect to the
standard Prony series; when
i
1, the fractional kernel (4.47) reduces to the standard one
(4.46), since E
1
(t) equals the exponential function.
The constitutive relation (4.45) is valid for every NLV material bearing no special char
acteristic. However carbon blacklled elastomers are usually treated as isotropic and incom
pressible materials. The rst property reects on the constitutive relation (4.45) requiring
all the functions,
(e)
, and , to be isotropic tensor functions of C. In particular, the
instantaneous stress, for an incompressible isotropic material, can be expressed as
(e)
(t) =
1
(I
1
(t), I
2
(t))F(t) +
2
(I
1
(t), I
2
(t))F(t)C(t) , (4.50)
where
1
and
2
are functions of the rst two invariants of C at time t, I
1
(t) = tr C(t) and
I
2
(t) = tr C
1
(t).
The incompressibility constraint, det F = 1, can be easily taken into account by adding
to the Cauchy stress a pressure eld p(t)I which is not constitutively assigned, but depends
upon the boundary and initial conditions (Ogden, 1997). By adding the unknown pressure in
the constitutive equation of the rst PiolaKirchho stress, it becomes
(t) = p(t)F
T
(t) +
ES
(t) ; (4.51)
therefore, the extra stress term
ES
(t) remains as the only part of the stressstrain relation
which is constitutively assigned:
ES
(t) =
(e)
ES
+F(t) (t)
_
t
0
k(t s)
(t s)
(s, t)ds ,
(e)
ES
(t) =
1
(t)F(t) +
2
(t)F(t)C(t) .
(4.52)
A suitable choice of the instantaneous stress
(e)
ES
, through the functions
1
and
2
, of the
functions , and of the viscoelastic kernel k in Eqs. (4.52) allows all the models under
investigation to be encompassed. Table 4.3 summarizes such choices for all the considered
constitutive equations.
4.3.1 Constitutive models under consideration
To the author knowledge all the major contributions involving hereditary integral viscoelastic
models have been here considered (refer to Table 4.3). With the intent of having the same
number of parameters for the instantaneous stress
(e)
ES
of each model, a MooneyRivlin rep
resentation has been chosen for the Models 1, 2, 3 and 6, in which the constitutive assumption
regarding
(e)
ES
were not specied by the respective authors. Moreover, Models 2, 4, 5, in their
original versions, had a single relaxation time (N = 1), which has revealed unable to t the
experimental data under investigation. A Pronys series with 3, or more, exponential terms
was alternatively used.
4.3 Joint Identication 97
In particular, we have labeled as number 1 Fungs model (Fung, 1972) which is commonly
referred to as the Quasi Linear Viscoelastic (QLV) model. While this represents one of the rst
major contributions to nite viscoelasticity, its main limit lies in the choice of the function
strictly related to the constitutive function
(e)
ES
, modeling the instantaneous material
response.
Recently, based on Fungs seminal work, many investigators have proposed new integral
type nonlinear viscoelastic constitutive relationships. Among them, the equations introduced
in (Fosdick & Yu, 1998; Hallquist, 1998; Hibbit et al., 2007; Shim et al., 2004; Yang et al.,
2000) are here considered and labeled as 2, 3, 4, 5 and 6 respectively.
In particular we have included several models which were intended to describe the vis
coelastic behavior of carbon blacklled rubber at high strain rates (Hoo Fatt & Ouyang,
2007; Shim et al., 2004; Yang et al., 2000). In all these models the time derivative of the
strain explicitly appears in the hereditary term. The model by Hoo Fatt & Ouyang (2007)
was not considered in the comparison since it has shown some shortcomings owing to a zero
Youngs modulus in the undeformed conguration. For this model the Cauchy stress arising
from a constant strain rate test, say
0
, such that = 1 +
0
t, is
=2
1
(I
1
3)
2
_
_
+
__
1
_
2
1
k(
0
)
_
2
+
2
3
__
1
2
__
d,
and, hence, the Youngs modulus around the undeformed conguration, that is [/]
=1
,
is equal to zero, in contrast with the experimental evidences (see Figs. 4.8 and 4.9).
The comparison includes also the constitutive equations used respectively in the FE com
mercial codes LSDyna (Hallquist, 1998) and Abaqus FEA (Hibbit et al., 2007). The latter
are two of the most used FE codes, which include nite viscoelastic eects. In particular,
the nite viscoelastic relationship (Hallquist, 1998), used in LSDyna, takes into account for
rate eects through a convolution integral which is linear with respect to the strain tensor.
The model used in Abaqus (Hibbit et al., 2007) is similar to a wellestablished model of nite
viscoelasticity, namely the PipkinRogers model (Pipkin & Rogers, 1968); this constitutive
equation, with an appropriate choice of the constitutive parameters, reduces to Fungs QLV
relation (Ciambella et al., 2009).
Finally, the fractional derivatives model introduced in (Lion & Kardelky, 2004) is also
considered and labeled as number 7 in Tab. 4.3. This model, introduced to describe the
amplitude dependence of the storage and loss moduli during dynamic motions, serves here to
investigate the eects of the fractional kernel (4.47) upon the prediction of hysteresis losses.
Since to dene the fractional kernel the additional parameters
i
must be introduced with
respect to the standard Pronys series, to compare properly Lion & Kardelkys model with
the others, the fractional coecients were initially kept xed,
i
= 1, leading to the standard
Pronys series. Thereafter, they were left free to evolve to further minimize the objective
function. The results between the standard and the dierential kernel are then compared in
terms of predicted energy dissipation.
4.3.2 Experimental Setup
For the sake of completeness in this section some of the experimental results already discussed
in Chap. 1 are here summarized.
The possibility of testing the dynamic behavior of the material through relaxation, creep
and fast loading/unloading cycles, that can be performed with a standard tensile machin
ery. While relaxation and creep experiments incorporate the longterm material response,
98 Chapter 4. Model Identication
Table 4.3 Material models under consideration: 1  Fung (1972); 2  Fosdick & Yu (1998); 3  Hallquist
(1998); 4  Yang et al. (2000); 5  Shim et al. (2004); 6  Hibbit et al. (2007); 7  Lion & Kardelky
(2004)
Model
(e)
ES
1
2
1 2 [
1
+
2
I
1
(t)] 2
2
I F
1
(s)
(e)
ES
(s)
2 2 [
1
+
2
I
1
(t)] 2
2
C
1
(t)
1
C(s)C
1
(t)
3 2 [
1
+
2
I
1
(t)] 2
2
I
1
C(s)
4
1
2
I [
1
+
2
I
2
(s)]
C(s)
5
1
2
1 +
1
I
2
(t)
I
1
(s)
I
1
(s)
C(s) + 2
2
C(s)
6 2 [
1
+
2
I
1
(t)] 2
2
I F
1
(s)
(e)
(s)C(s)C
1
(t)
7
1
I
1
C
1
(s)
C(s)C
1
(s)
4.3 Joint Identication 99
19 mm
25 mm lubricant
Figure 4.7 Specimens geometry for compression tests.
the stress arising from loading/unloading cycles at dierent strainrates involves the lowest
characteristic times, thus the highest characteristic frequencies (in our case of the order of 100
Hz). The identication results reported here, could be complemented with dynamic experi
mental results at higher frequencies (see Chap. 1). To this aim small harmonic deformations
superimposed on large static displacements are sometimes used; however the precision re
quested to the testing machinery, and its price, sensibly increase as the frequency of interest
increases. We have, therefore, limited our analysis to experimental setups that can be re
produced in any wellequipped laboratory which does not necessarily have the leadingedge
testing machineries.
Experiments were conducted with a Zwick/Roell z010. The specimen was cylindrically
shaped with a diametertoheight ratio D/H = 0.76, as shown in Fig. 4.7. During all the
experiment the plate were kept lubricated with graphite to guarantee uniform lateral displace
ment over the height and, consequently, avoid barrel deformation of the mantle. Bending and
torsional deformations, if present, were negligible; the strain eld could be reasonably assumed
uniform along the specimen. This turns out to be particularly important since the material
exhibits highly nonlinear behavior: constitutive nonlinearities coupled with nonuniform strain
would be dicult, even impossible, to analyze. All the tests were conducted at an average
temperature of 25. All the samples were subjected to loading/unloading cycles up to 25 %
to eliminate the Mullins eect.
Starting from the initial undeformed conguration, = 1, the specimen was compressed
up to the nal strain = 0.83 in t = 0.7 s with constant strain rate. Thereafter, the
deformation was held xed for 100 s, thus performing the stress relaxation test at a constant
stretch = 0.83. Figure 4.8 shows the recorded stretch and stress histories.
Theoretically the same deformation ramp between = 1 and = 0.83 with an innite
strain rate, t 0, would have allowed the direct measurement of the viscoelastic kernel k(t)
from the stress response, since, in this case, the nominal stress would have been proportional
to the kernel. However, in any actual lab test, to not account for the nite strain rate of
the initial ramp results in an underestimate of the material characteristic times (Antonakakis
et al., 2006).
Relaxation tests allow the capturing of the material behavior involving larger characteristic
times. Since in many engineering applications (e.g., tires, engine mounts, etc.), the shortest
intrinsic times are also signicant, loading/unloading cycles at high strain rate were also
performed. As shown in Fig. 4.9, the loading/unloading path was repeated at four dierent
velocity of the rising ramp in the range
min
= 0.14 s
1
to
max
= 1.09 s
1
. All the loading
paths from the undeformed conguration, = 1, to the maximum strain = 0.83 were
displacement controlled; all the unloading paths were force controlled to the zero force. This
has allowed us to perform after each cycle a 3 seconds creep test to recover the undeformed,
stressfree initial conguration. The timerate of the force controlled unloading paths were
proportional to those of the loading ramps.
Both relaxation and cyclic tests were repeated several times on dierent specimens, but
under the same environmental conditions. For all the repetitions, a good reproducibility has
100 Chapter 4. Model Identication
0 5 10 15 20 25 30
0.8
0.85
0.9
0.95
1
t (s)
0 5 10 15 20 25 30
1
0.5
0
t (s)
1
1
(
M
P
a
)
(a)
(b)
Figure 4.8 Relaxation test: (a) stretch versus time, (b) nominal stress versus time.
0 5 10 15 20 25 30 35
0.8
0.85
0.9
0.95
1
t (s)
0 5 10 15 20 25 30 35
1.5
1
0.5
0
t (s)
1
1
(
M
P
a
)
(a)
(b)
Figure 4.9 Loading/unloading/creep cyclic test at dierent strain rates: (a) strain versus time, (b) nominal
stress versus time.
[0.14, 1.09] s
1
.
4.3 Joint Identication 101
been obtained and the dierent experimental curves practically overlapped. For the sake of
clarity, only one curve is shown in Figs. 4.8 and 4.9.
The experimental results described above can properly be modeled using the standard
solution of the extension of an incompressible cylindrical body. Since for all the experiments
the strain can reasonably be assumed uniform along the specimen, the deformation is described
by
x
1
= (t)X
1
, x
2
= (t)
1/2
X
2
, x
3
= (t)
1/2
X
3
(4.53)
where (t) ( 1) is the stretch ratio in the direction of the uniaxial compression. The resulting
deformation gradient has a diagonal form
F(t) = Diag
_
(t), (t)
1/2
, (t)
1/2
_
. (4.54)
The condition of a vanishing contact force on the lateral mantle determines the Lagrangian
multiplier p(t); eliminating p(t) from Eq. (4.51) and neglecting the inertial forces yields
11
(t) =
(e)
ES 11
(t)
3/2
(t)
(e)
ES 33
(t)
+ (t)
11
(t)
_
t
0
k(t s)
(t s)
11
(s, t) ds
2
(t)
33
(t)
_
t
0
k(t s)
(t s)
33
(s, t) ds, (4.55)
which expresses the relation between the two observable quantities, force per reference area
11
and stretch . Equation (4.55) will be used in the next section to identify the constitutive
parameters for each of the models under investigation.
4.3.3 Identication of Material Parameters
The standard identication procedure of a linearly viscoelastic model is based on the solution
of the static equilibrium equation for the strain history under consideration (see, for instance,
Knauss & Zhao, 2007). This allows the separate identication of the longterm contribution,
associated with the stationary stress in a relaxation test, and thereafter of the shortterm
counterpart.
In the considered nonlinear case, however, the equilibrium equation arising from the con
stitutive equation (4.55) and the appropriate boundary conditions can only be solved numer
ically. Moreover, long and shortterm contributions can not be easily decoupled. Therefore,
the constitutive parameters are identied by minimizing an error functional which is based
on the entire stress history, namely:
min
pP
f(p) = min
pP
_
K
i=1
w
2
i
__
11
(p)
_
i
_
i
_
2
_
, (4.56)
where
is the vector of the experimental stress values recorded at the sampling times
(t
1
, ...t
K
), while
11
(p) is the vector of the model stress values. The latter depends upon
the vector of model parameters p which spans a subset P of IR
P
(P = A + B + 2 N), being
constituted by
p =
1
, ...,
A
,
1
, ...,
B
, k
1
, ..., k
N
,
1
, ...,
N
, (4.57)
with the constraints
i
0,
i
0, 0 k
i
< 1,
N
i=1
k
i
< 1,
i
0 . (4.58)
102 Chapter 4. Model Identication
The vector w = w
1
, w
2
, ..., w
K
contains the scalar weights used to compare properly the
dierent regions in the experimental curves, e.g., the short and the longterm responses. From
the experimental data, initial estimates for the material stinesses
i
have been obtained,
whereas the starting point for the material coecients
i
are chosen to be of the same order
of magnitude.
Due to the nonlinear dependence of
11
on the constitutive parameters
i
,
i
, k
i
and
i
, the resulting minimization problem (4.56) is generally nonconvex. A common choice to
overcome the numerical diculties due to the transcendental dependence of
11
upon the
i
is to x apriori the characteristic times, say
i
=
i
, inasmuch as one or two per decades of
the experimental time range (Knauss & Zhao, 2007). Therefore, the reduced minimization
problem becomes
min
pP
f(
1
, ...,
A
,
1
, ...,
B
, k
1
, ..., k
N
,
1
, ...,
N
) . (4.59)
with
p =
1
, ...
A
,
1
, ...,
B
, k
1
, ..., k
N
P IR
P
.
This apriori choice of the characteristic times, while often used by many authors (Antonakakis
et al., 2006; Knauss & Zhao, 2007), has led us to unsatisfactory identication results. As a
matter of fact, this choice does not account for the wellknown property of carbon blacklled
elastomers of having characteristic times very close to each other (Lion, 1997). Therefore, an
iterative scheme, which actually allows a better estimate of the characteristic times clusters,
is here considered. It consists of the following steps:
1. the minimum characteristic time
min
is chosen according to the acquisition sampling
rate, while the maximum characteristic time
max
is dictated by the duration of the
experiment;
2. within the range [
min
,
max
], N time constants are initially chosen to be equally spaced
in logarithmic scale, i.e.,
(0)
1
=
min
, ...,
(0)
i
=
min
10
(i1)
, ...,
(0)
N
=
max
(4.60)
where is the logarithmic time step, i.e., = [log
10
(
max
/
min
)]/(N 1);
3. once the reduced minimization problem
min
pP
f(
(h)
1
, ...,
(h)
A
,
(h)
B
, k
(h)
1
, ..., k
(h)
N
,
(h)
1
, ...,
(h)
N
) (4.61)
is solved, a new set
(h+1)
i
of characteristic times is considered. In particular every
(h)
i
is discarded, if the corresponding amplitude k
(h)
i
, minimizing (4.61), is below a given
threshold k, whereas if k
(h)
i
overcomes the threshold, the corresponding characteristic
time is rened and split into:
(h+1)
i
(h)
i
10
(
3
)
h
, 1, 10
(
3
)
h
. (4.62)
The starting point of the minimization (4.61) at step h is taken to be the solution at the
previous step (
(h1)
i
,
(h1)
i
, k
(h1)
i
). The initial value of the scalar coecients k
(h)
i
, corre
sponding to the newly added time constants, are set to zero. Step 3 is repeated either until
all the scalar coecients k
i
exceed k, or until the decrement of the objective function between
successive steps vanishes.
4.3 Joint Identication 103
10
1
10
0
10
1
1.2
1
0.8
0.6
0.4
0.2
0
t (s)
1
1
(
M
P
a
)
Experimental
Model 1
Model 2
Model 3
10
0
1.3
1.2
1.1
1
(a)
10
1
10
0
10
1
1.2
1
0.8
0.6
0.4
0.2
0
t (s)
1
1
(
M
P
a
)
Experimental
Model 4
Model 5
Model 6
Model 7
10
0
1.3
1.2
1.1
1
(b)
Figure 4.10 Identied nominal stresses
11
versus logarithmic time for the relaxation test: (a) Models 1, 2,
3 and (b) Models 4, 5, 6 and 7. In the nested graph, the tting models are shown in the range
t [0.4, 1].
To solve numerically the nonlinear minimization problem (4.61), the combination of the
functions patternsearch and fmincon in the Matlab Optimization Toolbox were used. In
particular, the rst is a Pattern Search (PS) algorithm (Audet & Dennis, 2003); since it
does not require the evaluation of the gradient of the objective function, it turns out to
be particularly useful when minimizing highgradient functions or functions with multiple
minima. The latter, instead, is based on a standard Interior Point (IP) algorithm. The
combination of PS and IP was needed because the former demonstrated robustness with
respect to the choice of the initial points, but poor convergence properties, the latter converged
faster, but showed high sensitivity with respect to the starting point. Therefore, at each step,
the PS was used to rene the starting point of the IP method.
4.3.4 Results and Discussion
The results of the iterative identication procedure are presented and discussed for all the
models under investigation.
We rstly report on the results concerning the relaxation test. In this case the time range
of the experimental data is [0.01 s, 30 s]; therefore, at iteration 0 the characteristic times were
chosen according to Eq. (4.60):
(0)
1
= 0.01 s,
(0)
2
= 0.547 s,
(0)
3
= 30 s. (4.63)
The smallest characteristic time is equal to the sampling rate of the acquisition channel, which
in turn was chosen to be 100 Hz.
The results of the optimization algorithm are shown in Fig. 4.10. For almost all the
models three steps of the outlined iterative procedure were generally sucient to reach the
stationarity of the objective function. All the models are able to describe the overall behavior
of the material also with a low number of initial time constants, N = 3 in Eq. (4.63). Fungs
104 Chapter 4. Model Identication
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
10
4
10
3
10
2
10
1
t (s)
e
r
r
(
%
)
1
st
iter
2
nd
iter
3
rd
iter
standard
Figure 4.11 Relative error versus time as identied by the standard procedure with N = 27 characteristic
times and by each of the three steps of the proposed iterative procedure starting with N = 3
(Model 4).
Table 4.4 Material parameters as identied through the relaxation test. Only the most signicant charac
teristic times
i
, associated with the highest amplitudes k
i
, are reported. In the last column, the
corresponding value of the objective function is shown.
Model
1
2
1
2
k
I
I
k
II
II
k
III
III
fval
[MPa] [MPa] [MPa] [MPa] [s] [s] [s]
1 1.787 0 1.787 0 0.2941 0.312 0.1984 0.252 0.0280 3.646 7.76
2 0.818 0 9.671 0.3523 0.038 0.2165 7.9 0.1557 0.005 19.04
3 0.674 0 89.73 0.9922 246.9 0.0043 3.646 0.0016 0.082 8.86
4 0 0.731 79.09 0 0.9924 246.9 0.0032 4.507 0.0025 3.646 7.52
5 0 0.865 0 1.946 0.2001 7.9 0.0301 2.08 0.0219 0.082 29.39
6 1.859 0.175 1.859 0.175 0.4777 0.083 0.0249 2.08 0.0156 246.9 9.82
7 1.591 7.036 0.8589 113.9 0.0486 13.84 0.0445 7.9 12.16
4.3 Joint Identication 105
model (Model 1) underestimates the slope and the curvature of the relaxation curve around
the kink; this fact is actually related to an underestimate of the lowest characteristic times.
Figure 4.11 shows the eects of the iterative procedure on the tting results of Model 4.
In order to avoid a division by small values of the force when the stretches are close to one,
the relative error is computed as proposed in (Ogden et al., 2004), i.e.,
err
i
:=
11
(p)
_
i
_
i
max
_
0.5,
_
_
i
_ , (4.64)
and is plotted versus time in Fig. 4.11 for four dierent cases. The nal outcome of the
standard procedure with N = 27 characteristic times, chosen as suggested in (Knauss &
Zhao, 2007), is reported as a solid black line: correspondingly the nal value of the objective
function was 85.87. The outcomes of all the three steps of the proposed iterative procedure,
starting with N = 3, are reported as dashed, dotteddashed and dotted lines, respectively; the
nal value of the objective function in the third step drops to 7.52. It is evident a monotone
reduction of the relative error from the rst to the third iteration. Similar results, not shown
here, has been obtained for all the other models, proving the eectiveness of the proposed
procedure to estimate the characteristic time constants.
Table 4.4 lists the material parameters identied from the relaxation data. On inspec
tion, it is evident that dierent models are able to t correctly the same data by dierent
modalities, which in turn call for dierent stress patterns. Concerning the characteristic times
and amplitudes, modeling the material fading memory, we observe that the identication of
Pronys series is sensibly illconditioned, since several set of characteristic constants are able
to produce a good agreement with the experimental data. We remark that for Models 3 and
4, the optimal constitutive parameters are characterized by a factor max(
i
)/ max(
i
) 100
between the  and stinesses. This happens despite we have chosen the same set of initial
data for all the models which is characterized by max(
i
)/ max(
i
) 1, as deduced from the
experimental curves. Moreover, the optimal values reported in Table 4.4 for Models 3 and
4 happen to be almost insensitive to variations of the initial point in P. We note that, in
these two models, the increase of the stiness values corresponds to
k
i
1. This cir
cumstance can perhaps be explained as in the Fungs model, the only one in which we could
obtain some closedform results, the limit stress value
11
:=
11
(t ) is proportional
to (
1
+
2
) (
k
i
1). Our guess is that the relaxation data are not sucient to avoid the
spurious behavior of Models 3 and 4, which tends to enlarge the stinesses
i
, while letting
k
i
1.
The results of the iterative tting for the loading/unloading cyclic tests at dierent strain
rates are shown in Figs. 4.12 and 4.13. While in Fig. 4.12 the whole stress responses are
plotted against time, the rst and fourth cycles of the same responses are plotted against the
stretch in Fig. 4.13a, b and Fig. 4.13c, d, respectively. The latter gure allows the visualization
of the error in the energy dissipation for all the considered models. Indeed, a small phase
shift in Fig. 4.12 might result in a bad prediction of the dissipated energy along the cycle.
We observe a relevant error for the energy dissipation prediction at low frequencies for
all the models; as soon as the strain rate of the loading path increases, and the highest time
constants become more signicant, a considerable improvement in the tting is obtained.
Figure 4.14 summarizes these results: while, for all the models, the percentage error in the
energy dissipation is about 50% for the rst cycle at the smallest strain rate, Models 3, 4,
6 and 7 are able to predict the exact dissipated energy within a 3% error at high strain
rate. The poor performance of Model 1 (Fungs model) could be due to the strict relation
between the elastic instantaneous stress
(e)
ES
and the integrand of the dissipative part . By
modifying Fungs constitutive equation to have dierent coecients in the two terms, a slight
106 Chapter 4. Model Identication
0 5 10 15 20 25 30 35
1.4
1.2
1
0.8
0.6
0.4
0.2
0
t (s)
1
1
(
M
P
a
)
Experimental
Model 1
Model 2
Model 3
(a)
0 5 10 15 20 25 30 35
1.4
1.2
1
0.8
0.6
0.4
0.2
0
t (s)
1
1
(
M
P
a
)
Experimental
Model 4
Model 5
Model 6
Model7
(b)
Figure 4.12 Identied nominal stresses
11
versus time t for the loading/unloading cycles: (a) Models 1, 2,
3 and (b) Models 4, 5, 6 and 7.
0.8 0.85 0.9 0.95 1
1.5
1
0.5
0
1
1
(
M
P
a
)
Experimental
Model 1
Model 2
Model 3
0.8 0.85 0.9 0.95 1
1.5
1
0.5
0
1
1
(
M
P
a
)
Experimental
Model 4
Model 5
Model 6
Model 7
(a) (b)
0.8 0.85 0.9 0.95 1
1.5
1
0.5
0
1
1
(
M
P
a
)
Experimental
Model 1
Model 2
Model 3
0.8 0.85 0.9 0.95 1
1.5
1
0.5
0
1
1
(
M
P
a
)
Experimental
Model 4
Model 5
Model 6
Model 7
(c) (d)
Figure 4.13 Identied nominal stresses
11
versus stretch the loading/unloading cycles: (a) Models 1, 2,
3 for
= 0.14 s
1
, (b) Models 4, 5, 6, 7 for
= 0.14 s
1
(rst cycle), (c) Models 1, 2, 3 for
= 1.09 s
1
, (d) Models 4, 5, 6, 7 for
= 1.09 s
1
(fourth cycle).
4.3 Joint Identication 107
1 2 3 4 5 6 7
20
0
20
40
60
80
100
Model
e
r
r
(
%
)
1
st
cycle
4
th
cycle
Total
Figure 4.14 Relative error on the dissipated energy for all the considered models in the loading/unloading
tests: black, gray and white bars concern the energy dissipation in the rst cycle, in the fourth
cycle and averaged in the all cycles respectively. A positive error means that the dissipation is
underestimated.
0 5 10 15
10
4
10
3
10
2
10
1
10
0
t (s)
e
r
r
(
%
)
1
st
iter
2
nd
iter
3
rd
iter
standard
Figure 4.15 Relative errors for the tting Model 4 through the loading/unloading cyclic tests.
improvement of the tting model has been obtained. We also recall that when
2
0, the
elastic part of Fungs model reduces to a neoHookean constitutive relation and, in this case,
it can be proved that its steadystate response is characterized by a vanishing dissipation.
Figure 4.15 shows, for the case of loading/unloading cycles, the advantages of the proposed
iterative technique with respect to the standard apriori choice of the material characteristic
times. Again, the solid black curve represents the nal outcome of the standard procedure
with N = 27 Pronys parameters (the nal value of the objective function is 2.60), while
the dashed, dotteddashed and dotted lines represent the rst three steps of the iterative
procedure (the nal value of the objective function drops from 3.03 to 2.44). Since the cyclic
tests at dierent strainrates involve a wide range of characteristic times, the advantages of
the iterative procedure are modest, yet its computational eort remains absolutely equivalent
to the standard apriori choice.
Table 4.5 lists the material constants identied from the loading/unloading cyclic tests. For
almost all the models, the identied numerical values of the constitutive parameters sensibly
diers from the ones listed in Tab. 4.4, which were identied from the relaxation experiment.
It is then interesting the following question: is each set of constitutive parameters, identied
from one experimental test, able to predict the material response of the other experimental
test?
To this aim, Fig. 4.16 shows the material response in loading/unloading cycles as extrap
olated by the constitutive parameters identied through the relaxation test. In particular
panels ab of Fig. 4.16 refer to the rst cycle at low strain rate, while panels cd of the same
gure refer to the fourth cycle at high strainrate. Viceversa, Fig. 4.17 shows the material re
108 Chapter 4. Model Identication
Table 4.5 Material parameters as identied through the loading/unloading cyclic tests. Only the most
signicant characteristic times
i
, associated with the highest amplitudes k
i
, are reported. In the
last column, the corresponding value of the objective function is shown.
Model
1
2
1
2
k
I
I
k
II
II
k
III
III
fval
[MPa] [MPa] [MPa] [MPa] [s] [s] [s]
1 4.106 0 4.106 0 0.1575 1.187 0.133 0.144 0.1299 2.08 10.11
2 0.881 0 13.08 0.3393 0.006 0.2702 0.001 0.0956 0.005 3.15
3 0.003 0 545.9 0.9156 246.9 0.0005 0.960 0.0002 0.082 2.48
4 0 0.909 0.934 0.002 0.5002 2.08 0.1485 1.187 0.1234 0.082 2.44
5 0 0.955 0 0.199 0.5569 2.08 0.2678 0.082 0.0624 0.960 4.51
6 9.778 0 9.778 0 0.1249 0.006 0.1151 0.005 0.0757 0.082 2.78
7 1.736 0.718 0.6025 2.08 0.103 0.960 0.0760 0.082 2.59
0.8 0.85 0.9 0.95 1
1.5
1
0.5
0
1
1
(
P
a
)
Experimental
Model 1
Model 2
Model 3
0.8 0.85 0.9 0.95 1
1.5
1
0.5
0
1
1
(
P
a
)
Experimental
Model 1
Model 2
Model 3
(a) (b)
0.8 0.85 0.9 0.95 1
1.5
1
0.5
0
1
1
(
M
P
a
)
Experimental
Model 4
Model 5
Model 6
Model 7
0.8 0.85 0.9 0.95 1
1.5
1
0.5
0
1
1
(
M
P
a
)
Experimental
Model 4
Model 5
Model 6
Model 7
(c) (d)
Figure 4.16 Panels (a) and (b) show the hysteresis loop at
= 0.03 s
1
and
= 0.3 s
1
, respectively as
predicted by the tting models inferred through the relaxation test.
4.3 Joint Identication 109
10
1
10
0
10
1
2.5
2
1.5
1
0.5
0
1
1
(
P
a
)
Experimental
Model 1
Model 2
Model 3
(a)
10
1
10
0
10
1
1.4
1.2
1
0.8
0.6
0.4
0.2
0
1
1
(
P
a
)
Experimental
Model 4
Model 5
Model 6
Model 7
(b)
Figure 4.17 Panels (a) and (b) show the relaxation curves as predicted by tting models inferred through
the cyclic test.
sponse in the relaxation test as extrapolated by the constitutive parameters identied through
the loading/unloading cycles.
Table 4.6, which refers to the data reported in Fig. 4.16, summarizes the percentage errors
in the crosspredictions of the energy dissipation when the material is identied through
the relaxation test. Models 1, 2, and 5 show relevant errors at both low and high strain
rates. In particular, we observe that i) Model 1 exhibits the already discussed tendency to
underestimate signicantly the dissipation, ii) Model 2 completely misses the correct material
behavior by predicting positive tension stresses under compression strains. On the other hand,
Models 3, 4, 6 and 7 perform extremely well by predicting the correct energy dissipation
within a percentage error of 67% at low strain rates and within a percentage error of 15% at
high strain rates. These two numbers must be compared with the direct prediction errors of
Fig. 4.11, i.e., about 50% and 3%, respectively. Models 4, 6 and 7 give also good estimates of
the loading and unloading tangent stinesses.
Similarly, Tab. 4.7, which refers to the data reported in Fig. 4.17, summarizes the percent
age errors in the crosspredictions of maximum and limit stress in the relaxation curve when
the material is identied through the cyclic tests. Apart from Model 1, which completely
misses the correct behavior, the other models are able to predict the main features of the
relaxation test within a 10% error. Model 2 distinguishes by predicting the maximum stress
with 1.5% error, while particularly good estimates of the limit stress are given by Models 3,
4, 6 and 7.
Finally, we discuss the fractional kernel (4.47) introduced in (Lion & Kardelky, 2004)
(Model 7). As discussed before, up to now the additional parameters
i
have been kept xed
to allow a fair comparison with the other models. Now these additional parameters are left
free to evolve to further minimize the objective function. We remark that in the case of the
fractional kernel (4.47), the computational eort in the minimization process considerably
increases. Indeed, if the standard viscoelastic kernel is used, since at each step the coecients
i
are kept xed, the resulting objective function is linear with respect to the unknowns
i
and bilinear in
i
and k
i
. Instead, when using the fractional model, the kernel function k(t)
110 Chapter 4. Model Identication
Table 4.6 Percentage errors in the crossprediction of energy dissipation when the material is identied
through the relaxation test. A positive value means an underestimate.
Model 1
st
cycle [%] 4
th
cycle [%] Total [%]
1 93 95 96
2 79 48 56
3 66 13 43
4 66 15 44
5 71 33 56
6 67 12 44
7 66 11 43
Table 4.7 Percentage errors in the crosspredictions of maximum and limit stress when the material is iden
tied through the cyclic tests.
Model max 
11
 [%]
11
(t = 30 s) [%]
1 78 16
2 1.5 7.9
3 11 1.5
4 10 5.5
5 8.4 10
6 10 5.7
7 9.9 6.4
4.3 Joint Identication 111
Table 4.8 Percentage errors of Model 7 with
i
= 1 and
i
=
i
. The rst two columns show the errors
obtained by inferring the material parameters through the cyclic tests, the latter columns refer to
the material parameters inferred through the relaxation test.
i
= 1
i
=
i
i
= 1
i
=
i
Direct Prediction CrossPrediction
1
st
Cicle [%] 54.09 54.03 65.53 65.51
4
th
Cicle [%] 0.67 0.57 11.4 11.38
Total [%] 29.2 29.08 42.6 42.5
CrossPrediction Direct Prediction
max 
11
 [%] 9.93 9.91 2.04 2.01
11
(t = 30 s) [%] 6.41 6.4 0.21 0.19
must be reevaluated for each tentative value of the parameters
i
, while the resulting objective
function depends transcendentally on the unknown coecients. The optimal set
i
=
i
for
the fractional kernel is characterized by 0.75 <
i
1; Tab. 4.8 shows the associated results.
We observe modest improvements in the direct and cross predictions of Model 7 for what
concern both the energy dissipation (at low and high strain rates) and the maximum and
limit stresses in relaxation.
4.3.5 Conclusions and perspectives
We have presented a literature survey of the main nonlinear viscoelastic constitutive models
and we have identied their parameters on the basis of two dierent compression tests in the
range /
0
[0.83, 1]. Relaxation and loading/unloading/creep cycles at dierent strain rates
have been performed on a carbon blacklled rubber compound. An iterative identication
procedure has been introduced and has led to a better estimate of the characteristic time
constants with respect to their apriori choice used by many authors. The presented results
allow us to draw the following conclusions:
Model 1 (Fungs model), which was rstly introduced to describe the behavior of soft
biological tissues, has revealed unable to describe the dissipated energy both at higher
and lower strain rates.
Models 3 (Hallquist, 1998), 4 (Yang et al., 2000), 6 (Hibbit et al., 2007) and 7 (Lion
& Kardelky, 2004) have shown the best predicting capabilities: these models are able
to predict the dissipated energy in high strainrate cycles, a crucial quantity in many
engineering and biomedical applications, within a 3% error. The same error raises to
1015% if their constitutive parameters are identied by the only mean of the relaxation
test.
The relaxation data are not sucient to resolve both the
i
stinesses and the sum of
characteristic amplitudes
k
i
for Models 3 and 4. Thus the limit stress value
11
is
112 Chapter 4. Model Identication
correctly predicted by the product of a spuriously large stiness
i
with a spuriously
small term (
k
i
1). The cyclic tests eliminate this indetermination for Model 4 but
not for Model 3 (cfr. Tabs. 4.4 and 4.5).
For the purpose of viscoelastic properties identication, the proposed loading/unloading
cycles revealed more eective with respect to the simple relaxation test; yet they require
the same hardware and the same amount of time to be performed.
All the considered models sensibly underestimate the energy dissipation in loading/unloading
cycles at low strainrates with 50% errors at least; we recall that the constitutive rela
tions 4 and 5 were specically aimed at modeling the lled rubber behaviour at high
strain rates.
At least for the presented experimental cases, the use of additional constitutive pa
rameters modeling a fractional derivative viscoelastic kernel does not introduce sensible
improvements in the predictions.
Moreover we have observed that dierent sets of material parameters (as listed in Tabs. 4.4
and 4.5) are able to cross predict the stress response within an acceptable tolerance. Being
pretty condent that the described minimization procedure has allowed us to bypass meaning
less local minima for the constitutive parameters, we conclude that the presented experimental
tests are not sucient to characterize unequivocally the threedimensional viscoelastic consti
tutive relation. In particular tests involving shear deformations, tensile and mixed deformation
patterns, yet dicult to perform, could usefully complement the presented data.
Finally we have enlightened a general illconditioning in the problem of the identication
of the Pronys series modeling the actual material memory. It could be worth curing this
diculty by reducing the description of the viscoelastic kernel to a minimal set of parameters.
Chapter 5
Numerical Applications
Chapter Outline. In this chapter the behavior of isotropic, almostincompressible, nonlinear elastic and
viscoelastic materials is simulated by means of the ABAQUS FEA code. Simple deformations are considered
and the numerical results are compared with the analytical solutions. Finally, shortcomings of the ABAQUS
nite viscoelasticity model are highlighted and discussed.
114 Chapter 5. Numerical Applications
5.1 Hyperelastic Material Simulation 115
5.1 Hyperelastic Material Simulation
When modeling a nonlinear hyperelastic material in ABAQUS, the program makes the fol
lowing assumptions:
1. the material behavior is elastic,
2. the material behavior is isotropic,
3. the simulation includes nonlinear geometrical eects.
The use of hybridtype elements is highly recommended when dealing with nearincompressible
materials. Such elements are based on mixedtype formulations, where independent interpo
lations for displacement and stress elds are assumed and two sets of governing equations,
both equilibrium and compatibility are enforced in weak form.
As wellknown, in the case of incompressible materials, the volumetric strain locking prob
lem appears, if standard displacementbased FE formulations are employed. In order to
overcome such numerical drawbacks, various mixed/hybrid FE have been proposed in the
literature based on distinct interpolations of displacement and pressure elds. It has been
widely demonstrated that such enhanced FE are able to avoid volumetric strain locking. On
the other hand, since most of the experimental tests are performed on 2D specimens, where
plane stress conditions are satised, the use of 2D displacementbased FE together with the
plane stress assumption is also justied. In this case the determination of the pressure eld
is not aected by the volumetric locking problem (see, e.g., Belytschko, 2000).
In ABAQUS, the identication of the hyperelastic material properties can be performed
on the basis of experimental stressstrain curves. Four dierent tests can be used to get an
accurate evaluation of the material parameters. These are uniaxial, planar (pure shear test),
equibiaxial and volumetric tests. Clearly, when already known, material parameters can be
directly specied in ABAQUS to describe hyperelastic material models. Since the problem of
the identication of material parameters has been already addressed in Chap. 4, the second
task will be used in this chapter.
In order to investigate ABAQUS performances comparison between analytical and nu
merical solutions will be carried out. With the aim of investigating how the adopted solution
scheme inuences the results, standard patch test on single element will be performed and
a comparison among dierent kinds of elements will be addressed (hybrid C3D8H and plane
stress CPS4R).
5.2 Implicit vs. Explicit Formulation
Finite Element Analysis (FEA) involving shorttime dynamical problems with large deforma
tion, quasistatic problems with large deformations and multiple nonlinearities, or complex
contact/impact problems requires the use of either implicit or explicit solution techniques. Ex
amples of these types of simulations are crashworthiness analysis, drop testing, deep drawing,
rolling, extruding, pipe whip, bird strike, fan containment and many more.
The ABAQUS FEA program includes the ability to address both implicit as well as explicit
solutions. Both the solution procedures are based on a numerical time integration scheme to
solve the discrete dynamical equilibrium equations in terms of displacements, velocities and
accelerations, then strains and stresses (Hibbit et al., 2007). Implicit integration schemes
(ABAQUS/Standard uses a HilberHughesTaylor algorithm for implicit integration) assume
a constant average acceleration over each time step, t = t
n+1
t
n
, where t
n
and t
n+1
are
the starting and ending points of the time interval t. The governing equations are solved
and the resulting accelerations and velocities at t
n+1
are calculated. Then the unknown
116 Chapter 5. Numerical Applications
displacements at t
n+1
are determined. Explicit integration schemes (ABAQUS/Explicit uses
a Central Dierence method) assume a linear change of the displacement in each time step.
The governing equations are evaluated and the resulting accelerations and velocities at t
n
are
calculated. Then, the unknown displacements at t
n+1
are determined.
There is one major dierence between the two techniques in the equations that are used
to solve for displacements at t
n+1
. The implicit solution method requires matrix inversion of
the structural stiness matrix, the explicit solution does not. However, unlike the implicit
solution scheme, which is unconditionally stable independently on the time step size, the
explicit scheme is stable only for time step size smaller than a critical size evaluated for the
analyzed structure. The undamped critical time step size is 2/
n
(where
n
is the largest
natural circular frequency), which is usually a very small value. This very small time step size
requirement for stability thereby makes explicit solutions useful only for very short transient
analyses. But, even though the number of time steps in an explicit solution may be orders
of magnitude greater than in an implicit solution, it is signicantly more ecient than an
implicit solution since no matrix inversion is required. Therefore the choice of the integration
scheme strongly relies on the problem under investigation.
5.3 Static Analysis
In this section some of the results obtained in Chap. 2 for the simple deformation of hyper
elastic materials are recalled. A comparison between analytical and numerical solutions are
presented and discussed. The FE analysis was performed using the C3D8R element, which is
a general purpose linear brick element, with reduced integration (1 integration point).
5.3.1 Uniaxial Extension
In an uniaxial extension test of an isotropic compressible solid, the (homogeneous) strain
tensor is described by
F = Diag
1
,
2
,
2
, (5.1)
where
1
and
2
are the longitudinal and lateral stretch, respectively. The stretch
1
is known,
since it depends on the imposed displacement eld, but
2
has to be determined from the
boundary condition. In particular, since
22
=
33
= 0, the stretch
2
can be derived from
the implicit relation:
22
(
1
,
2
) = 0 (5.2)
From equations (2.34)(2.35) the Cauchy stress arising from the deformation (5.1) is
= 2
0
I + 2
1
B
D
+ 2
2
_
BB
_
D
(5.3)
where
0
= I
1/2
3
V
I
3
,
1
= I
1/2
3
_
I
I
1
+ I
1
I
I
2
_
,
1
= I
1/2
3
I
I
2
, (5.4)
and ()
D
is the deviatoric part of the tensor.
The dependence of Psi
I
and Psi
V
, i.e., of the isochoric and deviatoric parts of the strain
energy function, is known once a constitutive equation is chosen. In the following three
dierent nonlinear models (NeoHooke, MooneyRivlin, Yeoh) are analyzed.
The NeoHookean constitutive law reads as:
I
(I
1
, I
2
) = c
10
_
I
1
3
_
,
V
(I
3
) =
1
D
_
I
1/2
3
1
_
2
,
5.3 Static Analysis 117
thus,
=
2
D
(I
1/2
3
1)I + 2 c
10
I
1/2
3
B
D
. (5.5)
and the initial shear and bulk moduli (moduli around the reference conguration) are given
by:
0
= 2 c
10
, k
0
= 2/D. (5.6)
From equation (5.5) one gets the following expression for the Cauchy stress components:
11
=
4
3
(
2
1
2
2
)(
1
2
2
)
5/3
+
2
D
(
1
2
2
1), (5.7)
22
=
33
=
2
3
(
2
2
2
1
)(
1
2
2
)
5/3
+
2
D
(
1
2
2
1), (5.8)
Using a numerical root nding method, the lateral strain
2
can be obtained from equation
(5.2) for a given longitudinal strain
1
. We only consider the socalled simple materials, where
the implicit relation (5.2) has only one root for a given
1
(Eihlers & Eppers, 1998).
In the case of linear elasticity, the Poisson ration is known to be a material constant. If a
solid is assumed to be incompressible, is equal to 1/2. In what follows, we dene a nonlinear
Poisson function via the longitudinal and lateral strain (see, e.g., Beatty & Stalnaker, 1986),
viz.
=
2
1
1
1
. (5.9)
Obviously, is not a constant in the case of nite elasticity.
The wellknown relation between the initial shear modulus
0
, the initial bulk modu
lus k
0
and the initial Poisson function
0
follows from a linearization around the reference
conguration of equation (5.3) (Eihlers & Eppers, 1998), i.e.,
0
=
3k
0
2
0
6k
0
+ 2
0
(5.10)
The values of the nonlinear Poisson function (5.9) have been plotted in gure (5.1) for a
NeoHookean material and for dierent values of D and increasing
1
.
The same results have been obtained for MooneyRivlin and Yeoh models, for which the
strain energy functions are:
I
(I
1
, I
2
) = c
10
_
I
1
3
_
+ c
01
(I
2
3), (5.11)
V
(I
3
) =
1
D
_
I
1/2
3
1
_
2
, (5.12)
and
I
(I
1
, I
2
) = c
10
_
I
1
3
_
+ c
20
(I
2
3), (5.13)
V
(I
3
) =
1
D
_
I
1/2
3
1
_
2
, (5.14)
respectively
From equation (5.3), the following stressstrain relations are valid:
=
2
D
(I
1/2
3
1)I + 2 I
1/2
3
_
c
10
+ c
01
I
1
_
B
D
2 c
01
I
1/2
3
(B
2
)
D
. (5.15)
for a MooneyRivlin material, and
=
2
D
(I
1/2
3
1)I + 2 I
1/2
3
_
c
10
+ c
20
(I
1
3)
B
D
. (5.16)
118 Chapter 5. Numerical Applications
10
6
10
7
10
8
5 10
8
5 10
9
10
10
Figure 5.1 Nonlinear Poisson function for a NeoHookean material with dierent values of the compressibility
coecient D for the analytical model (solid curves) and the ABAQUS FEA model (dotted curves).
The two curves perfectly overlap.
Yeoh
MR
NH
Figure 5.2 Nonlinear Poisson function in the case of dierent material models for the analytical model (solid
curves) and the ABAQUS FEA model (dotted curves) with D = 10
8
. The analytical and
numerical curves perfectly overlap.
5.3 Static Analysis 119
1
1
0
10
6
10
7
10
8
5 10
8
5 10
9
10
10
Figure 5.3 Normalized uniaxial stress
11
/
0
plotted against stretch
1
for a NeoHookean material with
dierent values of the compressibility coecient D for the analytical model (solid curves) and the
ABAQUS FEA model (dotted curves). The analytical and numerical curves perfectly overlap.
for a Yeoh material.
Moreover, the initial moduli are:
0
= 2 (c
10
+ c
01
), k
0
= 2/D (5.17)
from equations (5.15), and
0
= 2 c
10
, k
0
= 2/D (5.18)
from (5.16).
The results of the combined numerical/analytical simulations are shown in Figs. (5.1)
(5.4). The constitutive parameters were chosen as
c
10
= c
01
=
0
/4, (5.19)
for MooneyRivlin model and
c
20
=
0
/4 (5.20)
for Yeoh model. In all the cases the volumetric constant was D = 10
8
.
Figure 5.1 shows the nonlinear Poissons function in terms of the stretch
1
in the case
of a NeoHookean material model. For every value of the compressibility constraint D in the
range D
_
10
10
, 10
6
1
1
0
Yeoh
MR
NH
Figure 5.4 Normalized uniaxial stress
11
/
0
plotted against stretch
1
for dierent analytical models (solid
curves) and the corresponding ABAQUS FEA models (dotted curves) in the simple tension case.
The two curves perfectly overlap.
5.3.2 Simple Shear
A simple shear deformation has been described in Chap. 2. Recalling those results, the
deformation gradient is given by:
F =
_
_
1 0
0 1 0
0 0 1
_
_
(5.21)
thus, det F = 1.
The boundary conditions applied to subject a solid to a simple shear deformation state are such
that the deformation (5.21) is volumepreserving, independently on the assumed constitutive
equation.
We refer here and henceforth to the analytical solutions reported in table (2.2). The
comparison between numerical and analytical results are plotted in gures (5.5) and (5.6).
The material parameters are the same as in the uniaxial extension case. The volumetric
constant D has been set to D = 10
8
, which for
0
= 10
6
corresponds to = 0.495.
5.4 Dynamic Simulations
In this section we derive analytical solutions for the two boundaryvalue problems considered
also in the static case: simple shear and uniaxial extension.
We henceforth focus on an incompressible viscoelastic solid for which the instantaneous
response is modeled by a two terms Yeoh stressstrain relationship (a standard implementation
in ABAQUS):
e
= pI +
0
(1 3 + I
1
) B, (5.22)
where
0
and are positive constants (
0
is the shear modulus in the reference conguration)
and I
1
(t) = trB(t) . Also, for simplicity, the time relaxation of the solid is assumed to be
5.4 Dynamic Simulations 121
1
1
0
Yeoh
MR
NH
Figure 5.5 Normalized stress
11
/
0
plotted against stretch for dierent analytical models (solid curves)
and the corresponding ABAQUS FEA models (dotted curves) in the simple shear case.
1
2
0
Yeoh
MR
NH
Figure 5.6 Normalized shear stress
12
/
0
plotted against stretch for dierent analytical models (solid
curves) and the corresponding ABAQUS FEA models (dotted curves) in the simple shear case.
122 Chapter 5. Numerical Applications
governed by a oneterm Prony series expansion given by:
G(t) =
0
+
_
1
0
_
e
t/
, (5.23)
where
is the asymptotic value to which the shear modulus settles after an innite time
and is a characteristic time constant.
For the sake of clarity, let recall the ABAQUS FEA nite viscoelastic model already
introduced in Section 3.3.3; for an incompressible solid it reads as
(t) = p(t)I +
1
(t)B(t) +
2
(t)B(t)
2
+
2
i=1
SYM
_
F(t)
__
t
0
G(t s)
i
(s)C(s)
i
ds
_
F
1
(t)
_
,
(5.24)
where p(t) = p(t) +
_
t
0
0
)
SYM
_
F(t)
__
t
0
e
(ts)/
(1 3 + I
1
(s)) C(s)ds
_
F
1
(t)
_
.
(5.25)
By contrast, for the Fungs QLV model introduced in Sec. 3.3.1, the result is
0
(s, t s) =
0
[1 3 + I
1
(s)] G(t s),
1
(s, t s) = 0,
2
(s, t s) = 0, (5.26)
yielding
(t) =p(t)I +
0
[1 3 + I
1
(s)] B(t)
+
_
t
0
e
(ts)/
[1 3 + I
1
(s)] ds B(t).
(5.27)
In order to further highlight the dierences between Abaqus and QLV model, we will
compare the average rate of working in the case of the two simple motions considered here.
To this end, we dene the dissipated energy density E
d
over a period T referred to the current
conguration.
We recall that the internal rate of working of the stress per unit current volume is (t)D(t)
tr (t)D(t). Assuming for D(t) a periodic time law, the timeaveraged work done per unit
current volume over a period T starting at time T
0
, is given by
E
d
=
1
T
_
T
0
+T
T
0
(t) D(t)dt, (5.28)
where T
0
and T
0
+ T are the starting and ending time points.
Since the elastic part of the stress does not contribute to this expression (5.28) represents
the energy dissipated in a cycle (which depends on T
0
in general).
5.4.1 Uniaxial Extension
When an incompressible isotropic solid is subjected to an uniaxial extension, the (homoge
neous) motion is described by:
x
1
= (t)X
1
, x
2
= (t)
1/2
X
2
, x
3
= (t)
1/2
X
3
, (5.29)
5.4 Dynamic Simulations 123
where (t) ( 1) is the stretch ratio in the direction of the uniaxial tension
11
( 0). The
resulting deformation gradient has the following diagonal form:
F(t) = Diag
_
(t), (t)
1/2
, (t)
1/2
_
. (5.30)
Vanishing of the lateral stresses
22
=
33
= 0 determines the Lagrange multiplier p(t)
and elimination of p(t) yields
11
(t) =
0
_
(t)
2
(t)
_
2 + (1 3)(t) +
3
(t)
_
t
0
e
(ts)/
_
(s)
2
(s)
_
2 + (1 + 3)(s) +
3
(s)
ds,
(5.31)
for the Abaqus model and
11
(t) =
0
_
(t)
2
(t)
_
2 + (1 3)(t) +
3
(t)
2
(t)
1
(t)
_
t
0
e
(ts)/
1
(s)
_
2 + (1 + 3)(s) +
3
(s)
ds,
(5.32)
for the QLV model.
The dierence between these models is now clear, and it reects signicantly on the rate of
working, as we conrm numerically below. From an experimental point of view, it is common
practice to employ a dynamic displacement superimposed on a large static deformation. Here
we consider a NeoHookean solid deformed in tension to a stretch of 1.5 from time t = 0 to
time t = 1, and then made to oscillate with a superimposed amplitude such that the stretch
(t) ranges from 1. to 1.5. Thus,
(t) =
_
1 + 0.5t, 0 t 1,
1.5 + 0.2 sin (t 1), t 1.
(5.33)
For numerical purposes the stresses are nondimensionalized by dividing by
0
. For the
remaining parameters we set
/
0
= 0.5, = 0.01 s, = 48 s
1
, (5.34)
except when the frequency dependence of the energy dissipation rate is investigated where
varies.
The numerical results have been obtained with ABAQUS 6.71 using a single C3D8H
element (brick, 8 nodes, trilinear, hybrid with constant pressure) with an implicit solution
scheme. As a check, the same test was done with a 2D plane stress element CPS4 (4 nodes,
bilinear, plane stress) and the same results were obtained. Furthermore a comparison among
implicit solution scheme, explicit solution scheme and analytical solution has been carried out
and the results are shown in Fig. 5.7. In the case of Explicit time integration a C3D8 element
has been used.
The results settle very rapidly to the steady state so that the expression (5.28) becomes
essentially independent on T
0
. Figure 5.8, in which E
d
/
0
is plotted against frequency
in the steady state, shows clearly that the ABAQUS model overestimates the steady state
energy dissipation substantially (by a factor of 4 to 5) compared with the QLV model.
5.4.2 Simple Shear
We now consider a simple shear motion of the form:
x
1
= X
1
+ (t)X
2
, x
2
= X
2
, x
3
= X
3
, (5.35)
124 Chapter 5. Numerical Applications
1
1
0
Figure 5.7 Dependence of the dimensionless axial tension
11
on time t for the Abaqus FEA model (solid
curve), the ABAQUS output (dotted curve) and the QLV model (dashed curve).
0 100 200 300 400 500 600
0
1
2
3
4
s
1
E
d
Figure 5.8 Dependence of the rate of working per unit volume in dimensionless form E
d
/
0
on frequency
over a single period in the steady state for the ABAQUS model (solid curve), the ABAQUS
output (dotted curve) and the QLV model (dashed curve) in the case of uniaxial extension.
5.4 Dynamic Simulations 125
where (t) is the amount of shear strain. The (nonsymmetric) deformation gradient results
as:
F(t) =
_
_
1 (t) 0
0 1 0
0 0 1
_
_
. (5.36)
Simple calculations reveal that a sheared solid described by the ABAQUS model is in a
state of plane stress (
ij
,= 0 for i, j 1, 2;
3j
= 0 for j = 1, 2, 3) when the Lagrange
multiplier p(t) is taken as
p(t) =
0
)e
t/
+
0
2
(t)
_
t
0
e
(ts)/
2
(s)ds. (5.37)
Note that the combination of plane strain and plane stress (in the (1, 2) plane) is permissible
for an incompressible material provided p(t) (and hence
11
(t) and
22
(t)) is (are) adjusted
accordingly. Then, for the ABAQUS model we obtain
11
(t) =
0
2
(t)
_
1 +
2
(t)
(t)
_
t
0
e
(ts)/
(s)
_
1 +
2
(s)
ds,
12
(t) =
0
(t)
_
1 +
2
(t)
+
+
(
0
)
2
_
t
0
e
(ts)/
(s)
_
2 + (s)(t)
2
(t)
_
1 +
2
(s)
ds,
22
(t) =
(
0
)
_
t
0
e
(ts)/
(s) [(s) (t)]
_
1 +
2
(s)
ds,
(5.38)
and for the QLV model
11
(t) =
0
2
(t)
_
1 +
2
(t)
2
(t)
_
t
0
e
(ts)/
_
1 +
2
(s)
ds,
12
(t) =
0
(t)
_
1 +
2
(t)
(t)
_
t
0
e
(ts)/
_
1 +
2
(s)
ds,
22
(t) = 0.
(5.39)
Note that Rivlins universal relation
11
22
=
12
from isotropic elasticity holds also for
the present Yeohbased viscoelastic model. The expressions in (5.39) are expected intuitively
for the Yeoh model because for its instantaneous response at very short times the Cauchy stress
has components
e 11
=
0
2
(1 +
2
) ,
e 12
=
0
(1 +
2
),
e 22
= 0. By contrast, the
ABAQUS model gives rise to a
22
component generated purely by the viscoelastic eects,
in which case the relaxation process creates such a component ex nihilo! However, since
simple shear is a displacement controlled motion the stress components adjust automatically
to accommodate the geometry and they are therefore very much dependent on the form of
the constitutive law. The dierence in the shear stress, however, has serious consequences for
the rate of working, and hence the dissipation, because here
D(t) =
_
_
0 (t)/2 0
(t)/2 0 0
0 0 0
_
_
, (5.40)
and hence D =
12
(t) (t), which is obviously not the same for both models.
We conrm these ndings by testing the ABAQUS software against the formulas (5.38)
(5.40). We take the amount of shear (t) to vary as
(t) =
_
t, 0 t 1
1 + 0.2 sin [ (t 1)] , t 1,
(5.41)
126 Chapter 5. Numerical Applications
1
2
0
Figure 5.9 Dependence of
12
/
0
on time t in simple shear state for the ABAQUS analytical model (solid
curve), the ABAQUS output (dotted curve) and the QLV model (dashed curve).
with the other parameters given by (5.34). Figures 5.9 and 5.10 display the variations of
the
12
and
22
components computed from (5.38) and (5.39) in dimensionless form. As in
the case of simple tension there is not a great dierence in the active stress (in this case
12
)
between the two models, but the reactive stress
22
is very dierent. Finally, as also for simple
tension, we nd that the ABAQUS model overestimates the rate of working with respect to
the QLV model (by a factor of 2 to 3), as shown in Fig. 5.11, again in dimensionless form
with E
d
/
0
plotted against frequency in the steady state. Using the ABAQUS software,
we recovered the thick curves (in Figs. 5.9 and 5.10), which conrms that equation (5.24) is
actually implemented in the ABAQUS code.
A comparison between implicit and explicit solution scheme has been performed. In the
latter case a C3D8 element has been used (C3D8H is not available in explicit time integration).
We nd that the explicit solution scheme overestimates the rate of working with respect to
the implicit/theoretical model, as shown in Fig. 5.11, again in dimensionless form with E
d
/
0
plotted against frequency in the steady state.
5.4 Dynamic Simulations 127
2
2
0
Figure 5.10 Dependence of
22
(t)/
0
on time t in simple shear for the ABAQUS FEA model (solid curve),
the ABAQUS output (dotted curve) and the QLV model (dashed curve). The latter is zero for
all times.
0 100 200 300 400 500 600
0
1
2
3
4
5
s
1
E
d
s
.
1