Viscoelastic Models for Nearly Incompressible Materials

by Mike Scheidler

ARL-TR-4992 September 2009

Approved for public release; distribution is unlimited.
NOTICES

Disclaimers

The findings in this report are not to be construed as an official Department of the Army position unless so designated

by other authorized documents.

Citation of manufacturer’s or trade names does not constitute an official endorsement or approval of the use thereof.

Destroy this report when it is no longer needed. Do not return it to the originator.
Army Research Laboratory

Aberdeen Proving Ground, MD 21005-5066

ARL-TR-4992 September 2009

Viscoelastic Models for Nearly Incompressible Materials

Mike Scheidler

Weapons and Materials Research Directorate, ARL

Approved for public release; distribution is unlimited.
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering 

and maintaining the data needed, and completing and reviewing the collection information. Send comments regarding this burden estimate or any other aspect of this collection of information, 

including suggestions for reducing the burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson 

Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to 

comply with a collection of information if it does not display a currently valid OMB control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY) 2. REPORT TYPE 3. DATES COVERED (From - To) 

5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

4. TITLE AND SUBTITLE 

5c. PROGRAM ELEMENT NUMBER 

5d. PROJECT NUMBER 

5e. TASK NUMBER 

6. AUTHOR(S) 

5f. WORK UNIT NUMBER 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION 

 REPORT NUMBER 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR'S ACRONYM(S) 

11. SPONSOR/MONITOR'S REPORT 

 NUMBER(S) 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

15. SUBJECT TERMS

16. SECURITY CLASSIFICATION OF: 

19a. NAME OF RESPONSIBLE PERSON

a. REPORT b. ABSTRACT c. THIS PAGE 

17. LIMITATION 

OF ABSTRACT 

18. NUMBER 

OF PAGES 

19b. TELEPHONE NUMBER (Include area code)

Standard Form 298 (Rev. 8/98) 

Prescribed by ANSI Std. Z39.18

September 2009 Final

ARL-TR-4992

Approved for public release; distribution is unlimited.

August 2008-September 2009

Viscoelastic Models for Nearly Incompressible Materials

Mike Scheidler AH80

U.S. Army Research Laboratory

ATTN: RDRL-WMT-D

Aberdeen Proving Ground, MD 21005-5066

This report discusses linear and nonlinear constitutive models for the viscoelastic response of nearly incompressible solids

undergoing moderate to large deformations and arbitrarily large rotations. For both the linear and nonlinear cases, the general

theory is outlined first, then the Prony series approximation to the stress relaxation function is introduced, and this in turn is used

to derive various incremental relations for the stress. These incremental relations can be easily implemented in Lagrangean

finite element codes; they eliminate the costly evaluation of a hereditary integral at every time step. For several simple strain

histories, explicit closed-form solutions for the stress are derived. These solutions are useful for verifying the model

implementation.

nonlinear, viscoelastic, rate-dependence, nearly incompressible, Prony series, incremental relations

UU 62

Mike Scheidler

UNCLASSIFIED UNCLASSIFIED UNCLASSIFIED 410-278-5436

ii
Contents

List of Figures vi

Acknowledgments vii

1. Introduction 1

2. Stress and Deformation Tensors 2

3. A Linear Viscoelastic Model 5

3.1 General Remarks........................................................................... 5

3.2 Constitutive Assumptions ................................................................. 6

3.2.1 The Constitutive Relation for S .................................................. 6

3.2.2 The Constitutive Relation for the Pressure....................................... 7

3.2.3 Jump Discontinuities in Stress and Strain ....................................... 9

3.2.4 Rate Form of the Relation for S.................................................. 10

3.3 Response to Special Strain Histories...................................................... 12

3.3.1 Stress Relaxation Tests ........................................................... 12

3.3.2 Constant Strain-Rate and Strain-Rate Jump Tests ............................... 12

3.4 Prony Series Approximation and Incremental Relations................................. 14

3.4.1 Prony Series Approximation ..................................................... 15

3.4.2 Incremental Relations for Piecewise-Linear Strain Histories ................... 17

3.4.3 The functions f, α, and β ......................................................... 19

iii
3.4.4 Alternative Relations for Piecewise-Linear Strain Histories .................... 20

3.4.5 Special Piecewise-Linear Strain Histories ....................................... 23

3.4.6 Approximate Incremental Relations for General Strain Histories .............. 26

4. A Nonlinear Viscoelastic Model 26

4.1 The General Model ........................................................................ 27

4.1.1 Basic Assumptions ................................................................ 27

4.1.2 Stress Relaxation Tests ........................................................... 30

4.1.3 Rate Form of the Relation for S.................................................. 31

4.2 The Separable Case ........................................................................ 32

4.2.1 Constitutive Assumptions for the Separable Case ............................... 32

4.2.2 Alternative Notation for the Separable Case ..................................... 34

4.2.3 Prony Series and Incremental Relations for the Separable Case................ 35

4.3 Prony Series Approximation for the General Model ..................................... 35

4.3.1 Prony Series Approximation for the Stress Relaxation Function ............... 36

4.3.2 Relations for the Stress ........................................................... 37

4.3.3 Alternative Notation .............................................................. 38

4.3.4 Approximate Incremental Relations ............................................. 40

4.4 Special Cases of the Prony Series Approximation ....................................... 42

4.4.1 A Simple Non-Separable Case ................................................... 42

4.4.2 One-Parameter Family of Functions: Part I..................................... 44

4.4.3 The Mooney-Rivlin Elastic Model ............................................... 46

4.4.4 One-Parameter Family of Functions: Part II .................................... 47

iv
5. References 50

Distribution List 51

v
List of Figures

Figure 1. f(x), α(x), and β(x) for 0 ≤ x ≤ 20, with other functions shown for comparison. . 20

Figure 2. f(x), α(x), and β(x) for 0 ≤ x ≤ 5, with other functions shown for comparison. .. 21

vi
Acknowledgments

I am grateful to Steven Segletes for providing the LATEX style files that were used to convert my

original manuscript to the U.S. Army Research Laboratory report style and for the time he spent

assisting me in this process.

vii
INTENTIONALLY LEFT BLANK.

viii
1. Introduction

This report discusses several constitutive models for the viscoelastic response of nearly

incompressible solids. All of the models presented here are properly invariant, which implies

that they are valid for arbitrarily large material rotations. They should prove useful for moderate

to large deformations, depending on the particular model. For infinitesimal deformations, these

finite deformation models reduce to the classical theory of isotropic, linear viscoelasticity with a

purely elastic pressure-volume response. Thermodynamic effects are neglected here, although

we intend to incorporate temperature dependence in a subsequent report.

The term “nearly incompressible” refers to materials whose response in shear is much softer than

in (purely volumetric) compression. A more precise definition, at least for elastic materials, is

that the bulk modulus is several orders of magnitude larger than the shear modulus. For

viscoelastic materials, the shear modulus is a time-dependent function (the stress relaxation

function), but a similar criterion can be applied: the material is nearly incompressible provided

the bulk modulus is several orders of magnitude larger than the instantaneous elastic shear

modulus (see section 3.2). For such materials, volumetric and shear strains of the same

magnitude will generate a pressure that is several orders of magnitude larger than the shear stress;

conversely, pressure and shear stresses of the same magnitude will generate a volumetric strain

that is several orders of smaller larger than the shear strain. Examples of nearly incompressible

materials include rubber, biological materials with high water content, and tissue simulants such

as ballistic gelatin. For such materials, the pressure depends primarily on the current volumetric

strain; in particular, the viscoelastic contribution to the pressure may be neglected to a first

approximation.

Thus, the models described here focus on the viscoelastic nature of the shear stress. Any such

model requires that this stress depend not only on the current value of the strain but also on the

past values, that is, on the history of the strain up to the current time. One way of doing this is by

means of a hereditary integral, as in the linear theory of viscoelasticity for infinitesimal strains; cf.

Wineman and Rajogopal (1). This is the approach taken here. However, for finite deformations

some care must be taken to ensure a properly invariant formulation. This is most easily done by

working with the material (referential) description of stress and strain. Hence, the models

considered here are expressed in terms of the second Piola-Kirchhoff stress tensor S and the

Green-Lagrange strain tensor E. The deviatoric part of the Cauchy stress tensor, which measures

1
shear stress in the current configuration, may then be obtained from S by standard relations. The

necessary background material from continuum mechanics is summarized in section 2.

In section 3, we consider a linear hereditary integral analogous to that used in infinitesimal strain

theory. The general theory is laid out in sections 3.1–3.2, and the response to some simple strain

histories is discussed in section 3.3. In section 3.4, we consider the Prony series approximation

to the stress relaxation function. This approximation has two main advantages. It allows for

easier calibration of the model to experimental data (a topic not addressed here), and it allows one

to derive incremental relations for the stress. Several forms of these incremental relations are

derived. Some of these are analogous to relations in the literature; others appear to be new.

Their equivalence for piecewise-linear strain histories and for arbitrary strain histories is

discussed. These incremental relations are also used to derive explicit, closed-form solutions for

the stress for several simple strain histories.

The linear model in section 3 cannot be expected to adequately represent the material response for

arbitrarily large deformations. In section 4, we consider a nonlinear generalization, which retains

some of the useful features of the linear model. The general theory is outlined in section 4.1, and

a simple special case is discussed in section 4.2. In section 4.3, we introduce the Prony series

approximation to the general, nonlinear model and derive incremental relations analogous to

those for the linear model. Some simple special cases of this Prony series approximation are

discussed in section 4.4.

The incremental relations derived here (for either the linear or nonlinear models) can be easily

implemented in Lagrangean finite element codes. They eliminate the costly evaluation of a

hereditary integral at every time step. The explicit, closed-form solutions for the stress (for

simple strain histories) are useful for verifying the model implementation. Intended applications

of this work include terminal ballistics simulations involving biological materials or ballistic

gelatin, where the high strain-rates dictate that rate-dependence in the material response be taken

into account.

2. Stress and Deformation Tensors

We follow (for the most part) the notation and terminology in Truesdell and Noll (2), Gurtin (3),

Bowen (4), and Holzapfel (5). Let F denote the deformation gradient, and let J denote the

Jacobian of the deformation:

J = detF . (1)

2
Let R, U, and V denote the rotation tensor and the right and left stretch tensors in the polar

decomposition of F:

F = RU = V R . (2)

Let C denote the right Cauchy-Green deformation tensor, and let E denote the Green–Lagrange

strain tensor:

C := F

TF = U

2 = I + 2E , E :=

1

2

(C − I), (3)

where I denotes the identity tensor. Then

C

−1 = F−1F−T = U−2 = I − 2E + O(E2

), (4)

where F

−T

denotes the inverse transpose of F, and

J =

√

det C = 1 + tr E + O(E

2

). (5)

Let σ denote the Cauchy stress tensor, σRthe rotated (Cauchy) stress tensor, P the 1st

Piola-Kirchhoff stress tensor, and S the 2nd Piola-Kirchhoff stress tensor. They are related as

follows:

σR:= RTσR =

1

J

USU =

1

J

RTP U , (6)

S := JF

−1σF −T = J U−1σRU−1 = F−1P , (7)

P := JσF −T = F S = JRσRU

−1

, (8)

σ =

1

J

F SFT = RσRRT =

1

J

P FT. (9)

The deviatoric part of any 2nd order tensor T is denoted by dev T . In particular, the deviatoric

stress tensor, dev σ, is the deviatoric part of the Cauchy stress tensor and represents a tensorial

measure of shear stress. The pressure p and deviatoric stress tensor are given in terms of the

Cauchy stress tensor by

p := −

1

3

trσ and dev σ := σ −

1

3

tr(σ)I , (10)

which yields the standard decomposition of the Cauchy stress tensor into spherical and deviatoric

parts:

σ = −pI + dev σ . (11)

3
On substituting this into equation 61, we see that the rotated stress tensor σR

can be expressed in

terms of the pressure and deviatoric stress by

σR = −pI + dev σR, dev σR = RT(dev σ)R . (12)

Conversely, we may solve equation 12 for the pressure and deviatoric stress in terms of σR:

p = −

1

3

trσR, dev σ = R(dev σR)RT. (13)

Similarly, on substituting equation 11 into equation 71, we see that the 2nd Piola-Kirchhoff stress

tensor S can be expressed in terms of the pressure and deviatoric stress by

S = −Jp C

−1 + S∗

, S

∗

:= JF

−1

(dev σ)F

−T

. (14)

Conversely, on substituting equation 91 into equation 101 we see that the pressure can be

expressed in terms of S by

p = −

1

3J

tr ¡F SFT

¢

= −

1

3J

tr ¡USU¢= −

1

3J

tr ¡SC¢. (15)

And from equation 14, we see that the deviatoric stress can be expressed in terms of S

∗

by

dev σ =

1

J

F S∗F

T =

1

J

R(US∗U)RT. (16)

It follows that F S∗F

T

and US∗U are deviatoric tensors; this property depends on the fact that

the scalar p in the decomposition 141 of S is the pressure, i.e., that p = −

1

3

trσ.

Finally, we note that there are relations involving the 1st Piola-Kirchhoff stress tensor P

analogous to those for S above, but we do not list them here.

4
3. A Linear Viscoelastic Model

3.1 General Remarks

We consider a simple, isotropic, viscoelastic constitutive relation expressed in terms of the 2nd

Piola-Kirchhoff stress tensor S and the Green-Lagrange strain tensor E; thermodynamic effects

are neglected here. The advantage of working with these stress and strain tensors is that any

constitutive relation that expresses S(t) as a function of the history of E (or C or U) up to time t

is necessarily properly invariant; cf. (2).

We decompose S as

S = −Jp ̄C

−1 + S , (17)

and assume that S is governed by the isotropic, linear viscoelastic constitutive relation 25 below.

The decomposition 17 is formally similar to equation 141 but is not equivalent to it. Indeed, the

constitutive relation for S is such that FSFTand USU will generally not be deviatoric tensors,

so by the remark following equation 16 it follows that the scalar p ̄ in equation 17 will be related

to, but not necessarily equal to, the pressure p; cf. equations 23–24 below. The constitutive

relation for p ̄ will be chosen in such a way that the pressure p depends (either exactly or

approximately) on the current volumetric strain only.

The other stress tensors are easily computed in terms of p ̄ and S. On substituting equation 17

into equation 82, we see that

P = −Jp ̄F

−T + FS . (18)

Similarly, on substituting equation 17 into equations 91 and 62, we see that

σ = −p ̄I + σ , σ :=

1

J

FSFT

, (19)

σR = −p ̄I + σR

, σR:=

1

J

USU , (20)

where neither σ nor σRare deviatoric, in general.

5
The above relations imply

dev σ = dev σ =

1

J

dev ¡FSFT

¢

, (21)

dev σR = dev σR =

1

J

dev ¡USU¢. (22)

On substituting equations 19 and 20 into equations 101 and 131, respectively, or on substituting

equation 17 into equation 15, we see that the pressure p is given in terms of p ̄ and S by

p =  ̄p −

1

3

trσ , (23)

where

trσ = trσR =

1

J

tr ¡FSFT

¢

=

1

J

tr ¡USU¢=

1

J

tr ¡SC¢. (24)

Equations 21–24, together with the spherical-deviatoric decompositions 11 and 121, give

alternative expressions for the Cauchy stress tensor σ and its rotated counterpart σR.

3.2 Constitutive Assumptions

3.2.1 The Constitutive Relation for S

We assume that the value of S at time t is given in terms of the history of the strain tensor E up to

time t by the linear single integral law

S(t) = 2G(0)E(t) + Z t

−∞

2

.

G(t − s)E(s) ds

= 2G(0)E(t) + Z ∞

0

2

.

G(s)E(t − s) ds .

(25)

Here G : [0, ∞) → (0, ∞) is the stress relaxation function in shear, and

.

G denotes the derivative

of G:

.

G(s) := d

ds G(s), (26)

so that

.

G(t − s) = ∂

∂t G(t − s) = −

∂

∂s G(t − s). (27)

6
The stress relaxation function G is assumed to be smooth, convex, and strictly decreasing with a

positive limit at infinity:

G(0) > G(s) > G(∞) := lims→∞G(s) > 0 , ∀ s > 0 . (28)

G(0) and G(∞) are the instantaneous and equilibrium elastic shear shear moduli, respectively.

.

G is negative and strictly increasing and asymptotes to 0:1

.

G(0) <

.

G(s) <

.

G(∞) = 0 , ∀ s > 0 . (29)

For the remainder of section 3, we assume that the material is undeformed for times t < 0, so that

F(t) = R(t) = U(t) = C(t) = I and E(t) = 0 for t < 0. (30)

Then S(t) is zero for t < 0, and the constitutive relation 25 for S reduces to

S(t) = 2G(0)E(t) + Z t

0

2

.

G(t − s)E(s) ds

= 2G(0)E(t) + Z t

0

2

.

G(s)E(t − s) ds .

(31)

3.2.2 The Constitutive Relation for the Pressure

We are interested in materials for which the pressure p depends primarily on the current

volumetric strain. To impose this condition exactly, we would assume a (possibly nonlinear)

elastic constitutive relation for p of the form

p = p(J), p(1) = 0 , (32)

where the condition on the right ensures that the pressure is zero in the undeformed state, where

J = 1. Then

p = κ0

h

(1 − J) + O

¡

(1 − J)

2

¢

i

, (33)

1Other desirable properties of G, such as positive relaxation spectrum, are not discussed here since we eventually

assume that G is given by the Prony series in equation 64.

7
where the (initial) bulk modulus κ0 is given by

κ0 = −

dp

dJ

 ̄

 ̄

 ̄

 ̄

J=1

= −p

0

(1). (34)

Given equation 32, it follows from equations 23–24 that the scalar p ̄ in the decomposition 17 of S

must satisfy

p ̄ = p(J) + 1

3

trσ = p(J) + 1

3J

tr ¡SC¢. (35)

In particular, the expression on the right introduces a viscoelastic contribution to p ̄. With S

determined from equation 31 and p ̄ from equation 35, the stress tensors P , σ, and σRcan be

determined from equations 18–20.

If one is ultimately only interested in computing σ or σR, then there would be no need to actually

compute p ̄: simply use the spherical-deviatoric decomposition 11 or 121, with p determined from

equation 32, and dev σ or dev σRdetermined from S via equation 21 or 22.

2

Next, suppose we replace the constitutive assumption 32 for p by an analogous assumption for p ̄:

p ̄ = p(J), p(1) = 0 , (36)

so that

p ̄ = κ0

h

(1 − J) + O

¡

(1 − J)

2

¢

i

, (37)

with κ0 again given by −p

0

(1). Then by equations 23–24, the pressure p satisfies

p = p(J) −

1

3J

tr ¡SC¢. (38)

It is clear that in this case there will be a viscoelastic contribution to the pressure (via S). While

such an effect is not unreasonable for polymers, there is no reason to expect that it would be

determined by the stress relaxation function G in shear that characterizes S.

3 On the other hand,

for nearly incompressible materials, that is, materials for which the bulk modulus κ0 is several

orders of magnitude larger than the instantaneous elastic shear modulus G(0), the trace term in

equation 38 should be relatively small in many cases of interest, so that p ≈ p(J), i.e., we recover

2Of course, it is clear from equations 21 and 22 that the calculation of these deviatoric tensors from S involves the

same trace calculation as in equation 35, so this approach to computing σ or σR does not involve any less work than

the one in the preceding paragraph.

3One reason for considering equation 36 (in conjunction with equations 17 and 31) is that it is easy to construct

a viscoelastic strain energy function for the 2nd Piola-Kirchhoff stress tensor S in this case. It is not clear that this

conclusion necessarily holds when p ̄ is given by equation 35, i.e., when equation 32 holds.

8
equations 32–34 as an approximation. In this case p ≈ p ̄, so by equations 141 and 17, we see that

S ≈ S

∗

also.

The constitutive relation 17, with S given by equation 31 and p ̄ given by either equation 35 or

equation 36, is isotropic and properly invariant, with the reference configuration a natural state for

the material; it yields the elastic pressure-volume relation p = p(J), either exactly or

approximately. For terminal ballistics applications with nearly incompressible materials like

ballistic gelatin, for which the bulk modulus is four orders of magnitude larger than the shear

modulus, this modeling framework should be adequate for arbitrarily large rotations and for

moderate strains. Whether the linear constitutive relation 31 for S suffices to characterize the

material response for “large” shear strains would, of course, depend on the particular material.

Some nonlinear constitutive relations are discussed in section 4.

Finally, we note that for infinitesimal deformations, the finite deformation theory above reduces to

the classical theory of isotropic linear viscoelasticity with a purely elastic pressure-volume

response.

3.2.3 Jump Discontinuities in Stress and Strain

We assume that the strain history t 7→ E(t) is piecewise continuous. Then the (equivalent)

hereditary integrals in equation 31 are well-defined and are continuous functions of the time t.

The one-sided limits of E at time t are denoted by

E(t

+

) := lim

τ↓t

E(τ ) and E(t

−

) := lim

τ ↑t

E(τ ), (39)

and the jump in E at time t is defined by

[[E]](t) := E(t

+) − E(t−) ; (40)

similar conventions hold for the other variables. Of course, if E is continuous at the instant t

then E(t

+) = E(t−) = E(t) and [[E]](t) = 0. But if E suffers a jump discontinuity at time t,

then E(t) need not be defined.

If E suffers a jump discontinuity at time t, then by taking the one-sided limits of equation 31 and

using the continuity of the integrals, we see that equation 31 holds with the S(t) and E(t) terms

replaced by either of their one-sided limits. On taking the difference of these limiting relations,

9
we see that the integrals cancel, so that the jump in S at time t is given by

[[S]](t) = 2G(0) [[E]](t). (41)

In particular, past jump discontinuities in strain do not affect the current value of the stress.

The condition 30 implies that

E(0−) = 0 and S(0−) = 0 . (42)

We assume that E(0) equals its limiting value as t approaches 0 from above. Then by

equations 39–42 it follows that

E(0) = E(0+) = [[E]](0), (43)

S(0) = S(0+) = [[S]](0) = 2G(0)E(0). (44)

On taking the rate of equation 311, one finds that if E is continuous at the instant t but its rate

(i.e., material time derivative)

.

E suffers a jump discontinuity, then S is also continuous at time t,

but its rate suffers a jump discontinuity given by

[[

.

S]](t) = 2G(0) [[

.

E]](t). (45)

In view of equations 41 and 45, we see that jump discontinuities in S and its rate are governed by

the instantaneous elastic shear shear modulus G(0).

3.2.4 Rate Form of the Relation for S

To simplify things further, for the remainder of section 3 we will assume that the deformation

gradient is continuous, except possibly for a jump discontinuity at t = 0; then the same holds for

the strain tensor E. This assumption does not lead to any simplifications of the viscoelastic

constitutive relation 31, but it is used in deriving the alternative relation 47 below in terms of the

strain-rate, and it simplifies the incremental forms of the constitutive relation for S discussed in

10
section 3.4. Summarizing the assumptions on the strain history, we have

Continuity Assumptions on the Strain History:

E(t) = 0 for t < 0 ;

E(t) is a continuous function of t for t ≥ 0 ;

E(t) may have a jump discontinuity at t = 0 ;

.

E(t) is a piecewise continuous function of t .







(46)

On using the identity 27 in equation 31 and integrating by parts, we obtain a constitutive relation

for S in terms of the history of the strain-rate

.

E:

S(t) = 2G(t)E(0) + Z t

0

2G(t − s)

.

E(s) ds

= 2G(t)E(0) + Z t

0

2G(s)

.

E(t − s) ds .

(47)

For continuous strain histories, we must have E(0) = 0, in which case the first group of terms on

right in equation 47 drops out:

S(t) = Z t

0

2G(t − s)

.

E(s) ds =

Z t

0

2G(s)

.

E(t − s) ds . (48)

This should be appropriate for implementation in hydrocodes, where the use of artificial viscosity

yields continuous strain histories. However, for the discussion of certain theoretical results, such

as stress relaxation tests or shock waves treated as propagating jump discontinuities, we need to

retain the possibility of an initial jump in strain.4

4

In view of the assumptions 46 on the strain history, which are essential for the validity of equation 47, it follows

that this equation would only be valid for a shock wave arriving at time 0. Note that no such restriction applies to the

original relation 31. If strain jumps at times greater than zero are also allowed, then the relation 47 is no longer valid:

one must add terms involving all strain jumps up to time t; cf. Scheidler (6) for details.

11
3.3 Response to Special Strain Histories

3.3.1 Stress Relaxation Tests

Consider an ideal stress relaxation test:

F(t) =







I, t < 0 ,

F0, t ≥ 0 ,

(49)

where F0 is a constant tensor (with positive determinant J0); similarly for C, R, U, and V .

Then

E(t) =







0, t < 0 ,

E0, t ≥ 0;

(50)

so by equations 47 and 32 it follows that in a stress relaxation test

S(t) =







0, t < 0 ,

2G(t)E0, t ≥ 0;

(51)

lim

t→∞

S(t) = 2G(∞)E0 ; (52)

and

p(t) =







0, t < 0;

p(J0), t ≥ 0.

(53)

3.3.2 Constant Strain-Rate and Strain-Rate Jump Tests

Consider a continuous strain history which, for some time t∗ > 0, satisfies

.

E(t) =







.

E1, 0 < t < t∗ ,

.

E2, t > t∗ ,

(54)

where

.

E1 and

.

E2 are constant symmetric tensors. If

.

E2 6=

.

E1 then there is a strain-rate jump at

time t∗. The strain history corresponding to equation 54 is the continuous, piecewise-linear

function:

E(t) =







t

.

E1, 0 ≤ t ≤ t∗;

t∗

.

E1 + (t − t∗)

.

E2, t ≥ t∗.

(55)

12
When

.

E1 and

.

E2 are coaxial, the principal stretches (i.e., eigenvalues of U) are easily computed

from equation 55, but they are not piecewise-linear.

By equations 54 and 48, we have

S(t) =







·Z t

0

2G(s) ds ̧

.

E1, 0 ≤ t ≤ t∗;

·Z t−t∗

0

2G(s) ds ̧

.

E2 +

·Z t

t−t∗

2G(s) ds ̧

.

E1, t ≥ t∗ .

(56)

Hence the stress-rate is given by

.

S(t) =







2G(t)

.

E1, 0 < t < t∗;

2G(t − t∗)

.

E2 + 2[G(t) − G(t − t∗)]

.

E1, t > t∗ .

(57)

Since G(t) asymptotes to G(∞) as t → ∞,

.

S(t) → 2G(∞)

.

E2 , as t → ∞ . (58)

Note that for times t such that 0 < t < t∗ (i.e., the top lines in equations 54–57), we have a

constant strain-rate test. Of course, we have a constant strain-rate test for all times if

.

E2 =

.

E1.

For the special case where

.

E2 = 0, equation 55 reduces to a linear ramp in strain over the time

interval [0, t∗] to a fixed strain E∗ := E(t∗) = t∗

.

E1:

E(t) =







t

.

E1, 0 ≤ t ≤ t∗;

t∗

.

E1 = E∗, t ≥ t∗.

(59)

Since

.

E2 = 0, the relation 562 for S(t) for times t ≥ t∗ reduces to

S(t) = ·Z t

t−t∗

2G(s) ds ̧

.

E1 =

1

t∗

·Z t

t−t∗

2G(s) ds ̧E∗ , t ≥ t∗ . (60)

For real stress relaxation tests conducted in a laboratory (as opposed to the ideal case considered

in section 3.3.1), the strain ramps up continuously from 0 to its fixed value E∗ over some short

time interval [0, t∗]. The strain history 59 and the relation 60 for S(t) for times t ≥ t∗ provide

better approximations to such tests. Note that the coefficient of E∗ in the relation on the right in

13
equation 60 is the mean value of G over the time interval [t − t∗, t]; it approaches G(t) as t∗ → 0,

that is, as the ramp time approaches zero. Hence for fixed t,

S(t) → 2G(t)E∗ , as t∗ → 0 , (61)

which is consistent with the relation 512 for an ideal stress relaxation test with E0 = E∗.

3.4 Prony Series Approximation and Incremental Relations

In section 3.4.1, we consider the Prony series approximation to the stress relaxation function G ;

cf. equation 64. When this is substituted into the hereditary integrals, S(t) reduces to a sum of

certain simpler hereditary integrals. This sum can be written in several equivalent ways, three of

which are considered here. The corresponding simpler hereditary integrals are denoted by

An(t), An(t), and A

n

(t); cf. equations 68–75.

The strain tensor E(t) can be computed by integrating the relation

.

E = F

TDF , D = sym(∇v), (62)

where ∇v is the spatial velocity gradient. Of course, E(t) can also be computed directly from

F(t) using equation 3, with F obtained by integrating the relation for the material time derivative

of the deformation gradient:

.

F = (∇v)F . (63)

For the present discussion, the strain history E is regarded as given or computed at the times of

interest, and we focus attention on evaluating the integrals An(t), An(t), and A

n

(t). Recall that

the strain history satisfies the continuity assumptions 46. Sections 3.4.2 and 3.4.4 present exact

incremental relations for An, An, and A

n

for the case where the strain history is piecewise-linear.

In section 3.4.5, the results of section 3.4.4 are applied to some special piecewise-linear strain

histories. In section 3.4.6, we briefly discuss the corresponding approximate incremental

relations for arbitrary strain histories. When these incremental relations are substituted into the

corresponding expressions for S(t) (cf. equations 68, 71, or 74), we obtain incremental relations

for the stress. Section 3.4.3 discusses some special functions that appear in the incremental

relations.

14
3.4.1 Prony Series Approximation

Henceforth, we assume that the stress relaxation function in shear is given (or approximated) by

the Prony series

G(t) = μ∞ +

X

N

n=1

μne

−t/τn

, (64)

where the moduli μ∞ and μn and the relaxation times τn are positive. Then

G(0) = μ∞ +

X

N

n=1

μn =: μ0 , G(∞) = μ∞ , (65)

and

.

G(t) = −

X

N

n=1

μn

τn

e

−t/τn

. (66)

Since G(0) = μ0, by equation 44 we see that the value of S at time zero is given by

S(0) = 2μ0E(0). (67)

On substituting equation 66 (with t → t − s) into the constitutive relation 311 for S(t) and

replacing G(0) with μ0, we obtain the following relation for S(t) at any time t ≥ 0:

S(t) = 2μ0E(t) −

X

N

n=1

2μnAn(t), (68)

where the dimensionless, symmetric tensor An(t) is defined by

An(t) := Z t

0

1

τn

e

−(t−s)/τnE(s) ds , (69)

so that

An(0) = 0 . (70)

Note that on setting t = 0 in equation 68 and using equation 70, we recover equation 67.

Similarly, on substituting equation 64 into the constitutive relation 471 for S(t), we obtain

S(t) = 2μ∞E(t) +X

N

n=1

2μnAn(t) +X

N

n=1

2μne

−t/τnE(0), (71)

15
where the dimensionless, symmetric tensor An

(t) is defined by

An(t) := Z t

0

e

−(t−s)/τn

.

E(s) ds , (72)

so that

An(0) = 0 . (73)

Note that on setting t = 0 in equation 71 and and using equations 73 and 65, we recover

equation 67.

We can re-write equation 71 as

S(t) = 2μ∞E(t) +X

N

n=1

2μnA

n

(t), (74)

where the dimensionless, symmetric tensor A

n

(t) is defined by

A

n

(t) := An(t) + e

−t/τnE(0). (75)

Similarly, on substituting the expression for μ0 in equation 65 into equation 68 and re-arranging,

we obtain equation 74 with

A

n

(t) = E(t) − An(t). (76)

It is also easily verified directly that the relations 76 and 75 for A

n

(t) are equivalent; thus

An(t) + e

−t/τnE(0) = E(t) − An

(t). (77)

This may be used to express An(t) in terms of An(t) and visa versa; then equation 68 can be

obtained from equation 71 and visa versa.

Given that G is represented by the Prony series 64 and given the continuity assumptions 46 on the

strain history, the results above are exact, and An(t), An(t), and A

n

(t) are continuous functions

of t for t ≥ 0. The only constitutive property that they depend on is the nth relaxation time τn.

Note that

A

n

(0) = E(0) ; (78)

so unlike Anand An, the tensor A

n

(0) is not 0 unless E(0) = 0. The condition E(0) = 0 holds

when the strain history is continuous for all times, and in this case A

n

(t) = An(t) by equation 75,

the sum on the far right in equation 71 drops out, and the relations 71 and 74 are indistinguishable.

16
In the constitutive relations 68, 71, and 74 for S(t), the tensors An

, An, and A

n

can be regarded

as internal state variables that modify the elastic part of the stress (either 2μ0E(t) or 2μ∞E(t)).

One can easily derive evolution equations for Anand Anby differentiating the integral relations

69 and 72, then an evolution equation for A

n

can be obtained from equations 75 or 76. However,

for purposes of computing these internal state variables, it is simpler to work directly with the

integral relations, as discussed in section 3.4.2.

3.4.2 Incremental Relations for Piecewise-Linear Strain Histories

Here and in section 3.4.4, we assume that in addition to the conditions 46, E(t) is also a

piecewise-linear function of time t; equivalently,

.

E is a piecewise constant function. Thus, there

is some sequence of times tk (k = 0, 1, 2, . . .) satisfying

0 = t0 < t1 < t2 < . . . < tk−1 < tk < . . . (79)

and such that E is linear (more precisely, affine) over each time interval [tk−1, tk]. We set

∆tk = tk − tk−1 . (80)

Note that while any jump discontinuities in the strain-rate

.

E must occur at one of the discrete

times tk,

.

E need not suffer a jump discontinuity at every tk. In other words, E may be linear

(and hence,

.

E may be constant) over time intervals spanning more than one of the subintervals

[tk−1, tk]. For example, for the strain-rate jump test described in section 3.3.2, we need only

assume that tj = t∗ for some j ≥ 1 for the results of this section and section 3.4.4 to hold

exactly.5 On the other hand, for numerical simulations of this test (or of more general

piecewise-linear strain histories) with the tk such that ∆tk is the computational time step, the size

of ∆tk would be limited by a stability condition, so that many subintervals [tk−1, tk] might be

required within a larger time interval over which

.

E is constant.

We claim that An, An, and A

n

are given exactly at the times in equation 79 by the following

incremental relations (here k = 1, 2, . . ., and the functions α, β, and f below are defined in

5We will return to this special case in section 3.4.5.

17
section 3.4.3):

An(tk) = e

−∆tk/τnAn

(tk−1)

+ α(∆tk/τn)E(tk) + β(∆tk/τn)E(tk−1), (81)

An(t0) = An(0) = 0 ; (81a)

and

An(tk) = e

−∆tk/τnAn

(tk−1) + f(∆tk/τn)

£

E(tk) − E(tk−1)

¤

, (82)

An(t0) = An(0) = 0 ; (82a)

and

A

n

(tk) = e

−∆tk/τn An

(tk−1) + f(∆tk/τn)

£

E(tk) − E(tk−1)

¤

, (83)

A

n

(t0) = A

n

(0) = E(0). (83a)

Note that the incremental relations for Anand A

n have the same form but different initial

conditions.

The relations 81 and 82 are derived by the usual procedure6 of setting t = tk in equations 69 or 72

and breaking up the hereditary integral into two parts,

Z tk

0

. . . ds =

Z tk−1

0

. . . ds +

Z tk

tk−1

. . . ds . (84)

The details are straightforward and not included here.7 The relation 83 follows from

equations 82 and 75, or from equations 81 and 76. In fact, on using equations 75–77, we see that

any one of the incremental relations 81–83 implies the other two. Now set t = tk in the Prony

6Cf. Zocher et al. (7) and Clements (8). Both of these references consider the more general case where the

hereditary integral involves a reduced time variable.

7But the following observations are of interest. The piecewise-linearity assumption is not needed for the evaluation of the integral over [0, tk−1]; this integral yields first term on the right in equations 81 or 82 without any

assumptions on E or

.

E. However, the the linearity of E over [tk−1, tk] is crucial for the evaluation of the second

integral on the right, which yields the bottom right group of terms in equation 81 and the second group of terms on the

right in equation 82.

18
series relations 68, 71, and 74 for S(t). Then for k = 0, 1, 2, . . . we have

S(tk) = 2μ0E(tk) −

X

N

n=1

2μnAn(tk), (85)

S(tk) = 2μ∞E(tk) +X

N

n=1

2μnAn(tk) + "X

N

n=1

2μne

−tk/τn

#

E(0), (86)

S(tk) = 2μ∞E(tk) +X

N

n=1

2μnA

n

(tk). (87)

These relations, together with the incremental relations 81–83 for An(tk), An(tk), and A

n

(tk),

yield exact (and hence equivalent) relations for S(tk) for piecewise-linear strain histories.

Finally, we emphasize that that the continuity assumptions on the strain history (cf. equation 46)

have been used in the deriving the incremental relations above and in sections 3.4.4–3.4.6 below.

Appropriate modifications when these conditions are relaxed are discussed in Scheidler (6).

3.4.3 The functions f, α, and β

The function f in the incremental relations 82 and 83 for Anand A

n

is defined by

f(x) =







1 − e

−x

x

, x 6= 0 ;

1, x = 0 .

(88)

The functions α and β in the incremental relation 81 for Anare defined by

α(x) = 1 − f(x) and β(x) = f(x) − e

−x

. (89)

Then

α(x) =







1 +

e

−x − 1

x

=

e

−x + x − 1

x

, x 6= 0 ,

0 , x = 0 ;

(90)

and

β(x) =







1 − e

−x

x

− e

−x =

1 − (x + 1)e

−x

x

, x 6= 0 ,

0 , x = 0 .

(91)

See figures 1 and 2. f, α, and β are analytic functions on the real line, but we only need f(x),

α(x), and β(x) for x > 0, and these all lie between 0 and 1. f is strictly decreasing and

19
0 5 10 15 20

0.2

0.4

0.6

0.8

1.0

1

1-ã-x

ã

-x

ΒHxL

ΑHxL

fHxL

Figure 1. f(x), α(x), and β(x) for 0 ≤ x ≤ 20, with other functions shown for

comparison.

asymptotes to 0. α is strictly increasing and asymptotes to 1. β(x) increases from 0 at x = 0 to

an absolute maximum of 0.2984 . . . at x = 1.793 . . ., and then decreases asymptotically to 0.

Indeed,

β(x) ∼ f(x) ∼

1

x

, x → ∞ , (92)

and

α(x) ∼ 1 −

1

x

, x → ∞ . (93)

For small x, we have the approximations

f(x) = 1 −

1

2

x +

1

6

x

2 − · · · , (94a)

α(x) = 1

2

x −

1

6

x

2 + · · · , (94b)

β(x) = 1

2

x −

1

3

x

2 + · · · . (94c)

3.4.4 Alternative Relations for Piecewise-Linear Strain Histories

In this and the next paragraph, we assume that k = 1, 2, . . .. Given the assumptions in

section 3.4.2,

.

E is constant over any time interval (tk−1, tk), and we denoted this constant value

20
0 1 2 3 4 5

0.2

0.4

0.6

0.8

1.0

1x

1

1-1x

1-ã-x

ã

-x

ΒHxL

ΑHxL

fHxL

Figure 2. f(x), α(x), and β(x) for 0 ≤ x ≤ 5, with other functions shown for

comparison.

by

.

Ek. However,

.

E may suffer a jump discontinuity at any of the times tk, so

.

E(tk) and

.

E(tk−1) need not be defined, although the one-sided limits of

.

E(t) as t → tk from below or as

t → tk−1 from above necessarily exist and equal

.

Ek. Clearly, then, we have the following

expressions for

.

Ek:

.

E(t), t ∈ (tk−1, tk) =:

.

Ek =







.

E(t

−

k

) := lim

t ↑ tk

.

E(t),

.

E(t

+

k−1

) := lim

t ↓ tk−1

.

E(t),

.

E(t

−

k

) +

.

E(t

+

k−1

)

2

,

.

E

μ

tk + tk−1

2

¶

,

E(tk) − E(tk−1)

tk − tk−1

=

E(tk) − E(tk−1)

∆tk

.

(95)

21
The last relation above yields

E(tk) − E(tk−1) = ∆tk

.

Ek . (96)

When this is substituted into the incremental relation 82 for Anor the incremental relation 83 for

A

n

, we obtain alternative forms of these relations. In particular, equation 82 is equivalent to the

incremental relation

An(tk) = e

−∆tk/τnAn

(tk−1) + τn

¡

1 − e

−∆tk/τn

¢ .

Ek , (97)

An(t0) = An(0) = 0 , (97a)

where

.

Ek is given by any of the expressions in equation 95.

8

When equation 97 is substituted into equation 86, we obtain another (equivalent) incremental

relation for S(tk) in terms of E(tk) and An(tk). This relation may also be written (for

k = 0, 1, 2, . . .) as

S(tk) = 2μ∞E(tk) +X

N

n=1

Sn(tk) + "X

N

n=1

2μne

−tk/τn

#

E(0), (98)

where for k = 1, 2, . . . we have

Sn(tk) = e

−∆tk/τn Sn(tk−1) + 2μnτn

¡

1 − e

−∆tk/τn

¢ .

Ek , (99)

and

Sn(t0) = Sn(0) = 0 . (99a)

Since t0 = 0, on setting k = 0 in equation 98 and using equations 99a and 65, we recover

equation 67. When the strain is continuous for all times, so that E(0) = 0, the incremental

relation 98 reduces (for k = 0, 1, 2, . . .) to

S(tk) = 2μ∞E(tk) +X

N

n=1

Sn(tk), if E(0) = 0. (100)

8Of course, equation 97 with

.

Ek given by the bottom expression in equation 95 just gives back the relation 82.

22
3.4.5 Special Piecewise-Linear Strain Histories

Here we apply the results at the end of the preceding section to some special piecewise-linear

strain histories.

Constant Strain-Rate: We assume that the strain-rate has a constant value

.

E1 for t > 0, but

allow the possibility of a jump in strain at time t = 0, so that E(0) is not necessarily 0. Thus

.

E(t) =

.

E1 , t > 0 , (101)

and the strain is a linear function of time:

E(t) = E(0) + t

.

E1 , t ≥ 0 . (102)

Since E is linear over the time interval [0, t1] for any t1 > 0, and since ∆t1 = t1 − t0 = t1, on

setting k = 1 in equation 99 and using equation 99a, we obtain

Sn(t1) = 2μnτn

¡

1 − e

−t1/τn

¢ .

E1 . (103)

Then on setting k = 1 in equation 98 and using equations 103 and 102, and then replacing t1 with

t, we see that for any t > 0,

S(t) = "2μ∞t +

X

N

n=1

2μnτn

¡

1 − e

−t/τn

¢

#

.

E1 +

"

2μ∞ +

X

N

n=1

2μne

−t/τn

#

E(0). (104)

This also holds for t = 0, since it reduces to equation 67.

For purposes of checking the implementation of the incremental algorithms in codes, it is also

useful to have explicit relations for An(t) and An(t) in this special case. Proceeding as above

and using equations 81–82 and the the relations for f, α, and β in section 3.4.3, we find that

An(t) = α(t/τn)E(t) + β(t/τn)E(0)

= τn(e

−t/τn + t/τn − 1)

.

E1 + (1 − e

−t/τn

)E(0), (105)

and

An(t) = τn(1 − e

−t/τn

)

.

E1 . (106)

23
Of course, these relations can also be obtained directly from the definitions 69 and 72. Also, note

that on substituting equation 105 into equation 68 and using equations 102 and 65, we recover the

relation 104 for S(t). Similarly, we can recover this relation by substituting equation 106 into

equation 71 and using equation 102.

For the case where there is no initial jump in strain we have E(0) = 0, so the second group of

terms on the right in equations 104 and 105 drop out. The solution in this case can also be

obtained by substituting the Prony series 64 into the corresponding solution 561 for a general

stress relaxation function.

For a stress relaxation test we have E(0) = E0 6= 0 and

.

E1 = 0, so the first group of terms on

the right in equations 104 and 1052 drop out, and An(t) = 0. The solution in this case can also

be obtained by substituting the Prony series 64 into the corresponding solution 51 for a general

stress relaxation function.

Bilinear Strain (a single strain-rate jump): Here we assume that there is no initial jump in

strain, so that E(0) = 0, and that the strain-rate satisfies

.

E(t) =







.

E1, 0 < t < t∗ ,

.

E2, t > t∗ ,

(107)

so that there is a strain-rate jump at time t∗ if

.

E2 6=

.

E1. This is the case considered in

section 3.3.2 (for a general stress relaxation function G). Thus the piecewise-linear strain history

is given by equation 55.

For 0 < t < t∗, the strain-rate is constant and the results above apply. Thus for 0 ≤ t ≤ t∗, S(t)

is given by equation 104 with E(0) = 0. It remains to determine S(t) for t > t∗. We set

t1 = t∗. Then E is linear over the time interval [t1, t2] for any t2 > t∗, and ∆t2 = t2 − t∗. On

setting k = 2 in equations 100 and 99, and then replacing t2 with t, we see that for any t > t∗,

S(t) = 2μ∞E(t) +X

N

n=1

h

e

−(t−t∗)/τn Sn(t∗) + 2μnτn

¡

1 − e

−(t−t∗)/τn

¢ .

E2

i

, (108)

where Sn(t∗) is given by equation 103 with t1 = t∗. On substituting that relation into

24
equation 108 and using the relation 552 for E(t), we find that for any t > t∗,

S(t) = "2μ∞t∗ +

X

N

n=1

2μnτn

¡

1 − e

−t∗/τn

¢

e

−(t−t∗)/τn

#

.

E1

+

"

2μ∞(t − t∗) +X

N

n=1

2μnτn

¡

1 − e

−(t−t∗)/τn

¢

#

.

E2 .

(109)

This solution can also be obtained by substituting the Prony series 64 into the corresponding

solution 562 for a general stress relaxation function.

Similarly, we obtain the following explicit relations for An(t) and An(t) for t > t∗:

An(t) =£t∗ + τn

¡

e

−t∗/τn − 1

¢

e

−(t−t∗)/τn

¤ .

E1

+

£

(t − t∗) + τn

¡

e

−(t−t∗)/τn − 1

¢¤ .

E2 ,

(110)

and

An(t) = τn

¡

1 − e

−t∗/τn

¢

e

−(t−t∗)/τn

.

E1 + τn

¡

1 − e

−(t−t∗)/τn

¢ .

E2 . (111)

Note that on substituting equation 110 into equation 68 and using equations 552 and 65, we

recover the relation 109 for S(t). Similarly, we recover this relation by substituting equation 111

into equation 71 and using E(0) = 0 and equation 552.

For the special case where

.

E2 = 0, we have a linear ramp in strain over the time interval [0, t∗] to

a fixed strain E∗ := E(t∗) = t∗

.

E1. Then E(t) reduces to equation 59, and for t > t∗, S(t) is

given by the top line in equation 109, which may also be written as

S(t) = "2μ∞ +

1

t∗

X

N

n=1

2μnτn

¡

1 − e

−t∗/τn

¢

e

−(t−t∗)/τn

#

E∗ . (112)

These formulas can also be obtained by substituting the Prony series 64 into the corresponding

solution 60 for a general stress relaxation function.

Another special case where equation 109 simplifies occurs when

.

E2 and

.

E1 differ only in sign:

.

E2 = −

.

E1 . (113)

25
Then for t > t∗, equation 109 reduces to

S(t) = "2μ∞(2t∗ − t) +X

N

n=1

2μnτn

¡

2e

−(t−t∗)/τn − e−t/τn − 1

¢

#

.

E1 . (114)

3.4.6 Approximate Incremental Relations for General Strain Histories

The results in sections 3.4.2 and 3.4.4 are exact for piecewise-linear strain histories, and there is

no requirement that ∆tk be small. Now, we drop the assumption that the strain history is

piecewise-linear and retain only the continuity assumptions 46. The above results then give

approximate incremental relations for An, An, and A

n

, and consequently for S, provided ∆tk is

small enough that E is approximately linear over the time interval [tk−1, tk]. Here, we briefly

discuss these approximations (so keep in mind that the “=” is now replaced by “≈”, although we

won’t bother to re-write the equations).

The incremental relations 81–83 for An, An, and A

n

are equivalent approximations. That is, any

one of these approximations, together with the relations 75–77 and equation 89, implies the other

two. Consequently, they yield equivalent approximate incremental relations for S.

However, since

.

E is not necessarily constant over the time interval [tk−1, tk], the expressions on

the right-hand side of equation 95 generally yield different approximations for

.

E(t) on this time

interval. So equation 97, with

.

Ek given by any one of these expressions, generally yields a

different approximation for An(tk), and hence, for Sn(tk).

9 And it is only the bottom expression

for

.

Ek in equation 95 which yields an approximation incremental relation that is equivalent to

those in section 3.4.2.

4. A Nonlinear Viscoelastic Model

We continue to assume a decomposition of the 2nd Piola-Kirchhoff stress tensor S of the

form 17, with constitutive relations for p ̄ and the pressure p as discussed in section 3.2.2. But

now we consider a nonlinear generalization of the linear viscoelastic constitutive relation 25 for

S. As before, the Cauchy stress tensor σ, the rotated Cauchy stress tensor σR, and their

deviatoric parts can then be obtained from equations 19– 22.

9The relations 99–100, with the top expression for

.

Ek in equation 95, namely

.

Ek =

.

E(t

−

k

), are analogous to the

relations in Brad Clements’ notes (8) for the special case where the reduced time there is the ordinary time.

26
The nonlinear constitutive relation for S is introduced in section 4.1.1, assuming only that the

strain history is piecewise continuous. Stress relaxation tests are discussed in section 4.1.2. An

alternative relation for the stress in terms of the history of

.

E, which is only valid when the strain

history satisfies the continuity conditions in equation 46, is discussed in section 4.1.3. The

special case of separable strain and time dependence is considered in section 4.2.1, along with a

brief discussion of the Prony series approximation and incremental relations for this special case.

In section 4.3, we consider the Prony series approximation to the general nonlinear model, which

offers more flexibility than the special case considered in section 4.2.1, and discuss the

corresponding approximate incremental relations. Some special cases of this general Prony

series approximation are considered in section 4.4.

Prior to section 4.3.2, we assume only that the strain history is piecewise continuous unless

specifically stated otherwise. From that point on, we re-introduce the continuity and smoothness

conditions 46, since these are used in deriving the incremental relations.

4.1 The General Model

4.1.1 Basic Assumptions

Consider the viscoelastic constitutive relation:

S(t) = G

¡

E(t), 0

¢

+

Z t

−∞

.

G

¡

E(s), t − s

¢

ds

= G

¡

E(t), 0

¢

+

Z ∞

0

.

G

¡

E(t − s), s¢ds ,

(115)

where

G : Sym × [0,∞) → Sym (116)

is a smooth (and generally nonlinear) function of both arguments. Here Sym denotes the space of

symmetric tensors, and

.

G denotes the partial derivative of G with respect to its 2nd (temporal)

argument:

.

G(A, s) := ∂

∂s G(A, s) (117)

for any (fixed) symmetric tensor A. We assume that G satisfies

G(0, s) = 0 , ∀ s ≥ 0 ; (118)

27
then also

.

G(0, s) = 0 , ∀ s ≥ 0 . (119)

The conditions 118 and 119 guarantee that S(t) = 0 for t < t0 if E(t) = 0 for t < t0, consistent

with the assumption that the reference configuration is a natural state. Note that equations 115,

17, and 35 or 36 combine to give a properly invariant constitutive relation. A viscoelastic

constitutive model characterized by the above relations is often referred to as a Pipkin–Rogers

model10 ; cf. (9). The material is isotropic iff G is an isotropic function of its tensor argument.

The linear, isotropic model considered in section 3 is a special case of the above: if

G(E, s) = 2G(s)E , (120)

then

.

G(E, s) = 2

.

G(s)E , (121)

and equation 115 reduces to the constitutive relation 25. In fact, with certain qualifications,

2G(s)E is the first-order approximation to G(E, s). We have

G(E, s) = G(0, s) + DG(0, s)[E] + Os(E

2

), (122)

where Os(E

2

) denotes a function of E and s, which, for fixed s, is of order E

2

. Here DG(0, s)

denotes the derivative of G with respect to its first (tensor) argument, evaluated at 0; it is a

fourth-order tensor that acts linearly on E to produce a second-order tensor. If we assume that

the material is isotropic, then

DG(0, s)[E] = 2G(s)E + λ(s)(tr E)I

for some functions G and λ. But since tr E is approximately the volumetric strain (cf.

equation 5), and since the response to volumetric strain is to be modeled by an elastic relation for

p or p ̄ (cf. section 3.2.2), for the nearly incompressible materials considered here, it seems

reasonable to restrict attention to the case where λ ≡ 0. On using this assumption, together with

10Typically, this term is applied to models for which the entire 2nd Piola–Kirchhoff stress tensor S is given by a

relation analogous to equation 115, i.e., when equation 115 holds with S replaced by S and G replaced by some

function G: S(t) = G

¡

E(t), 0

¢

+

R t

−∞

.

G

¡

E(s), t − s

¢

ds. But the case where S is given by equation 17 with S

and p ̄ given by equations 115 and 36 can be obtained from this by taking G(E, s) = −Jp(J)C

−1 + G(E, s) and

noting that

.

G(E, s) =

.

G (E, s).

28
the condition 118, we see that equation 122 reduces to

G(E, s) = 2G(s)E + Os(E

2

). (123)

Then the linear constitutive model considered in section 3 can be regarded as a first-order

approximation to the nonlinear model considered in this section.

Depending on the particular function G, it may be simpler to express G(E, s) as a function of C

and s using the relations 3, say

G(E, s) = G(C, s). (124)

An example where this is the case is discussed in section 4.4.2. Then

.

G(E, s) =

.

G(C, s), (125)

and since C = I when E = 0, the conditions 118 and 119 imply that

G(I, s) =

.

G(I, s) = 0 , ∀ s ≥ 0 . (126)

Henceforth, as before (cf. equation 30), we assume that the material is undeformed for times

t < 0, so that E(t) = 0 for t < 0. Then S(t) = 0 for t < 0, and the constitutive relation 115

reduces to

S(t) = G0

¡

E(t)

¢

+

Z t

0

.

G

¡

E(s), t − s

¢

ds

= G0

¡

E(t)

¢

+

Z t

0

.

G

¡

E(t − s), s¢ds ,

(127)

where we have introduced the instantaneous elastic response function G0, defined by

G0(E) := G(E, 0) = G(C, 0) =: G0(C) ; (128)

in view of equations 118 and 126 it satisfies

G0(0) = G0(I) = 0 . (129)

The (equivalent) hereditary integrals above are well-defined for strain histories that are only

piecewise continuous, and they are continuous functions of t, so that past jump discontinuities in

29
strain do not affect the current stress. If E has a jump discontinuity at time t, then equation 127

holds with the terms S(t) and E(t) replaced by their one-sided limits. Then on taking the

difference of these limiting relations, we see that the integrals cancel, so that the jump in S at

time t is given by

[[S]](t) = G0

¡

E(t

+)

¢

− G0

¡

E(t

−)

¢

. (130)

As before, we assume that E(0) equals its limiting value as t approaches 0 from above. Then

equation 43 holds, and equation 44 is replaced with

S(0) = S(0+) = [[S]](0) = G0

¡

E(0)¢; (131)

that is, G0

¡

E(0)¢is the instantaneous step in the stress S corresponding to a step in strain of

amount E(0) from the undeformed state.

4.1.2 Stress Relaxation Tests

Now consider the stress relaxation test 50, i.e., E(t) = E0 and C(t) = C0 for t ≥ 0 . Then for

any time t ≥ 0,

S(t) = G(E0, 0) + Z t

0

.

G(E0, t − s) ds ;

and since

.

G(E0, t − s) = −

d

ds G(E0, t − s),

the integral above reduces to −

£

G(E0, 0) − G(E0, t)

¤

. Thus for any time t ≥ 0 , the stress S(t)

in a stress relaxation test is given by

S(t) = G(E0, t) = G(C0, t). (132)

In view of this result, G and G are called the stress relaxation functions (cf. (9)), although for the

linear case in equation 120 it is more common to refer to G as the stress relaxation function.

When we need to distinguish between these, we refer to G as the scalar stress relaxation function,

and to G and G as the tensor stress relaxation functions. Note that, in principle, the tensor stress

relaxation functions are completely determined by the response in ideal stress relaxation tests.11

By setting t = 0 in equation 132, we see that the stress at time zero in a stress relaxation test is

11This would no longer be the case for a more general single integral law in which the integrand in equation 127

also depends on the current value E(t) of the strain.

30
given by the instantaneous elastic response function:

S(0) = G(E0, 0) = G0(E0) = G0(C0), (133)

consistent with equation 131.

It is assumed that the stress S(t) in the stress relaxation test 50 approaches a limiting value

S(∞), called the equilibrium stress, as t → ∞. Then by equation 132 it follows that

S(∞) = limt→∞G(E0, t) = G(E0, ∞) =: G∞(E0)

= G(C0,∞) =: G∞(C0).

(134)

The functions G∞ and G∞ defined above are called the equilibrium elastic response functions.

By equations 118 and 126, we have

G∞(0) = G∞(I) = 0 , (135)

which implies that the equilibrium stress corresponding to zero strain is zero, consistent with the

assumption that the reference configuration is a natural state. Note that the equilibrium and

instantaneous elastic response functions are related by

G∞(E0) = G0(E0) + Z ∞

0

.

G(E0, s) ds ; (136)

an analogous relation holds for G∞, G0, and

.

G.

4.1.3 Rate Form of the Relation for S

The constitutive relations in section 4.1.1 are valid for strain histories that are only piecewise

continuous. For any instant s at which

.

E(s) exists, the chain rule gives

d

ds G

¡

E(s), t − s

¢

= DG

¡

E(s), t − s

¢

[

.

E(s)] −

.

G

¡

E(s), t − s

¢

. (137)

Recall that DG denotes the derivative of G with respect to its first (tensor) argument, so that

DG

¡

E(s), t − s

¢

is a fourth-order tensor that acts linearly on the strain rate

.

E(s) to produce a

second-order tensor.

On solving equation 137 for

.

G

¡

E(s), t − s

¢

and substituting the result into equation 127, we find

31
that

S(t) = G

¡

E(0), t¢+

Z t

0

DG

¡

E(s), t − s

¢

[

.

E(s)] ds , (138)

provided the strain history satisfies the stronger conditions in equation 46. In particular, for the

stress relaxation test 50, for which E(0) = E0 and

.

E(s) = 0 for s > 0, equation 138 reduces to

equation 132. If the strain history is continuous for all times (including t = 0), then we must

have E(0) = 0, so on using equation 118, we see that equation 138 reduces to

S(t) = Z t

0

DG

¡

E(s), t − s

¢

[

.

E(s)] ds . (139)

The presence of DG, rather than

.

G, in the integrands above may result in a more complicated

form of the constitutive relation as compared with equation 127. An exception to this statement

occurs when G is given by the linear relation 120; then

DG(E, s)[A] = 2G(s)A (140)

for any symmetric tensor A, and equations 138 and 139 reduce to equations 47 and 48,

respectively. A similar simplification in equations 138 and 139 holds for the separable case

considered in the section 4.2.

4.2 The Separable Case

As a simple example of the general theory in section 4.1, we consider the special case where the

stress relaxation function G(E, s) is a separable function of E and s.

4.2.1 Constitutive Assumptions for the Separable Case

Assume that G satisfies

G(E, s) = 2G(s)Φ(E), (141)

with G as in section 3, and Φ : Sym → Sym a smooth (and generally nonlinear), dimensionless

function satisfying

Φ(0) = 0 . (142)

The material is isotropic iff Φ is an isotropic function. In this case, if equation 123 also holds,

then Φ must satisfy

Φ(E) = E + O(E

2

). (143)

32
The simple, linear, isotropic model considered in section 3 is a special case of the above with the

remainder term in equation 143 equal to 0.

From equation 141 we have

.

G(E, s) = 2

.

G(s)Φ(E). (144)

And from equations 141, 128 and 134,

G0(E) = 2G(0)Φ(E), G∞(E) = 2G(∞)Φ(E). (145)

Then the general constitutive relation 127 reduces to

S(t) = 2G(0)Φ

¡

E(t)

¢

+

Z t

0

2

.

G(t − s)Φ

¡

E(s)

¢

ds

= 2G(0)Φ

¡

E(t)

¢

+

Z t

0

2

.

G(s)Φ

¡

E(t − s)

¢

ds .

(146)

And integration by parts yields the alternative form

S(t) = 2G(t)Φ

¡

E(0)¢+

Z t

0

2G(t − s)

d

ds Φ

¡

E(s)

¢

ds , (147)

provided the strain history satisfies the stronger conditions in equation 46. When E(0) = 0, this

reduces to

S(t) = Z t

0

2G(t − s)

d

ds Φ

¡

E(s)

¢

ds . (148)

The relation 148, or something analogous to it, is often encountered in the literature on nonlinear

viscoelasticity.

Note that equation 139 reduces to equation 148 when G is given by equation 141. Indeed, in this

case, we have

DG(E, s)[A] = 2G(s)DΦ(E)[A] , (149)

for any symmetric tensor A, so

DG

¡

E(s), t − s

¢

[

.

E(s)] = 2G(t − s)DΦ

¡

E(s)

¢

[

.

E(s)]

= 2G(t − s)

d

ds Φ

¡

E(s)

¢

.

(150)

33
4.2.2 Alternative Notation for the Separable Case

Depending on the particular form of the function Φ, it may be simpler to express Φ(E) as a

function of C using the relations 3, say Φ(E) = Ψ(C). Then by the condition 142, we must

have

Ψ(I) = 0 . (151)

Also, for use in the relations above, as well as for the discussion below, it is convenient to

introduce the notation

E := Ψ(C) = Φ(E) = E + O(E

2

), (152)

where the approximation on the right is just equation 143. Then we may regard E as a

generalized strain tensor. Since

E(s) = Φ

¡

E(s)

¢

= Ψ

¡

C(s)

¢

, (153)

we have

.

E(s) = d

ds Φ

¡

E(s)

¢

= DΦ

¡

E(s)

¢

[

.

E(s)] , (154)

and similarly in terms of Ψ. With this notation, equation 146 can be written as

S(t) = 2G(0)E(t) + Z t

0

2

.

G(t − s)E(s) ds

= 2G(0)E(t) + Z t

0

2

.

G(s)E(t − s) ds ;

(155)

and equation 147 can be written as

S(t) = 2G(t)E(0) + Z t

0

2G(t − s)

.

E(s) ds

= 2G(t)E(0) + Z t

0

2G(s)

.

E(t − s) ds ;

(156)

and then equation 148 is just equation 156 with E(0) = 0. These relations are completely

analogous to those for the linear theory in section 3; cf. equations 31, 47, and 48, respectively. In

particular, when Φ(E) = E, so that E = E, we recover the relations in section 3.

34
4.2.3 Prony Series and Incremental Relations for the Separable Case

Now assume the Prony series approximation 64 for G in equation 141. Then

G(E, t) = Ã2μ∞ +

X

N

n=1

2μne

−t/τn

!

Φ(E)

= 2μ∞Φ(E) +X

N

n=1

2μne

−t/τnΦ(E).

(157)

Also, assume the continuity conditions 46, so that both equations 155 and 156 hold. Then all of

the relations in section 3.4.1 hold with E replaced by E everywhere. Because of this close

correspondence, we will not bother to write these new relations here, but simply refer to them as

the “modified” relations in the discussion below.

Regarding the incremental forms of the modified relations for An, An, and A

n

, first note that

even if E(t) is a piecewise linear function of t, E(t) generally will not be, due to the possible

nonlinearity of Φ. However, if E is approximately linear over the time interval [tk−1, tk], then so

is E, provided ∆tk is sufficiently small. When this is the case, the incremental relations 81–83

with E → E yield approximate incremental relations for the modified tensors An, An, and A

n

.

These incremental relations are equivalent: any one of them, together with equations 75–77

(with E → E) and equation 89, implies the other two. Consequently, they yield equivalent

approximate incremental relations for S(t) via the modified form of the relations 85–87 for

S(tk). Other approximate incremental relations can be obtained from those in section 3.4.4 via

the replacements E → E and

.

E →

.

E.

4.3 Prony Series Approximation for the General Model

For the special case considered in section 4.2.3, where G has the separable form 141 with G given

by the Prony series 64, the stress relaxation function is given exactly by equation 157. This can

also be regarded as an approximation to a general stress relaxation function G in the

Pipkin-Rogers model. Clearly, a better approximation could be obtained if we allowed the

function Φ in the nth term of equation 157 to vary with n. This is the case considered here.

35
4.3.1 Prony Series Approximation for the Stress Relaxation Function

Assume that the stress relaxation function G in section 4.1.1 is given (or approximated) by the

Prony series

G(E, t) = 2μ∞Φ∞(E) +X

N

n=1

2μne

−t/τnΦn(E), (158)

where each Φn : Sym → Sym (n = 1, . . . , N, ∞) is a smooth (and generally nonlinear),

dimensionless function, and the moduli μ∞ and μn and the relaxation times τn are positive. Then

.

G(E, t) = −

X

N

n=1

2μn

τn

e

−t/τnΦn(E). (159)

The material is isotropic if each Φn is an isotropic function, which will be assumed here.

Depending on the particular form of the function Φn, it may be simpler to express Φn(E) as a

function Ψn of C using the relations 3; an example where this is the case is discussed in

section 4.4.2. Thus, for n = 1, . . . , N, ∞, we have12

Ψn(C) := Φn(E) = E + On(E

2

), (160)

where the approximation assumed on the right guarantees that the approximation 123 holds with

G(s) given by the Prony series 64:

G(E, t) = 2G(t)E + Ot(E

2

), G(t) = μ∞ +

X

N

n=1

μne

−t/τn

. (161)

In particular,

Φn(0) = Ψn(I) = 0 . (162)

Thus, as the strain E becomes smaller, the Prony series approximation 158 for the general stress

relaxation function G in section 4.1.1 approaches the Prony series approximation for the linear

theory considered in section 3. It follows that the moduli μ∞ and μn and the relaxation times τn

can be determined from the viscoelastic response to small deformations. Departure from the

linear viscoelastic theory enters only through the time-independent functions Φn. This is an

advantage (as well as a limitation) of the approximation 158.

12On(E

2

) denotes a function of E and n which, for fixed n, is of order E

2

.

36
From equation 158 it follows that the equilibrium elastic response function is given by

G∞(E) := G(E,∞) = 2μ∞Φ∞(E), (163)

with μ∞ again representing the equilibrium elastic shear modulus. We call Φ∞ the normalized

equilibrium elastic response function. The instantaneous elastic response function is given by

G0(E) := G(E, 0) = 2μ∞Φ∞(E) +X

N

n=1

2μnΦn(E). (164)

This may also be written as

G0(E) = 2μ0Φ0(E), (165)

where the instantaneous elastic shear modulus μ0 is defined as in equation 65, i.e.,

μ0 := μ∞ +

X

N

n=1

μn , (166)

and the dimensionless function Φ0 is defined by

Φ0(E) := μ∞

μ0

Φ∞(E) +X

N

n=1

μn

μ0

Φn(E). (167)

We call Φ0 the normalized instantaneous elastic response function. Note that the sum of the

dimensionless coefficients in equation 167 is unity, so on using the small strain approximations

from equation 160 we see that Φ0(E) = E + O0(E

2

), i.e., equation 160 also holds for n = 0.

Finally, note that in view of equation 1601, the relations 158, 159, 163–165 and 167 hold with the

replacements G → G, E → C and Φn → Ψn (n = 0, 1, . . . , N, ∞).

4.3.2 Relations for the Stress

On substituting equation 159 (with E → E(s) and t → t − s) and equation 165 into the

constitutive relation 1271 for S(t), we obtain

S(t) = 2μ0Φ0

¡

E(t)

¢

−

X

N

n=1

2μnAn(t), (168)

where the dimensionless, symmetric tensor An(t) is defined by

An(t) := Z t

0

1

τn

e

−(t−s)/τnΦn

¡

E(s)

¢

ds , (169)

37
so that

An(0) = 0 . (170)

The relations 168–170 generalize the relations 68–70 for the linear case.

Henceforth, we again assume the continuity and smoothness conditions 46 on the strain history.

Then on integrating equation 169 by parts, we obtain

An(t) = Φn

¡

E(t)

¢

− e

−t/τnΦn

¡

E(0)¢− An(t), (171)

where the dimensionless, symmetric tensor An(t) is defined by

An(t) := Z t

0

e

−(t−s)/τn

d

ds Φn

¡

E(s)

¢

ds , (172)

so that

An(0) = 0 . (173)

On substituting equation 171 into equation 168 and using equation 167, we obtain

S(t) = 2μ∞Φ∞

¡

E(t)

¢

+

X

N

n=1

2μnAn(t) +X

N

n=1

2μne

−t/τnΦn

¡

E(0)¢. (174)

The relations 174, 172, and 173 generalize the relations 71–73 for the linear case.

4.3.3 Alternative Notation

The expressions above for the tensors S(t), An(t), and An(t), as well as the approximate

incremental relations for them, which are given in the next section, can be simplified somewhat

by using introducing the following notation:13

E

n

(s) := Φn

¡

E(s)

¢

= Ψn

¡

C(s)

¢

, n = 0, 1, . . . , N, ∞. (175)

On using equation 175 in equations 168 and 169, we obtain

S(t) = 2μ0E

0

(t) −

X

N

n=1

2μnAn(t), (176)

An(t) = Z t

0

1

τn

e

−(t−s)/τn E

n

(s) ds , (177)

13Note that E

n

does not denote E raised to the nth power.

38
which more closely resemble the relations 68–69 for the linear case.

From equation 175, the material time derivative of E

n

is

.

E

n

(s) := d

ds E

n

(s) = d

ds Φn

¡

E(s)

¢

= DΦn

¡

E(s)

¢

[

.

E(s)] ,

(178)

and similarly in terms of Ψn. On using equations 178 and 175 in equations 174 and 172, we

obtain

S(t) = 2μ∞E

∞(t) +X

N

n=1

2μnAn(t) +X

N

n=1

2μne

−t/τn E

n

(0), (179)

An(t) = Z t

0

e

−(t−s)/τn

.

E

n

(s) ds , (180)

which more closely resemble the relations 71–72 for the linear case.

Also, by equation 175 the relation 171 between An(t) and An(t) can be written as

E

n

(t) − An(t) = An(t) + e

−t/τn E

n

(0) =: A

n

(t), (181)

which more closely resembles the relations 75–77 for the linear case. On using the relation on

the right above, we see that equation 179 can be written as

S(t) = 2μ∞E

∞(t) +X

N

n=1

2μnA

n

(t), (182)

analogous to the relation 74 for the linear case. And by equations 181, 173 or 170, and 175, we

see that

A

n

(0) = E

n

(0) = Φn

¡

E(0)¢= Ψn

¡

C(0)¢. (183)

The condition E(0) = 0 holds when the strain history is continuous for all times, and in this case

Φn

¡

E(0)¢= 0 by equation 162, so A

n

(0) = E

n

(0) = 0 by equation 183; then A

n

(t) = An(t) by

equation 181, the sum on the far right in equation 179 drops out, and the relations 179 and 182 are

indistinguishable.

39
Finally, note that the definition 175 and the condition 160 imply

E

n

(t) = E(t) + On

¡

E(t)

2

¢

, n = 0, 1, . . . , N, ∞. (184)

Thus we may think of the E

n

as generalized strain tensors. And substitution of equation 184 into

equations 176–177, 179–180, and 181–182 yields the relations for the linear theory, namely

equations 68–69, 71–72, and 74–77, respectively, as small strain approximations.

4.3.4 Approximate Incremental Relations

Even if E(t) is a piecewise linear function of t, the tensors E

n

(t) generally will not be, due to the

possible nonlinearity of the functions Φn in equation 175. On the other hand, given the

continuity conditions 46, we can assume that the sequence of times tk in equation 79 are chosen

so that not only E but each E

n

is approximately linear over each time interval [tk−1, tk]. Then,

on using the relations for An, An, and A

n

in section 4.3.3, and proceeding as in section 3.4.2, we

can derive approximate incremental relations for An, An, and A

n

. By equation 175, we have

E

n

(tk) = Φn

¡

E(tk)

¢

= Ψn

¡

C(tk)

¢

(185)

for n = 0, 1, . . . , N, ∞. And as before, we let ∆tk = tk − tk−1. Then, we find that

An(tk) ≈ e

−∆tk/τnAn

(tk−1)

+ α(∆tk/τn)E

n

(tk) + β(∆tk/τn)E

n

(tk−1), (186)

An(t0) = An(0) = 0 ; (186a)

and

An(tk) ≈ e

−∆tk/τnAn

(tk−1) + f(∆tk/τn)

£

E

n

(tk) − E

n

(tk−1)

¤

, (187)

An(t0) = An(0) = 0 ; (187a)

and

A

n

(tk) ≈ e

−∆tk/τn An

(tk−1) + f(∆tk/τn)

£

E

n

(tk) − E

n

(tk−1)

¤

, (188)

A

n

(t0) = A

n

(0) = E

n

(0). (188a)

40
These are analogous to the incremental relations 81–83a for the linear case.

Note that the incremental relations for Anand A

n have the same form but different initial

conditions. Furthermore, all three incremental relations are equivalent: any one of them,

together with the relations 181 and 89, implies the other two. Consequently, they yield

equivalent approximate incremental relations for S via equations 176, 179 and 182.

Other approximate relations may be obtained as in section 3.4.4. Since, by assumption, E

n

is

approximately linear over each time interval [tk−1, tk], its rate

.

E

n will be approximately constant

over that time interval. Let (

.

E

n

)k denote any approximation for

.

E

n over [tk−1, tk]. Five

(generally different) approximations, analogous to those in equation 95, are given below:

(

.

E

n

)k =







.

E

n

(t

−

k

) := lim

t ↑ tk

.

E

n

(t),

.

E

n

(t

+

k−1

) := lim

t ↓ tk−1

.

E

n

(t),

.

E

n

(t

−

k

) +

.

E

n

(t

+

k−1

)

2

,

.

E

n

μ

tk + tk−1

2

¶

,

E

n

(tk) − E

n

(tk−1)

tk − tk−1

=

E

n

(tk) − E

n

(tk−1)

∆tk

.

(189)

As before, we have allowed for the possibility that

.

E

n may suffer a jump discontinuity at tk

and/or tk−1. The last choice for (

.

E

n

)k above is equivalent to the relation

E

n

(tk) − E

n

(tk−1) = ∆tk(

.

E

n

)k . (190)

When this is substituted into the incremental relation 187 for Anor the incremental relation 188

for A

n

, we obtain alternative forms of these relations. For example, equation 187 is equivalent to

the incremental relation

An(tk) ≈ e

−∆tk/τnAn

(tk−1) + τn

¡

1 − e

−∆tk/τn

¢

(

.

E

n

)k (191)

when (

.

E

n

)k is given by the bottom expression in equation 95; this is analogous to equation 97.

Then equation 191 with any of the other choices for (

.

E

n

)k in equation 95 yields a (generally)

different approximate incremental relation for An. These relations and their analogs for A

n yield

41
approximate incremental relations for S via equations 179 and 182. Finally, one can proceed as

at the end of section 3.4.4 and obtain approximate incremental relations for S in a form analogous

to equations 98–100.

4.4 Special Cases of the Prony Series Approximation

Regarding the Prony series approximation 158 for the stress relaxation function G, we have not

placed any restrictions on the functions Φn (n = 1, . . . , N, ∞) other than those discussed in the

first paragraph in section 4.3.1. In particular, the relations in section 4.3 do not require that these

functions be related to each other in any simple way. Here we consider some special cases of this

general theory.

4.4.1 A Simple Non-Separable Case

When Φ1 = · · · = ΦN = Φ∞ =: Φ, the Prony series 158 reduces to equation 157, and we

recover the Prony series approximation to the separable stress relaxation function considered in

section 4.2. A special case of equation 158 that is not necessarily separable occurs when

Φ1 = · · · = ΦN =: Φ , (192)

with Φ possibly distinct from the normalized equilibrium elastic response function Φ∞. Then

the stress relaxation function in equation 158 reduces to

G(E, t) = 2μ∞Φ∞(E) + ÃX

N

n=1

2μne

−t/τn

!

Φ(E). (193)

As before, the equilibrium elastic response function is given by equation 163:

G∞(E) := G(E,∞) = 2μ∞Φ∞(E). (194)

And on setting t = 0 in equation 193 and using equation 166, or by using equations 192 and

164–167, we see that the instantaneous elastic response function is given by

G0(E) = 2μ∞Φ∞(E) + 2(μ0 − μ∞)Φ(E) = 2μ0Φ0(E), (195)

where the normalized instantaneous elastic response function is given by

Φ0(E) = μ∞

μ0

Φ∞(E) + μ1 −

μ∞

μ0

¶

Φ(E). (196)

42
In view of equation 192, the tensors E

n

defined in equation 175 are identical for n = 1, . . . , N:

E

n

(t) = E(t) := Φ

¡

E(t)

¢

= Ψ

¡

C(t)

¢

, n = 1, . . . , N . (197)

However, the tensors E

0

(t) and E

∞(t), which appear in the relations 176, 179, and 182 for S(t),

generally do not coincide with E(t): by equations 175 and 196, we have

E

∞(t) := Φ∞

¡

E(t)

¢

= Ψ∞

¡

C(t)

¢

, (198)

and

E

0

(t) := Φ0

¡

E(t)

¢

= Ψ0

¡

C(t)

¢

=

μ∞

μ0

E

∞(t) + μ

1 −

μ∞

μ0

¶

E(t).

(199)

When Φ∞ = Φ we recover the separable case in section 4.2. And when Φ∞(E) = Φ(E) = E,

this reduces to the linear case in section 3..

The Simplest Non-Separable Case: The simplest special case of the above, which does not

reduce to the separable case, occurs when the normalized equilibrium elastic response function

Φ∞ is allowed to be nonlinear but

Φ(E) = E , (200)

in which case

E

n

(t) = E(t) = E(t), n = 1, . . . , N . (201)

In other words, the equilibrium elastic response may be nonlinear, but the viscoelastic part of the

response is linear in E.

In this case the relations 177, 180, and 181 for An, An, and A

n

reduce to the relations 69, 72, and

75–76 of the linear theory. Hence the incremental relations for An, An, and A

n

in section 4.3.4

reduce to those of the linear theory as well. And just as in the linear theory, these incremental

relations are exact for continuous, piecewise-linear strain histories (with a possible jump

discontinuity at t = 0 only), provided that any jump discontinuities in the strain rate

.

E occur at

one of the discrete times tk.

We may use the conclusions in the preceding paragraph to derive exact relations for S(t) for

simple piecewise-linear strain histories. The results are similar to those in section 3.4.5, except

43
that any terms there involving μ∞ require appropriate modification to account for the fact that the

equilibrium elastic response function G∞ is no longer given by G∞(E) = 2μ∞E. For example,

for the constant strain-rate case with a possible initial jump in strain (cf. equations 101–102), the

relation 104 for S(t) needs to be changed to

S(t) = 2μ∞Φ∞

¡

E(t)

¢

+

"X

N

n=1

2μnτn

¡

1 − e

−t/τn

¢

#

.

E1 +

"X

N

n=1

2μne

−t/τn

#

E(0). (202)

4.4.2 One-Parameter Family of Functions: Part I

Next, we consider a stress relaxation function G that is intermediate in generality between the one

in the the previous section, equation 193, and the stress relaxation function for the general Prony

series 158 with the functions Φn unrelated. While the flexibility inherent in the latter case allows

for a wider range of materials to be approximated than would be the case, say, with the separable

stress relaxation function in equation 157 or the slightly more general stress relaxation function in

equation 193, when it comes to calibrating the model for a particular material we are eventually

forced to choose specific functional forms for the Φn or Ψn (recall that Ψn(C) := Φn(E)). This

task is simplified if each of these functions is assumed to belong to some specific subclass of the

set of all isotropic functions satisfying equation 160, and in particular if the functions in this

subclass are generated by a single real parameter.

Thus let Ψ and Φ be isotropic functions satisfying

Ψ(ω, C) = Φ(ω, E) = E + Oω(E

2

), (203)

where ω is a non-dimensional real parameter. And for each for n = 1, . . . , N, assume that

Φn(E) is given by Φ(ω, E) and Ψn(C) by Ψ(ω, C) for some particular choice of ω which may

depend on n, say ω = ωn. Then for n = 1, . . . , N we have

Φn(E) Φ(ωn, E)

°

°

°

°

°

°

Ψn(C) Ψ(ωn, C),

(204)

and by equation 203 the condition 160 is satisfied. We do not necessarily impose the

assumption 204 for n = ∞, although as a special case we might do so, as in the example that

follows. Note that if condition 204 is not imposed for n = ∞, then it generally will not hold for

n = 0 either, since Φ0 and Φ∞ are related by equation 167.

44
As an illustration of the above, we take

Ψ(ω, C) = ω

2

(I − C

−1

) + 1 − ω

2

(C

−1 − C−2

), 0 < ω ≤ 1 . (205)

Motivation for this choice is discussed below. By equation 4, we have

C

−2 = (C−1

)

2 = I − 4E + O(E

2

), (206)

which together with equation 4 implies that

1

2

(I − C

−1

) & 1

2

(C

−1 − C−2

) = E + O(E

2

). (207)

Then on substituting equation 207 into equation 205, we see that this choice for Ψ(ω, C) indeed

satisfies equation 203. We will assume that equation 204 holds not only for n = 1, . . . , N but

also for n = ∞. And each of the ωn in equation 204 must satisfy the inequality in equation 205:

0 < ωn ≤ 1 (208)

for n = 1, . . . , N, ∞. Then on using equations 204–206, 208, 167, and 166, we find that the

normalized instantaneous elastic response function is also of the form 205. More precisely,

Φ0(E) = Ψ0(C) = Ψ(ω0, C) (209)

for

ω0 :=

μ∞

μ0

ω∞ +

X

N

n=1

μn

μ0

ωn , (210)

and ω0 satisfies the inequality 0 < ω0 ≤ 1. Thus equations 208 and 204 hold for n = 0 also.

From the relations above, equation 165, and equation 163, we see that the instantaneous and

equilibrium elastic response functions are given by

G0(E) = G0(C) = 2μ0Ψ(ω0, C) (211)

and

G∞(E) = G∞(C) = 2μ∞Ψ(ω∞, C). (212)

This example is motivated by the Mooney-Rivlin model for incompressible elastic materials. We

45
give a brief discussion of that model in section 4.4.3 and then return to the discussion of the

example above.

4.4.3 The Mooney-Rivlin Elastic Model

The Mooney-Rivlin constitutive model for an isotropic, incompressible, elastic solid is given by14

σ =  ̃pI + μ

+B − μ−B−1

, (213)

where p ̃ is indeterminate and the coefficients μ

+ and μ− are elastic moduli satisfying

μ

+ > 0 , μ− ≥ 0 . (214)

B is the left Cauchy-Green deformation tensor:

B := F FT = V

2 = RCRT

; (215)

and

B−1 = F

−T F−1 = V−2

. (216)

On taking the small strain limit of equation 213, one finds that the moduli μ

+ and μ− are related

to the linear elastic shear modulus μ by

μ

+ + μ− = μ . (217)

If we define the non-dimensional parameter ω by

ω =

μ

+

μ+ + μ−

=

μ

+

μ

= 1 −

μ

−

μ

, (218)

then

ω > 0 and 0 ≤

μ

−

μ

= 1 − ω , (219)

so

μ

+ = μ ω , μ− = μ (1 − ω) and 0 < ω ≤ 1 . (220)

14Cf. Truesdell and Noll (2, §95). The indeterminate scalar p ̃ is sometimes referred to as the “pressure” and often

denoted by p, although it is not equal to the pressure p as we have defined it in equation 101. The latter is given by

p = −p ̃−

1

3

(μ

+tr B − μ− tr B

−1

).

46
Thus we may write the Mooney-Rivlin constitutive relation 213 as

σ = −p ̃I + μ

£

ωB − (1 − ω)B−1

¤

, 0 < ω ≤ 1 . (221)

The special case ω = 1 (equivalently, μ

− = 0) is called a neo-Hookean material.

Next, note that by equations 215, 216, and 41, we have

F

−1F−T = C−1

, F

−1BF −T = I , F−1B−1F−T = C−2

. (222)

Then by equations 213, 221, 222, and 71 and the fact that J ≡ 1 for an incompressible material,

we see that the Mooney-Rivlin constitutive relation is equivalent to

S = −p ̃C

−1 + μ

+

I − μ

−C

−2

= −p ̃C

−1 + μ

£

ωI − (1 − ω)C

−2

¤

.

(223)

This may be re-written in the equivalent forms

S = −p∗C

−1 + μ

+

¡

I − C

−1

¢

+ μ

−

¡

C

−1 − C−2

¢

= −p∗C

−1 + 2μΨ(ω, C),

(224)

where p∗ is also indeterminate but different than p ̃, and Ψ(ω, C) is given by equation 205.

Though the relations 224 are slightly more complicated than the equivalent relations 223, they are

more useful for our purposes. Indeed, as noted in the previous section, the function Ψ(ω, C) in

equation 2242 satisfies the requirement 203. On the other hand, if we rewrite equation 2232 as

S = −p ̃C

−1 + 2μ ·

1

2

£

ωI − (1 − ω)C

−2

¤

, (225)

we find that

1

2

£

ωI − (1 − ω)C

−2

¤

6= E + O(E

2

). (226)

4.4.4 One-Parameter Family of Functions: Part II

Now return to the discussion in section 4.4.2. When the assumptions 204 hold for

n = 1, . . . , N, ∞, the Prony series approximation 158 to the stress relaxation function G(E, t),

47
when expressed in terms of C, takes the form

G(C, t) = 2μ∞Ψ(ω∞, C) +X

N

n=1

e

−t/τn 2μnΨ(ωn, C). (227)

As mentioned previously, the Mooney-Rivlin model provides the motivation for choosing the

one-parameter family of functions Ψ(ω, C) to have the particular form in equation 205 . Each

term in this Prony series involves (the determinate part of) a particular Mooney-Rivlin

constitutive function, as is clear from equation 2242.

In particular, on comparing equations 211–212 with equation 2242, we see that the instantaneous

and equilibrium elastic response functions for the viscoelastic part S of the model are given by

(different) Mooney-Rivlin constitutive functions. If we recall the decomposition 17,

S = −Jp ̄C

−1 + S , (228)

and the constitutive assumptions for p ̄ in section 3.2.2, and compare with equation 2242, we see

that the instantaneous and equilibrium elastic response functions for the total 2nd Piola-Kirchhoff

stress tensor S are given by (different) compressible versions of the Mooney-Rivlin model.

Next, we examine more closely the form of the stress relaxation function in equation 227 for

arbitrary t when Ψ(ω, C) is given by equation 205. For this choice, the term 2μnΨ(ωn, C) in

equation 227 can be written as

2μnΨ(ωn, C) = μ

+

n

(I − C

−1

) + μ

−

n

(C

−1 − C−2

), (229)

where, analogous to equation 220, we have (for n = 1, . . . , N, ∞)

μ

+

n

:= μn ωn > 0 , μ−

n

:= μn (1 − ωn) ≥ 0 . (230)

Then on substituting equation 229 into equation 227, we see that the stress relaxation function is

given by

G(C, t) = μ

+(t)(I − C

−1

) + μ

−(t)(C

−1 − C−2

), (231)

where the functions μ

+ and μ− are given by the Prony series

μ

±(t) := μ±

∞ +

X

N

n=1

μ

±

n

e

−t/τn

. (232)

48
On comparing equation 231 with equation 2241, we see that at each instant t the function

C 7→ G(C, t) is (the determinate part of) a particular Mooney-Rivlin constitutive function. It

follows that the total stress in a stress relaxation test is given by a compressible and

“time-dependent” version of the Mooney-Rivlin model.

On setting t = ∞ and t = 0 in equation 232, we obtain

μ

±

(∞) = μ

±

∞ , (233)

and

μ

±

(0) = μ

±

∞ +

X

N

n=1

μ

±

n =: μ

±

0

. (234)

Then on using equation 234, the definition 166 of μ0, and the definition 210 of ω0, we find that

equation 230 also holds for n = 0. And on setting t = 0 and t = ∞ in equation 231 and using

equations 234, 233, 230 (for n = 0 and n = ∞), and equation 205, we recover equations 211 and

212.

The stress relaxation function in equations 230–232 represents a substantial simplification of the

general stress relaxation function G(E, s) = G(C, s) in the Pipkin-Rogers model as well as a

substantial simplification of the general Prony series approximation 158 to it. However, use of

equations 230–232 does not simplify the form of the approximate incremental relations in

section 4.3.4. On the other hand, the evaluation of E

n

(tk) at each of the discrete times tk is fairly

simple once C

−1

(tk) has been determined. Indeed, by equations 185, 204, and 205,

E

n

(tk) = Ψ(ωn, C(tk))

=

ωn

2

I +

μ

1

2

− ωn

¶

C

−1

(tk) −

1 − ωn

2

£

C

−1

(tk)

¤2

.

(235)

As observed in section 4.3.1, the moduli μ∞ and μn and the relaxation times τn for the general

Prony series approximation 158 (and hence, in particular, for the special approximation

considered here) can be determined from the viscoelastic response to small strain deformations,

for which the linear theory in section 3. is adequate. For the particular model considered in the

present section, the departure from the linear theory is completely characterized by the additional

material parameters ω1, . . . , ωN , ω∞. Thus, one could determine the moduli μn and relaxation

times τn from small strain tests first, and then optimize the choice of the ωn to fit the large strain

response.

49
5. References

1. Wineman, A. S.; Rajagopal, K. R. Mechanical Response of Polymers; Cambridge University

Press: New York, 2000.

2. Truesdell, C.; Noll, W. The Non-Linear Field Theories of Mechanics. In Handbuch der

Physik; Vol. III, No. 3; Flügge, S., Ed.; Springer-Verlag: Berlin, 1965.

3. Gurtin, M. E. An Intoduction to Continuum Mechanics; Academic Press: New York, 1981.

4. Bowen, R. M. Introduction to Continuum Mechanics for Engineers; Plenum Press: New

York, 1989.

5. Holzapfel, G. H. Nonlinear Solid Mechanics; Wiley: New York, 2000.

6. Scheidler, M. Notes on Approximate Incremental Relations for Viscoelasticity, July 2008.

7. Zocher, M. A.; Groves, S. E.; Allen, D. H. A three-dimensional finite element formulation

for thermoviscoelastic orthotropic media. Int. J. Numer. Methods Eng. 1997, 40, 2267-2288.

8. Clements, B. Notes on PC Model, Oct. 2005 and Feb. 2006.

9. Pipkin, A. C.; Rogers, T. G. A non-linear integral representation for viscoelastic behaviour.

J. Mech. Phys. Solids 1968, 16, 59-72.

50
NO. OF

COPIES ORGANIZATION

1

(PDF

ONLY)

DEFENSE TECHNICAL

INFORMATION CTR

DTIC OCA

8725 JOHN J KINGMAN RD

STE 0944

FORT BELVOIR VA 22060-6218

1 DIRECTOR

US ARMY RESEARCH LAB

IMNE ALC HRR

2800 POWDER MILL RD

ADELPHI MD 20783-1197

1 DIRECTOR

US ARMY RESEARCH LAB

RDRL CIM P

2800 POWDER MILL RD

ADELPHI MD 20783-1197

1 DIRECTOR

US ARMY RESEARCH LAB

RDRL CIM L

2800 POWDER MILL RD

ADELPHI MD 20783-1197

ABERDEEN PROVING GROUND

1 DIR USARL

RDRL CIM G (BLDG 4600)

51
NO. OF

COPIES ORGANIZATION

ABERDEEN PROVING GROUND

25 DIR USARL

RDRL WMM D

G GAZONAS

C RANDOW

RDRL WMT

P BAKER

RDRL WMT A

S SCHOENFELD

RDRL WMT C

T W BJERKE

M FERMEN-COKER

S SEGLETES

RDRL WMT D

N GNIAZDOWSKI

S R BILYK

D CASEM

J CLAYTON

D DANDEKAR

M GREENFIELD

Y HUANG

R KRAFT

B LOVE

M RAFTENBERG

E RAPACKI

M SCHEIDLER (6 copies)

T WEERASOORIYA

52