Mech Time-Depend Mater (2018) 22:421–433

DOI 10.1007/s11043-017-9364-x

An anisotropic linear thermo-viscoelastic constitutive law

Elastic relaxation and thermal expansion creep in the time domain

Heinz E. Pettermann1 · Antonio DeSimone2

Received: 17 October 2016 / Accepted: 1 September 2017 / Published online: 19 September 2017
© The Author(s) 2017. This article is published with open access at Springerlink.com

Abstract A constitutive material law for linear thermo-viscoelasticity in the time domain is
presented. The time-dependent relaxation formulation is given for full anisotropy, i.e., both
the elastic and the viscous properties are anisotropic. Thereby, each element of the relax-
ation tensor is described by its own and independent Prony series expansion. Exceeding
common viscoelasticity, time-dependent thermal expansion relaxation/creep is treated as in-
herent material behavior. The pertinent equations are derived and an incremental, implicit
time integration scheme is presented.
The developments are implemented into an implicit FEM software for orthotropic mate-
rial symmetry under plane stress assumption. Even if this is a reduced problem, all essential
features are present and allow for the entire verification and validation of the approach. Var-
ious simulations on isotropic and orthotropic problems are carried out to demonstrate the
material behavior under investigation.

Keywords Viscoelastic · Anisotropic · Constitutive laws · Thermal expansion creep ·

Finite element method implementation

1 Introduction
The apparent effects of viscoelastic material behavior manifest themselves in a time-
dependent response to loading and accompanying energy dissipation. Relaxation or creep
occurs when a viscoelastic material is exposed to quasi-static loads and load changes. Un-
der cyclic excitation damping is exhibited. These phenomena are treated typically in the
time and in the frequency domain, respectively.
Such viscoelastic effects are widespread in natural as well as in engineering materials.
Among them are almost all biological tissues and most polymers, in particular thermoplastic

B H.E. Pettermann
pettermann@ilsb.tuwien.ac.at
1 Institute of Lightweight Design and Structural Biomechanics, Vienna University of Technology,
Vienna, Austria
2 Scuola Internazionale Superiore di Studi Avanzati – SISSA, Trieste, Italy
422 Mech Time-Depend Mater (2018) 22:421–433

materials. Biological materials are composites “by nature”, whereas engineering polymers
are often mixed with other constituents to improve their performance. These materials are
likely to exhibit thermal expansion relaxation/creep behavior, i.e. time-dependent changes
of unconstrained thermal strains as a response to temperature changes. Such composites,
natural and man made ones, often show direction dependent properties and the consideration
of anisotropy becomes inevitable.
A general introduction to viscoelasticity can be found, e.g. in Lakes (2009), Gurtin and
Sternberg (1962), Schapery (2000), which focus predominantly on isotropic behavior and
treat the time as well as the frequency domain. Transversely isotropic (and orthotropic)
viscoelastic models have been presented, e.g. in Holzapfel (2000), Bergström and Boyce
(2001), Pandolfi and Manganiello (2006), Puso and Weiss (1998), Taylor et al. (2009), Ned-
jar (2007), based on invariants of the strain representation, which are commonly used for
biological soft tissues. Typically, they combine non-linear orthotropic elasticity with lin-
ear isotropic viscosity. Linear viscoelasticity of orthotropic media is presented in Carcione
(1995), Dong and McMechan (1995) in the context of geo-materials. Linear and non-linear
orthotropic viscoelasticity adopting not the full set of orthotropic viscous effects is presented
in Melo and Radford (2004), Kaliske (2000) and applied to composites and foam materials,
respectively. Full elastic anisotropy is treated in Lévesque et al. (2008), but with one single
relaxation function operating on the entire tensor. In Schapery (1964) full anisotropy is mod-
eled with various relaxation functions applied to elasticity tensor contributions. A review on
damping of composites is given in Chandra et al. (1999). Anisotropic non-linear viscoelas-
tic constitutive laws have been introduced in Schapery (1969), Zocher et al. (1997), Sawant
and Muliana (2008), Poon and Ahmad (1998) which found widespread application to study
composite materials. Implementation aspects are discussed in Crochon et al. (2010), Sorvari
and Hämäläinen (2010).
Relaxation type effects occurring for the thermal expansion behavior are mentioned in
Lakes (2009), Zocher et al. (1997), Schapery (1964), and observed experimentally for poly-
mers, e.g. Spencer and Boyer (1946). A hygro-thermal expansion relaxation function is
presented in Bažant (1970) in the context of concrete behavior. Otherwise, to the knowledge
of the authors, there are no modeling attempts in the literature for thermal expansion creep.
It is noted that thermal expansion relaxation/creep, here, does not refer to secondary effects
arising, e.g. from transient thermal fields, but solely denotes the time-dependent response as
inherent material behavior.
In the present work linear thermo-viscoelasticity of anisotropic materials is treated in the
time domain. Thereby the thermo-elastic as well as the viscous responses are considered
to behave anisotropically. The temperature dependence for thermorheologically simple ma-
terials is accounted for by a time–temperature shift function. Time-dependent, anisotropic
thermal expansion relaxation/creep is treated.
The constitutive formulations are developed and the pertinent equations are given, includ-
ing the consistent material tangent operator. The state of the material is handled by internal
variables. The algorithms are implemented into an implicit Finite Element Method (FEM)
software package. This is done under the assumption of plane stress states and orthotropic
materials, but the approach is general and the implementation can be extended to tri-axial
problems and general anisotropy in a straightforward manner just by adding additional terms
of the same form.
For examples of applications, homogenized orthotropic material data is computed from
a simple unit cell approach which serves as input to the developed constitutive laws. This
allows for a consistent verification and validation of the developments. Not only the principal
load cases are simulated to extract the input data, but also general loading scenarios and
loading histories.
Mech Time-Depend Mater (2018) 22:421–433 423

2 Linear viscoelastic constitutive model

The constitutive model for considerations in time domain has to give the stress response as
a function of the strain (rate) history and the appropriate instantaneous material behavior.
For isotropic material symmetry the stress and strain tensors can be split into their volu-
metric and deviatoric part. The corresponding scalar valued material properties are the bulk
and shear relaxation functions. A similar decomposition can be performed for cubic mate-
rial symmetry (e.g. used in Pettermann and Hüsing 2012) and for the transversely isotropic
case; see e.g. Suquet and Bornert (2001). Alternatively, one can deal with the tensorial equa-
tions, as is pursed in the present work which allows for the treatment of the most general
anisotropic case. Both the elastic and the viscous contributions to the response are allowed
to be of general anisotropy.
Adopting Voigt notation, the hereditary integral in tensorial form reads
 t
σi (t) = Rij (t − s)ε̇j (s) ds. (1)
0

Stress and strain are given by

σi = (σ11 σ22 σ33 σ12 σ13 σ23 )T , εi = (ε11 ε22 ε33 γ12 γ13 γ23 )T , (2)

with γij = 2εij . Note that in Voigt notation, coordinate transformation operators must be
formulated correspondingly. The time-dependent material tensor for anisotropic materials is
given by the symmetric 6×6 matrix, Rij (t) = Rj i (t). In the most general case it is composed
of 21 individual relaxation functions of the form,
  
  
Rij (t) = Rij 0 1 − rij k 1 − exp −t/τ rij k (no sum on ij ), (3)
k

which are given by Prony series expansions with Rij 0 being the instantaneous values,
i.e., the elasticity tensor elements giving the short term behavior. The sums over k Prony
terms contain the relative relaxation elements, rij k , and the corresponding characteris-
tic times, τ rij k . Note that the relaxation tensor, Rij (t), possesses off-diagonal terms with
i, j ≤ 3, implicitly containing some Poisson relaxation which is not necessarily monotonic;
see e.g. Lakes and Wineman (2006). The stress components can be expressed by combining
Eqs. (1) and (3) as,

   rij k  t 
6
 
σi (t) = Rij 0 εj − exp −s/τ rij k
ε j (t − s) ds , (4)
j =1 k
τ rij k 0

where the sum over j gives the contributions of the individual strain components. The “rel-
ative creep contribution” for each Prony term,

j 1 t  
ei k = exp −s/τ rij k εj (t − s) ds, (5)
τ rij k 0

is introduced and summed to the “total creep contributions” by

j
 j
i = rij k ei k , (6)
k
424 Mech Time-Depend Mater (2018) 22:421–433

where the superscript indicates the cause (i.e., the applied strain) and the subscript denotes
the stress component in Eq. (4) which is affected. This unusual super/subscript notation is
chosen to discern different types of “strain like” quantities unambiguously. Finally, the stress
components read,

6
 j
σi = Rij 0 εj − i . (7)
j =1

j j j
Note that the expressions i and ei k cannot be interpreted as creep strains, and i = ji .
Consequently, they are rather property type quantities carrying information on how much
j
of the relaxation potential has been consumed until time t . Nevertheless, the ei k for every
Prony term k represent the internal state variables.
The proposed approach may be viewed as an ad hoc extension of common linear vis-
coelastic models. Its thermodynamic consistency, however, is not ensured and could be vio-
lated by a unsuited set of material parameters.
The temperature dependence of viscoelastic thermorheologically simple materials is
widely treated by introducing a reduced time; see e.g. Lakes (2009),
t
tred = (8)
A(T )

with the shift function, A(T ), depending on the temperature, T . In the present work the
empirical Williams–Landel–Ferry (WLF) function, see e.g. Lakes (2009),

C1 (T − Tref )
log A(T ) = − , (9)
C2 + (T − Tref )
is used with Tref being the reference temperature and C1 and C2 are constants pertaining to
the particular material under consideration. In general cases, each relaxation element, Rij (t),
can have its own time–temperature shift function, Aij (T ).

3 Thermal expansion relaxation/creep constitutive model

In linear thermo-(visco)elasticity assuming small strains, additive decomposition can be ap-

plied, which gives the total strain (in Voigt notation),

εitot = εive + εith , (10)

as the sum of the viscoelastic ve and the thermal th contributions.

Now, materials may show thermal expansion which exhibits time dependence as an in-
herent material behavior (e.g. composites). Following Lakes (2009), Zocher et al. (1997),
the time-dependent thermal expansion relaxation is given by
 t
εith (t) = αi (t − s)ϑ̇(s) ds, (11)
0

with αi being the time-dependent coefficients of thermal expansion and, here, temperature
independent behavior is assumed. ϑ is the temperature change from some (stress free) start-
ing temperature. For the time independent case, the usual thermal expansion results from
integration of Eq. (11) as εith = αi ϑ . Note that Eq. (11) exactly resembles the hereditary
Mech Time-Depend Mater (2018) 22:421–433 425

integral in Eq. (1) for the time-dependent stresses under strain loading. Thus, the further
treatment goes in perfect analogy to linear viscoelasticity as studied before.
For the most general anisotropic case, the time-dependent thermal expansion tensor (in
correspondence to the strain in Eq. (2)) reads
 T
αi (t) = α11 (t) α22 (t) α33 (t) α12 (t) α13 (t) α23 (t) . (12)

The components are given by the six individual relaxation functions,

   
 
αi (t) = αi 0 1 − ai k 1 − exp −t/τ ai k (no sum on i), (13)
k

with the instantaneous thermal expansion coefficients αi 0 and the relative relaxation values
ai k as well as characteristic times τ ai k for the Prony series with k terms. For common ma-
terials, the relative relaxation values, ai k , are negative which give rise to creep type thermal
expansion behavior.
Applying Eqs. (4) to (7) analogously, one obtains the thermal strain:
  
εith = αi 0 ϑ − ai k θi k (no sum on i). (14)
k

The “relative creep temperatures”, θi k , have no direct physical meaning, however, they carry
the information on the amount of relaxation and, consequently, are used as internal state vari-
ables. The concept of reduced time and a shift function, Eq. (8), can be applied equivalently.
The viscoelasticity and the thermal expansion relaxation have to be solve simultaneously,
fulfilling Eq. (10).

4 Time integration and FEM implementation

The scheme of the implicit incremental time solution for arbitrary strain histories in the con-
text of FEM can be found e.g. in Crochon et al. (2010), Sorvari and Hämäläinen (2010), Poon
and Ahmad (1998), Muliana and Khan (2008). The pertinent equations are given briefly in
the following. Key to the stress update and the consistent Jacobian formulation is the time
integration of Eq. (1) with the material properties given in Eq. (3). Besides the strain rate,
the only time-dependent contributions are the exponential functions in Eq. (3). They appear
in all elements of Rij (t), just with different characteristic times. With focus on these terms
one can extract an integral of the form,
 t
j   dεj
ei = 1 − exp s − t/τ rij k ds. (15)
0 ds
j
where dεj /ds is the appropriate strain rate. Here, ei is the “relative creep contribution” as
introduced in Eq. (5). The above equation is solved for a finite increment of time, t , added
to time t n at which the state is completely known and under the assumption that the strain
increases linearly with time. This yields the increment of the relative creep contribution as,

j τ rij k t t t  j n
ei = rij k + exp − rij k −1 εj + 1 − exp − rij k εjn − ei , (16)
t τ τ τ
426 Mech Time-Depend Mater (2018) 22:421–433

with the superscript n denoting the state at the beginning of the increment. With this expres-
sion the stress update can be given as


6   
j
σi = Rij 0 εj − rij k ei k . (17)
j =1 k

j
The role of the solution dependent state variables is taken by ei k , accumulated over the time
increments. The equation for computing the consistent tangent moduli elements depends on
the time increment only,
  
tangent τ rij k t t
Rij = Rij 0 1 − rij k + exp − − 1 , (18)
k
t τ rij k τ rij k

and does not include any matrix summation. In the implementation of this solution scheme,
j j
“symmetry” on the creep contributions ei k (and i ) must not be enforced; see Eq. (7). How-
ever, the material input requires Rij 0 = Rj i 0 , rij k = rj i k , and τ rij k = τ rj i k , to maintain
symmetry in the relaxation tensor at any time.
The treatment of thermal expansion relaxation/creep is equivalent to the stress update,
because of the identical structure of the governing equations, and it yields the increment of
thermal strain components,
  
εith = αi 0 ϑ − ai k θi k , (no sum on i). (19)
k

Similar to the consistent tangent moduli, the coefficients of thermal expansion at the end
of the time increment are required for effective convergence. They are computed in perfect
analogy with Eq. (18).

4.1 Orthotropic plane stress implementation

The implementation of the viscoelastic model is carried out with the UMAT option, the
thermal expansion with the UEXPAN option, for the implicit FEM software package
ABAQUS/Standard v6.14 (Dassault Systèmes Simulia Corp., Providence, RI, USA). The
present implementation allows for 13 Prony terms maximum for any of the relaxation type
properties. It is carried out for the case of orthotropic material symmetry under plane stress
conditions. Even if this poses a reduced problem, all essential features of the proposed con-
stitutive law are included, and the model formulation and the software realization can be
tested and verified entirely. Extensions to more general cases just need the consideration of
additional terms in the summation loops of exactly the same type.
For constant temperatures which are not equal to the WLF-reference temperature,
Eqs. (8) and (9) can be applied directly to compute the reduced time and likewise the reduced
time increment. If the temperature changes within a time increment, an approximation for
small temperature increments is utilized,
2
tred = t, (20)
A(T n+1 ) + A(T n )
assuming one unique time–temperature shift function for all Rij (t) elements. Extension to
individual shift functions is straightforward by modifying Eqs. (16) and (18). The time–
temperature shift option is implemented for the viscoelastic behavior only, but not (yet) for
the time-dependent thermal expansion relaxation/creep.
Mech Time-Depend Mater (2018) 22:421–433 427

Table 1 Isotropic linear thermo-viscoelastic plane stress material data as input to the ABAQUS material
law; instantaneous elastic moduli, bulk and shear relaxation by one Prony term each, WLF temperature shift
data, and coefficient of thermal expansion

E0 = 2500 MPa ν0 = 0.25 (G0 = 1000 MPa)

b = 0 (τ = 1.0 s) g = 0.5 τ g = 1.0 s
b

Tref = 0◦ C C1 = 17 C2 = 50
α0 = 1 × 10−4 /◦ C

To control the accuracy of the solution, i.e. controlling t in Eq. (16), an automatic time
stepping procedure is implemented.

4.2 Material parameters

For the measurement of relaxation data which can serve as input to a viscoelastic constitu-
tive law there are standard procedures for simple loading scenarios. For intricate material
response as given by Rij , but in particular by R12 , a direct measurement can be very difficult
or impossible. So, to calibrate the input from a series of (feasible) experiments some kind of
optimization procedure might be useful.
In the present work the input data is generated by micromechanical type unit cell predic-
tions which directly yields the viscoelastic relaxation matrix elements and thermal expansion
relaxation coefficients. These functions of time are fitted by Prony series expansions.

5 Verification and examples

For clarity of the following, the stress and strain components will be given with their full
indication. For elements of the plane stress relaxation tensors (elastic and thermal expan-
sion), the compact Voigt notation is maintained. Unless stated otherwise, all simulations
are carried out at the WLF-reference temperature. All simulations employ four noded plane
stress elements, using reduced integration for single element tests and full integration for the
structural simulations. The element thickness is 1 mm.

5.1 Isotropic material—single element tests

Single element tests with isotropic linear viscoelastic material are carried out first. The ele-
ment size is 1 × 1 mm. The ABAQUS built in material is used to obtain reference solutions,
the generic material data (of some glassy polymer) is given in Table 1, which show shear
relaxation only, plotted in Fig. 1 (left) as “R33”. Simulations are performed to compute the
material parameters pertaining to the developed material model, i.e. the response to “uni-
axial strain” loading is sought for. Here, uni-axial strain denotes the in-plane state, in the
perpendicular direction plane stress is maintained throughout this work. A Heaviside step
function is applied at t = 0 and kept constant afterwards as ε11 (t > 0) = 1. The other in-
plane strain components are ε22 = 0 and γ12 = 0 at any time. From this load case the Ri1 (t)
components of the relaxation matrix can be deduced directly, see Fig. 1 (thin solid lines).
Note that the response which contains the Poisson effect, R12 (t), is non-monotonic. For
completeness, the shear behavior R33 (t) is added which exactly resembles the input, G(t).
Based on this data, the material parameters for the UMAT are calibrated. The (sum of the)
relative relaxation elements can be computed directly from the long term and short term
428 Mech Time-Depend Mater (2018) 22:421–433

Fig. 1 Single element predictions of the relaxation functions Rij (t) as response to Heaviside step strains by
the ABAQUS built in material law (“abamat”) and the developed constitutive material law (“umat”) calibrated
to the former one

Table 2 Isotropic linear viscoelastic plane stress material data as input to the UMAT; instantaneous elasticity
matrix elements, relaxation matrix elements by two Prony terms each, and WLF temperature shift data

R11 0 = 2666 MPa R22 0 = 2666 MPa R12 0 = R21 0 = 666 MPa R33 0 = 1000 MPa

k=1 r11 = 0.41 r22 = 0.41 r12 = r21 = 0.25 r33 = 0.25
τ r11 = 1.25 s τ r22 = 1.25 s τ r12 = τ r21 = 2.0 s τ r33 = 1.0 s

k=2 r11 = 0 r22 = 0 r12 = r21 = −0.11 r33 = 0.25

τ r11 = 1.25 s τ r22 = 1.25 s τ r12 = τ r21 = 0.6 s τ r33 = 1.0 s

Tref = 0◦ C C1 = 17 C2 = 50

response. The characteristic times can be deduced from the initial slope. If two Prony terms
are required, a manual “trail and error” procedure is used.
The diagonal element R11 (t) is well captured by one Prony term, for the off-diagonal
element, R21 (t) = R12 (t), at least two Prony terms are required to fit the non-monotonic
behavior. The resulting material parameters for the UMAT are listed in Table 2. Now, these
parameters and the developed UMAT are employed to simulate the same load case as before;
see Fig. 1 (bold dashed lines). Comparison to the input shows some deviations but the result-
ing material behavior is considered sufficiently accurate. Note the difference in relaxation
for various Rij term.
A series of additional tests and verifications are run based on single elements. Among
them are classical relaxation test (under uni-axial stress), bi-axial loading, and arbitrary load
cases. Besides the response to step functions, arbitrary ramp loading and load reversals are
realized. All comparisons give satisfactory results. Also, the time–temperature shift function
for temperature dependent viscoelastic material data is tested successfully as well as the
automatic time stepping algorithm.

5.2 Isotropic material—structural simulations

Next, simulations of structural problems with isotropic material are carried out which al-
low direct comparison with reference solutions. For this purpose a quadratic patch of unit
length is discretized with 20 × 20 plane stress elements. The patch contains two parallel
rectangular voids with 0.2 × 0.6 mm; see Fig. 2. The boundary conditions are “unit cell”
Mech Time-Depend Mater (2018) 22:421–433 429

Fig. 2 Structure with voids

(light gray), also representing a
unit cell of an orthotropic
material

Fig. 3 Orthotropic relaxation functions Rij (t) as response to Heaviside step strains; homogenization of
a voided structure (“unit cell”) and single element predictions by the developed constitutive material law
(“umat”) calibrated to the homogenized behavior

like to be used later for homogenization; see Pettermann and Suresh (2000), Böhm (2004).
The deformations of the nodes at the edges are coupled to the adjacent corner nodes of the
model to remain on a straight line. The deformation of the model is controlled by the corner
nodes, whereas the upper right corner is coupled to the others to maintain a parallelepipedic
shape. This way, a unit cell type model is set up with three master nodes to apply loads
and to read the response. Note that this does not give correct periodic boundary conditions
(see e.g. Pettermann and Suresh 2000; Böhm 2004) under shear loads. However, it gives
some homogenized behavior which sufficiently serves for the purpose of comparison in the
present study.
The master nodes are used to apply a variety of loading scenarios as discussed above.
The material model employed is either the ABAQUS built in one or the developed UMAT,
whereas the voids are realized by setting the Young’s modulus to Evoid = R33 0 × 10−3 . The
predicted master node responses and the local, time-dependent stress fields are compared,
and for all studied load cases almost one to one agreement is found.
The implemented automatic time stepping as well as the consistent tangent matrix,
i.e., the Jacobian, work correctly and the performance of the implemented material law is
very satisfactory.

5.3 Orthotropic porous material

Now, the structure investigated in the previous section can be interpreted as porous material
with orthotropic properties. By homogenization, the effective linear viscoelastic orthotropic
plane stress behavior can be predicted. Applying the load cases of Heaviside step uni-axial
strain in the two principal directions and for shear, the entire constitutive behavior can be
predicted (as introduced in Sect. 5.1). The relaxation response is presented in Fig. 3 (de-
noted by “unit cell”), with R12 (t) = R21 (t). All elements of Rij (t) are approximated by one
430 Mech Time-Depend Mater (2018) 22:421–433

Table 3 Orthotropic linear viscoelastic plane stress material data as input to the UMAT; instantaneous elas-
ticity matrix elements, relaxation matrix elements by one Prony term each, and WLF temperature shift data
(taken to be equal to that of the input material)

R11 0 = 1693 MPa R22 0 = 1124 MPa R12 0 = R21 0 = 227 MPa R33 0 = 242 MPa

k=1 r11 = 0.440 r22 = 0.439 r12 = r21 = 0.234 r33 = 0.453
τ r11 = 1.09 s τ r22 = 1.09 s τ r12 = τ r21 = 1.98 s τ r33 = 1.08 s

Tref = 0◦ C C1 = 17 C2 = 50

Prony term, the material parameters are listed in Table 3. These material parameters together
with the developed UMAT are now employed to run single element simulations. Compari-
son to the unit cell response is shown in Fig. 3 and excellent agreement can be seen for
all elements. In addition to the direct comparison of the fitted response, a series of general
load cases are investigated and compared successfully. The assumption of taking the same
time–temperature shift as for the base material is confirmed to hold.
Note the orthotropic behavior of the relaxation, too, which is fully captured by the present
approach.

5.4 Thermo-viscoelastic orthotropic composite

In this section an example is presented to study the thermal expansion relaxation. To this end
the geometrical model of the structure from the previous examples is used, but the regions of
the voids are now filled with rigid material to obtain a model composite. The reinforcement
material is taken to be quasi-rigid with a Young’s modulus of 1 × 106 MPa, the Poisson ratio
is set to 0.25. The isotropic coefficient of thermal expansion is set to 1 × 10−6 /◦ C. For the
matrix material the data from Table 1 is used, disregarding temperature dependence (i.e., no
WLF shift is applied).
As in the previous examples, the unit cell is employed to perform material characteri-
zation, i.e., to compute the effective linear thermo-viscoelastic properties of the composite.
The necessary (mechanical) load cases are simulated and the input data to the UMAT is fitted.
Again, the off-diagonal term, R12 (t), exhibits some curvature change around t = 2.0 s and
is fitted by two Prony terms. The calibrated homogenized material data is listed in Table 4,
graphical presentation is omitted.
In addition, the effective thermal expansion behavior of the composite can be predicted
by the unit cell simulations. A Heaviside temperature step, ϑ(t > 0) = 1, from an unde-
formed, stress free configuration is applied and the thermal expansion response is evaluated.
The constituents of the model composite posses different coefficients of thermal expansion,
thus, a uniform temperature change of the composite induces stresses and strains in the con-
stituents (the stresses, of course, are self-equilibrated). Since these stresses and strains show
relaxation and creep, the effective thermal expansion shows relaxation/creep, too. Note that
no thermal diffusion is involved here. Due to the orthotropic material symmetry there are
two independent thermal expansion relaxation/creep functions. The corresponding unit cell
predictions are presented in Fig. 4 (thin lines).
To obtain input data to the UEXPAN subroutine, the thermal expansion response is fitted,
and the material data is listed in Table 4. It comprises the two instantaneous coefficients of
thermal expansion and the corresponding relaxation/creep functions with one Prony term
each. The fitted behavior is presented in Fig. 4 (thick dashed lines) for comparison with the
unit cell predictions.
Mech Time-Depend Mater (2018) 22:421–433 431

Table 4 Orthotropic linear thermo-viscoelastic plane stress material data of a rigid inclusion reinforced
model composite as input to the UMAT and the UEXPAN; instantaneous elasticity matrix elements, relaxation
matrix elements, instantaneous coefficients of thermal expansion, and thermal expansion relaxation data

R11 0 = 5301 MPa R22 0 = 3857 MPa R12 0 = R21 0 = 817 MPa R33 0 = 714 MPa

k=1 r11 = 0.419 r22 = 0.413 r12 = r21 = 0.203 r33 = 0.434
τ r11 = 1.20 s τ r22 = 1.20 s τ r12 = τ r21 = 2.00 s τ r33 = 1.00 s

k=2 r11 = 0.000 r22 = 0.000 r12 = r21 = −0.034 r33 = 0.000
τ r11 = 1.20 s τ r22 = 1.20 s τ r12 = τ r21 = 0.35 s τ r33 = 1.00 s

α1 0 = 5.16 × 10−5 /◦ C α2 0 = 7.58 × 10−5 /◦ C

k = 1 a1 = −0.019 a2 = −0.045
τ a1 = 1.00 s τ a2 = 1.30 s

Fig. 4 Orthotropic thermal

expansion relaxation/creep
functions αi (t) as response to
unit temperature Heaviside step;
homogenization of a rigid
inclusion reinforced model
composite (“unit cell”), and
single element predictions by the
developed thermal expansion law
(“uexpan”) calibrated to the
homogenized behavior

So far, the material data input is calibrated and other thermo–viscoelastic load cases are
studied. Temperature changes under fully and uni-axial in-plane constraint are simulated
which give rise to thermally induced stresses. For both cases excellent agreement is found
for the unit cell response and the UMAT/UEXPAN prediction. For comparison, the same load
cases are simulated with coefficients of thermal expansion taken to be time independent,
i.e. suppressing thermal expansion relaxation. As expected, the stress responses deviate for
increasing time.

6 Summary

A constitutive material law for linear thermo-viscoelasticity in the time domain with full
anisotropy is presented. For the linear viscoelastic material behavior the anisotropic mate-
rial symmetry is not only considered for the elastic part. Also for the relaxation response,
full anisotropy is realized with the appropriate number of independent material values. The
formulation is based on a time-dependent elasticity tensor for which each element pos-
sesses its own relaxation function. This way, the mutual coupling of the normal components,
i.e., Poisson type effects, is accounted for. The relaxation function for each tensor element
is prescribed by its individual Prony series expansions.
An incremental algorithm is set up which gives the time-dependent stress response to
a prescribed incremental strain history. State dependent internal variables are utilized to
432 Mech Time-Depend Mater (2018) 22:421–433

account for the loading history, and the consistent material Jacobian is formulated. Temper-
ature dependence is realized by a time–temperature shift function of the Williams–Landel–
Ferry type. It is implemented for orthotropic material symmetry under plane stress assump-
tion into a FEM package. The extension to tri-axial stress and strain states and to general
anisotropy is straightforward.
The thermal expansion relaxation/creep is considered as inherent material behavior. It
is treated in analogy to the mechanical, i.e., viscoelastic, model. The material behavior is
described by instantaneous coefficients of thermal expansion and their individual relaxation
functions. The time-dependent thermal strain is modeled as a function of the applied temper-
ature change. Again, internal state variables are used to account for the temperature history.
Various tests on isotropic and orthotropic plane stress problems are carried out success-
fully for verification and validation. The performance is found to be very satisfactory. Or-
thotropic test cases are set up by considering heterogeneous model composites and employ-
ing a periodic unit cell type approach. Thereof, homogenized material data is computed,
serving as input to the developed material law. Also the thermal expansion relaxation be-
havior is investigated.
Finally, the developed constitutive laws provide tools for structural analyses of compo-
nents made from orthotropic linear thermo–viscoelastic materials.

Acknowledgements Open access funding provided by TU Wien (TUW). The first author thanks Florian
Toth, TU Wien, for the fruitful and inspiring discussions on viscoelasticity.

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 Inter-
national License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution,
and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source,
provide a link to the Creative Commons license, and indicate if changes were made.

References
Bažant, Z.: Delayed thermal dilatations of cement paste and concrete due to mass transport. Nucl. Eng. Des.
14(2), 308–318 (1970)
Bergström, J., Boyce, M.: Constitutive modeling of the time-dependent and cyclic loading of elastomers and
application to soft biological tissues. Mech. Mater. 33(9), 523–530 (2001)
Böhm, H.J.: A short introduction to continuum micromechanics. In: Mechanics of Microstructured Materials,
pp. 1–40. Springer, Berlin (2004)
Carcione, J.M.: Constitutive model and wave equations for linear, viscoelastic, anisotropic media. Geophysics
60(2), 537–548 (1995)
Chandra, R., Singh, S., Gupta, K.: Damping studies in fiber-reinforced composites—a review. Compos.
Struct. 46(1), 41–51 (1999)
Crochon, T., Schönherr, T., Li, C., Lévesque, M.: On finite-element implementation strategies of Schapery-
type constitutive theories. Mech. Time-Depend. Mater. 14(4), 359–387 (2010)
Dong, Z., McMechan, G.A.: 3-d viscoelastic anisotropic modeling of data from a multicomponent, multiaz-
imuth seismic experiment in northeast Texas. Geophysics 60(4), 1128–1138 (1995)
Gurtin, M.E., Sternberg, E.: On the linear theory of viscoelasticity. Arch. Ration. Mech. Anal. 11(1), 291–356
(1962)
Holzapfel, G.A.: Nonlinear Solid Mechanics, vol. 24. Wiley, Chichester (2000)
Kaliske, M.: A formulation of elasticity and viscoelasticity for fibre reinforced material at small and finite
strains. Comput. Methods Appl. Mech. Eng. 185(2), 225–243 (2000)
Lakes, R.S.: Viscoelastic Materials. Cambridge University Press, Cambridge (2009)
Lakes, R.S., Wineman, A.: On Poisson’s ratio in linearly viscoelastic solids. J. Elast. 85(1), 45–63 (2006)
Lévesque, M., Derrien, K., Baptiste, D., Gilchrist, M.D.: On the development and parameter identification of
Schapery-type constitutive theories. Mech. Time-Depend. Mater. 12(2), 95–127 (2008)
Melo, J.D.D., Radford, D.W.: Time and temperature dependence of the viscoelastic properties of peek/im7.
J. Compos. Mater. 38(20), 1815–1830 (2004)
Mech Time-Depend Mater (2018) 22:421–433 433

Muliana, A., Khan, K.A.: A time-integration algorithm for thermo-rheologically complex polymers. Comput.
Mater. Sci. 41(4), 576–588 (2008)
Nedjar, B.: An anisotropic viscoelastic fibre–matrix model at finite strains: continuum formulation and com-
putational aspects. Comput. Methods Appl. Mech. Eng. 196(6), 1745–1756 (2007)
Pandolfi, A., Manganiello, F.: A model for the human cornea: constitutive formulation and numerical analysis.
Biomech. Model. Mechanobiol. 5(4), 237–246 (2006)
Pettermann, H., Hüsing, J.: Modeling and simulation of relaxation in viscoelastic open cell materials and
structures. Int. J. Solids Struct. 49(19), 2848–2853 (2012)
Pettermann, H.E., Suresh, S.: A comprehensive unit cell model: a study of coupled effects in piezoelectric
1–3 composites. Int. J. Solids Struct. 37(39), 5447–5464 (2000)
Poon, H., Ahmad, M.: A material point time integration procedure for anisotropic, thermo rheologically
simple, viscoelastic solids. Comput. Mech. 21(3), 236–242 (1998)
Puso, M., Weiss, J.: Finite element implementation of anisotropic quasi-linear viscoelasticity using a discrete
spectrum approximation. J. Biomech. Eng. 120(1), 62–70 (1998)
Sawant, S., Muliana, A.: A thermo-mechanical viscoelastic analysis of orthotropic materials. Compos. Struct.
83(1), 61–72 (2008)
Schapery, R.A.: Application of thermodynamics to thermomechanical, fracture, and birefringent phenomena
in viscoelastic media. J. Appl. Phys. 35(5), 1451–1465 (1964)
Schapery, R.A.: On the characterization of nonlinear viscoelastic materials. Polym. Eng. Sci. 9(4), 295–310
(1969)
Schapery, R.A.: Nonlinear viscoelastic solids. Int. J. Solids Struct. 37(1), 359–366 (2000)
Sorvari, J., Hämäläinen, J.: Time integration in linear viscoelasticity—a comparative study. Mech. Time-
Depend. Mater. 14(3), 307–328 (2010)
Spencer, R., Boyer, R.: Thermal expansion and second-order transition effects in high polymers: Iii. time
effects. J. Appl. Phys. 17(5), 398–404 (1946)
Suquet, P., Bornert, M.: Rappels de calcul tensoriel et d’élasticité. Homogénéisation en mécanique des matéri-
aux 2, 171–202 (2001)
Taylor, Z.A., Comas, O., Cheng, M., Passenger, J., Hawkes, D.J., Atkinson, D., Ourselin, S.: On modelling
of anisotropic viscoelasticity for soft tissue simulation: numerical solution and GPU execution. Med.
Image Anal. 13(2), 234–244 (2009)
Zocher, M., Groves, S., Allen, D.: A three-dimensional finite element formulation for thermoviscoelastic
orthotropic media. Int. J. Numer. Methods Eng. 40(12), 2267–2288 (1997)