Jtam.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

47, 2, pp. 397-409, Warsaw 2009

A NEW INCREMENTAL FORMULATION FOR LINEAR
VISCOELASTIC ANALYSIS: CREEP DIFFERENTIAL

APPROACH

Claude Chazal

Rostand Moutu Pitti

GEMH-GCD Laboratory, Limoges University, Centre Universitaire Génie Civil, Egletons, France

e-mail: chazal@unilim.fr; rostand.moutou-pitti@etu.unilim.fr

This paper presents a new incremental formulation in the time domain
for linear, non-aging viscoelastic materials undergoing mechanical de-
formation. The formulation is derived from linear differential equations
based on a discrete spectrum representation for the creep tensor. The
incremental constitutive equations are then obtained by finite difference
integration.Thus thedifficulty of retaining the stresshistory in computer
solutions is avoided.Acomplete general formulationof linear viscoelastic
stress analysis is developed in terms of increments of strains and stres-
ses. Numerical results of good accuracy are achieved.The analytical and
numerical solutions are compared using the energy release rate in pure
mode I and pure mode II.

Key words: incremental constitutive law, strain history, discrete creep
spectrum

1. Introduction

The important use of viscoelastic materials in civil engineering structures re-
quires understanding of the behaviour of time dependent mechanical fields
which can lead to collapse of such structures. Themain problem in computa-
tion mechanics is to know the response of a viscoelastic material taking into
account its complete past history of stress and strain.Most of analytical solu-
tions proposed in the literature can not deal with real and complex problems
because these methods require the retaining of the complete past history of
stress and strain in the memory of a digital computer.
In this context, a number of theories have been proposed in the past in

order to formulate incremental constitutive equations for the linear viscoelastic



398 C. Chazal, R. Moutou Pitti

behaviour. Among them, Kim and Sung Lee (2007), Ghazlan et al. (1995),
Chazal andDubois (2001), Klasztorny (2008), Dubois et al. (1999a) proposed
the incremental formulation and constitutive equations in the finite element
context. In fracture of viscoelastic mechanics, Dubois et al. (1999b, 2002),
Dubois and Petit (2005) and Moutou Pitti et al. (2007, 2008) applied the
incremental formulation in order to evaluate the crack growthprocess inwood.
However, the formulation usedwas based on the spectral decomposition using
a generalized Kelvin Voigt model.

To avoid the use of the generalized Kelvin Voigt model, we develop in this
paper a new incremental formulation based on a discrete creep spectrum and
the finite differencemethodusing generalized differential equations in the time
domain. The incremental stress-strain constitutive equation is not restricted
to isotropicmaterials and canbeused to resolve complex viscoelastic problems
without retaining the past history of the material.

The first section recalls the discrete creep spectrum representation and
its use in Boltzmann’s superposition principle (Boltzmann, 1878). The
one-dimensional linear viscoelastic behaviour is used to reduce the three-
dimensional response.

The second section contains the development of the generalized differential
equations in terms of one-dimensional stress and strain components.

In the third section, the solution of the differential equations is propo-
sed using the finite difference method and the new constitutive stress-strain
relations are then obtained.

Finally, the constitutive law is implemented in finite element software
CASTEM(Charvet-Quemin et al. 1986) and the numerical results are compa-
red to the analytical solution given by Moutou Pitti et al. (2007).

2. Creep spectrum representation

In this work, we consider only small strains. According to the results obtained
by Mandel (1978), Ghazlan et al. (1995), Chazal and Dubois (2001), Moutou
Pitti et al. (2007, 2008) and Dubois and Petit (2005), the components of the
creep tensor J(t) can be represented in terms of an exponential series

Jijkℓ(t)=
[
J
(0)
ijkℓ
+
M∑

m=1

J
(m)
ijkℓ

(
1−e

−tλ
(m)
ijkℓ

)]
H(t) (2.1)



A new incremental formulation for linear... 399

Where λ
(m)
ijkℓ
, m = 1, . . . ,M, are strictly positive scalars, and the repeated

indices do not imply summation convention. J
(0)
ijkℓ
and J

(m)
ijkℓ
represent the

instantaneous and the differed creep tensor, respectively, and H(t) is the He-
aviside unit step function.

According to Boltzmann’s principle (Boltzmann, 1878), the constitutive
equations between the components σij(t) of the stress tensor and the com-
ponents of the strain tensor εij(t) for linear viscoelastic materials can be
expressed in the time domain by the hereditary integral

εij(t)=
∑

k

∑

ℓ

t∫

−∞

Jijkℓ(t− τ)
∂σkℓ(τ)

∂τ
dτ (2.2)

Let us consider the fourth order tensor Π(t) of the components Πijkℓ(t) de-
fined by

Πijkℓ(t)=

t∫

−∞

Jijkℓ(t− τ)
∂σkℓ(τ)

∂τ
dτ ∀i,j,k,ℓ∈ [1,2,3], ∀t∈ℜ (2.3)

The components Πijkℓ(t) can be interpreted as the contribution of the stress
history {σkℓ(τ),τ ¬ t} of the components σkℓ(t) of the stress tensor σ(t) to
the strain components εij(t).

Introducing equation (2.3) into (2.2), we obtain

εij(t)=
∑

k

∑

ℓ

Πijkℓ(t) ∀i,j,k,ℓ∈ [1,2,3], ∀t∈ℜ (2.4)

Each equation of relation (2.4) represents a one-dimensional non-aging linear
viscoelastic material defined by its creep function J(t).

3. Formulation of differential equations

When we apply the mechanical stress defined by the stress history
{σkℓ(τ),τ ∈ ℜ}, the response of the material is then given by the history
of strains {Πijkℓ(t), t ∈ℜ} defined by behaviour equation (2.3) in which the
creep function is given by equation (2.1).

We note by σkℓ(t) the stress applied to the material at the time t and by
Πijkℓ(t) the total strain at the same time t. Then the response in strains can



400 C. Chazal, R. Moutou Pitti

be obtained using the finite creep spectrum representation given by equation
(2.1)

Πijkℓ(t)=

t∫

−∞

[
J
(0)
ijkℓ
+
M∑

m=1

J
(m)
ijkℓ

(
1− e

−λ
(m)
ijkℓ
(t−τ)
)]∂σkℓ(τ)
∂τ

dτ (3.1)

This equation can be rewritten in the following form

Πijkℓ(t)=Π
(0)
ijkℓ
(t)+

M∑

m=1

Π
(m)
ijkℓ
(t) (3.2)

with

Π
(0)
ijkℓ
(t)=

t∫

−∞

J
(0)
ijkℓ

∂σkℓ(τ)

∂τ
dτ = J

(0)
ijkℓ
σkℓ(t)

(3.3)

Π
(m)
ijkℓ
(t)=

t∫

−∞

J
(m)
ijkℓ

(
1− e

−λ
(m)
ijkℓ
(t−τ)
)∂σkℓ(τ)
∂τ

dτ

In these equations, Π
(0)
ijkℓ
(t) and Π

(m)
ijkℓ
(t) represent the instantaneous and the

differed part of the one-dimensional strain of the material.
Using equation (2.4), the rate of the total strain is determined by

∂εij(t)

∂t
=
∑

k

∑

ℓ

∂Πijkℓ(t)

∂t
=
∑

k

∑

ℓ

(∂Π(0)
ijkℓ
(t)

∂t
+
M∑

m=1

∂Π
(m)
ijkℓ
(t)

∂t

)
(3.4)

According to equation (3.3)1, the rate of the instantaneous part of the one-

dimensional strain Π
(0)
ijkℓ
(t) is given by

∂Π
(0)
ijkℓ
(t)

∂t
= J

(0)
ijkℓ

∂σkℓ(t)

∂t
(3.5)

However, the rate of the differed part of the one-dimensional strain Π
(m)
ijkℓ
(t)

is more complicated to be determined. Using equation (3.3)2, we can write

∂Π
(m)
ijkℓ
(t)

∂t
= J

(m)
ijkℓ

(
1−e

−λ
(m)
ijkℓ
(t−t)
)∂σkℓ(τ)
∂τ

+

(3.6)

+

t∫

−∞

J
(m)
ijkℓ

(
0+λ

(m)
ijkℓ
e
−λ
(m)
ijkℓ
(t−τ)
)∂σkℓ(τ)
∂τ

dτ



A new incremental formulation for linear... 401

or

∂Π
(m)
ijkℓ
(t)

∂t
= J

(m)
ijkℓ
λ
(m)
ijkℓ

t∫

−∞

e
−λ
(m)
ijkℓ
(t−τ)∂σkℓ(τ)

∂τ
dτ (3.7)

knowing that

λ
(m)
ijkℓ
Π
(m)
ijkℓ
(t)= J

(m)
ijkℓ
λ
(m)
ijkℓ

t∫

−∞

(
1−e

−λ
(m)
ijkℓ
(t−τ)
)∂σkℓ(τ)
∂τ

dτ =

= J
(m)
ijkℓ
λ
(m)
ijkℓ
σkℓ(t)−λ

(m)
ijkℓ
J
(m)
ijkℓ

t∫

−∞

e
−λ
(m)
ijkℓ
(t−τ)∂σkℓ(τ)

∂τ
dτ = (3.8)

= J
(m)
ijkℓ
λ
(m)
ijkℓ
σkℓ(t)−

∂Π
(m)
ijkℓ
(t)

∂t

The last equation can be written as a linear differential equation and can be
integrated analytically

∂Π
(m)
ijkℓ
(t)

∂t
+λ
(m)
ijkℓ
Π
(m)
ijkℓ
(t)= J

(m)
ijkℓ
λ
(m)
ijkℓ
σkℓ(t) (3.9)

The solution to this differential equation gives the rate of the one-dimensional

strain Π
(m)
ijkℓ
(t).

Finally, the general differential equations governing the non-aging linear
viscoelastic behaviour can be obtained from equation (3.4) after summation
on k and ℓ indices. One finds

∂εij(t)

∂t
=
∑

k

∑

ℓ

J
(0)
ijkℓ

∂σkℓ(t)

∂t
+
M∑

m=1

∂Λ
(m)
ij (t)

∂t
(3.10)

where Λ
(m)
ij (t), i,j ∈ {1,2,3}, m∈ {1, . . . ,M} are the solutions to the follo-

wing equations

Λ
(m)
ij (t)=

3∑

k=1

3∑

ℓ=1

Π
(m)
ijkℓ
(t) (3.11)

with
∂Π
(m)
ijkℓ
(t)

∂t
+λ
(m)
ijkℓ
Π
(m)
ijkℓ
(t)= J

(m)
ijkℓ
λ
(m)
ijkℓ
σkℓ(t) (3.12)

The non-aging linear viscoelastic behaviour is completely defined by differen-
tial equations (3.10), (3.11) and (3.12).
We note that this formulation, written in terms of strain and stress rates,

is easily adapted to temporal discretisation methods such as finite difference
ones.



402 C. Chazal, R. Moutou Pitti

4. Finite difference integration

Here we describe the solution process of a step-by-step nature in which the
loads are applied stepwise at various time intervals. Let us consider the time
step ∆tn = tn − tn−1. The subscript n− 1 and n refer to the values at the
beginning and end of the time step, respectively. We assume that the time
derivative during each time increment is constant and is expressed by

∂ζij

∂t
=
ζij(tn)− ζij(tn−1)

∆tn
=
∆(ζij)n
∆tn

(4.1)

where ζij represent strains or stresses. The following expressions can be then
written

∂Λ
(m)
ij (tn)

∂t
=
Λ
(m)
ij (tn+1)−Λ

(m)
ij (tn)

∆tn
=
∆Λ
(m)
ij (tn)

∆tn
(4.2)

∂σij(tn)

∂t
=
σij(tn+1)−σij(tn)

∆tn
=
∆σij(tn)

∆tn

A linear approximation is used for stresses, and is expressed by

σkℓ(τ)=σkℓ(tn)+
τ− tn
∆tn

[σkℓ(tn+1)−σkℓ(tn)]H(τ − tn) (4.3)

By integrating equation (3.10) between tn and tn+1, it can be written in the
form

∆εij(tn)=
∑

k

∑

ℓ

J
(0)
ijkℓ
∆σkℓ(tn)+

M∑

m=1

∆Λ
(m)
ij (tn) (4.4)

In order to determine the strain increments from this equation, we have to

determine the strain increments ∆Λ
(m)
ij (tn).

First, let us consider differential equation (3.12). The analytical solution
to this differential equation can be expressed as

Π
(m)
ijkℓ
(tn+1)−Π

(m)
ijkℓ
(tn)=

(
e
−λ
(m)
ijkℓ
∆tn −1

)
Π
(m)
ijkℓ
(tn)+

(4.5)

+J
(m)
ijkℓ

{(
1− e

−λ
(m)
ijkℓ
∆tn
)
σkℓ(tn)+∆σkℓ(tn)

[
1−

1

∆tnλ
(m)
ijkℓ

(
1−e

−λ
(m)
ijkℓ
∆tn
)]}

Consequently, when we substitute equation (4.5) into equation (3.11), we ob-

tain the strain increments ∆Λ
(m)
ij (tn)

M∑

m=1

∆Λ
(m)
ij (tn)=

M∑

m=1

3∑

k=1

3∑

ℓ=1

[Π
(m)
ijkℓ
(tn+1)−Π

(m)
ijkℓ
(tn)] (4.6)



A new incremental formulation for linear... 403

or

M∑

m=1

∆Λ
(m)
ij =

M∑

m=1

3∑

k=1

3∑

ℓ=1

(
e
−λ
(m)
ijkℓ
∆tn −1

)
Π
(m)
ijkℓ
(tn)+

(4.7)

+J
(m)
ijkℓ

{(
1− e

−λ
(m)
ijkℓ
∆tn
)
σkℓ(tn)+∆σkℓ(tn)

[
1−

1

∆tnλ
(m)
ijkℓ

(
1−e

−λ
(m)
ijkℓ
∆tn
)]}

5. Incremental viscoelastic constitutive equations

In this section, the incremental constitutive equations can now be obtained
from equation (4.4). Substituting equation (4.7) into (4.4), we find

∆εij(tn)=
∑

k

∑

ℓ

Dijkℓ(∆tn)∆σkℓ(tn)+ ε̃ij(tn) (5.1)

where Dijkℓ(∆tn) is a fourth order tensor which can be interpreted as a com-
pliance tensor, it is given by

Dijkℓ(∆tn)= J
(0)
ijkℓ
+
M∑

m=1

J
(m)
ijkℓ

[
1−

1

∆tnλ
(m)
ijkℓ

(
1− e

−λ
(m)
ijkℓ
∆tn
)]

(5.2)

and ε̃ij(tn) is a pseudo-strain tensor which represents the influence of the
complete past history of stresses. It is given by

ε̃ij(tn)=−
3∑

k=1

3∑

ℓ=1

M∑

m=1

(
1− e

−λ
(m)
ijkℓ
∆tn
)
Π
(m)
ijkℓ
(tn)+

(5.3)

+
3∑

k=1

3∑

ℓ=1

σkℓ(tn)
[ M∑

m=1

J
(m)
ijkℓ

(
1− e

−λ
(m)
ijkℓ
∆tn
)]

Finally, the incremental constitutive law given by equation (5.1) can now be
inverted to obtain

∆σij(tn)=
∑

k

∑

ℓ

Cijkℓ(∆tn)∆εkℓ(tn)− σ̃ij(tn) (5.4)

where Cijkℓ = (Dijkℓ)
−1 is the inverse of the compliance tensor and σ̃ij(tn)

is a pseudo-stress tensor which represents the influence of the complete past
history of strain. It is given by



404 C. Chazal, R. Moutou Pitti

σ̃ij(tn)=
3∑

k=1

3∑

ℓ=1

Cijkℓ(∆tn)ε̃ij(tn) (5.5)

The incremental constitutive law represented by equation (5.4) can be intro-
duced in a finite element discretisation in order to obtain solutions to complex
viscoelastic problems.

6. Numerical results

The finite element computation is compared with an analytical solution. The
incremental constitutive viscoelastic law given by equation (5.4) is implemen-
ted in Finite Element software CASTEM (Charvet-Quemin et al., 1986). In
order to validate our method, we employ a timber specimen of side 50mm.
The crack length chosen is 25mm.The external load is a unit loading applied
to steel Arcan as seen in Fig.1 (Moutou Pitti et al., 2008).

Fig. 1. CTS specimen (Moutou Pitti et al., 2008)

This specimen has similar properties of CTS (Compact Tension Shear)
specimens used by Zhang et al. (2006), Ma et al. (2006) and developed by
Richard and Benitz (1983) for an isotropic material. The points Aα and Bα
with α∈ (1, . . . ,7) are holes where unspecified forces can be applied with the
angle β directed according to the crack in the trigonometrical direction. Pure
mode I (opening mode) is obtained by using opposite forces in A1 and B1
with β=0◦. A loading with β=90◦, in A7 and B7, corresponds to the case



A new incremental formulation for linear... 405

of pure mode II (shear mode). Intermediary positions induce different mixed
mode configurations. The timber element is framed with perfectly rigid steel
Arcan.

In order to simplify the analytic development, a time proportionality for
the creep tensor is chosen

J(t)=
1

E(t)
C0 (6.1)

in which C0 is a constant compliance tensor composed by a unit elastic mo-
dulus and a constant Poisson’s coefficient ν = 0.4, and E(t) represents the
tangentmodulus for the longitudinal direction. In this context, the creep pro-
perties are given in terms of the creep function as given in equation (2.1)

1

E(t)
=
1

EX

[
1+

1

74.3

(
1− e−

74.3
3.37
t
)
+
1

74.4

(
1− e−

74.4
33.37

t
)
+

+
1

22.9

(
1− e−

22.9
104.09

t
)
+
1

27.6

(
1− e−

27.6
1251
t
)
+ (6.2)

+
1

7.83

(
1− e−

7.83
3554
t
)
+
1

3.23

(
1− e−

3.23
14660

t
)]

where EX is the longitudinal modulus and is equal to 15000MPa (Guitard,
1987). In this context, C0 admits the following definition for plane configura-
tions

C0 =




1 −ν 0

−ν
EX

EY
0

0 0
EX

GXY




(6.3)

where EY and GXY are the transverse and shear moduli, respectively. Their
values are fixed to: EY =600MPa and GXY =700MPa (elastic pine spruce
properties, Guitard, 1987).

In this test, the numerical results are compared to the analytical solu-
tion given by the isothermal Helmholtz free energy (Staverman and Schwarzl,
1952). According to the last creep tensor form, the viscoelastic compliance
takes the following form in puremode I and pure mode II, respectively

C1(t)=C
0
1f(t)=7.35 ·10

−3
f(t)

(6.4)

C2(t)=C
0
2f(t)=1.47 ·10

−3f(t)



406 C. Chazal, R. Moutou Pitti

in which C01 and C
0
2 are the reduced elastic compliances and

f(t)=
[
1+

1

74.3

(
1− e−

74.3
3.37
t
)
+
1

74.4

(
1− e−

74.4
33.37

t
)
+
1

22.9

(
1− e−

22.9
104.09

t
)
+

(6.5)

+
1

27.6

(
1− e−

27.6
1251
t
)
+
1

7.83

(
1− e−

7.83
3554
t
)
+
1

3.23

(
1− e−

3.23
14660

t
)]

In bi-dimensional analysis, we can express the energy release rate by the
expression

1
Gv(t)=

1

8
[2C1(t)−C1(2t)](

u
K
0
1)
2

(6.6)

2
Gv(t)=

1

8
[2C2(t)−C2(2t)](

u
K
0
2)
2

where uK01 and
uK02 are the instantaneous stress intensity factors for mode I

and mode II, respectively, computed with a classical elastic finite element
process. 1Gv and

2Gv are viscoelastic energy release rates in mode I and
mode II, respectively. Now, we present the comparison between the numerical
results, given by incremental formulation (5.4), and the analytical solution
given by expressions (6.6). The results are presented in Fig.2 and Fig.3 for
puremode I and puremode II versus time. The average error observed in the
numerical solution is less than 1% in puremodes I and II.

Fig. 2. Analytical and numerical solution in pure mode I for
energy release rate 1G

v

7. Conclusions

A new linear incremental formulation in the time domain for non-aging vi-
scoelastic materials undergoingmechanical deformation have been presented.
The formulation is based on a differential approach using a discrete spectrum



A new incremental formulation for linear... 407

Fig. 3. Analytical and numerical solution in pure mode II for
energy release rate 2G

v

representation for the creep tensor. Thegoverning equations are then obtained
using a discretised form of Boltzmann’s principle. The analytical solution to
differential equations is then obtained using a finite difference discretisation
in the time domain. In this way, the incremental constitutive equations for
linear viscoelastic material use a pseudo fourth order rigidity tensor. The in-
fluence of the whole past history on the behaviour at the time t is given by a
pseudo second order tensor. This formulation is introduced in a finite element
discretisation. The numerical results obtained are compared with the analy-
tical solution in terms of the energy release rate. The method can be easily
extended to deal with ageing boundary viscoelastic problems.

References

1. Boltzmann L., 1878, Zur Theorie der elastischen Nachwirkung Sitzungsber,
Mat. Naturwiss. Kl. Kaiser. Akad. Wiss., 70, 275

2. Charvet-Quemin F., Combescure A., Ebersol L., Charras T., Mil-
lard A., 1986, Méthode de calcul du taux de restitution de l’énergie en
élastique et en non linéaire matériau,Report DEMT, 86/438

3. Chazal C., Dubois F., 2001, A new incremental formulation in the time do-
main of crack initiation in an orthotropic linearly viscoelastic solid,Mechanics
of Time-Dependent Materials, 5, 229-253

4. Dubois F., Chazal C., Petit C., 1999a, A finite element analysis of creep-
crack growth in viscoelastic media, Mechanics of Time-Dependent Materials,
2, 269-286



408 C. Chazal, R. Moutou Pitti

5. DuboisF., ChazalC., PetitC., 1999b,Modelling of crackgrowth initiation
in a linear viscoelasticmaterial, Journal of Theoretical and Applied Mechanics,
37, 207-222

6. Dubois F., Chazal C., Petit C., 2002,Viscoelastic crack growth process in
wood timbers: An approach by the finite element method for mode I fracture,
International Journal of Fracture, 113, 367-388

7. Dubois F., Petit C., 2005,Modelling of the crack growth initiation in visco-
elastic media by the Gθ

v
-integral, Engineering Fracture Mechanics, 72, 2821-

2836

8. Ghazlan G., Caperaa S., Petit C., 1995, An incremental formulation for
the linear analysis of thin viscoelastic structures using generalized variables,
International Journal of Numeric Methods Engineering, 38, 3315-3333

9. Guitard D., 1987, Mécanique du matériau bois et composites, Edition
Cépaudes

10. KimK.S., SingLeeH., 2007,An incremental formulation for thepredictionof
two-dimensional fatigue crack growthwith curved paths, International Journal
for Numerical Methods in Engineering, 72, 697-721

11. KlasztornyM., 2008,Coupledanduncoupled constitutive equations of linear
elasticityandviscoelasticityof orthotropicmaterials,Journal ofTheoretical and
Applied Mechanics, 46, 1, 3-20

12. Ma S., Zhang X.B., Recho N., Li J., 2006, The mixed-mode investigation
of the fatigue crack inCTSmetallic specimen, International Journal of Fatigue,
28, 1780-1790

13. Mandel J., 1978, Dissipativité normale et variables cachés, Mechanics Rese-
arch Communications, 5, 4, 225-229

14. Moutou Pitti R., Dubois F., Petit C., Sauvat N., 2007, Mixed mode
fracture separation in viscoelastic orthotropic media: numerical and analytical
approach by theMθv-integral, International Journal of Fracture,145, 181-193

15. MoutouPitti R., Dubois F., Petit C., SauvatN., PopO., 2008,A new
M integral parameter for mixedmode crack growth in orthotropic viscoelastic
material,Engineering Fracture Mechanics, 75, 4450-4465

16. Richard H., Benitz K., 1983, A loading device for the creation of mixed
mode in fracture mechanics, International Journal of Fracture, 22, R55-58

17. StavermanA.J., SchwarzlP., 1952,Thermodynamics of viscoelastic beha-
vior,Proceeding Academic Science, 55, 474-492

18. Zhang X.B., Ma S., Recho N., Li J., 2006, Bifurcation and propagation of
amixed-mode crack in a ductilematerial,Engineering FractureMechanics, 73,
1925-1939



A new incremental formulation for linear... 409

Nowe sformułowanie przyrostowe w liniowej analizie lepkosprężystości:
różniczkowy opis pełzania

Streszczenie

Przedmiotem pracy jest prezentacja nowego przyrostowego opisu niestarzejących
się materiałów lepkosprężystych poddanych deformacji w dziedzinie czasu. Sformuło-
wanie wyprowadzono z równań różniczkowych opartych na dyskretnej reprezentacji
widma tensora pełzania. Następnie, przyrostowe równania konstytutywne otrzymano
w drodze całkowania różnicowego. W ten sposób uniknięto konieczności zachowy-
wania w pamięci komputera historii naprężenia. Kompletna i ogólna liniowa analiza
naprężeń lepkosprężystych została przedstawiona za pomocą przyrostów odkształceń
i naprężeń. Otrzymane wyniki symulacji numerycznych uzyskano z dobrą dokładno-
ścią. Analityczne i numeryczne rozwiązania porównano poprzez zestawienie tempa
uwalnianej energii dla czystej postaci I i II.

Manuscript received June 12, 2008; accepted for print November 18, 2008