Finite element analysis of viscoelastic composite structures based

on a micromechanical material model

B. Kurnatowski, A. Matzenmiller
*
University of Kassel, Institute of Mechanics, D-34109 Kassel, Germany
a r t i c l e i n f o
Article history:
Received 11 January 2008
Accepted 18 February 2008
Available online xxxx
Keywords:
Generalized method of cells
Homogenization
Two-scale analysis
Composites
Viscoelasticity
a b s t r a c t
The stress and creep analysis of structures made of micro-heterogeneous composite materials is treated
as a two-scale problem, dened as a mechanical investigation on different length scales. Reinforced com-
posites show by denition a heterogeneous texture on the microlevel, determined by the constitutive
behaviour of the matrix material and the embedded bres as well as the characteristics of the bonding
properties in the interphase. All these heterogeneities are neglected by the nite element analysis of
structural elements on the macroscale, since a cticious and homogeneous continuum with averaged
properties is assumed. Therefore, the constitutive equations of the substitute material should well reect
the mechanical behaviour of the existing micro-heterogeneous composite in an average sense.
The paper at hand starts with the brief outline of a micromechanical model, named generalized method
of cells (GMC), which provides the macrostress responses due to macrostrain processes as well as the
homogenised constitutive tensor of the substitute material. The macroscopic stresses and strains are
obtained as volume averages of the corresponding microelds within a representative volume element.
The effective material tensor constitutes the mapping between the macro-strains and the macro-stresses.
The cells method is used for the homogenisation of the unidirectionally reinforced single layers of lam-
inates made of viscoelastic resins and exibly embedded elastic bres. The algorithm for the homogeni-
sation of the constitutive properties runs simultaneously to the nite element analysis at each point of
numerical integration and provides the macro-stresses and the homogenised constitutive properties.
The validity of the proposed two-scale simulation is investigated by solving boundary value problems
and comparing the numerical results for the structures to the experimental data of creep and relaxation
tests or analytical solutions.
2008 Elsevier B.V. All rights reserved.
1. Introduction
A diversity of composite materials is nowadays available and in
use for the design of structural elements. The effective and macro-
scopically observable material behaviour, the knowledge of which
is required for the mechanical analyses, depends on the mechani-
cal properties of the constituting phase materials, the quality of
the bonding as well as the geometrical arrangements of the bres
and the volume fractions. It is the aim in this paper to predict the
effective properties not on a purely phenomenological basis for
each and every composite build-up, but to compute them from
the generally well-known properties of the constituents, the bond-
ing characteristics and the geometrical setup of the reinforcing
bres.
The paper is focused on the two-scale simulation of viscoelastic
composite structures with viscoelastically embedded bres. While
the so-called macrolevel of the structural component is modeled
by the nite element method, the microscale behaviour of the het-
erogeneous bre-reinforced material is analysed by means of the
generalized method of cells (GMC), see, e.g. Paley and Aboudi [19].
The homogenisation of composite materials with a viscoelastic re-
sin has already been discussed by Aboudi [2] for nite deforma-
tions and for rate dependent phases with a compliant bre
matrix bond by Matzenmiller and Gerlach [14] for small strains
both on the basis of the GMC. The present paper uses the efciently
reformulated cells method of Pindera and Bednarcyk [20], which
provides numerically the same results as the original version, but
reduces the computational effort noteworthily. Especially the
mathematical formulation of the viscoelastic bonding simplies
in comparison to Matzenmiller and Gerlach [14], since the stresses
are taken as the primary unknowns in this new approach. Orozco
and Pindera [18] applied the reformulated GMC to composites
made of viscoelastic matrix materials with rigid brematrix
bonding. The identication of the viscoelastic model parameters
in conjunction with the cells method is the subject in Matzenmiller
and Gerlach [15] and Gerlach [6]. Recently, a rened version of the
method of cells has been developed, named as high-delity-GMC
0927-0256/$ - see front matter 2008 Elsevier B.V. All rights reserved.
doi:10.1016/j.commatsci.2008.02.026
* Corresponding author.
E-mail address: amat@ifm.maschinenbau.uni-kassel.de (A. Matzenmiller).
Computational Materials Science xxx (2008) xxxxxx
Contents lists available at ScienceDirect
Computational Materials Science
j our nal homepage: www. el sevi er . com/ l ocat e/ commat sci
ARTICLE IN PRESS
Please cite this article in press as: B. Kurnatowski, A. Matzenmiller, Comput. Mater. Sci. (2008), doi:10.1016/j.commatsci.2008.02.026
(HFGMC), for details see, e.g. Bansal and Pindera [3]. While the
GMC applies a linear expansion of the microscopic displacement
eld, the more sophisticated HFGMC approach uses higher order
polynomials instead and is furthermore rested upon the idea of
perturbation methods, see Nayfeh [17].
The numerical efciency of the reformulated GMC approach to-
gether with its sufcient accuracy recommends the method of cells
to be used as the constitutive model within the framework of a
two-scale simulation of composite structures. The notorious in-
crease of the computational effort, which lies in the nature of mul-
tiscale computations, is kept at a decent level, compared to the
nite element based micromechanical modelling. Hence, the
simultaneous two-scale analysis of engineering structures, requir-
ing even larger nite element meshes, becomes feasible.
In this paper, the interphase as the transition zone from the -
bre to the pure matrix phase is assumed to be viscoelastic and is
described mathematically by an interface model as it was similarly
done by Hashin [9]. While the compliance in normal direction is
assumed to be elastic, the tangential exibility shows viscoelastic
behaviour. The time-dependent shear modulus of the matrix mate-
rial is taken as a function of time and is approximated by a nite
DirichletProny series. The bulk modulus is considered to be
time-independent. The bres possess elastically isotropic or trans-
versely isotropic properties. The equations of the micromechanical
model are numerically integrated by the help of a recursive inte-
gration scheme presented by Taylor et al. [24].
The paper starts with the discussion of the reformulated cells
method for viscoelastic composites. For the macrostructural analy-
ses the nite element analysis program (FEAP
1
) is used. The link of
the GMC with the nite element method is described briey. Finally,
the validity of the model is investigated by comparing numerical re-
sults with experimental data or analytical solutions for four
problems.
2. Micromechanical constitutive model
The micromechanical analysis is rested upon the concept of the
so-called representative volume element (RVE), see Hashin [8]. By
this is meant an adequately chosen sample of the heterogeneous
material representing the microstructure as a whole. The solution
of the appropriately formulated boundary value problem on the
domain X
RVE
provides the microelds of stresses and strains within
the RVE. The volume averaged eld quantities then dene the mac-
roscopic stress and strain tensors at a material point of the macro-
structure. The effective material tensor serves as the mapping
between the averages of stress and strain. The method of cells dis-
cretises the representative volume element X
RVE
by subdividing it
into the total of N
b
N
c
rectangular subdomains X
(bc)
X
RVE
which
are referred to as the subcells, see Fig. 1. The thickness of the RVE in
bre direction x
1
is set to unit length, while the transversal x
2
; x
3
-
plane has the cross-section h l. The microscopic displacement
eld u(x), x X
RVE
inside the RVE is approximated piecewisely
by linear functions of the local coordinate systems x
(bc)
i
located at
the centre of each X
(bc)
:
u
(bc)
i
= w
(bc)
i
U
(bc)
i
x
(bc)
2
W
(bc)
i
x
(bc)
3
; i = 1; 2; 3 (1)
with the displacements w
(bc)
i
at the origin O
(bc)
and the microvari-
ables U
(bc)
i
and W
(bc)
i
. The average innitesimal strain tensors
(bc)
)
of the subcells are computed from the subcell displacement elds:

(bc)
ij
) =
1
V
X
(bc)
_
V
X
(bc)

(bc)
ij
dV =
1
V
X
(bc)
_
V
X
(bc)
1
2
(o
i
u
(bc)
j
o
j
u
(bc)
i
) dV (2)
Two subcells X
(bc)
and X
(
^
bc)
with
^
b = b 1 are attached to each
other in global x
2
-direction at the subcell surface
(2)
oX
(bc)
with the
outward unit-normal vector
(2)
e
(bc)
n
|e
2
. The differences of the sur-
face averaged displacement elds u
(
^
bc)
i
and u
(bc)
i
dene the coordi-
nates of the displacement discontinuity vector
(2)
sut
(bc)
at the
interface
(2)
oX
(bc)
for j = 1; 2; 3:
(2)
su
j
t
(bc)
=
_
lc=2
lc=2
_
1=2
1=2
u
(
^
bc)
j
[
x
(
^
b)
2
=h^
b
=2
dx
1
dx
(
^
bc)
3

_
lc=2
lc=2
_
1=2
1=2
u
(bc)
j
[
x
(b)
2
=h
b
=2
dx
1
dx
(bc)
3
(3)
Analogously, it holds for interfaces
(3)
oX
(bc)
with the positive out-
ward normal
(3)
e
(bc)
|e
3
and ^c = c 1:
(3)
su
j
t
(bc)
=
_
h
b
=2
h
b
=2
_
1=2
1=2
u
(b^c)
j
[
x
(^c)
3
=l^c
=2
dx
1
dx
(b^c)
2

_
h
b
=2
h
b
=2
_
1=2
1=2
u
(bc)
j
[
x
(c)
3
=lc=2
dx
1
dx
(bc)
2
(4)
The boundary conditions imposed on the discrete cells model are
given in terms of the RVE-surface displacements u
i
(x) =
ij
)x
j
. They
lead together with the displacement jump denitions (3) and (4) to
the following system of equations:

(bc)
11
) =
11
); b = 1; N
b
; c = 1; N
c
(5)

N
b
b=1
(h
b

(bc)
22
)
(2)
su
2
t
(bc)
) = h
22
); c = 1; N
c
(6)

Nc
c=1
(l
c

(bc)
33
)
(3)
su
3
t
(bc)
) = l
33
); b = 1; N
b
(7)

N
b
b=1
(h
b

(bc)
12
)
(2)
su
1
t
(bc)
) = h
12
); c = 1; N
c
(8)

Nc
c=1
(l
c

(bc)
13
)
(3)
su
1
t
(bc)
) = l
13
); b = 1; N
b
(9)

N
b
b=1

Nc
c=1
(h
b
l
c

(bc)
23
) l
c
(2)
su
3
t
(bc)
h
b
(3)
su
2
t
(bc)
) = hl
23
) (10)
The ranges of the subcell indices are b = 1; 2; . . . ; N
b
=: 1; N
b
and
c = 1; 2; . . . ; N
c
=: 1; N
c
. Pindera and Bednarcyk [20] dene the fol-
lowing components T
(:)
ij
of tractions on the subcell planes with the
normal n|e
i
on the basis of the original stress continuity conditions
given by Aboudi [1]:
T
(bc)
11
:= r
(bc)
11
); b = 1; N
b
; c = 1; N
c
(11)
T
(c)
21
:= r
(1c)
21
) = r
(2c)
21
) = . . . = r
(N
b
c)
21
); c = 1; N
c
(12)
T
(c)
22
:= r
(1c)
22
) = r
(2c)
22
) = . . . = r
(N
b
c)
22
); c = 1; N
c
(13)
T
(c)
23
:= r
(1c)
23
) = r
(2c)
23
) = . . . = r
(N
b
c)
23
); c = 1; N
c
(14)
T
(b)
33
:= r
(b1)
33
) = r
(b2)
33
) = . . . = r
(bNc)
33
); b = 1; N
b
(15)
T
(b)
31
:= r
(b1)
31
) = r
(b2)
31
) = . . . = r
(bNc)
31
); b = 1; N
b
(16)
T
(b)
32
:= r
(b1)
32
) = r
(b2)
32
) = . . . = r
(bNc)
32
); b = 1; N
b
(17)
where r
(bc)
) denotes the volume average of the subcell stress tensor
r
(bc)
. Due to the complementary property of shear the symmetry of
the stress tensor r) holds in classical continuum mechanics and
Eqs. (14) and (17) provide: T
(b)
23
equals T
(c)
32
for all arbitrary pairings
of b and c, i.e. the transversal shear stress eld r
23
(x), x X
RVE
is
homogeneous and quantied by T
23
:= T
(b)
23
= T
(c)
32
as spatially con-
stant for all subcells in the RVE. For the sake of a unied notation,
the normal stresses r
(bc)
11
) are renamed as T
(bc)
11
. The stresses T
(:)
ij
might be interpreted as the tractions on the surface of the discre-
tised volume element X
RVE
. Up to this point, the model equations
of the GMC are independent from the mechanical properties of 1
URL: http://www.ce.berkeley.edu/~rlt/feap/.
2 B. Kurnatowski, A. Matzenmiller / Computational Materials Science xxx (2008) xxxxxx
ARTICLE IN PRESS
Please cite this article in press as: B. Kurnatowski, A. Matzenmiller, Comput. Mater. Sci. (2008), doi:10.1016/j.commatsci.2008.02.026
the phase materials. In order to solve for the displacement Eqs. (5)
(10) together with the stress continuity conditions (11)(17), the
constitutive equations for the phase materials and the brematrix
bond must be considered additionally as discussed next.
2.1. Linear viscoelasticity for phase materials
The original method of cells is a displacement-based microme-
chanical approach, where the set of equilibrium, kinematical and
material equations are solved for the strains. Thus, the relaxation
integrals are used for linear rate-dependent composites as the con-
stitutive equations. However, for the efciently reformulated cells
method see Bednarcyk and Pindera [4] the micromechanical
equations are solved for the stresses as the primary variables.
Therefore, the strains are replaced in favour of the stresses by
means of the creep integrals in the case of time-dependent mate-
rial behaviour. This is in contrast to the approach in Matzenmiller
and Gerlach [14]. In the theory of linear viscoelasticity, the consti-
tutive equations can be written as convolution integrals. The visco-
elastic strain response (t) at present time t is then given by
(t) =
_
t
s=
S(t s) :
dr(s)
ds
ds (18)
Hence, the actual strain tensor (t) is a linear functional of the stress
history r(s)
6s6t
in past times s and the time dependent compli-
ance tensor S(t) of fourth order contracted twice (:) with the rate
of the stress tensor. The stress and strain tensors can be decom-
posed into the volumetric parts e
V
=
1
3
tr()1 and s
V
=
1
3
tr(r)1 and
the deviatoric tensors e
D
= e
V
and s
D
= r s
V
, respectively.
Assuming isotropic symmetry and a rate-independent bulk modu-
lus K, the convolution integral (18) merges with tr =
1
3K
trr into:
(t) =
1
9K
trr1
_
t

J(t s)
ds
D
(s)
ds
ds (19)
where J(t) is the creep function, describing the materials response
due to a unit stress jump. On the basis of the generalized Kelvin-
model, the function J(t) can be specied by the nite Dirichlet
Prony-series, factoring it into the retardation strengths J
k
and the
exponentials of the creep times ^s
k
:
J(t) = J

N
k=1
J
k
e

t
^s
k ; J

:= J
0

N
k=1
J
k
(20)
The (2N 1)-parameter model (20) is uniquely dened by the dis-
crete creep spectrum L(s) : J
k
; ^s
k
; k = 1; . . . ; N. The parameter
J

= 1=(2l

) is given by the long time shear modulus G

= l

at
equilibrium. It is further on assumed that the stress history is con-
stantly equal to zero for times s < 0. With the ansatz of the kernel
function in (20) and the contribution factors ^m
k
= J
k
=J

the solution
of the convolution integral in Eq. (19) can be approximated by the
numerical time integration algorithm of Taylor et al. [24]. Applying
the recursive integration scheme at t = t
n1
= t
n
Dt to the convo-
lution of the rate of the average stress deviator
_
s
(bc)
D
) with the kernel
J(t)
(bc)
for the subcell X
(bc)
it holds:

n1
e
(bc)
D
) =
1
2l
(bc)

_
tn
0
1

N
k=1
^m
(bc)
k
e

tnDts
^s
(bc)
k
_ _
ds
(bc)
D
(s))
ds
ds
_

_
tnDt
tn
1

N
k=1
^m
(bc)
k
e

tnDts
^s
(bc)
k
_ _
ds
(bc)
D
(s))
ds
ds
_
(21)
Hence, the original integral
_
tnDt
0
() is decomposed into a rst part,
representing the stress deviator
n
s
D
) =
_
tn
0
() at the previous time
step t = t
n
, and an incremental part D
n1
s
D
) =
_
tnDt
tn
(). After some
elementary operations, Eq. (21) becomes:

n1
e
(bc)
D
) =
1
2l
(bc)

n1
s
(bc)
D
)

N
k=1
e

tnDt
^s
(bc)
k ^m
(bc)
k
_
tn
0
e
s
^s
(bc)
k
ds
(bc)
D
(s))
ds
ds
_

N
k=1
e

tnDt
^s
(bc)
k ^m
(bc)
k
_
tnDt
tn
e
s
^s
(bc)
k
ds
(bc)
D
(s))
ds
ds
_
(22)
With the help of the denitions
n
d
(bc)
k
:= e

tn
^s
(bc)
k
_
tn
0
e
s
^s
(bc)
k
ds
(bc)
D
(s))
ds
ds (23)
Fig. 1. Representative volume element of a periodical microstructure. Discrete GMC-model of the RVE with subcell X
(bc)
.
B. Kurnatowski, A. Matzenmiller / Computational Materials Science xxx (2008) xxxxxx 3
ARTICLE IN PRESS
Please cite this article in press as: B. Kurnatowski, A. Matzenmiller, Comput. Mater. Sci. (2008), doi:10.1016/j.commatsci.2008.02.026
D
n1
d
(bc)
k
:= e

tnDt
^s
(bc)
k
_
tnDt
tn
e
s
^s
(bc)
k
ds
(bc)
D
(s))
ds
ds (24)
and the initial condition
0
d
(bc)
k
=
0
s
(bc)
D
) (25)
Eq. (22) merges to:

n1
e
(bc)
D
) =
1
2l
(bc)

n1
s
(bc)
D
)

N
k=1
^m
(bc)
k
e

Dt
^s
(bc)
k
n
d
(bc)
k
D
n1
d
(bc)
k
_ _ _ _
(26)
If the time derivative of the stress deviator is approximately re-
placed by the nite difference quotient, i.e.:
ds
D
(t))
dt
~
1
Dt
(
n1
s
D
)
n
s
D
)); (27)
the integral in Eq. (24) can be solved analytically to get:
D
n1
d
(bc)
k
=
^s
(bc)
k
Dt
1 e

Dt
^s
(bc)
k
_ _
(
n1
s
(bc)
D
)
n
s
(bc)
D
)) (28)
Hence, the actual subcell tensors
n1
d
(bc)
k
at time
n1
t are found by
their values
n
d
(bc)
k
at the previous point in time
n
t and their incre-
ments in Eq. (28) with the help of a highly efcient recursion
formula:
n1
d
(bc)
k
= e

Dt
^s
(bc)
k
n
d
(bc)
k
D
n1
d
(bc)
k
(29)
with which the total present strain tensor amounts to:

n1

(bc)
) =
1
9K
(bc)
tr(
n1
r
(bc)
))1

1
2l
(bc)

n1
s
(bc)
D
)

N
k=1
^m
(bc)
k
n1
d
(bc)
k
_ _
(30)
Eq. (30) might also be written in terms of the total stress increment
D
n1
r
(bc)
) and the average strain history tensor
n
y
(bc)
):

n1

(bc)
) =

S
(bc)
ve
: D
n1
r
(bc)
)
n
y
(bc)
) (31)
Herein, the instantaneous viscoelastic compliance tensor

S
(bc)
ve
is de-
ned as

S
(bc)
ve
=
1
9K
(bc)
1 1

1
2l
(bc)

N
k=1
^m
(bc)
k
^s
(bc)
k
Dt
1 e

Dt
^s
(bc)
k
_ _ _ _
II
1
3
1 1
_ _
; (32)
by using the fourth order tensor II =
1
2
(d
ik
d
jl
d
il
d
kj
)e
i
e
j
e
k
e
l
.
The strain history tensors
n
y
(bc)
) are computed from:

n
y
(bc)
) =
1
9K
(bc)
tr(
n
r
(bc)
))1

1
2l
(bc)

n
s
(bc)
D
)

N
k=1
^m
(bc)
k
e

Dt
^s
(bc)
k
n
d
(bc)
k
_ _
(33)
The constitutive description of linearly elastic material behaviour,
which is assumed for the bres, is also covered by Eqs. (31) and
(33) as the special case with ^m
k
= 0 for all k = 1; . . . ; N and
l

= l
0
. In order to write the constitutive Eq. (31) in vector-matrix
notation, the following stress and strain vectors are stipulated:
T
(bc)
:=
T
(bc)
11
T
(c)
22
T
(b)
33
T
23
T
(b)
31
T
(c)
12
_

_
_

_
=
r
(bc)
11
)
r
(bc)
22
)
r
(bc)
33
)
r
(bc)
23
)
r
(bc)
31
)
r
(bc)
12
)
_

_
_

_
= r
(bc)
);
(bc)
) :=

(bc)
11
)

(bc)
22
)

(bc)
33
)

(bc)
23
)

(bc)
31
)

(bc)
12
)
_

_
_

_
(34)
The denition (34)
1
makes use of the stress continuity conditions
given in Eqs. (11)(17). Now, the constitutive Eq. (31) read in matrix
notation as

n1

(bc)
) =

S
(bc)
ve
D
n1
T
(bc)

n
y
(bc)
) (35)
with
n
y
(bc)
) analogously dened as
(bc)
) in Eq. (34)
2
and

S
(bc)
ve
as the
matrix form of the constitutive tensor in Eq. (32). It is mentioned
that by writing down the constitutive equations for the subcells
in the manner of Eqs. (35), the stress continuity conditions (11)
(17) of the cells model will be inevitably fullled by the mathemat-
ical solution. The only remaining equations of the model are the dis-
placement continuity conditions (5)(10) as well as the constitutive
equations of the viscoelastic interface model, which will be given
briey in the following Section 2.2. In Eq. (35) above the matrix

S
(bc)
ve
=
1
9K
(bc)
1 1 1 0 0 0
1 1 1 0 0 0
1 1 1 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
_

_
_

1
2G
(bc)
ve
2=3 1=3 1=3 0 0 0
1=3 2=3 1=3 0 0 0
1=3 1=3 2=3 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
_

_
_

_
(36)
stores the components of the instantaneous viscoelastic compliance
tensor for the phases with the time-dependent isotropic shear com-
pliance factor
1
G
(bc)
ve
=
1
l
(bc)

N
k=1
^m
(bc)
k
^s
(bc)
k
Dt
1 e

Dt
^s
(bc)
k
_ _ _ _
: (37)
The strain history vector
n
y
(bc)
) follows accordingly to Eq. (33):

n
y
(bc)
) =
1
3K
(bc)
n
T
(bc)
V

1
2l
(bc)

n
T
(bc)
D

N
k=1
^m
(bc)
k
e

Dt
^s
(bc)
k
n
d
(bc)
k
_ _
(38)
The herein used vector notations of the stress deviator and the vol-
umetric stress tensor are dened with T
(bc)
kk
:= T
(bc)
11
T
(c)
22
T
(b)
33
as
s
(bc)
D
) T
(bc)
D
=
T
(bc)
11
T
(bc)
kk
=3
T
(c)
22
T
(bc)
kk
=3
T
(b)
33
T
(bc)
kk
=3
T
23
T
(b)
31
T
(c)
12
_

_
_

_
and s
(bc)
V
) T
V
=
T
(bc)
kk
=3
T
(bc)
kk
=3
T
(bc)
kk
=3
0
0
0
_

_
_

_
(39)
Eqs. (25), (28) and (29) read in vector notation as
0
d
(bc)
k
=
0
T
(bc)
D
(40)
D
n1
d
(bc)
k
=
^s
(bc)
k
Dt
1 e

Dt
^s
(bc)
k
_ _
(
n1
T
(bc)
D

n
T
(bc)
D
) (41)
n1
d
(bc)
k
= e

Dt
^s
(bc)
k
n
d
(bc)
k
D
n1
d
(bc)
k
(42)
2.2. Viscoelastic interface model
Similar to the constitutive relations for the strains in Eq. (18)
the actual interface displacement jump vector su(t)t is determined
from a linear functional over the interface creep function R(t) and
the time history of the traction vector t on the interface see Mat-
zenmiller and Gerlach [14]:
4 B. Kurnatowski, A. Matzenmiller / Computational Materials Science xxx (2008) xxxxxx
ARTICLE IN PRESS
Please cite this article in press as: B. Kurnatowski, A. Matzenmiller, Comput. Mater. Sci. (2008), doi:10.1016/j.commatsci.2008.02.026
su(t)t =
_
t

R(t s)
dt
ds
ds (43)
As before, the second order tensor of the creep functions R(t) is
specied by a nite DirichletProny-series. Assuming that only
the tangential components of sut are explicitly time-dependent
and the normal one is elastic, the elements of R(t) might be given as
R
ij
(t) = d
i1
R
n
d
ij
(1 d
i1
)R

(d
i2
d
i3
)

N
k=1
R
k
e
t=n
k
_ _
(44)
with the creep times n
k
, the spectrum strengths R
k
and R
n
as the
elastic compliance in normal direction. Hence, the initial tangential
compliance R
t
follows as R
t
= R

N
k=1
R
k
. Inserting the functions
(44) into Eq. (43) and applying the integration scheme of Taylor
et al. [24] to the resulting convolution integral, the vector of the dis-
placement discontinuity
n1(i)
u
(bc)
follows similarly to Eq. (35) for
each brematrix interface
(i)
oX
(bc)
of the cells model:
n1(i)
u
(bc)
=
(i)
R
(bc)
D
n1(i)
t
(bc)

n(i)
v
(bc)
(45)
Herein,
(i)
R
(bc)
means the instantaneous interface compliance ma-
trix, assumed as
(i)
R
(bc)
=
(i)
R
(bc)
nn
0 0
0
(i)
R
(bc)
tt
0
0 0
(i)
R
(bc)
bb
_

_
_

_
(46)
with
(i)
R
(bc)
tt
=
(i)
R
(bc)
bb
=
(i)
R
(bc)

N
k=1
(i)
R
(bc)
k
(i)
n
(bc)
k
Dt
(1 e

Dt
(i)
n
(bc)
k ) (47)
and the elastic compliance
(i)
R
(bc)
nn
= R
n
in normal direction of the
interface
(i)
oX
(bc)
. The local n; t; b-coordinates in normal, tangential
and binormal direction to the interface
(i)
u
(bc)
k
=
(i)
su
i
t
(bc)
e
i

(i)
e
(bc)
k
of the jump vectors
(i)
sut
(bc)
are compiled as
(2)
u
(bc)
=
(2)
su
2
t
(bc)
(2)
su
3
t
(bc)
(2)
su
1
t
(bc)
_

_
_

_ and
(3)
u
(bc)
=
(3)
su
3
t
(bc)
(3)
su
1
t
(bc)
(3)
su
2
t
(bc)
_

_
_

_ (48)
for the bonding zone with the unit normal vectors
(2)
e
(bc)
n
|e
2
and
(3)
e
(bc)
n
|e
3
, respectively see Fig. 1. The coordinates
D
(i)
t
(bc)
k
= D
(i)
t
(bc)

(i)
e
(bc)
k
of the incremental interfacial stress vectors
D
(i)
t
(bc)
= Dr
(bc)
)
(i)
e
(bc)
n
are obtained from the increments of the sub-
cell stress tensors Dr
(bc)
) respecting the denition (34)
1
:
D
(2)
t
(bc)
=
DT
(c)
22
DT
23
DT
(c)
21
_

_
_

_ and D
(3)
t
(bc)
=
DT
(b)
33
DT
(b)
31
DT
23
_

_
_

_: (49)
The time history of the displacement jump vector
n(i)
v
(bc)
in Eq. (45)
is computed to
n
(i)
v
(bc)
= R
n
n(i)
t
(bc)
n
0
0
_

_
_

_ R

0
n(i)
t
(bc)
t
n(i)
t
(bc)
b
_

_
_

N
k=1
R
k
e

Dt
(i)
n
(bc)
k
n(i)
q
(bc)
k
(50)
with the initial condition
0(i)
q
(bc)
k
= [0
0(i)
t
(bc)
t
0
(i)
t
(bc)
b
[
T
; (51)
the updated value
n1(i)
q
(bc)
k
= e

Dt
(i)
n
(bc)
k
n(i)
q
(bc)
k
D
n1(i)
q
(bc)
k
(52)
and the increment from time t
n
to t
n1
.
D
n1(i)
q
(bc)
k
=
(i)
n
(bc)
k
Dt
1 e

Dt
(i)
n
(bc)
k
_ _
0 0 0
0 1 0
0 0 1
_

_
_

_D
n1(i)
t
(bc)
(53)
The terms in Eqs. (51)(53) are analogously derived as those in Eqs.
(40)(42).
2.3. Constitutive equation of homogenised viscoelastic composite
As the GMC has to provide the constitutive relation for the nite
element analysis of composite structures, the effective constitutive
equation of the homogenised continuum should be derived in the
following incremental vector-matrix representation for the step-
wise time integration of the viscoelastic model:

n1
r) =

C
+
D
n1
)
n
v
r
) (54)
Here, the total macrostress
n1
r) at a point in time, say t = t
n1
,
depends on the effective viscoelastic stiffness matrix

C
+
, the
macrostrain increment D
n1
) and on the not as yet specied vol-
ume averaged history vector similar to the notation in Eq. (34)
1

of the macrostress tensor

n
v
r
). The incremental constitutive Eq.
(54) of the homogeneous substitute material is derived in three
steps. First, the primary displacement variables
(bc)
ij
) and
(i)
sut
(bc)
on the left hand sides of the continuity conditions (5)(10) are
replaced by means of the right hand sides of the expressions in
Eq. (35):

n1

(bc)
) =

S
(bc)
ve
D
n1
T
(bc)

n
y
(bc)
) (55)
and (45):
n1(i)
u
(bc)
=
(i)
R
(bc)
D
n1(i)
t
(bc)

n(i)
v
(bc)
(56)
considering Eq. (49) additionally. Eqs. (55) and (56) were obtained
above as the incremental constitutive Eqs. (35) and (45) of all
subcells and interfaces. Thus, the stress increments DT
(bc)
become
the unknowns and will be computed below from the continuity
conditions, containing the instantaneous subcell compliance ten-
sors in matrix form

S
(bc)
ve
, the subcell strain history vector
n
y
(bc)
),
the instantaneous interface compliance
(i)
R
(bc)
, and the history in
vector notation
n(i)
v
(bc)
of the displacement jumps across the
interfaces.
Next, the reformulated displacement continuity conditions will
be solved for the microstress increments D
n1
T
(bc)
. Since the micro-
scopic normal and shear stress elds inside the RVE are decoupled
from each other in the framework of the GMC, the normal and
shear stress components of D
n1
T
(bc)
can be computed separately
as it is carried out in Sections 2.3.1 and 2.3.2, respectively. The ex-
plicit form of the solution for the stresses renders a very effective
numerical scheme for the two-scale analysis of composite struc-
tures in contrast to the elaborate matrix inversion for a general,
coupled, linear system of algebraic equations for all subcell
stresses.
At last, by averaging the increments D
n1
T
(bc)
over the RVE and
numerical time-integration of the resulting macrostress incre-
ments, the effective total stresses and the homogenised instanta-
neous stiffness tensor is obtained. Due to the notorious lack of
normal-shear coupling of the GMC-approach, the homogenised
material stiffness tensor in matrix form will possess the following
simple structure:

C
+
=

C
+
N
0
0

C
+
S
_
_
_
_
(57)
derived in the following subchapters.
Only Eqs. (5)(10), describing the longitudinal deformations
with the three normal strain components, lead to a coupled system
of algebraic equations, the solution of which yields the submatrix

C
+
N
. The non-zero elements of the diagonal submatrix

C
+
S
are com-
puted individually from the analytical solutions of each shear
deformation continuity condition (8)(10), see Pindera and
B. Kurnatowski, A. Matzenmiller / Computational Materials Science xxx (2008) xxxxxx 5
ARTICLE IN PRESS
Please cite this article in press as: B. Kurnatowski, A. Matzenmiller, Comput. Mater. Sci. (2008), doi:10.1016/j.commatsci.2008.02.026
Bednarcyk [20]. The high numerical efciency of the elaborate two-
scale approach, used in this paper, is caused by the shift from the
strain based original formulation to the stress based version of
the GMC. It is obvious from the stress continuity conditions in
Eqs. (11)(17) that the number of unknown independent subcell
stress components T
(:)
ij
is signicantly less than the total number
of all subcell stresses r
(bc)
) which equals the count of unknown
subcell strains
(bc)
), occurring in the original equations (5)(10).
2.3.1. Solution of the micromechanical equations for normal stresses
To begin with, the displacement conditions 6 in bre direction
x
1
are written together with the rst of Eqs. (55) for unidirection-
ally reinforced composites as

n1

11
) =
n1

(bc)
11
)
=

S
(bc)
11
D
n1
T
(bc)
11

S
(bc)
12
D
n1
T
(c)
22

S
(bc)
13
D
n1
T
(b)
33

n
y
(bc)
1
) (58)
The (N
b
N
c
) Eqs. (58) allow for the computation of the normal
stress increments D
n1
T
(bc)
11
of each subcell X
(bc)
:
D
n1
T
(bc)
11
=

n1

11
)
n
y
(bc)
1
)

S
(bc)
11

S
(bc)
12

S
(bc)
11
D
n1
T
(c)
22

S
(bc)
13

S
(bc)
11
D
n1
T
(b)
33
(59)
Hence, the normal stress increments D
n1
T
(bc)
11
in longitudinal bre
direction are given by the total macrostrain
n1

11
) at time
n1
t,
the strain history variable
n
y
(bc)
1
) for the 11-component at the pre-
vious time step
n
t and the normal stress increments D
n1
T
(c)
22
and
D
n1
T
(b)
33
in the transversal x
2
; x
3
-plane. Eqs. (59) are summarised
in matrix notation for the solution of the longitudinal stress compo-
nents of all subcells in the RVE as
D
n1
T
11
= K
11

n1

11
) S
12
D
n1
T
22
S
13
D
n1
T
33

n
v
11
(60)
The matrices used to write down Eq. (60) are dened in the
Appendix.
Next, the continuity equations (5) in x
2
-direction can be refor-
mulated such that no kinematical variables do appear explicitly
by replacing the strains
n1

(bc)
22
) with the help of Eq. (55) and the
displacement jumps in normal direction
(2)
s
n1
u
2
t
(bc)
on the basis
of Eqs. (56) together with Eq. (49). The N
c
reformulated continuity
Eqs. (5) read then:

N
b
b=1
h
b
[

S
(bc)
21
D
n1
T
(bc)
11

S
(bc)
22
D
n1
T
(c)
22

S
(bc)
23
D
n1
T
(b)
33

n
y
(bc)
2
)[

N
b
b=1
[
(2)
R
(bc)
nn
D
n1
T
(c)
22

n(2)
v
(bc)
n
[
= h
n1

22
); c = 1; N
c
(61)
Combining the increments D
n1
T
(bc)
11
of the normal stress in direction
of the bres in Eqs. (59) with Eq. (61) leads to

N
b
b=1
h
b

S
(bc)
21

n1

11
)
n
y
(bc)
1
)

S
(bc)
11

S
(bc)
12

S
(bc)
11
D
n1
T
(c)
22

S
(bc)
13

S
(bc)
11
D
n1
T
(b)
33
_ _ _

S
(bc)
22
D
n1
T
(c)
22

S
(bc)
23
D
n1
T
(b)
33

n
y
(bc)
2
)
_

N
b
b=1
(
(2)
R
(bc)
nn
D
n1
T
(c)
22

n(2)
v
(bc)
n
) = h
n1

22
); c = 1; N
c
(62)
By gathering in Eq. (62) all coefcients in front of D
n1
T
(c)
22
and
D
n1
T
(b)
33
on the left hand side and dividing it by the total hight h,
Eq. (62) simplies to:
S
(c)
22
D
n1
T
(c)
22

N
b
b=1
S
(bc)
23
D
n1
T
(b)
33
=
n1

22
)
n
v
(c)
22
S
(c)
21

n1

11
)
n
v
(c)
21
;
c = 1; N
c
(63)
Herein, the premultipliers, superscribed with a bar, are dened as
S
(c)
22
:=
1
h

N
b
b=1
h
b

S
(bc)
22

S
(bc)
21

S
(bc)
12

S
(bc)
11
_ _

(2)
R
(bc)
nn
_ _
S
(bc)
23
:=
h
b
h

S
(bc)
23

S
(bc)
21

S
(bc)
13

S
(bc)
11
_ _
; S
(c)
21
:=
1
h

N
b
b=1
h
b

S
(bc)
21

S
(bc)
11
n
v
(c)
22
:=
1
h

N
b
b=1
[h
b

n
y
(bc)
2
)
n(2)
v
(bc)
n
[;
n
v
(c)
21
:=
1
h

N
b
b=1
h
b

S
(bc)
21

S
(bc)
11

n
y
(bc)
1
)
(64)
Analogous considerations lead to the reformulation of the continu-
ity Eqs. (7) in x
3
-direction with b = 1; N
b
:
S
(b)
33
D
n1
T
(b)
33

Nc
c=1
S
(bc)
32
D
n1
T
(c)
22
=
n1

33
)
n
v
(b)
33
S
(b)
31

n1

11
)
n
v
(b)
31
;
(65)
where the auxiliary terms are dened as
S
(b)
33
:=
1
l

Nc
c=1
l
c

S
(bc)
33

S
(bc)
31

S
(bc)
13

S
(bc)
11
_ _

(3)
R
(bc)
nn
_ _
S
(bc)
32
:=
l
c
l

S
(bc)
32

S
(bc)
31

S
(bc)
12

S
(bc)
11
_ _
; S
(b)
31
:=
1
l

Nc
c=1
l
c

S
(bc)
31

S
(bc)
11
n
v
(b)
33
:=
1
l

Nc
c=1
[l
c

n
y
(bc)
3
)
n(3)
v
(bc)
n
[;
n
v
(b)
31
:=
1
l

Nc
c=1
l
c

S
(bc)
31

S
(bc)
11

n
y
(bc)
1
)
(66)
The (N
b
N
c
) Eqs. (63) and (65) are collected in matrix notation:
S
22
S
23
S
32
S
33
_ _
D
n1
T
22
D
n1
T
33
_ _
=
K
21
K
31
_ _

n1

11
)

n1

22
)

n1

33
)
_

_
_

n
v

22
n
v

33
_ _
(67)
The matrices, used to assemble Eq. (67), are given explicitly in the
Appendix. The solution of the linear system of equations (67) for
the stress increments D
n1
T
22
and D
n1
T
33
is combined with Eq.
(60) to get the nal result in explicit form:
D
n1
T
11
D
n1
T
22
D
n1
T
33
_

_
_

_
=

C
11

C
12

C
13

C
21

C
22

C
23

C
31

C
32

C
33
_

_
_

n1

11
)

n1

22
)

n1

33
)
_

_
_

n
^v
T
11
n
^v
T
22
n
^v
T
33
_

_
_

_
(68)
The stiffness coefcients provide the submatrix

C
+
N
in Eq. (57). The
increments of the normal stresses in Eq. (68) are obtained after sev-
eral matrix operations briey outlined in the Appendix. The vectors

C
ij
= C
()
ij
contain the average stiffness components C
()
ij
of all col-
umns c and all rows b of the discrete RVE.
2.3.2. Solution of micromechanical equations for shear stresses
Since the three shear-strain elds are decoupled from each
other and from the normal strains, the shear-stress components
of the macrostress tensor can be derived analytically from the
shear continuity equations. Combining the constitutive Eq. (35)
for the axial shear strains
(bc)
12
) and Eqs. (45) for the discontinu-
ities
(2)
su
1
t
(bc)
with the continuity Eq. (8) leads to:

N
b
b=1
h
b
(

S
(bc)
66
D
n1
T
(c)
12

n
y
(bc)
6
))

N
b
b=1
(
(2)
R
(bc)
bb
D
n1
T
(c)
12

n(2)
v
(bc)
b
)
= h
n1

21
); c = 1; N
c
(69)
Eqs. (69) are immediately solved in favour for the independent
microscopical stress increments of axial shear:
D
n1
T
(c)
12
=
h
n1

21
)

N
b
b=1
[h
b

n
y
(bc)
6
)
n(2)
v
(bc)
b
[

N
b
b=1
[h
b

S
(bc)
66

(2)
R
(bc)
bb
[
; c = 1; N
c
(70)
6 B. Kurnatowski, A. Matzenmiller / Computational Materials Science xxx (2008) xxxxxx
ARTICLE IN PRESS
Please cite this article in press as: B. Kurnatowski, A. Matzenmiller, Comput. Mater. Sci. (2008), doi:10.1016/j.commatsci.2008.02.026
With the help of the independent average shear stiffnesses
C
(c)
66
:=
h

N
b
b=1
[h
b

S
(bc)
66

(2)
R
(bc)
bb
[
; c = 1; N
c
(71)
of the subcells, beaded in column c, and the associated history
terms:
n
^v
(c)
T
12
:=

N
b
b=1
[h
b

n
y
(bc)
6
)
n(2)
v
(bc)
b
[

N
b
b=1
[h
b

S
(bc)
66

(2)
R
(bc)
bb
[
; c = 1; N
c
(72)
Eq. (71) is written in the short form:
D
n1
T
(c)
12
= C
(c)
66

n1

21
)
n
^v
(c)
T12
; c = 1; N
c
(73)
Eqs. (73) are summarised by the matrix equation:
D
n1
T
12
=

C
66

n1

21
)
n
^v
T
12
(74)
with the column vectors D
n1
T
12
,

C
66
and
n
^v
T
12
which contain the
stresses D
n1
T
(c)
12
, stiffnesses C
(c)
66
and history variables
n
^v
(c)
T
12
ordered
by the index c over all columns in the RVE.
Analogously to Eq. (70), the stress increments D
n1
T
(b)
31
follow
from the continuity Eqs. (9) together with the constitutive Eqs.
(35) and (45) as
D
n1
T
(b)
13
=
l
n1

13
)

Nc
c=1
[l
c

n
y
(bc)
5
)
n(3)
v
(bc)
t
[

Nc
c=1
[l
c

S
(bc)
55

(3)
R
(bc)
tt
[
; b = 1; N
b
: (75)
With the help of the average shear stiffness C
(b)
55
of columns b
C
(b)
55
:=
l

Nc
c=1
[l
c

S
(bc)
55

(3)
R
(bc)
tt
[
; b = 1; N
b
(76)
and the stress history variables dened as
n
^v
(b)
T
13
:=

Nc
c=1
[l
c

n
y
(bc)
5
)
n(3)
v
(bc)
t
[

Nc
c=1
[l
c

S
(bc)
55

(3)
R
(bc)
tt
[
; b = 1; N
b
(77)
Eqs. (75) are notated as
D
n1
T
(b)
13
= C
(b)
55

n1

13
)
n
^v
(b)
T13
; b = 1; N
b
(78)
The matrix notation of Eqs. (78) reads:
D
n1
T
13
=

C
55

n1

31
)
n
^v
T
13
(79)
Finally, the transversal shear stress D
n1
T
23
has to be computed.
Inserting the constitutive equations (35) for
(bc)
23
) as well as Eqs.
(45) for the displacement jumps at the bre matrix interface into
the continuity condition (10) leads to the single equation for the
shear stress eld D
n1
T
23
being homogeneous within each RVE:

N
b
b=1

Nc
c=1
(h
b
l
c
[

S
(bc)
44
D
n1
T
23

n
y
(bc)
4
[ l
c
[
(2)
R
(bc)
tt
D
n1
T
23

n(2)
v
(bc)
t
[
h
b
[
(3)
R
(bc)
bb
D
n1
T
23

n(3)
v
(bc)
b
[) = hl
23
) (80)
Eq. (80) can be solved analytically for D
n1
T
23
:
D
n1
T
23
=
hl
n1

23
)

N
b
b=1

Nc
c=1
[h
b
l
c
n
y
(bc)
4
l
c
n(2)
v
(bc)
t
h
b
n(3)
v
(bc)
b
[

N
b
b=1

Nc
c=1
[h
b
l
c

S
(bc)
44
l
(2)
c

R
(bc)
tt
h
(3)
b

R
(bc)
bb
[
(81)
By means of the average stiffness term
C
44
:=
hl

N
b
b=1

Nc
c=1
[h
b
l
c

S
(bc)
44
l
c
(2)
R
(bc)
tt
h
b
(3)
R
(bc)
bb
[
(82)
and the history variable of the transversal shear stress
n
^v
T
23
:=

N
b
b=1

Nc
c=1
[h
b
l
c
n
y
(bc)
4
l
c
n(2)
v
(bc)
t
h
b
n(3)
v
(bc)
b
[

N
b
b=1

Nc
c=1
[h
b
l
c

S
(bc)
44
l
c
(2)
R
(bc)
tt
h
b
(3)
R
(bc)
bb
[
(83)
Eq. (81) simplies to
D
n1
T
23
= C
44

n1

23
)
n
^v
T
23
(84)
With the intention of a uniform representation of all incremental
stressstrain relations, the single Eq. (84) is written in matrix nota-
tion as
D
n1
T
23
=

C
44

n1

23
)
n
^v
T
23
(85)
Subsumming Eqs. (74), (79) and (85) the increments of the indepen-
dent shear stresses D
n1
T
ij
with i ,= j in all subcells are given by
D
n1
T
23
D
n1
T
13
D
n1
T
12
_

_
_

_
=

C
44
0 0
0

C
55
0
0 0

C
66
_

_
_

n1

23
)

n1

13
)

n1

12
)
_

_
_

n
^v
T
23
n
^v
T
13
n
^v
T
12
_

_
_

_
; (86)
where the coefcient matrix denes the submatrix

C
+
S
in Eq. (57).
2.3.3. Macrostresses of composite material
The normal stress components D
n1
r)
N
of the incremental mac-
roscopic stress tensor in vector form D
n1
r) are obtained by aver-
aging the increments D
n1
T
ij
of the independent microstresses of
the subcells in Eq. (68):
D
n1
r)
N
=
D
n1
r
11
)
D
n1
r
22
)
D
n1
r
33
)
_

_
_

_
=
1
hl

N
b
b=1

Nc
c=1
h
b
l
c
D
n1
T
(bc)
11
1
l

Nc
c=1
l
c
D
n1
T
(c)
22
1
h

N
b
b=1
h
b
D
n1
T
(b)
33
_

_
_

_
(87)
The rst three components of the average stress history tensor
stored in vector form
n
^v
r
), are computed from the associated ele-
ments of the vectors
n
^v
T
ij
in Eq. (68) for the subcells:

n
^v
r
)
N
=

n
^v
r
11
)

n
^v
r
22
)

n
^v
r
33
)
_

_
_

_
=
1
hl

N
b
b=1

Nc
c=1
h
b
l
n
c
^v
(bc)
T
11
1
l

Nc
c=1
l
n
c
^v
(c)
T
22
1
h

N
b
b=1
h
n
b
^v
(b)
T
33
_

_
_

_
(88)
The total normal stresses at the current time step
n1
t amount to

n1
r)
N
=
n
r)
N
D
n1
r)
N
(89)
with the stresses
n
r)
N
at the preceeding point in time
n
t and the
increment D
n1
r)
N
as in Eq. (87). The macroscopic shear stresses
are derived as the weighted sums of the components of the column
vectors D
n1
T
ij
= D
n1
T
()
ij
in Eq. (86):
D
n1
r)
S
=
D
n1
r
23
)
D
n1
r
13
)
D
n1
r
12
)
_

_
_

_
=
D
n1
T
23
1
h

N
b
b=1
h
b
D
n1
T
(b)
13
1
l

Nc
c=1
l
c
D
n1
T
(c)
12
_

_
_

_
(90)
The total shear stresses at the current time step
n1
t amount to

n1
r)
S
=
n
r)
S
D
n1
r)
S
(91)
with the stresses
n
r)
S
at the preceeding point in time
n
t and their
increments D
n1
r)
S
. Likewise averaging of the history variables
n
^v
T
ij
in Eq. (86) leads to the macroscopic history variables of the
shear stresses:

n
^v
r
)
S
=

n
^v
r
23
)

n
^v
r
13
)

n
^v
r
12
)
_

_
_

_
=
n
^v
T
23
1
h

N
b
b=1
h
b
n
^v
(b)
T
13
1
l

Nc
c=1
l
c
n
^v
(c)
T
12
_

_
_

_
(92)
B. Kurnatowski, A. Matzenmiller / Computational Materials Science xxx (2008) xxxxxx 7
ARTICLE IN PRESS
Please cite this article in press as: B. Kurnatowski, A. Matzenmiller, Comput. Mater. Sci. (2008), doi:10.1016/j.commatsci.2008.02.026
2.3.4. Homogenised instantaneous stiffness tensor and incremental
constitutive equation
Similar to the computation of the macrostresses D
n1
r) and the
macroscopic history variables
n
^v
r
) in Eqs. (87), (88), (90) and (92),
the homogenised instantaneous stiffness tensor in matrix form

C
+
of the composite is derived from the incremental subcell stresses
as a relation of the macroscopic strains and the history variables
on the basis of Eqs. (68) and (86) recapitulated as
D
n1
T
11
D
n1
T
22
D
n1
T
33
D
n1
T
23
D
n1
T
13
D
n1
T
12
_

_
_

_
=

C
11

C
12

C
13
0 0 0

C
21

C
22

C
23
0 0 0

C
31

C
32

C
33
0 0 0
0 0 0

C
44
0 0
0 0 0 0

C
55
0
0 0 0 0 0

C
66
_

_
_

n1

11
)

n1

22
)

n1

33
)

n1

23
)

n1

13
)

n1

12
)
_

_
_

n
^v
T
11
n
^v
T
22
n
^v
T
33
n
^v
T
23
n
^v
T
13
n
^v
T
12
_

_
_

_
(93)
and written symbolically with the help of the hypermatrix
^
C of size
(N
b
N
c
2(N
b
N
c
) 1) 6:
D
n1
T =
^
C
n1
)
n
^v
T
(94)
Note that all elements in the C-matrix of Eq. (93) are vectors of the
length as specied above. The elements of the upper left submatrix

C
+
N
in the homogenised instantaneous stiffness tensor of Eq. (57) re-
sult from averaging the entries of the column vectors

C
ij
= C
()
ij

with i; j 6 3 in the hypermatrix

C of Eqs. (93):

C
+
N
=
1
hl

N
b
b=1

Nc
c=1
h
b
l
c
C
(bc)
11
1
hl

N
b
b=1

Nc
c=1
h
b
l
c
C
(bc)
12
1
hl

N
b
b=1

Nc
c=1
h
b
l
c
C
(bc)
13
1
l

Nc
c=1
l
c
C
(c)
21
1
l

Nc
c=1
l
c
C
(c)
22
1
l

Nc
c=1
l
c
C
(c)
23
1
h

N
b
b=1
h
b
C
(b)
31
1
h

N
b
b=1
h
b
C
(b)
32
1
h

N
b
b=1
h
b
C
(b)
33
_

_
_

_
(95)
The three effective shear stiffness components

C
+
44
,

C
+
55
and

C
+
66
align with the diagonal terms of

C
+
S
in Eq. (57). They are obtained
by averaging the vectors

C
ij
= C
()
ij
with indices i = j _ 4 over the
columns c and rows b of the RVE-model, respectively:

C
+
S
=

C
+
44
0 0
0

C
+
55
0
0 0

C
+
66
_

_
_

_
=
C
44
0 0
0
1
h

N
b
b=1
h
b
C
(b)
55
0
0 0
1
l

Nc
c=1
l
c
C
(c)
66
_

_
_

_
(96)
By using the matrices given in Eqs. (96) and (95), the macroscopic
stress increments may be assembled in vector notation and written
as
D
n1
r)
N
D
n1
r)
S
_ _
=

C
+
N
0
0

C
+
S
_
_
_
_

n1
)

n
v
r
)
N

n
v
r
)
S
_ _
(97)
which reads in symbolic representation:
D
n1
r) =

C
+

n1
)
n
v
r
) (98)
Hence, the total stress
n1
r) =
n
r) D
n1
r) of the homogenised
substitute continuum is computed with the help of the increments
in Eq. (98) to:

n1
r) =
n
r)

C
+

n1
)
n
^v
r
) (99)
Here, the total stresses depend formally on the total strains, the ten-
sor of time dependent moduli

C
+
, the stresses of the previous time
step and the stress history variables. By splitting the total macro-
strain into
n1
) =
n
) D
n1
), Eq. (99) can be transformed in or-
der to represent the macrostresses as a function of the macrostrain
increments:

n1
r) =

C
+
D
n1
)
n
v
r
) (100)
wherein the newly introduced average macrostress history tensor
in vector notation
n
v
r
) is dened as

n
v
r
) :=
n
r)

C
+

n
)
n
^v
r
) (101)
2.4. Derivation of strain concentration tensors
Within the theory of homogenisation for linear elastic compos-
ites the average strain tensor
(i)
) of the phase (i) inside the RVE
can be described by making use of Hills strain concentration ten-
sor A see Hill [10]:

(i)
) = A
(i)
: ) (102)
Eq. (102) states a linear mapping of the macrostrain onto the phase
strain. In the case of time-dependent material behaviour, the phase
strains
(i)
) become a tensor functional of the macrostrain history:

(i)
(t)) = A); s [ [
s=t
s>
(103)
The functional A for the viscoelastic composite is computed
from the solution of the micromechanical model in Eq. (94). There-
fore, the incidence matrices B
(bc)
explicitly given in Appendix A.1
are introduced in order to extract the average subcell stress tensor
in matrix notation r
(bc)
) from the hypervector of all microstresses

T. The matrices B
(bc)
must not be confused with the stress concen-
tration tensor which maps the macrostresses onto the average
stresses r
(i)
) of phase (i). The increments of the subcell stresses
read then:
D
n1
r
(bc)
) := B
(bc)
D
n1
T (104)
The stress increments in Eq. (104) are combined with the constitu-
tive Eqs. (35) and (94) for the subcell X
(bc)
to get:

n1

(bc)
) =

S
(bc)
ve
[B
(bc)
D
n1
T[
n
y
(bc)
)
=

S
(bc)
ve
B
(bc)
(

C
n1
)
n
^v
T
)
_ _

n
y
(bc)
) (105)
Hence, the subcell strains become a linear functional of the
macrostrains:

n1

(bc)
) = A[); s[
s=t
s>
= A
(bc)

n1
) a
(bc)
h
(106)
with the actual concentration matrix
A
(bc)
:=

S
(bc)
ve
B
(bc)
C (107)
and the history variable
a
(bc)
h
:=

S
(bc)
ve
B
(bc)n
^v
T

n
y
(bc)
) (108)
The 6 6 matrices A
(bc)
are determined by the RVE-geometry as
well as by the material properties of the subcells. Choosing a
constant time stepping Dt, the elements of the matrices A
(bc)
are
constant. The time-dependency of the subcell strains is governed
by the history vector a
(bc)
h
solely, if a constant time-stepping Dt is
applied. The average strain concentration matrix of phase (i) is
computed as the weighted sum of the concentration matrices
of all subcells X
(i)
= X
(bc)
which are lled with the material phase
(i):
8 B. Kurnatowski, A. Matzenmiller / Computational Materials Science xxx (2008) xxxxxx
ARTICLE IN PRESS
Please cite this article in press as: B. Kurnatowski, A. Matzenmiller, Comput. Mater. Sci. (2008), doi:10.1016/j.commatsci.2008.02.026
A
(i)
:=
1
hl

X
(i)
h
b
l
c

S
(bc)
ve
B
(bc)
C (109)
The strain history vector a
(i)
h
of the phase (i) reads analogeously:
a
(i)
h
:=
1
hl

X
(i)
h
b
l
c

S
(bc)
ve
B
(bc)n
^v
T

n
y
(bc)
) (110)
3. Twoscale analysis with FEM/GMC of composite structures
The nite element code FEAP of Taylor [23], used for the subse-
quent macroscopic analyses, allows for the implementation of con-
stitutive models for elements with a three-dimensional material
formulation. The GMC-model was implemented for the 8-node vol-
ume element SOLId8 which is based on linear shape functions N.
The discrete displacement eld u
h
= Nu
k
depends on the nodal
values u
k
. The macroscopic strain eld is given as
h
= LNu
k
=
Bu
k
with the well-known kinematic operator L of the small defor-
mation theory, see Taylor [23]. It is assigned one RVE to each
integration point x
p
of the nite element mesh. The macroscopic
strain
h
(x
p
) = )
p
denes the boundary conditions of the GMC-
analyses. The volume averaged stress tensor r) and the tangent
stiffness matrix

C
+
p
are returned from the micromechanical subpro-
gram and enter the vector of internal forces f
in
and the stiffness
matrix K of the macrostructure respectively. The solution scheme
for the micromechanically based two-scale analyses is sketched
in Fig. 2.
If the time increment Dt is not changed during the analysis, the
repeated computation of K and its inverse is obsolete. This is due to
the fact that the stiffness matrices of the phases and interfaces are
constant for a xed time increment, see Eqs. (36), (46) and (47).
Hence, the computation of the displacement increment Du
k
can
be performed with the initial stiffness matrix
0
K =
n
K[
Dtfixed
so that
the GMC is only used repeatedly to calculate the history depen-
dend macroscopic stress tensor r)
p
, needed at all Gauss points p
in each time step n for the assembly of the vector of internal forces
f
in
. For this the time-dependent state variables d
(bc)
k
and
(i)
q
(bc)
k
of all
subcells and interfaces need to be stored for the next time step as
well as the components of the hyper vector of stresses

T. The kine-
matical variables like
(bc)
ij
) and
(i)
su
i
t
(bc)
do not appear explicitly in
Fig. 2. Solution scheme for quasi-static two-scale analysis of composite structures.
B. Kurnatowski, A. Matzenmiller / Computational Materials Science xxx (2008) xxxxxx 9
ARTICLE IN PRESS
Please cite this article in press as: B. Kurnatowski, A. Matzenmiller, Comput. Mater. Sci. (2008), doi:10.1016/j.commatsci.2008.02.026
the stress based formulation of the GMC and thus do not have to be
memorised.
4. Numerical analyses and validation of experimental tests
The nal section of the publication is concerned with the veri-
cation and validation of the proposed micromechanical material
model. For this purpose the continuums theory of heterogeneous
materials is applied for the solution of boundary value problems
as given by structures made of bre reinforced composites.
4.1. Verication example: thick-walled cylinder under internal
pressure
The thick-walled tube under internal pressure shown in Fig. 3
serves as the rst example in order to verify the implementation
of the two-scale model. The exact analytical solution of the plane
strain problem can be found in Lekhnitskii [12]. It is remarked that
the distribution of stresses in axial, tangential and radial direction
of the tube is temporally invariant for this boundary problem
although the material behaviour is time-dependent, since the
AIRY-stress-function of the plane problem does not depend on
any material parameters see Findley et al. [5]. The nite element
model for the radially symmetric eld problem consists of the cir-
cular segment with the angle a = 5

as sketched in Fig. 3. The inner

and outer diameters of the tube are r
i
= 2 cm and r
o
= 4 cm,
respectively. The segment is discretised with forty 8-node solid
elements in radial direction. There is only one element in circum-
ferential and axial direction, which is sufcient since the tube un-
der internal pressure means a mathematically 1-dimensional
problemdepending on the radius r as the only coordinate. The tube
is exposed to an internal pressure p = 1 MPa.
The nodal degrees of freedom perpendicular to the cut edges are
locked. There is also no displacement parallel to the tubes axis of
polar symmetry. Hence, the system is in the state of plane strain.
The creep behaviour of the tube is investigated for two different
types of viscoelastic constitutive modelling. First, the linear visco-
elastic material model of FEAP is used. There, the relaxation mod-
ulus function of the homogeneous material is dened as
G(t) = 2G

N
i=1
G
i
e
t
s
i (111)
The bulk modulus K is time-independent. Afterwards, the GMC-
model is applied with all subcells being perfectly bonded and lled
with the same viscoelastic material namely epoxy resin, leading to
an isotropic constitutive behaviour. The numerical results for stres-
ses and displacements of both approaches are compared. Hereby, it
is intended to assure that the implementation of the GMC-algo-
rithm has been done correctly and to appraise the increase of com-
putational time caused by the latter two-scale modelling. The
creep-function, used for the GMC-model, reads:
J(t) = J
0

N
i=1
J
k
1 e
t
^s
k
_ _
(112)
Obviously, the two approaches are not immediately comparable
since the homogeneous model is based on the relaxation function
G(t) and the micromodel is premised on the creep function J(t). In
order to answer the purpose of verication the creep and relaxation
spectra have to meet specic requirements. For the Laplace trans-
formed J(s) and G(s) of the creep and relaxation functions J(t) and
G(t) it holds with the complex Laplace variable s:
G(s) =
1
s
2
J(s)
(113)
On the basis of Eq. (113) Gross [7] computed the discrete creep
spectrum from the given relaxation spectrum. For the special case
of rheological models with N = 1 in Eqs. (111) and (112) the Eq.
(113) leads to:
^s
1
= 2(J
0
J
1
)(G

G
1
)s
1
(114)
Setting J

= J
0
J
1
= 1=(2G

) as in Eq. (20), the relaxation and

retardation times s
1
and ^s
1
have to meet the condition:
^s
1
= 1
G
1
G

_ _
s
1
(115)
Fig. 3. Thick-walled tube under internal pressure p = 1 MPa. FE-model of the 5 cut-out. Material properties are given in MPa, MPa
1
and hours, respectively.
Fig. 4. Verication of displacement predictions in radial direction. FEM/GMC versus
FEM/homogeneous for isotropic material.
10 B. Kurnatowski, A. Matzenmiller / Computational Materials Science xxx (2008) xxxxxx
ARTICLE IN PRESS
Please cite this article in press as: B. Kurnatowski, A. Matzenmiller, Comput. Mater. Sci. (2008), doi:10.1016/j.commatsci.2008.02.026
Choosing the initial value as J(0) = J
0
= 1=(G(0)) =
1
2G2G
1
, the creep
strength J
1
is calculated to:
J
1
= J

J
0
=
G
1
G

G(0)
=
1
2G

1
G
G
1
_ _ (116)
The numerical values chosen for the parameters of the rheological
models are tabulated in Fig. 3. The evolution of the displacements
u(r; t) in radial direction at the inner (r = r
i
) and outer (r = r
o
)
periphery of the thick-walled specimen is shown in Fig. 4.
Obviously, the GMC-model with a regular pattern of four
equally sized subcells renders exactly the same numerical results
as the homogeneous viscoelastic model of FEAP.
The diagrams in Fig. 5 show the absolute value of the radial and
hoop stresses [r
r
=p[ and [r
h
=p[ related to the internal pressure p at
time t = 0 h and t = 200 h as a function of the radius r. Due to the
afore mentioned statical determination of the tube, the stresses
should not vary with time. The outcome of the homogeneous and
the GMC-based constitutive model are again identical and almost
equal to the analytical solution for an elastic isotropic material.
The analytical solution for the distribution of radial and hoop
stresses through the thickness of the tube is given by the formula:
r
r
r
h
_ _
= p
c
k1
1c
2k
(q
k1
q
(k1)
)
kc
k1
1c
2k
(q
k1
q
(k1)
)
_ _
(117)
which is taken fromLekhinitskii [12]. The abbreviations are dened as
q =
r
r
o
c =
r
i
r
o
k =

b
11
b
22

(118)
The values of b
ij
in the parameter k depend on the components S
ij
of
the compliance tensor S:
b
ij
= S
ij

S
i3
S
j3
S
33
(119)
The compliances S
ij
are calculated from the engineering constants
to:
S
11
=
1
E
11
; S
22
=
1
E
22
; S
13
=
m
13
E
11
; S
23
=
m
23
E
33
(120)
Now that the results of the GMC approach are veried for a homo-
geneous, isotropic material, the attention is focused next on rein-
forced composites with bres in circumferential direction. The
bonding of the investigated glass bres and the epoxy ller is con-
sidered as rigid in the rst example. The elastic and rheological
properties of both phases are summarised in the table of Fig. 3.
The volume fraction of the bres is set to c
f
= 30%. In order to verify
the implementation of the numerical time integration scheme, the
solution for the radial and circumferential normal stresses is com-
pared with the analytical solution for the non-isotropic material gi-
ven in Eq. (117). The homogenisation algorithm leads to the
constitutive tensor in matrix notation:

C
+
=
8036 3387 3477 0 0 0
3387 25050 3387 0 0 0
3477 3387 8036 0 0 0
0 0 0 2015 0 0
0 0 0 0 2015 0
0 0 0 0 0 1764
_

_
_

_
(121)
Fig. 5. Verication of radial and circumferential stresses. FEM/hom., FEM/GMC and analytical solution for isotropic case of absolute values for rr and r#.
Fig. 6. Verication of radial and circumferential stresses. FEM/GMC versus analytical solution for bre-reinforced resin.
B. Kurnatowski, A. Matzenmiller / Computational Materials Science xxx (2008) xxxxxx 11
ARTICLE IN PRESS
Please cite this article in press as: B. Kurnatowski, A. Matzenmiller, Comput. Mater. Sci. (2008), doi:10.1016/j.commatsci.2008.02.026
which is given here with reference to the cylindrical base system of
the tube, see Fig. 3. The factor k for the analytical solution in Eq.
(117) is computed to k = 1:7656 on the basis of the homogenised
material properties. The diagrams in Fig. 6 show the distributions
of the normal and hoop stresses [r
r
=p[ and [r
h
=p[ in radial direction.
The results obtained by the FEM/GMC approach t well the analyt-
ical solution with some small deviation at the inner edge of the tube
where the pressure is imposed.
Various generic cells with a different number of subcells are
investigated with respect to the required run time of the FEM/
GMC approach. The increase of the relative computational time
as a function of the number of subcells can be read from Table 1
where the relative run time is given with respect to the computa-
tion time of the homogeneous model. The growth is a little more
than linear for the efciently reformulated GMC used in the pres-
ent study. If the discretisation of the RVE is chosen such that
N = N
b
= N
c
is the number of subcells along each edge of the
RVE, the number of independent unknowns grows linearly with
N in the reformulated version of the GMC. The amount of unknown
variables increases with N
2
in the original formulation.
4.2. Verication example: perforated tension strip
Fig. 7 shows the meshing of one quarter piece of a perforated
tension strip cut out at the lines of symmetry. The bres are ori-
ented under an angle a relative to the longer side of the thin plate.
The thickness of the plate is discretised with one solid element. The
unidirectionally bre reinforced sample is stressed all at once with
the uniaxial tension load p = 1 MPa at time t = 0. The loading is
kept constant for the period of 100 hours. For all nodes located on
the axes of symmetry with x
1
= const: = 0 and x
2
= const: = 0,
the degrees of freedom perpendicular to the cut edge are locked.
The tension strip is assumed to be in a state of plane stress.
The interface parameters are derived from the viscoelastic ma-
trix properties E, G

, the bre diameter d

f
and the ratio of inter-
phase thickness to bre thickness r
fi
:
R
n
=
r
fi
d
f
E
R
1
=
r
fi
d
f
G
1
R

=
r
fi
d
f
G

The relation of height to length is h=l = 1=2. The radius of the circu-
lar hole amounts to R = 1=5h.
The computation is based on the material parameters given in
Table 2.
While the bres behaviour is elastic, the matrix phase behaves
viscoelastic. The bonding of the phases is assumed to depend on
time in tangential direction and to be elastic in normal direction.
The bre volume content is c
f
= 30%.
The numerical analysis is done for the bre angles a = 0 and
a = 90

and with the three discretisations of the RVE given in

Fig. 8. The displacement-time history of the nodal point P see
Fig. 7 is depicted in the left diagram of Fig. 9 for the case of the
bre angle a = 0

. The creeping of the perforated tension strip is

unincisive since the orientation of the stiff elastic bres coincides
with the loading direction. The normal stress distribution along
the vertical axis of symmetry is plotted in the right diagram of
Fig. 9. It is obvious that the stresses do not strongly depend on
the level of discretisation. The spatial distribution of the stresseld
for r
11
slightly changes with time as it can be seen also from
Fig. 10. The stress r
11
concentrates around the circular hole.
The results for the creep displacement u
1
at point P and the
normal stress r
11
along the vertical axis with x
1
= 0 are depicted
in Fig. 11. The initial displacement for a = 90

is approximately
three times higher than the displacement for a = 0

. This is expect-
edly due to the fact that the reinforcing elastic and stiff bres are
perpendicular to the loading direction. Hence, the effective stiff-
ness is dominated by the viscoelastic matrix phase and the compli-
ant bonding in the rst case.
The deformation for the coarse micromodel Q2 with just four
subcells leads to a more compliant effective material than the
models Q5 and Q7 with 25 and 49 cells, respectively. This is not
so obvious for the rst case with a = 0

where the bres are aligned

along the loading direction.
The stress distribution, shown in the right diagram of Fig. 11, is
less concentrated in the vicinity of the perforation in contrast to
the one in Fig. 9. The change of the stress intensity in time is small
but recognisable again at the rim of the hole. The macro-stresses
do not depend signicantly on the RVE-Type, i.e. the number of
subcells used to discretise the microstructure.
Table 1
Evolution of computation time for increased number of subcells
Nb Nc: 1 1 2 2 4 4 5 5 6 6 7 7 8 8
Time
tGMC
thom
: 1.25 1.85 4.08 5.73 7.61 10.08 13.00
Fig. 7. Tension strip under uniaxial loading p.
Table 2
Parameters of bre, matrix and interface
Material Elastic parameters Viscoelastic parameters
Eq. Shear
modulus, G
(MPa)
Youngs
modulus, E
(MPa)
Creep time,
^s
1
[h]
Contribution
factor, ^m
1
= J
1
=J

Glass 28725 68940 0

Epoxy 893.3 3419.5 14.83 0.3
Interface parameters
n (h) Rn R R
1
r
fi
14.286 1:462 10
3
5:597 10
3
1:679 10
3
1/30
Fig. 8. RVE-types assume a quadratic packing of the bres.
12 B. Kurnatowski, A. Matzenmiller / Computational Materials Science xxx (2008) xxxxxx
ARTICLE IN PRESS
Please cite this article in press as: B. Kurnatowski, A. Matzenmiller, Comput. Mater. Sci. (2008), doi:10.1016/j.commatsci.2008.02.026
4.3. Validation example: uniaxial relaxation test
Masuko and Kawai [13] investigated the time-dependent
behaviour of a symmetric angle-ply as well as an UD-laminate
made of T800H/3631 carbon/epoxy composite see Table 3 for
constitutive properties of the isotropic epoxy resin and the trans-
versly isotropic carbon bres at an elevated temperature of
T = 100 C. The experimental data were obtained from uniaxial
strain-controlled relaxation tests. The initial elastic properties of
the single UD-layer were also given by Masuko and Kawai [13]
and are compared to the GMC predictions in Table 4. The four
RVE-models, applied to the analyses, are depicted in Fig. 12. The
rst three sketches of the RVEs assume a quadratic packing of
the bres, which leads to orthotropic properties. In order to obtain
the elasticity parameters of the transversely isotropic macroscopic
continuum, the direction dependent stiffness tensor is integrated
about the bre axis and averaged according to:
C
Ti
ijkl
=
1
p
_
p
0
C
Ortho
ijkl
d/ (122)
The fourth RVE-model aims directly at the analysis with trans-
versely isotropic parameters by supposing a hexagonal bre pack-
ing. The initial elastic properties show minor dependence on the
rened approximation of the circular bre shape as it can be seen
from Table 3. The differences between the quadratic and hexagonal
modelling are at least partly due to the smaller bre content of the
latter. While the upper row and the left column of the quadratic
models allow for the continuous change of the relative bre con-
tent, the bre volume fraction of the hexagonal model is altered
by increasing the diameter of the bres. This is restricted to a step-
wise enlargement because of the discrete subcell structure.
The conformity of the results in Table 4 recommends the use of
the Q2 model for the subsequent two-scale simulations with re-
gard to the balanced efciency of modelling errors and computa-
tional effort.
The test specimens with dimensions l=w=t = 200=20=2 mm are
clamped at both ends and have a gauge length of 100 mm. The an-
gle between the embedded bres and the direction of the applied
loading is given by h. For the nite element analyses, the free gauge
length of the UD-specimen is discretised with a total of 50 10 1
eight-node isoparametric volume elements. The nodal displace-
ments u at both clamped ends of the nite element mesh are pre-
scribed in loading direction in order to obtain the average
longitudinal strain
a
= Du=l = 0:44%. While the inplane displace-
ments perpendicular to the longer side of the strip are inhibited
at the clamped ends, there is no obstruction of any thickness con-
traction (see Fig. 13). The corresponding mean stress r
a
= P=A is
computed from the resulting longitudinal force P. The creep func-
tion of the epoxy matrix is assumed as
J(t) =4:316 10
4
0:270 10
4
1e

t
0:1h
_ _
5:935 10
4
(1e

t
0:75h)
0:270 10
4
e

t
10:0h
_ _
[MPa
1
[ (123)
with the initial value J(0) = 1=2G
0
and G
0
= 1:16 (GPa). The creep
spectrum L(s) is chosen by trial and error to t the relaxation curve
for h = 60

qualitatively well.
The tangential interface creep function R(t) is taken as j times
the matrix creep function J(t) premultiplied by the interphase
thickness t
(i)
~
fibre
=30 = 0:2 lm:
R
t
(t) = t
(i)
jJ(t); R
n
= t
(i)
j
E
m
(124)
The tangential stiffness R
t
(t) of the bonding tends to innity
as j 0 approaches zero. Analogously, the normal interface
Fig. 9. Stress and displacement plots for a = 0

. Left: Displacement u1 versus time at point P. Right: Nodal projections of the normal stress r11 along the edge with
x1 = const = 0.
Time = 1.00E+00
2.40E-01
5.00E-01
7.60E-01
1.02E+00
1.28E+00
1.54E+00
1.80E+00
2.06E+00
2.32E+00
2.58E+00
2.84E+00
-1.97E-02
3.10E+00
S T R E S S 1
Time = 1.00E+00
Time = 1.00E+02
2.41E-01
5.14E-01
7.86E-01
1.06E+00
1.33E+00
1.60E+00
1.88E+00
2.15E+00
2.42E+00
2.69E+00
2.96E+00
-3.11E-02
3.24E+00
S T R E S S 1
Time = 1.00E+02
Fig. 10. Normal stress distribution r11 at time t = 0 h and t = 100 h for bre dire-
ction a = 0

.
B. Kurnatowski, A. Matzenmiller / Computational Materials Science xxx (2008) xxxxxx 13
ARTICLE IN PRESS
Please cite this article in press as: B. Kurnatowski, A. Matzenmiller, Comput. Mater. Sci. (2008), doi:10.1016/j.commatsci.2008.02.026
compliance R
n
is chosen proportionally to the Youngs modulus
E
m
= E
a
= E
t
= 3:2 of the isotropic matrix material. The numerical
predictions of the two-scale computations are compared to the test
data in Fig. 14. The length of the time step is set to Dt = 0:01 h. The
two-scale model is solved with 500 time-steps each of which
requiring 5 s of computational time. The relative computational
time t
GMC
=t
hom
is given as a function of the RVE-model in Table 5.
By comparing Table 5 with Table 1 it becomes obvious that the size
of the nite element mesh inuences the increase of the relative
computational time: the larger the nite element mesh the smaller
the relative increase of calculation time due to the number of
subcells.
The numerical results show that the assumption of perfect
bonding (j = 0) overestimates the macroscopic stiffness of the
specimen and, hence, leads to too high stresses r
a
in the axial
direction (loading direction) of the test specimen. By taking the
Fig. 11. Stress and displacement plots for a = 90

. Left: Displacement u1 versus time at point P. Right: Nodal projections of the normal stress r11 along the edge with
x1 = const = 0.
Table 3
Elastic properties of the phases
Ea (GPa) Ga (GPa) Et (GPa) Gt (GPa) ma []
Epoxy-matrix 3.20 1.16 3.20 1.16 0.35
Carbon-bre 288.0 24.0 16.0 6.7 0.35
Table 4
Comparison of initial elastic properties from experiments with GMC predictions
v
f
Ea (GPa) Ga(GPa) Et(GPa) Gt (GPa) ma []
UD-Exp. ~0.52 151.0 3.5 7.3
GMC Q2 0.52 151.30 2.98 8.05 2.90 0.36
GMC Q10 0.52 151.30 3.01 7.98 2.86 0.36
GMC Q50 0.52 151.30 3.02 7.99 2.87 0.36
GMC H30 0.51 149.31 2.36 6.53 2.22 0.36
Fig. 12. Discrete RVEs for quadratic and hexagonal bre packing.
Fig. 13. Specimens geometry in mm and nite element mesh of free gauge length.
14 B. Kurnatowski, A. Matzenmiller / Computational Materials Science xxx (2008) xxxxxx
ARTICLE IN PRESS
Please cite this article in press as: B. Kurnatowski, A. Matzenmiller, Comput. Mater. Sci. (2008), doi:10.1016/j.commatsci.2008.02.026
exible interface with j = 1. . . 2 into account, it is possible to t
the model response to the relaxation data well. The best result
for h = 60

is obtained by setting the interface parameter j = 2.

In that case, the interface behaves twice as soft as the matrix phase.
For h = 45

the t with j = 2 is not as good, since the computed

stress is to small compared to the experimental data. This is sup-
posed to be due to the non-linearity of the viscoelastic response
which is experimentally detected by Masuko and Kawai [13]. The
linear viscoelastic parameters are identied at the higher strain
of 0.56%. Hence, the stiffnesses may be to small to obtain good re-
sults at the lower level of strain of 0.44%.
4.4. Validation example: thin-walled tube under axial loading
Robinson et al. [21] stressed thin-walled tube specimens of a
composite material with the axial tension loading p = 28 MPa at
room temperature. The single-layer test specimens consist of a
glass-bre reinforced epoxy matrix material with the bre volume
fraction c
f
= 60% and a bre diameter of 6 lm. The inner and outer
diameters of the tubes are
i
= 33 mm and
o
= 37 mm, respec-
tively. The axial length amounts to l = 100 mm. The winding angle
of the bres relative to the longitudinal axis x
3
is h = 45

. The tubes
are clamped at their ends to the load controlling test device. At the
one end of the nite element model the nodal degrees of freedom
are locked in longitudinal as well as in radial direction. At the other
end only the displacements of the radial degrees of freedom are
prescribed equal to zero. The displacements in the circumferential
direction are not constraint. In longitudinal direction the force
boundary conditions are given in terms of a constant normal stress
of 28 MPa. The compliance of the interface model is given as j-
times the compliance of the matrix phase taking into account the
interphase thickness t
(i)
= 0:2 lm as it was done in the previous
numerical example (see Fig. 15).
The experimental results of the creep tests are given in terms
of the longitudinal and shear creep strains (t) =
total
(t)
elastic
versus time plots. The diagrams in Fig. 16 show the creep re-
sponse for (a) rigid bond (j = 0), (b) an interface with properties
equivalent to the matrix phase (j = 1) and (c) a rather soft bond-
ing (j = 2). As expected, an increasing interface compliance
intensies the axial creep. The shear creep strain, which is not
caused by any torsional loading but by the bres angle H = 45

versus the load direction, does not show a discernible dependence

on the interface modelling.
5. Summary and concluding remarks
The present article is concerned with the two-scale simulation
of viscoelastic composite structures. While the macrostructure is
modelled with standard nite element tools, the simultaneous
microscale analysis is performed with the numerically efcient
generalized method of cells. The micromechanically based ap-
Fig. 14. Experimental and simulated stress-relaxation-curves of carbon/epoxy composite for the two off-axis angles h = 60

and h = 45

at moderate constant axial strains.

Table 5
Computational time t
rel
= tGMC=t
hom
GMC-model Q2 Q10 Q50 H30
Time
tGMC
thom
1.1 1.5 75.1 34.8
Fig. 15. Elastic and viscoelastic model parameters.
B. Kurnatowski, A. Matzenmiller / Computational Materials Science xxx (2008) xxxxxx 15
ARTICLE IN PRESS
Please cite this article in press as: B. Kurnatowski, A. Matzenmiller, Comput. Mater. Sci. (2008), doi:10.1016/j.commatsci.2008.02.026
proach provides the effective macroscopic stress response to strain
processes as well as the homogenised constitutive tangent tensor
from the phase and interphase properties as well as the bre archi-
tecture. The numerical examples show that two-scale simulations
can be put into practice with a reasonable expenditure of compu-
tational time. Nevertheless, material testings evince a distinct non-
linear behaviour of epoxy materials especially under shear loading
conditions see among others Hu et al. [11] or Megnis and Varna
[16]. In order to cope with nonlinear phenomena at higher stress
strain levels, it is adviseable to replace the linear model with, e.g.
the nonlinear viscoelastic approach dating from Schapery [22],
the implementation of which is feasible.
Acknowledgement
The nancial support of the Deutsche Forschungsgemeinschaft,
DFG, under Contract MA1186/4 is gratefully acknowledged.
Appendix A
A.1. Cell-stress components stored in hyper vectors
The normal stresses T
(:)
ij
are compiled by the vectors:

T
11
:=
T
(11)
11
T
(21)
11
.
.
.
T
(N
b
Nc)
11
_

_
_

_
(N
b
Nc)
;

T
22
:=
T
(1)
22
T
(2)
22
.
.
.
T
(Nc)
22
_

_
_

_
Nc
;

T
33
:=
T
(1)
33
T
(2)
33
.
.
.
T
(N
b
)
33
_

_
_

_
N
b
(125)
and the column matrices of all shear stresses read:

T
23
:= T
23
;

T
13
:=
T
(1)
13
T
(2)
13
.
.
.
T
(N
b
)
13
_

_
_

_
N
b
;

T
12
:=
T
(1)
12
T
(2)
12
.
.
.
T
(N
b
)
12
_

_
_

_
Nc
(126)
The hyper vector

T assembles the vectors just dened in Eqs. (125)
and (126):

T :=

T
T
11

T
T
22

T
T
33

T
T
23

T
T
13

T
T
12
_ _T
(127)
In order to read out the components of subcell stress tensor
in matrix notation r
(bc)
) see Eq. (34) from the hyper vector

T,
the 6 (N
b
N
c
2(N
b
N
c
) 1) incidence matrices B
(bc)
of the
subcells X
(bc)
are dened with the help of Kroneckers symbol d
ij
as
B
(bc)
= b
ij
_
:=
b
1j
= d
jn
1
_ _
b
2j
= d
jn
2
_ _
b
3j
= d
jn
3
_ _
b
4j
= d
jn
4
_ _
b
5j
= d
jn
5
_ _
b
6j
= d
jn
6
_ _
_

_
_

_
with
n
1
= N
b
(c 1) b
n
2
= N
b
N
c
c
n
3
= N
b
N
c
N
c
b
n
4
= N
b
N
c
N
c
N
b
1
n
5
= N
b
N
c
N
c
N
b
1 b
n
6
= N
b
N
c
N
c
2N
b
1 c
(128)
The range of the column index j is j =1; 2; . . . ; (N
b
N
c
2(N
b
N
c
) 1)
for all subcells X
(bc)
. Hence, it holds for the subcell stresses:
r
(bc)
) = T
(bc)
= B
(bc)
T (129)
Fig. 16. Creep strains versus time for different interface parameters j.
16 B. Kurnatowski, A. Matzenmiller / Computational Materials Science xxx (2008) xxxxxx
ARTICLE IN PRESS
Please cite this article in press as: B. Kurnatowski, A. Matzenmiller, Comput. Mater. Sci. (2008), doi:10.1016/j.commatsci.2008.02.026
A.2. Matrices for the coupled system of equations
The vectors in Eq. (60) are dened as
K
11
:=
1=S
(11)
11
1=S
(21)
11
.
.
.
1=S
(N
b
Nc)
11
_

_
_

_
(N
b
Nc)
;
n
v
11
:=

n
y
(11)
1
)=S
(11)
11

n
y
(21)
1
)=S
(21)
11
.
.
.

n
y
(N
b
Nc)
1
)=S
(N
b
Nc)
11
_

_
_

_
(N
b
Nc)
(130)
and the non-quadratic matrices
S
12
:=

S
(11)
12
=

S
(11)
11
.
.
.

S
(N
b
1)
12
=

S
(N
b
1)
11

S
(12)
12
=

S
(12)
11
.
.
.

S
(N
b
2)
12
=

S
(N
b
2)
11
.
.
.

S
(1Nc)
12
=

S
(1Nc)
11
.
.
.

S
(N
b
Nc)
12
=

S
(N
b
Nc)
11
_

_
_

_
(N
b
Nc)Nc
(131)
and
S
13
:=

S
(11)
13
=

S
(11)
11

S
(21)
13
=

S
(21)
11
.
.
.

S
(N
b
1)
13
=

S
(N
b
1)
11

S
(12)
13
=

S
(12)
11

S
(22)
13
=

S
(22)
11
.
.
.

S
(N
b
2)
13
=

S
(N
b
2)
11
.
.
.

S
(1Nc)
13
=

S
(1Nc)
11

S
(2Nc)
13
=

S
(2Nc)
11
.
.
.

S
(N
b
Nc)
13
=

S
(N
b
Nc)
11
_

_
_

_
(N
b
Nc)N
b
(132)
The entries of the following matrices, used to assemble Eq. (67), are
given by the denitions in Eqs. (64) and (66). The compliance ma-
trix
^
S with the submatrices S
22
, S
23
, S
32
and S
33
reads:

S :=
S
22
S
23
S
32
S
33
_ _
=
S
(1)
22
S
(11)
23
S
(21)
23
. . . S
(N
b
1)
23
S
(2)
22
S
(12)
23
S
(22)
23
. . . S
(N
b
2)
23
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
S
(Nc)
22
S
(1Nc)
23
S
(2Nc)
23
. . . S
(N
b
Nc)
23
S
(11)
32
S
(12)
32
. . . S
(1Nc)
32
S
(1)
33
S
(21)
32
S
(22)
32
. . . S
(2Nc)
32
S
(2)
33
.
.
.
S
(N
b
1)
32
S
(N
b
2)
32
. . . S
(N
b
Nc)
32
S
(N
b
)
33
_

_
_

_
(N
b
Nc)(N
b
Nc)
(133)
The right hand side of Eq. (67) is built with:
K
21
:=
S
(1)
21
1 0
.
.
.
S
(Nc)
21
1 0
_

_
_

_
; K
31
:=
S
(1)
31
0 1
.
.
.
S
(N
b
)
31
0 1
_

_
_

_
(134)
and
n
v

22
:=
n
v
(1)
22

n
v
(1)
21
.
.
.
n
v
(Nc)
22

n
v
(Nc)
21
_

_
_

_
;
n
v

33
:=
n
v
(1)
33

n
v
(1)
31
.
.
.
n
v
(N
b
)
33

n
v
(N
b
)
31
_

_
_

_
(135)
The solution of Eq. (67) given by Eq. (68) depends upon the six col-
umn vectors C
ij
collected in the following matrix:
C
21
C
22
C
23
C
31
C
32
C
33
_ _
=
S
22
S
23
S
32
S
33
_ _
1
K
21
K
31
_ _
(136)
with the size (N
c
N
b
) 3 and in addition on the three vectors
C
11
= K
11
S
12
C
21
S
13
C
31
C
12
= S
12
C
22
S
13
C
32
C
13
= S
12
C
23
S
13
C
33
(137)
The history vector in Eq. (67) is given by
n
v
T
22
n
v
T
33
_ _
=
S
22
S
23
S
32
S
33
_ _
1
n
v

22
n
v

33
_ _
(138)
and
n
v
T
11
=
n
v
11
S
12
n
v
T
22
S
13
n
v
T
33
(139)
References
[1] J. Aboudi, Mechanics of Composite Materials. A Unied Micromechanical
Approach, Studies in Applied Mechanics, vol. 29, Elsevier, 1991.
[2] J. Aboudi, Mathematik und Physik 51 (2000) 114134.
[3] Y. Bansal, M.-J. Pindera, Journal of Applied Mechanics (2005) 177195.
[4] B.A. Bednarcyk, M.-J. Pindera, Journal of Composite Materials 34 (2000) 299
331.
[5] W. Findley, J. Lai, K. Onaran, Creep and Relaxation of Nonlinear Viscoelastic
Materials, Dover Publications Inc., New York, 1981.
[6] S. Gerlach, Modellbildung und Parameteridentikation viskoelastischer
Faserverbundstrukturen, PhD thesis, Institut fr Mechanik, Fachbereich
Maschinenbau, Universitt Kassel, 2003. <http://www.ifm.maschinenbau.uni-
kassel.de/ifm/index.htm>.
[7] B. Gross, Mathematical Structure of the Theories of Viscoelasticity, Hermann,
Publisher in Arts and Science, 1968.
[8] Z. Hashin, Journal of Applied Mechanics 50 (1983) 481504.
[9] Z. Hashin, Mechanics of Materials 11 (1991).
[10] R. Hill, Journal of Mechanical and Physical Solids 11 (1963) 357372.
[11] Y. Hu, F. Ellyin, Z. Xia, Polymer Engineering and Science 41 (11) (2001) 2047
2060.
[12] S. Lekhnitskii, Theory of Elasticity of an Anisotropic Body, Mir Publisher,
Moscow, 1981 (English translation).
[13] Y. Masuko, M. Kawai, Composites Part A 35 (2004) 817826.
[14] A. Matzenmiller, S. Gerlach, Computational Materials Science 29 (3) (2004)
283300.
[15] A. Matzenmiller, S. Gerlach, International Journal for Numerical Methods in
Engineering 63 (2005) 428454.
[16] M. Megnis, J. Varna, Composite Science and Technology 63 (2003) 1931.
[17] A.H. Nayfeh, Perturbation Methods, Pure & Applied Mathematics, A Wiley-
Interscience Series of Texts, Monographs & Tracts, John Wiley & Sons, 1973.
[18] C. Orozco, M.-J. Pindera, American Institute of Aeronautics and Astronautics
Journal 40 (8) (2002) 16191626.
[19] M. Paley, J. Aboudi, Mechanics of Materials 14 (1992) 127139.
[20] M.-J. Pindera, B. Bednarcyk, Composites Part B: Engineering 30 (1999) 87105.
[21] D.N. Robinson, W.K. Binienda, M.B. Ruggles, Journal of Engineering Mechanics
129 (3) (2003) 310317.
[22] R.A. Schapery, Polymer Engineering and Science 9 (4) (1969) 295310.
[23] R.L. Taylor, FEAP A Finite Element Analysis Program, Theory Manual, 8.1 ed.,
Department of Civil and Environmental Engineering, University of California at
Berkeley, Berkeley, USA, 2007. <http://www.ce.berkeley.edu/rlt/feap>/.
[24] R.L. Taylor, K.S. Pister, G.L. Goudreau, International Journal of Numerical and
Methods in Engineering 2 (1970) 4559.
B. Kurnatowski, A. Matzenmiller / Computational Materials Science xxx (2008) xxxxxx 17
ARTICLE IN PRESS
Please cite this article in press as: B. Kurnatowski, A. Matzenmiller, Comput. Mater. Sci. (2008), doi:10.1016/j.commatsci.2008.02.026