Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

53, 3, pp. 593-604, Warsaw 2015
DOI: 10.15632/jtam-pl.53.3.593

AN INTERNAL-STATE-VARIABLE BASED VISCOELASTIC-PLASTIC

MODEL FOR POLYMERS

Cyprian Suchocki

Warsaw University of Technology, Department of Mechanics and Armament Technology, Warsaw, Poland

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

In this study, a new viscoelastic-plastic constitutive model which has been formulated by
utilizing the formalism of stress-like internal state variables is introduced. The developed
constitutive equation allows for a good description of the inelastic material response of
polymericmaterials over awide rangeof strain rates.An algorithm for numerical integration
of themodel equationshasbeenderived.TheFE implementationof the constitutive equation
is widely discussed and the results of solving several exemplary problems are presented.

Keywords: polymers, rheology, constitutive equation, finite element method

Nomenclature

b – Knowles material parameter
B,B – left Cauchy-Green (C-G) and isochoric left C-G deformation tensor
C,C – right and isochoric right C-G deformation tensor
C,Cve−p – elasticity tensor and viscoelastic-plastic material tensor
C
τc – material tensor related to convected stress rate

C
MJ – material Jacobian tensor used by Abaqus
D – strain rate tensor
D1 – inverse of bulkmodulus
ek – unit vector of Cartesian base, k=1,2,3,
F,F – deformation gradient and isochoric deformation gradient tensor
Hk – viscoelastic overstress tensor k=1,2, . . . ,N

H̃k – endochronic overstress tensor k=1,2, . . . ,P
I – fourth order identity tensor
IC−1 – fourth order identity tensor in reference configuration
I – auxiliary fourth order tensor
Ik,Ik – algebraic invariants of right and of isochoric right C-G deformation tensor,

k=1,2,3
Jk – algebraic invariants of material time derivative of isochoric right C-G defor-

mation tensor, k=1,2,3
J – Jacobian determinant
L – velocity gradient tensor
Q – orthogonal tensor
S,S – second Piola-Kirchhoff (P-K) and auxiliary second P-K total stress tensor
S0 – second P-K elastic stress tensor
T – Lagrange total stress tensor
t – time since loading
U – volumetric stored elastic energy potential



594 C. Suchocki

W – spin tensor
W,W – stored and isochoric stored elastic energy potential
z,ζ – fictitious time variable at time t and at time τ
τ,τ ′ – time and auxiliary time variable
λk – stretch ratio in k-th direction, k=1,2,3
τ,σ – Kirchhoff and Cauchy stress tensors
µ,κ – Knowles material shear modulus and Knowles material hardening parameter
Γk – k-th relaxation coefficient, k=1,2, . . . ,N
γk – k-th endochronic coefficient, k=1,2, . . . ,P
τk – k-th relaxation time, k=1,2, . . . ,N

D̃k – k-th endochronic parameter, k=1,2, . . . ,N
1 – second order identity tensor
DEV[•] – operator extracting deviatoric part of a tensor in reference configuration
⊗ – dyadic product operator
(•)∇ – Zaremba-Jaumann (Z-J) objective rate operator

1. Introduction

In recent years, various constitutive equations of viscoelasticity and viscoplasticity have been
developed in order to describe the mechanical properties of thermoplastic polymers and resins.
Ayoub et al. (2011) presented a viscoplastic model aimed at capturing the large strain response
of polyethylene over a wide range of strain rates. A better approximation of the stress-strain
relation of polyethylene, however for a narrower range of the strain rates, was achieved by
Bergström and Bischoff (2010) who developed the so-called three network viscoplastic model.
Abdul Hameed et al. (2014) formulated a constitutive model based on the multiplicative split
of the deformation gradient. It was further used to simulate the monotonic stress-strain curves
of polyethylene. Bardella (2001) developed a nonlinear viscoelastic constitutive model oriented
for the description of the epoxy resin behavior during cyclic loadings and creep.

Instead ofmodeling themechanical behavior of a polymer over the entire range of physically
possible strains, some researchers focus on improving the constitutive equations for the spectrum
of small and moderate deformations which are usually experienced by the structural elements.
Drozdov and Christiansen (2007) developed a viscoplastic constitutive model aimed at predic-
ting the stresses in highdensity polyethylene (HDPE) subjected to the cyclic loadings in tension.
Hassan et al. (2011) utilized an overstress-based viscoplasticity theory to model the hysteresis
loops of ultra high molecular weigth polyethylene (UHMWPE) for the maximum strains equal
to 10%.BenHadjHamouda et al. (2007) developed a viscoplasticmodel to capture the nonlinear
stress-strain relation, creep and relaxation of the medium density polyethylene (MDPE). Mo-
notonic stress-strain curves with themaximum strain of 12% were analyzed. Krairi and Doghri
(2014) proposed a model utilizing the additive split of the strain tensor into viscoelastic and
viscoplastic components. Monotonic and cyclic stress-strain curves were approximated for both
polyamide (PA) and HDPE with themaximum axial strain of 14%.

Constitutive theories based on the endochronic plasticity theory enjoy recently an increasin-
gly high interest. This theory was originally developed by Pipkin and Rivlin (1965). Kästner et
al. (2012) proposed a viscoplastic constitutive equation for polyprophylene (PP) based on the
generalized Maxwell mechanistic model. The polymer viscosity was assumed to depend on the
non-equilibriumoverstress while endochronic plasticity was used to helpmodeling the hysteresis
loop. Strains up to 5% in tension were analyzed and infinitesimal strain tensor was used. This
approach was later generalized to large displacements and large rotations by Alkas Yonan et al.
(2013) for the purpose of modeling the inelastic response of polyvinylchloride (PVC).



An internal-state-variable based viscoelastic-plastic... 595

All thementioned above constitutive theories requiredetermination of about twentymaterial
parameters. The numerical algorithms used for integrating themodels are usually complicated.
In the current study, a new viscoelastic-plastic constitutive model for thermoplastics and

resins is presented. The model is an extension of the popular nonlinear viscoelasticity theory
proposed byHolzapfel (2010) andwidely used in FE codes. The formalism of internal stress-like
state variables previously used by Suchocki (2013) to develop a nonlinear viscoelastic model
for polymers is extended by taking advantage of the endochronic plasticity framework (Pipkin
andRivlin, 1965). The common assumptions of isothermal deformations and treating the effects
of volumetric viscoelasticity as negligible have been applied. It was experimentally observed
that irreversible volume changes are usually small (e.g. Wu, 2005). Thus, the postulate of pla-
stic incompressibility is well established and utilized by the classical theory of plasticity. This
assumption is adopted in the constitutive model formulated in the present study for sake of
simplicity. The constitutive equation has a modular structure which enables adjusting it for a
specific polymeric material. It is applied to capture the inelastic behavior of polyethylene. For
that purpose the specific values of the material parameters have been determined. It is found
that the proposed constitutive model accurately reproduces the stress relaxation curve and the
hysteresis loop for the deformation ratemagnitude from the interval 0.0005-0.05s−1. The consti-
tutive equation is discretized and an algorithm for efficient numerical computations is presented.
The implementation of the model into the FE software Abaqus is discussed.

2. Decoupled finite elasticity

In this section, thebasicnotions of hyperelasticity arepresented.Amaterial is called hyperelastic
if it posseses a stored-energy function W = W(C) (Holzapfel, 2010). In the case of material
isotropy, this function must be invariant with respect to a rotation Q, i.e.

W(C)=W(QCQT) (2.1)

This requirement is fulfilled if W is a function of algebraic invariants ofC.
In terms of FEM, it is profitable to decouple the volumetric and isochoric responses within

the constitutivemodel. This is facilitated by utilizing propermultiplicative decomposition of the
deformation gradient tensor F

F=FvolF Fvol = J
1

31 F= J−
1

3F C=F
T
F= J−

2

3C (2.2)

where Fvol and F represent purely volumetric and isochoric deformations, respectively. Conse-
quently, the algebraic invariants ofC

I1 = trC I2 =
1

2
[(trC)2− trC

2
] I3 =detC=1 (2.3)

The stored-energy function is formulated in a decoupled form

W(C)=U(J)+W(C) (2.4)

whereU andW are the stored-energy components associated with the volumetric and isochoric
deformations, respectively. The stored-energy function as given by (2.4) results in

S= JpC−1+J−
2

3 DEV[S] (2.5)

where

p=
∂U

∂J
DEV[S] =S−

1

3
(S ·C)C

−1
S=2

∂W

∂C
(2.6)

and DEV[•] = [•]− 1
3
([•]·C)C

−1
is an operator extracting a deviator of a tensor in the reference

configuration, whereas S is an auxiliary stress tensor.



596 C. Suchocki

3. Internal-state-variable based viscoelastic-plastic model

The total stress is assumed in the form of a sum

S=S0+
P∑

k=1

H̃k+
N∑

j=1

Hj (3.1)

where the components S0 and H̃k (k=1,2, . . . ,P) describe the equilibrium material response,
whereas Hj (j = 1,2, . . . ,N) account for the viscoelastic effects. S0 is taken in the decoupled
form, i.e.

S0 =S
vol
0 +S

iso
0 S

vol
0 = J

∂U

∂J
C
−1

S
iso
0 =2J

−2
3 DEV

[∂W
∂C

]
(3.2)

The evolution of overstresses H̃k is governed by the equations of the form

˙̃
Hk+

1

D̃kM
(∣∣Ċ
∣∣)
H̃k = γkṠ

iso
0 k=1,2, . . . ,P (3.3)

where

M
(∣∣Ċ
∣∣)= J

−1
2

2 J2 = trĊ
2

(3.4)

This specific choice of M enables one to eliminate the time from Eq. (3.3), thus making the
evolution of H̃k rate-independent. By integrating Eq. (3.3), it is found that

H̃k(t)=

z(t)∫

0

γke
−z−ζ
D̃k
∂Siso0 (ζ)

∂ζ
dζ (3.5)

where the following equalities are used

ζ̇ =
1

M
(∣∣Ċ
∣∣)

ζ(τ)=

τ∫

−∞

dτ ′

M(τ ′)
ζ(t)= z (3.6)

Substituting Eqs (3.4) into Eq. (3.6)2 leads to the result

ζ(τ)=

τ∫

−∞

[
Ċ(τ ′) · Ċ(τ ′)

]1
2 dτ ′ ˙(•)=

d

dτ ′
(•) (3.7)

which in turn gives

ζ(τ)=

τ∫

−∞

[
dC(τ ′) ·dC(τ ′)

]1
2 (3.8)

The variable ζ is a fictious internal time associated with the material. For this reason the
overstresses defined by Eq. (3.5) are called endochronic. Equation (3.5) is a special case of
the stress-strain relation proposed by Pipkin and Rivlin (1965). The stress tensor S0 together
with the endochronic overstresses H̃k (k = 1,2, . . . ,P) form up the elasto-plastic part of the
constitutive equation which is responsible for reproducing the hysteresis loop of the modeled



An internal-state-variable based viscoelastic-plastic... 597

polymeric material. The strain rate effects are captured by the viscoelastic overstresses Hj
(j=1,2, . . . ,N) which evolve in time according to the following differential equation

Ḣj +
1

τj
Hj =ΓjṠ

iso
0 j=1,2, . . . ,N (3.9)

After time integrating Eq. (3.9), one obtains

Hj(t)=

t∫

−∞

Γje
−t−τ
τj
∂Siso0 (τ)

∂τ
dτ (3.10)

For the purpose of modeling the mechanical properties of thermoplastics and resins, the
stored-energy functionW is adopted in the form proposed by Knowles (1977). The advantages
of using the isochoricKnowles functionwere discussedbySuchocki (2011). Thevolumetric strain
energy is assumed in the standard form, thus

W =
µ

2b

{[
1+
b

κ
(I1−3)

]κ
−1
}

U =
1

D1
(J−1)2 (3.11)

Summing up, the developed constitutive equation utilizes four elasticity constants, 2P en-
dochronic plasticity constants and 2N constants of viscoelasticity all of which are independent.

4. Discretization of constitutive equation

Below a numerical algorithm for the incremental integration of the developed model is intro-
duced. For the time increment n+1, the total stress is, according to Eq. (3.1), given by the
equation

Sn+1 =S0n+1+
P∑

k=1

H̃kn+1+
N∑

j=1

Hjn+1 (4.1)

The reccurence-update formula for the integration of Eq. (3.3) is obtained using the central
difference method, thus the following approximations are utilized

ẏ
∣∣∣
n+1
2

≈
yn+1−yn

∆t
=
∆yn+1

∆t
y
∣∣∣
n+1
2

≈ yn+
1

2
∆yn+1 (4.2)

Applying the rules given by Eqs (4.2) to Eq. (3.3) one obtains the following formula

H̃kn+1 =

(
1− 1

D̃k

∆zn+1
2

)
H̃kn+γk(S

iso
0n+1−S

iso
0n)

1+ 1
D̃k

∆zn+1
2

(4.3)

where

∆zn+1 =(∆Cn+1 ·∆Cn+1)
1

2 (4.4)

The integration of Eq. (3.9) is performed bymeans of the reccurence-update formula developed
by Taylor et al. (1970), i.e.

Hjn+1 =e
−∆t
τj Hjn+Γj

1− e
−∆t
τj

∆t
τj

(Siso0n+1−S
iso
0n) (4.5)

Equations (4.1), (4.3) and (4.5) form up a discretized set of equations which can be utilized for
numerical simulations. For the use of FEM, amaterial tangent stiffness tensor has to be defined
in addition. The derivation of a tangent stiffness is discussed further in the text.



598 C. Suchocki

5. Exemplary application

In order to check the ability of the developed viscoelastic-plastic constitutive equation to fit the
experimental data, it has been employed to capture the mechanical properties of UHMWPE
which is a popular biomaterial with numerous applications in orthopaedics, such as hip joint
implants, elbow implants or artificial intervertebral discs for instance.

For the purpose of modeling the inelastic behavior of UHMWPE, a version of the develo-
ped constitutive equation has been adopted with N = 3 viscoelastic overstresses and P = 1
endochronic overstress. The graphic interpretation of themodel can be seen in Fig. 1.

Fig. 1. Mechanical scheme of the rheological model assumed for UHMWPE

Anumber of mechanical tests have been performed on cylindrical specimensmachined from
a UHMWPE rod. The dimensions of the specimens are ∅17 and l0 = 21mm. A medical grade
polyethyleneCHIRULEN1050 has been used.All the experiments have been conducted onMTS
Bionix 2500 testing machine at a constant room temperature of 20◦C.
A number of loading-unloading compression tests at a constant deformation rate λ̇ has been

performed. The used values of the deformation rate are 0.0005, 0.005 and 0.05s−1 with λ equal
to 0.93 as the minimum stretch ratio. Furthermore, a relaxation test in compression has been
carried out. The specimen has been loaded with a constant deformation rate of 0.003s−1 and
subsequently allowed to relax for 30 minutes.

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

The specimen stretch ratios in both axial an perpendicular directions have been measured
using a video extensometer. The axial deformation has been at the same timemeasured with a
strain gauge extensometer in order to verify the measured values. A sketch of the experimental
setup can be seen inFig. 2.APoisson ratio of 0.46 has been determined,which justifies adopting
an assumption of material incompressibility.



An internal-state-variable based viscoelastic-plastic... 599

The material parameters of the constitutive equation have been evaluated using the least
squares method. Thus, the following error function has beenminimized

F(p) =
P∑

k=1

[(T11(p))k − (T̃11)k]
2 (5.1)

where T11(p) is the theoretical and T̃11 is the experimental Lagrange stress in the axial
direction of the specimen at the time instant tk (k = 1,2, . . . ,P). The column matrix
p= [µ,b,κ,γ1,D̃1,Γ1,Γ2,Γ3,τ1,τ2,τ3]

T contains the material parameters being optimized. The
experimental data from three loading-unloading tests at various deformation rates together with
the data from the relaxation test have been utilized to formup theminimized objective function
given by Eq. (5.1). The results of the curve fitting are shown in Fig. 3, whereas the determined
values of thematerial constants of elasticity, viscoelasticity and plasticity have been collected in
Table 1. The inverese of the bulk modulus D1 does not participate in the optimization. In FE
implementation,D1 is set to a value corresponding to the bulkmodulus high enough to account
for the material near incompressibility.

Fig. 3. Fitting of the viscoelastic-plastic constitutive model to the experimental data.
Loading-unloading compression tests with the deformation rate of: (a) 0.0005s−1, (b) 0.005s−1,
(c) 0.05s−1 and stress relaxation tests with the relaxation time of: (d) 10s, (e) 200s, (f) 1700s

Table 1.Model parameters for UHMWPE

Elastic constants Viscoelastic constants Plastic constants

µ 52.56MPa Γ1 2.89 τ1 0.16 s γ1 3.25 D̃1 0.029

b 209.28 Γ2 0.93 τ2 8.803 s

κ 0.81 Γ3 0.62 τ3 279.16 s

D1 0.00033 MPa
−1



600 C. Suchocki

6. Finite element implementation

The linearized constitutive equation used for the FE implementation is obtained by taking a
directional derivative of Eq. (4.1) with respect toCn+1. Thus, for the n+1 increment

∆Sn+1 =C
ve−p
n+1 ·

1

2
∆Cn+1 C

ve−p
n+1 =2

∂Sn+1

∂Cn+1
(6.1)

where the approximate material tangent tensor takes the form

C
ve−p
n+1 =C

vol
0n+1+

[
1+

P∑

k=1

γk

(
1+
∆zn+1

2D̃k

)−1
+
N∑

j=1

Γj
1−e

−∆t
τj

∆t
τj

]
C
iso
0n+1 (6.2)

The incremental constitutive rate equation given by Eq. (6.1) can be expressed using the
Zaremba-Jaumann (Z-J) objective rate of the Kirchhoff stress τ, thus

τ
∇
n+1 = Jn+1C

MJ
n+1 ·∆Dn+1 (6.3)

where the incremental Z-J objective rate of the Kirchhoff stress

τ
∇
n+1 =∆τn+1−∆Wn+1τn+1−τn+1∆W

T
n+1 (6.4)

and

∆Wn+1 =
1

2

[
∆Fn+1F

−1
n+1− (∆Fn+1F

−1
n+1)

T
]

∆Dn+1 =
1

2

[
∆Fn+1F

−1
n+1+(∆Fn+1F

−1
n+1)

T
]

∆Fn+1 =Fn+1F
−1
n

(6.5)

Thematerial stiffness tensor takes the corresponding form

C
MJ
n+1 =

1

Jn+1
(Cτcn+1+In+1) (6.6)

where

C
τc
n+1 =(FiPFjQFkRFlSC

ve−p
PQRS)n+1ei⊗ej ⊗ek⊗el

In+1 =
1

2
(δikτjl+ τikδjl+ δilτjk+ τilδjk)n+1ei⊗ej ⊗ek⊗el

(6.7)

and ek (k=1,2,3) are the unit vectors of the Cartesian basis. TheKirchhoff stress in Eq. (6.4)
is determined from the stress transformation law, i.e.

τn+1 =Fn+1Sn+1F
T
n+1 (6.8)

where the second Piola-Kirchhoff stress, Sn+1 is calculated utilizing Eq. (4.1).
Substitution of Eq. (6.2) into Eq. (6.7)1 results in the following form of thematerial stiffness

in the current configuration

C
τc
n+1 =C

vol
n+1+

[
1+

P∑

k=1

γk

(
1+
∆zn+1

2D̃k

)−1
+
N∑

j=1

Γj
1−e

−∆t
τj

∆t
τj

]
C
iso
n+1 (6.9)

where

C
vol =

2

D1
J(J−1)(1⊗1−2I)+J2

2

D1
1⊗1 (6.10)



An internal-state-variable based viscoelastic-plastic... 601

and

C
iso =

2

3
µ
[
1+
b

κ
(I1−3)

]κ−1
I1

(
I+
1

3
1⊗1

)
−
2

3
J−

2

3µ
[
1+
b

κ
(I1−3)

]κ−1
(B⊗1+1⊗B)

+2J−
4

3µ
b(κ−1)

κ

[
1+
b

κ
(I1−3)

]κ−2
B⊗B

−
2

3
J−

2

3µ
b(κ−1)

κ

[
1+
b

κ
(I1−3)

]κ−2
I1(B⊗1+1⊗B)

+
2

9
µ
b(κ−1)

κ

[
1+
b

κ
(I1−3)

]κ−2
I
2
11⊗1

(6.11)

are, respectively, the volumetric and isochoric components of the elasticity tensor corresponding
to the stored-energy functions as given in Eqs (3.11), cf Suchocki (2011).
The viscoelastic-plastic model has been implemented into FE software Abaqus by taking

advantage of the user subroutineUMAT(UserMATerial), which is called by theFE solver during
every iteration of the Newton-Raphson numerical procedure (Hibbit et al. 2008). The written
subroutine calculates Cauchy stress tensor and material Jacobian defined in Eq. (6.6) for each
finite element.These quantities are subsequentlyusedbyAbaqus to formupthe element stiffness
matrix. Finally, the global stiffness matrix is assembled by Abaqus using the element stiffness
matrices. The calculations done by the written subroutine UMAT are listed further in the text.
The user subroutines used in other FE programs to define custom constitutive equations have a
similar structure.
In order to verify the performance of the developed UMAT code, several exemplary si-

mulations have been performed using the Abaqus FEM program. The simulations involved a
15mm×15mm×15mmUHMWPE block undergoing ramp tension test, ramp compression test
and sinusoidal deformation with variable amplitude. In the first approach, the polymeric block
is meshed using a single finite element C3D8H1. Subsequently, all the simulations are conduc-
ted using the polyethylene block meshed with 125 elements. In each of the simulated processes
the excitation was kinematic, i.e. the frontal face of the block performs a designated displace-
ment program.The displacement distribution and the used boundary conditions are depicted in
Fig. 4. Due to the assumed incompressibility, the deformation gradient tensor has the following
components

F3×3 =



λ1 0 0
0 1√

λ1
0

0 0 1√
λ1


 (6.12)

where λ1 denotes the stretch ratio in the direction x (Fig. 4).

Fig. 4. Homogeneous deformation of a single finite element. (a) Distribution of the displacement
component in the tension/compression direction. (b) Applied boundary conditions

All the simulations have been repeated utilizing programs written in MATLAB in which
the discretized equations presented in Section 4 are used to calculate the stresses for the given
deformation histories.
1Cubic, three-dimensional, 8 nodes, hybrid.



602 C. Suchocki

Fig. 5. Polyethylene block undergoing (a) ramp tension and (b) ramp compression tests with the
deformation rate of 0.003s−1 and (c) cyclic deformation history – a sinusoidal function with the

increasing amplitudeA

Algorithm for the implementation in Abaqus

Input:Fn+1, no. of direct and shear stress components

1. Calculate strain measures from the current increment

Cn+1 =F
T
n+1Fn+1 Fn+1 = J

−1
3

n+1Fn+1 Bn+1 =Fn+1F
T
n+1

2. Extract variables from the previous increment

3. Calculate elastic stresses from the current increment

S
vol
0n+1 = Jn+1pn+1C

−1
n+1 pn+1 = ∂Jn+1U(Jn+1)

S
iso
0n+1 = J

−2
3

n+1DEV[Sn+1] Sn+1 =2∂Cn+1W(Cn+1)

S0n+1 =S
vol
0n+1+S

iso
0n+1

4. Update endochronic and viscoelastic overstresses (k=1,2, . . . ,P),

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

∆zn+1 =(∆Cn+1 ·∆Cn+1)
1

2 ∆Cn+1 =Cn+1−Cn

H̃kn+1 =

(
1− 1

D̃k

∆zn+1
2

)
H̃kn+γk(S

iso
0n+1−S

iso
0n)

1+ 1
D̃k

∆zn+1
2

Hjn+1 =e
−∆t
τj Hjn+Γj

1− e
−∆t
τj

∆t
τj

(Siso0n+1−S
iso
0n)

5. Calculate total stress from the current increment

Sn+1 =S0n+1+
P∑

k=1

H̃kn+1+
N∑

j=1

Hjn+1 σn+1 =
1

Jn+1
Fn+1Sn+1F

T
n+1

6. Calculate viscoelastic-plastic stiffness from the current increment

7. Store stresses and isochoric C-G tensor from the current increment



An internal-state-variable based viscoelastic-plastic... 603

In Fig. 5, the results of ramp tension (Fig. 5a) and compression (Fig. 5b) are shown. In the
case of ramptension simulation, thematerial is deformedataconstant rate λ̇1 =−0.003s

−1until
the stretch ratio λ1 = 1.1 is reached. Subsequently, unloading is performed at the deformation
rate of λ̇1 =0.003s

−1. The compression process was simulated analogously.

In Fig. 5c, the results of sinusoidal loading simulation are presented. The simulated process
comprises of four periods, each performed at different amplitude.

The predictions of the FEM simulations utilizing UMAT are in a good agreement with the
results produced by MATLAB programs (see Fig. 5). All the simulations have been performed
utilizing the material parameter values collected in Table 1.

7. Conclusions

In this study, anewviscoelastic-plasticmodel for thermoplastic polymers and resins is presented.
Themodel is formulated by utilizing the formalism of stress-like internal state variables. It has
been implemented into the FE software Abaqus and applied to model the inelastic behavior of
polyethylene with very good results.

References

1. Abdul-Hameed H., Messager T., Zäiri F., Näit-Abdelaziz M., 2014, Large-strain
viscoelastic-viscoplastic constitutive modeling of semi-crystalline polymers and model identifica-
tion by deterministic/evolutionary approach,Computational Materials Science, 90, 241-252

2. AlkasYonanS., SoyarslanC.,HauptP.,Kwiatkowski L.,TekkayaA.E., 2013,Asimple
finite strainnon-linearvisco-plasticmodel for thermoplasticsand its application to the simulationof
incremental cold formingofpolyvinylchloride (PVC), International Journal ofMechanical Sciences,
66, 192-201

3. Ayoub G., Zäiri F., Fréderix C., Gloaguen J.M., Näit-Abdelaziz M., Seguela R.,
Lefebvre J.M., 2011, Effects of crystal content on the mechanical behaviour of polyethylene
under finite strains: experiments and constitutive modelling, International Journal of Plasticity,
27, 492-511

4. Bardella L., 2001, A phenomenological constitutive law for the nonlinear viscoelastic behaviour
of epoxy resins in the glassy state,European Journal of Mechanics – A/Solids, 20, 907-924

5. Ben Hadj Hamouda H., Laiarinandrasana L., Piques R., 2007, Viscoplastic behaviour of a
medium density polyethylene (MDPE): constitutive equations based on double nonlinear deforma-
tion model, International Journal of Plasticity, 23, 1307-1327

6. Bergström J.S., Bischoff J.E., 2010, An advanced thermomechanical constitutive model for
UHMWPE, International Journal of Structural Changes in Solids, 2, 1, 31-39

7. Drozdov A.D., Christiansen J.deC., 2007, Cyclic viscoplasticity of high-density polyethyle-
ne/montmorillonite clay nanocomposite,European Polymer Journal, 43, 1, 10-25

8. Hassan T., Colak O.U., Clayton P.M., 2011, Uniaxial strain and stress-controlled cyclic
responses of ultrahighmolecular weigth polyethylene: experiments andmodel simulations, Journal
of Engineering Materials and Technology, 133, 021010-1–021010-9

9. Hibbit B., Karlsson B., Sorensen P., 2008, ABAQUS Theory Manual, Hibbit, Karlsson &
Sorensen Inc.

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

11. 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, 5, 611-639



604 C. Suchocki

12. Krairi A., Doghri I., 2014, A thermodynamically-based constitutive model for thermoplastic
polymers coupling viscoelasticity, viscoplasticity and ductile damage, International Journal of Pla-
sticity, 60, 163-181

13. Kästner M., Obst M., Brummund J., Thielsch K., Ulbricht V., 2012, Inelastic behavior
of polymers – experimental characterization, formulation and implementation of amaterialmodel,
Mechanics of Materials, 52, 40-57

14. PipkinA.C., Rivlin R.S., 1965,Mechanics of rate-independentmaterials,ZAMP, 16, 3, 313-327

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

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

17. 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

18. Wu H.C., 2005,Continuum Mechanics and Plasticity, CRCPress, NewYork

Manuscript received October 26, 2014; accepted for print January 7, 2015