Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

51, 3, pp. 569-580, Warsaw 2013

DETERMINATION OF MATERIAL PARAMETERS OF QUASI-LINEAR
VISCOELASTIC RHEOLOGICAL MODEL FOR THERMOPLASTICS

AND RESINS

Cyprian Suchocki, Marek Pawlikowski, Konstanty Skalski

Warsaw University of Technology, Institute of Mechanics and Printing, Warsaw, Poland

e-mail: c.suchocki@imik.wip.pw.edu.pl

Cezary Jasiński, Łukasz Morawiński

Warsaw University of Technology, Institute of Manufacturing Technology, Warsaw, Poland

In this study, an algorithm for identification of elastic and viscoelastic constants of a recen-
tly developed constitutive equation for thermoplastics and resins is presented. The equation
has been applied tomodel the viscoelastic response of ultra highmolecularweigth polyethy-
lene. In order to determine the material parameters, a series of rheological tests has been
performed. A set of equations describing one-dimensional processes has been derived and is
utilized to determine the material constants. Thematerial parameter identification leads to
a set of 9 constants, i.e. 3 constants of elasticity and 6 constants of viscoelasticity. For the
determined values of parameters, several validation tests have been performed.

Key words: rheology, constitutive equation, quasi-linear viscoelasticity

1. Introduction

Although the determination of constants of viscoelasticity is a problemwidely discussed in nu-
merous papers, e.g.Wilczyński (1968), Bradshaw andBrinson (1997), no universal identification
algorithmhas been proposed so far. There aremany opinions regarding the sort of experimental
data and the technique that should be used for material parameter identification.

Both in the case of linear and nonlinear viscoelasticity, the data utilized for determining
material constants are usually obtained from creep or relaxation experiments, cf Christensen
(1971). It is often assumed that the loading (force-induced or kinematic) is given by the Heavi-
side function, thus the loading ramp is simply not taken into account during the identification
of the material parameters. In linear theory, this assumption is admissible if the ramp time
is short compared to the time of the examined rheological process. For nonlinear viscoelastic
models, the described simplification may lead to incorrect values of the parameters describing
the instantaneous material response. Thus, in the case of nonlinear viscoelastic models, it is
usually recommended that a complete deformation history should be used for the parameter
identification, e.g. Goh et al. (2004).

An increasingly large number of researchers recommend determination ofmaterial constants
from experiments involving complex loading histories. For instance, Garbarski (1988) propo-
sed a rheological test which includes loading alternation from force-induced to kinematic. The
complete deformation paths were used for the parameter identification. Usually, complex load
histories are limited to the cases when a closed form solution of the hereditary integral can be
found. However, Goh et al. (2004) demonstrated that for the class of quasi-linear viscoelastic
(QLV) constitutive equations, replacing a hereditary integral with Taylor’s recurrence-update
numerical formula (Taylor et al., 1970) and using it in the least squares optimization procedure,



570 C. Suchocki et al.

does not influence the determined values of material constants seriously. Thus, for this particu-
lar class of viscoelastic materials, even the complex loading histories for which the hereditary
integral has no analytical solution are admissible for the parameter identification.
In the case of linear viscoelasticity, a mixed analytical-numerical algorithm for the para-

meter identification is often postulated, e.g. Wilczyński (1968). In this method, some of the
parameter values are calculated analytically utilizing the experimental data, whereas the re-
maining constants, usually the characteristic times, are determined numerically using the least
squares optimization procedure. The identification of material constants is usually a strongly
nonlinear optimization problem. This is the case even for the simplest generic function, i.e. a
series of exponentials, also called the Prony series. Thus, in order to simplify the optimization
problem, it is often postulated that the values of characteristic times should be predetermined,
e.g. Bradshaw and Brinson (1997), Laksari et al. (2012), Pawlikowski (2012), and should not
participate in the optimization. Usually, it is assumed that thematerial characteristic times are
distributed uniformly in a logarithmic scale of time of the experiment. The major drawback of
this approach is that the final set of material parameters does not include the optimal values of
the characteristic times.

Ciambella et al. (2010) proposedanalgorithm for thematerial parameter identificationwhere
the optimization process is launched iteratively with a predetermined set of characteristic times
being updated for each successive iteration according to the rule proposed by the authors. This
algorithm allows for a better estimation of the material characteristic times without involving
them in the least squares optimization.
In this paper, the concept presentedbyCiambella et al. (2010) is used to develop a computer-

-aided material parameter identification algorithm utilizing Taylor’s recurrence-update formula
for discretization of the constitutive equation.This combination enables the algorithm tobeused
for arbitrary deformation histories. AMatlab programhas beenwritten based on the developed
algorithm. It is applied to the identification of material parameters of the recently developed
QLV rheological model (see Suchocki, 2013) corresponding to themechanical properties of ultra
high molecular weigth polyethylene (UHMWPE). An appropriate set of equations describing
one-dimensional processes has been derived from theQLV constitutive equation. The developed
algorithm is used to determine two separate sets ofmaterial constants. The first set of constants
is determined by utilizing test data as recommended in the literature (e.g. Goh et al., 2004).
The second set of parameters is determined by using a combination of relaxation and cyclic test
data. The obtained results are compared.

2. Quasi-linear viscoelastic constitutive model

This point summarizes the basic notions of the QLV model for polymeric materials developed
by Suchocki (2013).

The following assumptions are adopted:

• The polymer is a homogeneous and isotropic, nonlinear viscoelastic body.

• The considered deformation processes are isothermal.

• The deformation is restricted to the processes where no plastic strains occur or the plastic
strains are negligible.

• The nonlinear viscoelastic material properties are described by amodification of the Bolt-
zmann superposition principle, cf Fung (1981).

• The instantaneous material elastic response follows from the modified Knowles stored-
-energy potential, cf Knowles (1977), Suchocki (2011).

• Thementioned Prony series is used as the relaxation function.



Determination of material parameters of quasi-linear viscoelastic rheological model ... 571

The formulation of the constitutive equation involves uncoupling of the volumetric and de-
viatoric stress tensor components.
The total stress is given by the following equation (see Suchocki, 2013)

S= g∞S
e+

N∑

j=1

Hj (2.1)

where

Se = J
∂U

∂J
C−1+2J−

2

3 DEV
[∂W
∂C

]

Ḣj +
1

τj
Hj = gjṠ

e j=1,2, . . . ,N

DEV[•] = [•]−
1

3
([•] ·C)C

−1

(2.2)

The utilized symbols have the followingmeaning: S – second Piola-Kirchhoff (P-K) total stress
tensor, Se – instantaneous second P-K elastic stress tensor, Hj – viscoelastic overstress tensor
(j=1,2, . . . ,N), U –volumetric stored-energypotential, W – isochoric stored-energypotential,
C – right Cauchy-Green (C-G) deformation tensor, C – isochoric C-G deformation tensor,
J – Jacobian determinant, τj – relaxation time (j = 1,2, . . . ,N), g∞ and gj – relaxation
coefficients (j=1,2, . . . ,N), N – assumed number of overstress tensors. The position and time
arguments are omitted for convenience.
After time-integrating, Eq. (2.2)2 takes the form of a functional

Hj(t)=

t∫

0

gje
−
t−τ
τj Ṡe dτ j=1,2, . . . ,N (2.3)

Thus, by substituting Eq. (2.3) into Eq. (2.1), one can obtain

S(t)=

t∫

0

G(t− τ)
∂Se[E(τ)]

∂τ
dτ (2.4)

where

G(t) = g∞+
N∑

j=1

gje
−
t
τj (2.5)

is the reduced relaxation function, cf Fung (1981).
Following the postulate mentioned above, the isochoric stored-energy potential W(C) is

taken in the form of theKnowlesmodel, whereas the volumetric stored-energy U(J) is assumed
in the simplest possible form, i.e.

U(J)=
1

D1
(J−1)2 W(C)=

µ

2b

{[
1+
b

κ
(Ī1−3)

]κ
−1
}

(2.6)

Inserting Eqs (2.6) into Eq. (2.2)1, yields

Se =
2

D1
J(J−1)C−1+µ

[
1+
b

κ
(Ī1−3)

]κ−1(
J−

2

31−
1

3
Ī1C

−1
)

(2.7)

The symbols used inEqs (2.6) and (2.7) denote: µ – shearmodulus, κ – “hardening” parameter,
b – auxiliary constant, D1 – inverse of the bulk modulus, Ī1 – first invariant of the isochoric
right C-G tensor, 1 – identity tensor.



572 C. Suchocki et al.

Fig. 1. Mechanical scheme of the rheological model

Summing up, the developed model is described by four elasticity constants and N + 1
viscoelasticity constants, of which N are independent.

It follows from Eqs (2.1) and (2.2)2 that the presented constitutive equation has an inter-
pretation of the generalized Maxwell model (Fig. 1).

Equations (2.1) and (2.3) can be discretized by taking advantage of the algorithm developed
by Taylor et al. (1970), i.e.

Sn+1 = g∞S
e
n+1+

N∑

j=1

Hjn+1

Hjn+1 = e
−
∆t
τj Hj n+gj

1−e
−
∆t
τj

∆t
τj

(
Sen+1−S

e
n

)
(2.8)

The recurrence-update formula requires stresses at the time increment n to be stored inmemory
for computation of stresses in the successive time increment n+1.

3. Algorithm for identification of elasticity and viscoelasticity constants

In order to develop an algorithm for determination of material parameters, the neccessary equ-
ations describing the selected deformation process have to be derived. The uniaxial compression
experiments are the easiest to perform bymeans of the standardmaterial testing toolset. Thus,
an appropriate set of scalar equations has beenderived,which describes a one-dimensional stress
state and a three-dimensional strain state in a QLV body subjected to a compressive loading.
The obtained equations are valid for arbitrary uniaxial stress loading histories.

By using Ogden’s (1997) approach, the most general equation describing a homogeneous
deformation of a continuum takes the form

x=FX+c (3.1)

where the following symbols are used: x – position vector of a particle in the current configura-
tion, X – position vector of a particle in the reference configuration, F – deformation gradient,
c – vector representing rigid translation of the whole body.

In the case of a cylindrical specimen subjected to axial compression it is convenient to use
one Cartesian basis for analyzing both reference and current configurations. It is assumed that
one axis is alignedwith the axis of the specimenand the other pair lying in aplane perpendicular
to the specimen axis, so as to form a rectangular coordinate system.



Determination of material parameters of quasi-linear viscoelastic rheological model ... 573

Since the specimendoes not undergo rigid translationalmotion, the vector c is equal to zero.
Thus, Eq. (3.1) has the following component form

x1 =λ1X1 x2 =λ2X2 x3 =λ3X3 (3.2)

where λk (k=1,2,3) are the principal stretches.
The component matrix of the deformation gradient tensor takes the form

F3×3 =



λ1 0 0
0 λ2 0
0 0 λ3


 (3.3)

where the stretch ratio λ1 is chosen to correspond to the axial direction. Thus λ1 = l/l0, where
l0 and l denote initial height of the specimen and height at a certain time instant, accordingly.
Equation (2.4), related to the axial compression of a cylindrical specimen, gives the following

relations for the second P-K stress components

S11(t)=

t∫

0

G(t− τ)
∂Se11(τ)

∂τ
dτ S22(t)=

t∫

0

G(t− τ)
∂Se22(τ)

∂τ
dτ

S33(t)=

t∫

0

G(t− τ)
∂Se33(τ)

∂τ
dτ S12(t)=S23(t)=S31(t)= 0

(3.4)

where the relations for instantaneous elastic stress components follow from Eq. (2.7) and the
diagonal form of the deformation gradient, i.e.

Se11 =
2

D1
J(J−1)C−111 +µ

[
1+
b

κ
(Ī1−3)

]κ−1(
J−

2

3 −
1

3
Ī1C

−1
11

)

Se22 =
2

D1
J(J−1)C−122 +µ

[
1+
b

κ
(Ī1−3)

]κ−1(
J−

2

3 −
1

3
Ī1C

−1
22

)

Se33 =
2

D1
J(J−1)C−133 +µ

[
1+
b

κ
(Ī1−3)

]κ−1(
J−

2

3 −
1

3
Ī1C

−1
33

)

Se12 =S
e
23 =S

e
31 =0

(3.5)

with

Ī1 = J
−
2

3

(
λ21+λ

2
2+λ

2
3

)
J =λ1λ2λ3 (3.6)

For convenience, the time and position arguments are omitted.
By inserting into Eqs (3.5), the right C-G deformation tensor components expressed by

means of the principal stretches, i.e.

C−111 =
1

λ21
C−122 =

1

λ22
C−133 =

1

λ23
(3.7)

and substituting p=2(J−1)/D1, one obtains

Se11 = Jp
1

λ21
+µ
[
1+
b

κ
(Ī1−3)

]κ−1(
J−

2

3 −
1

3
Ī1
1

λ21

)

Se22 = Jp
1

λ22
+µ
[
1+
b

κ
(Ī1−3)

]κ−1(
J−

2

3 −
1

3
Ī1
1

λ22

)

Se33 = Jp
1

λ23
+µ
[
1+
b

κ
(Ī1−3)

]κ−1(
J−

2

3 −
1

3
Ī1
1

λ23

)
(3.8)



574 C. Suchocki et al.

Due to the symmetry of the specimen and material isotropy, the principal stretches in the
directions perpendicular to the axis are equal, i.e., λ2 = λ3 = λ. The boundary condition of
the stress-free lateral surface requires S22 =S33 =0. In addition, the reference configuration is
assumed tobe stress-free. It follows fromEqs (3.4)2 and (3.4)3 that these conditions are satisfied
only if the instantenous stresses Se22 =S

e
33 =0.Hence, Eqs (3.8)2 and (3.8)3 are equivalent and

allow for determining the term Jp, i.e.

Jp=−µ
[
1+
b

κ
(Ī1−3)

]κ−1(
J−

2

3λ2−
1

3
Ī1
)

(3.9)

where

Ī1 = J
−
2

3(λ21+2λ
2) J =λ1λ

2 (3.10)

Substituting Eq. (3.9) into Eq. (3.8)1 yields

Se11 = J
−
2

3µ
[
1+
b

κ
(Ī1−3)

]κ−1[
1−
( λ
λ1

)2]
(3.11)

which is valid for both compressible and incompressible materials.
In the case of material incompressibility, it follows from the constraint J = 1, cf Holzapfel

(2010), that

λ2 =
1

λ1
Ī1 = I1 =

(
λ21+2

1

λ1

)
(3.12)

Thus, by inserting Eqs (3.12) into Eq. (3.11) one finally obtains

Se11 =µ
[
1+
b

κ
(I1−3)

]κ−1(
1−
1

λ31

)
(3.13)

The Taylor numerical integration algorithm (see Taylor et al., 1970) is utilized to discretize
Eqs (3.4)1 and (3.13), thus

S11(tn+1)= g∞S
e
11(tn+1)+

N∑

j=1

Hj(tn+1)

Hj(tn+1)= e
−
∆t
τj Hj(tn)+gj

1−e
−
∆t
τj

∆t
τj

[Se11(tn+1)−S
e
11(tn)]

Se11(tn+1)=µ
[
1+
b

κ

(
I1(tn+1)−3

)]κ−1(
1−

1

λ31(tn+1)

)

T11(tn+1)=λ1(tn+1)S11(tn+1)

(3.14)

where I1(tn+1) = λ
2
1(tn+1)+2/λ1(tn+1), T11 is the first P-K stress, and the axial components

of the overstress tensors are denoted as Hj (j=1,2, . . . ,N).
Adiscretized time history of λ1 and corresponding ∆t are the input arguments of Eqs (3.14)

allowing one to determine the material response. In each time increment, Eq. (3.14)3 is used to
calculate the instantenous elastic stress, subsequently Eq. (3.14)2 for j =1,2, . . . ,N is used to
update the viscoelastic overstresses. The total stress is updated according to Eq. (3.14)1.
Equations (3.14) are utilized within the least squares optimization procedure to calculate

the theoretical values of stresses for a given deformation process. The material parameters are
found byminimization of the following stress error function

F(p) =
M∑

k=1

[(
T11(p)

)
k
− (T̃11)k

]2
(3.15)



Determination of material parameters of quasi-linear viscoelastic rheological model ... 575

where T11(p) is the theoretical and T̃11 is the experimentalLagrange stress at the time instant tk
(k=1,2, . . . ,M), M being the number of collocation points, and

p= [µ,b,κ,g1,g2, . . . ,gN,g∞]
T (3.16)

is a vector of thematerial parameters involved in the least squares optimization. The subsequent
constraints must be satisfied

µ­ 0 b­ 0 κ­ 0

0¬ gj ¬ 1 0¬ g∞ ¬ 1 j=1,2, . . . ,N
(3.17)

with

g∞+
N∑

j=1

gj =1 (3.18)

The reduced relaxation function given by Eq. (2.5) must bemonotonically decreasing, thus

τj ­ 0 j=1,2, . . . ,N (3.19)

The concept presented in Ciambella et al. (2010) is utilized to define a recurrence-update
scheme in which the error function given by Eq. (3.15) is iteratively minimized for a predeter-
mined set of relaxation times τj (j = 1,2, . . . ,N). The relaxation times are updated for each
successive iteration according to a rule specified below.

The initial values of the relaxation times, used in the iteration i = 1, are distributed uni-
formly in a logarithmic scale within the range 〈τmin,τmax〉, i.e.

τ
(i=1)
1 = τmin, . . . ,τ

(i=1)
j = τmin10

(j−1)∆, . . . ,τ
(i=1)
N = τmax (3.20)

where τmin is set equal to the sampling interval, τmax is equal to the total timeof the experiment
and ∆ is the logarithimic time step.

The logarithmic step and the initial values of the viscoelastic coefficients are defined as
follows

∆=
1

N−1
log
τmax
τmin

g
(i=1)
j = g

(i=1)
∞
=
1

N+1
j=1,2, . . . ,N (3.21)

After minimization of F(p) at a least squares optimization step, all the viscoelastic coeffi-
cients gj (j =1,2, . . . ,N) are checked upon overcoming a predefined threshold. If a coefficient
is below the threshold, then both the given coefficient and the corresponding relaxation time
are discarded. If the viscoelastic coefficient overcomes the threshold, then the corresponding
relaxation time is split according to the rule

τ
(i)
j ⇒ τ

(i)
j

{
10
−

(
∆
3

)i

,1,10

(
∆
3

)i
}

j=1,2, . . . ,N (3.22)

where i denotes the number of the optimization step.

The initial values of the viscoelastic coefficients associated with the new relaxation times
are set to zero. Subsequently, the least squares optimization procedure is launched again. This
process is repeated iteratively in a loop until the stop condition is satisfied.



576 C. Suchocki et al.

4. Experimental set-up and rheological tests

Thematerial chosen for the rheological tests wasCHIRULEN1050medical ultra highmolecular
weight polyethylene (UHMWPE). UHMWPE is a popular biomaterial with numerous applica-
tions as hip joint implants, elbow implants or artificial intervertebral discs for instance. The
usual role of the elementsmade of UHMWPE is dissipating themechanical energy, thusmaking
an implant less vulnerable to a failure. The nonlinear viscoelastic properties of UHMWPE are
representative for a large number of polymers, aspecially thermoplastics and resins.

Cylindrical specimens of dimensions ∅17 and l0 =21 mm were machined from a rod made
of CHIRULEN 1050. The tests were performed on MTS Bionix 2500 testing machine at the
constant room temperature of 20◦C. Two types of tests were performed, i.e. loading/unloading
compression experiments and relaxation tests.The lowest deformation rate usedwas 0.00023s−1

whereas the maximumwas 0.03s−1. During the relaxation tests the specimens were allowed to
relax for 30 minutes.

Thestrainmeasurementmethodutilizes avideoextensometer.Theelements of themeasuring
system are a CCD camera (Charge Coupled Device) with a resolution of 1600×1200 pixels, a
camera lens with focal length of 35mm. The axial deformation was additionally measured with
a strain gauge extensometer in order to verify themeasured values. A sketch of the experimental
setup can be seen in Fig. 2.

Fig. 2. (a) Experimental set-up, (b) view of undeformed specimen, (c) view of deformed specimen

Fig. 3. Exemplary deformation histories recorded by video extensometer; (a) axial stretch ratio versus
time, (b) transverse stretch ratio versus time

In Fig. 3, the principal stretches recorded by the video extensometer during a compression
loading/unloading experiment are presented. The experimental values of the principal stretch
in the transverse direction were compared to the theoretical values calculated fromEq. (3.12)1.
The comparison shows good agreement, thus the assumption of incompressibility is justified in
the case of UHMWPE.



Determination of material parameters of quasi-linear viscoelastic rheological model ... 577

5. Identification of material parameters for uniaxial compression and relaxation
tests data

Theprocedure proposed byGoh et al. (2004) was applied to determine the parameters of incom-
pressible QLVmodel for UHMWPE. The initial values of the constants were evaluated using a
Matlab program based on the algorithm discussed in Section 3. The relaxation test data were
utilized. Subsequently, the elasticity constants were recalculated to fit the test data from com-
pression experiments performed at different deformation rates. The results of curve fitting can
be seen in Fig. 4.

Fig. 4. Fitting of quasi-linear viscoelastic material model to experimental data; (a) compression ramp
tests with deformation rate of: 0.0003, 0.003 and 0.03s−1, (b) stress relaxation test

The identification led to a set of 15 independent parameters, i.e., 3 elasticity constants and
12 viscoelasticity constants. The parameters have been gathered in Table 1.

Table 1. Incompressible QLVmaterial parameters of UHMWPE

Elastic Viscoelastic Relaxation
constants coefficients time [s]

µ=296.13MPa g1 =0.096 τ1 =0.49
b=314.85 g2 =0.12 τ2 =0.97
κ=0.67 g3 =0.034 τ3 =3.73

g4 =0.10 τ4 =7.34
g5 =0.16 τ5 =55.81
g6 =0.15 τ6 =834.49

The constant D1 did not participate in the optimization. In FE implementation, it was
assumed that D1 = 33E-9MPa

−1, which corresponds to the bulk modulus high enough to
account for material near incompressibility.

In order to verify the performance of the model, a complex multistep rheological test was
performed. The comparison of the experimental results and model predictions can be seen in
Fig. 5.

6. Identification of material parameters for uniaxial compressive
loading/unloading and relaxation tests data

In the second approach, the material parameters were evaluated from compressive lo-
ading/unloading and relaxation test data.
Firstly, the algorithm discussed in Section 3 was utilized in curve-fitting of the relaxation

data. Subsequently, both elasticity constants and viscoelastic coefficients were recalculated by
fitting the loading/unloading data. The results of the optimization are shown in Fig. 6.



578 C. Suchocki et al.

Fig. 5. Model validation; (a) experimental deformation history, (b) material response, comparison of
experimental data andmodel predictions

Fig. 6. Fitting of quasi-linear viscoelastic material model to experimental data. Cyclic compression tests
with deformation rate of: (a) 0.00023, (b) 0.00038 and (c) 0.0005 s−1, (d) stress relaxation test

The identification led to a set of 9 independent parameters, i.e., 3 elasticity constants and
6 viscoelasticity constants. The parameters have been gathered in Table 2. Again, they are
supplemented by D1 =33E-9MPa

−1.

Table 2. Incompressible QLVmaterial parameters of UHMWPE

Elastic Viscoelastic Relaxation
constants coefficients time [s]

µ=309.85MPa g1 =0.14 τ1 =3.73
b=247.93 g2 =0.19 τ2 =7.34
κ=0.79 g3 =0.41 τ3 =55.81

In order to verify the performance of the model, several validation tests were performed
including a complex, multistep rheological test and compression loading/unloading experiments
withvarieddeformation rates.Thecomparison of the experimental results andmodel predictions
can be seen in Fig. 7 and Fig. 8.



Determination of material parameters of quasi-linear viscoelastic rheological model ... 579

Fig. 7. Validation of quasi-linear viscoelastic material model. Experimental data andmodel predictions
for cyclic compression tests with deformation rate of: (a) 0.00028, (b) 0.00046 and (c) 0.0022s−1

Fig. 8. Model validation; (a) experimental deformation history, (b) material response, comparison of
experimental data andmodel predictions

7. Conclusions

In the paper, an algorithm allowing one to identify the material constants of the QLV consti-
tutive equation is presented. The algorithm combines the concept of iterative determination of
relaxation timeswithdiscretization ofQLVmodel byTaylor’s recurrence-update formula (Taylor
et al., 1970). This enables evaluating the constants from arbitrary deformation histories.

Two different experimental data sets were used for the curve-fitting. In the first approach,
the data from compression tests performed at different deformation rates were used along with
the relaxation test data, as proposed in Goh et al. (2004). For the determined set of constants,
themodel predictions are in a good agreementwith the experimental data, whichwas confirmed
by the complex validation test (see Figs 4 and 5). However the number of material parameters
is large (15 independent constants). What is more, it was found that for loading/unloading
simulations the error of predicted energy dissipation is significant. Thus in the second approach
the loading/unloading data were used for parameter evaluation. The spectrum of deformation



580 C. Suchocki et al.

rates was narrowed. This resulted in the number of constants decreasing to 9 and obtaining
much better prediction of the energy dissipation, which was confirmed by the validation tests
(see Figs 6, 7 and 8).

References

1. Bradshaw R.D., Brinson L.C., 1997, A sign control method for fitting and interconverting
material functions for linearly viscoelastic solids, Mechanics of Time-Dependent Materials, 1,
85-108

2. Christensen R.M., 1971,Theory of Viscoelasticity, Academic Press

3. Ciambella J., Paolone A., Vidoli S., 2010, A comparison of nonlinear integral-based visco-
elastic models through compression tests on filled rubber,Mechanics of Materials, 42, 932-944

4. Fung Y.C., 1981,Biomechanics: Mechanical Properties of Living Tissues, Springer-Verlag

5. Garbarski J., 1988, Application of the exponential function to the description of viscoelasticity
in some solid polymers, International Journal of Mechanical Sciences, 3, 165-178

6. Goh S.M., Charalambides M.N., Williams J.G., 2004, Determination of the constitutive
constants of non-linear viscoelasticmaterials,Mechanics of Time-Dependent Materials, 8, 255-268

7. Holzapfel G.A., 2010,Nonlinear Solid Mechanics, JohnWiley & Sons Ltd., NewYork

8. Knowles J.K., 1977, The finite anti-plane shear field near the tip of a crack for a class of incom-
pressible elastic solids, International Journal of Fracture, 13, 611-639

9. Laksari K., Shafieian M., Darvish K., 2012, Constitutive model for brain tissue under finite
compression, Journal of Biomechanics, 45, 642-346

10. Ogden R.W., 1997, Non-Linear Elastic Deformations, Dover Publications, Inc., Mineola, New
York

11. Pawlikowski M., 2012, Cortical bone tissue viscoelastic properties and its constitutive equation
– preliminary studies,The Archive of Mechanical Engineering, 59, 31-52

12. Suchocki C., 2011, A finite element implementation of knowles stored-energy function: theory,
coding and applications,The Archive of Mechanical Engineering, 58, 319-346

13. Suchocki C., 2013, A quasi-linear viscoelastic rheological model for thermoplastics and resins,
Journal of Theoretical and Applied Mechanics, 51, 1, 117-129

14. Taylor R.L., Pister K.S., Goudreau G.L., 1970, Thermomechanical analysis of viscoelastic
solids, International Journal for Numerical Methods in Engineering, 2, 45-59

15. WilczyńskiA., 1968, Selected problems of testing ofmechanical properties of linearly viscoelastic
solids, Scientific Surveys Warsaw University of Technology, Mechanics, 1 [in Polish]

Manuscript received June 28, 2012; accepted for print September 27, 2012