Lecture notes: Composite Mechanics

Martin Fagerstrom

Division of Material and Computational Mechanics, Department of Applied Mechanics

Chalmers University of Technology, Goteborg, Sweden
e-mail: martin.fagerstrom@chalmers.se
2
Chapter 1

Elastic anisotropy

1.1 Basics of tensors and index notation

A 3-dimensional vector a in the 3D space spanned by the base vectors e1 , e2 , e3 can be represented in
the following way:

a = a1 e1 + a2 e2 + a3 e3 (1.1)
where ai , i = 1, 2, 3 are the vector components. This can of course also be written as a sum according
to:
3
X
a= ai ei (1.2)
i=1

Now, in order to simplify things, we can drop the summation sign whenever two indices are the same (i
in the vector form above). Thus, we use the short notation according to:
3
X
a= ai ei = ai ei (1.3)
i=1

meaning that any index (e.g. i) can take the numbers 1-3 and when two identical indices are present,
this implicitly means a summation over that index from 1 to 3.
We can also introduce the matrix representation of the vector a according to:

a1
[a]e1 ,e2 ,e3 = a2 (1.4)
a3
where the subscripts indicate that the representation is made with respect to the coordinate system with
base vectors ei .
In the same way, we can write a second order tensor T as

T = T11 e1 e1 + T12 e1 e2 + T13 e1 e3

+ T21 e2 e1 + T22 e2 e2 + T23 e2 e3
+ T31 e3 e1 + T32 e3 e2 + T33 e3 e3 (1.5)
3 X
X 3
= Tij ei ej = { index notation } = Tij ei ej
i=1 j=1

where is the so-called open product between two vectors (resulting in a second order tensor) defined as

A = a b, Aij = ai bj . (1.6)

3
Please note that on matrix form can be visualised as the multiplication of a column vector and a row
vector into a matrix. Furthermore, in Eq. (1.5) there are two indices i and j present precisely two times
meaning that a summation should be conducted both for i = 1, 2, 3 and j = 1, 2, 3. And consequently,
the matrix representation of the second order tensor T is given by:

T11 T12 T13
[T]e1 ,e2 ,e3 = T21 T22 T23 (1.7)
T31 T32 T33

1.2 Repetition of strains

Figure 1.1: Deformation sketch (from Agarwal et al., figure A2-2).

A point P is subjected to the displacement field u with components in the x-direction u(x, y, z),
the y-direction v(x, y, z) and the z-direction w(x, y, z). To identify the normal (longitudinal) and shear
strains at this point, we study the change of length and angle of two infinitesimal and perpendicular line
segments P A and P B, cf. Figure 1.1.
The normal strain in the x-direction at point P is obtained as the ratio between the increase of
length of P A in the x-direction due to deformation and the original length of P A. Thus, we obtain the
longitudinal strain in the x-direction, x as:
0 0
|P A |x |P A|x x + u(x + x, y, z) u(x, y, z) x u(x, y, z)
x = lim = lim = (1.8)
x0 |P A|x x0 x x
In the same manner, we obtain the longitudinal strain in the y-direction, y , as the ratio between the
increase of length of P B in the y-direction due to deformation and the original length of P B:
0 0
|P B |y |P B|y y + v(x, y + y, z) v(x, y, z) y v(x, y, z)
y = lim = lim = (1.9)
y0 |P B|y y0 y y
Finally, we obtain the shear strain (or shear angle) as the sum of change of angles of P A and P B due
to deformation as:
u
v+ v
x x v u+ y y u v u
xy = yx = + = lim + = + (1.10)
x0, y0 x y x y
where it was used that tan , tan for small strains.
Now, generalising these results to three dimensions, the remaining strains are found as:
w(x, y, z)
z = (1.11)
z
w u
xz = zx = + (1.12)
x z
w v
yz = zy = + (1.13)
y z

4
 Figure 1.2: Three-dimensional stress components (from Agarwal et al. figure A2-5).

Generally speaking, the strain is represented by a symmetric second order tensor = ij ei ej where
ei is the ith base vector. However, in this course we will restrict mainly to cartesian x,y,z-coordinate
systems thereby limiting ourselves to the fact that the strain tensor can be represented (with respect to
this coordinate system) by its cartesian matrix representation []xyz according to:
1 1

x xy xz x 2 xy 2 xz
1 1
[]xyz = xy y yz = 2 xy y 2 yz
(1.14)
1 1
xz yz z 2 xz 2 yz z

Please note that we in the following will omit the subscript xyz for convenience.

1.3 Repetition of stress

Generally, the stress at a point in a body is described by the nine components in the stress tensor ij
according to
= ij ei ej (1.15)
where (again) ei is the ith base vector. The stress tensor may be represented by a matrix representation
considering a particular coordinate system. If we again limit ourselves to the cartesian x, y, z-system (as
in Figure 1.2), the matrix representation of the stress tensor may be expressed as

xx xy xz xx xy xz
[] = yx yy yz = {ij = ji } = xy yy yz (1.16)
zx zy zz xz yz zz

Please note that it in the last equality was used that the (Cauchy) stress tensor is symmetric, which
can be realised by studying the moment equilibrium of an infinitesimal parallelepiped element (at point
P (x, y, z)) with sides x, y and z respectively, cf. any basic course in solid mechanics or Appendix
A-2 in the course book.

1.4 Coordinate transformation of stress and strain components

Later, we will consider the transformation of the stress and strain components between a fibre oriented
coordinate system and a global reference coordinate system. For that purpose, we will need the trans-
formation law for a 2D plane stress case for which the two coordinate systems differ by a rotation about
the z axis, cf. Figure 1.3. Here, we will show two ways to derive the expressions for the coordinate
transformations based on i) coordinate transformation of the stress components by consideration of force
equilibrium and ii) coordinate transformation by change of basis for tensors.

5
1.4.1 Coordinate transformation of stress components by consideration of
force equilibrium

0 0 0
Figure 1.3: Definition of coordinate systems xyz and x y z where the latter has been rotated + degrees
0
around the common z axis (z axis) (from Agarwal et al., Figure A2-3).

Figure 1.4: (from Agarwal et al., Figure A2-8).

In order to express the stress components in the x0 y 0 z 0 coordinate system (which is rotated + degrees
around the z axis relative to the xyz coordinate system, cf. Figure 1.3) in terms of the components with
respect to the xyz system, we consider the equilibrium of two infinitesimal triangular elements as in
Figure 1.4. Considering the left part of the figure, with the diagonal area A and Ax = A cos and
Ay = A sin we can via force equilibrium in the x0 direction conclude that:

0
x A = x cos Ax + y sin Ay + xy cos Ay + xy sin Ax
(1.17)
= A x cos2 + y sin2 + 2xy cos sin

0
x = x cos2 + y sin2 + 2xy cos sin (1.18)
0
In the same way, if we consider the equilibrium in the y direction we obtain:

0
xy A = x sin Ax + y cos Ay + xy cos Ax xy sin Ay

(1.19)
= A x sin cos + y cos sin + xy cos2 xy sin2
0
xy = (y x ) sin cos + xy cos2 sin2


(1.20)
Finally, if we do the same for the right part of the figure we obtain after some derivations:

0
y = x sin2 + y cos2 2xy sin cos (1.21)
0
(y x ) sin cos + xy cos2 sin2


xy = (1.22)

6
 We see that Eq. (1.20) and Eq. (1.22) give the same result. Combining Eqs. (1.18),(1.21) and (1.22)
into Voigt matrix form we finally obtain:
0
x x
0
= [T1 ] y (1.23)
0y
xy

xy
where the stress transformation matrix [T1 ] is defined as:

cos2 sin2

2 sin cos
2
[T1 ] = sin cos2 2 sin cos (1.24)
sin cos sin cos cos2 sin2
It should be remarked that the coordinate transformation could just as well have been determined by
using the transformation laws of a second order tensor under change of basis, cf. Appendix A-1, which
yields the same transformation matrix as in Eqs. (1.24).

1.4.2 Coordinate transformation of stress components by by change of basis

for tensors
In order to derive the expression for a coordinate transformation of the stress tensor, we first note that
the new base vectors of the rotated coordinate system, given a rotation about the z-axis as in Figure 1.3,
are obtained as:

e0x = R ex ex = Rt e0x (1.25)

e0y = R ey ey = R ey t
(1.26)
e0z = R ez = ez (1.27)

with the given inter-relationships (e0i ej = Rij ):

e0x ex = cos() (1.28)

e0x ey = sin() (1.29)
e0x ez = 0 (1.30)
e0y ex = sin() (1.31)
e0y ey = cos() (1.32)
e0y ez = 0 (1.33)
e0z ex = 0 (1.34)
e0z ey = 0 (1.35)
e0z ez = 0 (1.36)

whereby the transformation tensor R becomes:

cos() sin() 0
R = sin() cos() 0 . (1.37)
0 0 1
Now, if we consider a change of base for a tensor, say the stress tensor, we have that:
t 0 t 0
= ij ei ej = ij (Rik ek ) (Rjl t
el ) = Rki ij Rjl e0k e0l (1.38)
| {z }
0
kl

By restricting to plane stress, i.e. only accounting for components in the x- and y-directions, Eq. (1.38)
can be written on matrix form as:
 0 0
    
xx xy cos() sin() xx xy cos() sin()
0 0 = (1.39)
yx yy sin() cos() yx yy sin() cos()

which, if in turn rewritten on contracted Voigt form, becomes identical to Eqs. (1.23)-(1.24) (show this!).

7
1.4.3 Coordinate transformation of strain components
Since the strain tensor is an ordinary second order tensor exactly as the stress tensor, they both follow
the same transformation laws. Hence, we can start by concluding that:

0
x = x cos2 + y sin2 + 2xy cos sin (1.40)
0
2
y = x sin + y cos2 2xy sin cos (1.41)
0
2 2


xy = (y x ) sin cos + xy cos sin (1.42)
0 0 0 0
Now, since we have that xy = 21 xy , xz = 12 xz and yz = 12 yz we also have that xy = 21 xy , xz = 12 xz
0 0
and yz = 21 yz which finally yields:

0
x = x cos2 + y sin2 + xy cos sin (1.43)
0
2 2
y = x sin + y cos xy sin cos (1.44)
0
2 2


xy = 2 (y x ) sin cos + xy cos sin (1.45)
or
0
x x
0
= [T2 ] y (1.46)
0y
xy

xy
with

cos2 sin2

sin cos
[T2 ] = sin2 cos2 sin cos . (1.47)
2 sin cos 2 sin cos cos2 sin2
Please note that it is often more convenient to use numbers rather than characters when indicating
the components of a vector or tensor. Hence, we will in the following use that the x-direction cor-
responds to the 1-direction (or the direction of base vector x1 ), the y-directions corresponds
to the 2-direction and the z-direction corresponds to the 3-direction. Consequently, we get
xx 11 , yy 22 , xy 12 etc.

1.5 Hookes law and stiffness and compliance matrices

The most general linear relationship that relates the (Cauchy) stress with the linear strain under the
assumption of small strains can on component form be expressed as the generalised Hookes law as
ij = Eijkl kl (1.48)
which for e.g 11 yields (i=1; j=1; summation over k=1,2, 3 and l= 1,2,3):

11 = E1111 11 + E1112 12 + E1113 13 + E1121 21 + E1122 22 + E1123 23 + E1131 31 + E1132 32 + E1133 33

(1.49)
This can for all components of stress be written on Voigt matrix form as:


11 E1111 E1122 E1133 E1123 E1131 E1112 E1132 E1113 E1121 11
22

E2211 E2222 E2233 E2223 E2231 E2212 E2232 E2213 E2221
22



E E E E E E E E E

33 3311 3322 3333 3323 3331 3312 3332 3313 3321 33





23 E2311 E2322 E2333 E2323 E2331 E2312 E2332 E2313 E2321 23





31 = E3111 E3122 E3133 E3123 E3131 E3112 E3132 E3113 E3121
31 (1.50)
E1211 E1222 E1233 E1223 E1231 E1212 E1232 E1213 E1221

12 12







32 E3211 E3222 E3233 E3223 E3231 E3212 E3232 E3213 E3221 32







E1311 E1322 E1333 E1323 E1331 E1312 E1332 E1313 E1321 13

13





21 E2111 E2122 E2133 E2123 E2131 E2112 E2132 E2113 E2121 21

8
 Please distinguish between the matrix form of a tensor [] resulting in a 3x3 matrix and the Voigt
matrix form of the same tensor {} resulting in a 9x1 (or rather 6x1 due to symmetry) vector.
Please also note that the first row of Eq. (1.50) is the same as the expression for 11 in Eq. (1.49),
but with the terms in a different order. The index ordering of Eq. (1.50) is often used in the literature
and hence also adopted in the current course.
At a first glance it appears as if there in the most general elastic case are 81 elastic constants describing
the relation between stress and strain. Fortunately, due to symmetry arguments it can, cf. below, be
shown that in the most general elastically anisotropic case, 21 independent components are necessary to
describe the relation between and .

1.5.1 Reduction of 27 constants due to symmetry of the strain tensor

Since is symmetric (kl = lk ) we also have that

Eijkl = Eijlk . (1.51)

In order to see how this reduces the amount of independent constants of E, we first conclude that the two
first indices i and j can be combined in 9 different ways. Secondly, for each combination of i and j, we
can find three pairs of identical components of E (Eij12 = Eij21 , Eij13 = Eij31 , Eij23 = Eij32 ) whereby
three components (for each combination of i and j) can be replaced by their counterpart, leading to a
reduction of 9 x 3 = 27 constants (54 left).

1.5.2 Reduction of additionally 18 constants due to symmetry of the stress

tensor
Since is symmetric (ij = ji ) we also have that

Eijkl = Ejikl . (1.52)

With the respect to the first reduction we note that to two last indices (k and l) can be combined in six
independent ways. And for each of these combinations, there are as in the previous case three pairs of
equal components (E12kl = E21kl , E13kl = E31kl , E23kl = E32kl ). This leads to a reduction with another
18 constants (36 left).

1.5.3 Further reduction of 15 constants based on thermodynamical argu-

ments
Assume the existence of a strain energy density function U = U (ij ) with the property:

U
= ij = Eijkl kl (1.53)
ij
 
U
= Eijkl (1.54)
kl ij
Interchanging order of taking the derivative results in:
 
U
= Eklij (1.55)
ij kl
which finally yields:

Ejikl = Eklij (1.56)

This results in a reduction with additionally 15 components. Thus, finally we end up with 21
independent components in the most general anisotropic case for linear elasticity.

9
1.6 Specially orthotropic material
Orthotropic materials have properties that exhibit symmetry with respect to certain planes. Or, as
an alternative interpretation, the elastic constants do not change when the direction perpendicular (or
normal) to the plane of symmetry is reversed.
An example is a composite laminate reinforced with fibres. Thus, we let x1 and x2 be two base vectors
in the lamina plane and x3 to be the base vector pointing out of plane. Then, it can be realised that
the lamina experience the same properties in the x1 -direction meaning that one plane of symmetry
is the x2 x3 -plane. In the same way it can be realised that the properties are also the same in x2 -
direction (x1 x3 -plane is plane of symmetry) and in the x3 -direction (x1 x2 -plane is plane of symmetry)
respectively. Based on these symmetry arguments, it can be shown that a total number of 9 elastic
constants are enough to describe the constitutive relation between stresses and strains, cf. the course
book. Hence, E can be expressed (with due consideration to the symmetry arguments) on Voigt form as:

E1111 E1122 E1133 0 0 0
E1122 E2222 E3322 0 0 0

E1133 E2233 E3333 0 0 0
[E] =
(1.57)
0 0 0 E2323 0 0
0 0 0 0 E1313 0
0 0 0 0 0 E1212
or equivalently with a more contracted notation as i = Cij j , i, j = 1, 2, 3, 4, 5, 6

1
C11 C12 C13 0 0 0
1
C12 C22 C23 0 0 0 2

2






3 C13 C23 C33 0 0 0

3

= (1.58)

23

0 0 0 C 44 0 0 23

13
0 0 0 0 C55 0 13




12 0 0 0 0 0 C66 12

Please note that the forms of Eqs. (1.57) and (1.58) are obtained when the coordinate axes are placed
along the normals to the symmetry planes. In general, if this is not the case, the [E] and [C] matrices
are full. But in this case, the components are not independent and there is still only necessary to have
nine independent material constants to describe the constitutive response.

1.7 Transversely isotropic material

A transversely isotropic material has one plane of isotropy. As an example consider a unidirectional fibre
reinforced composite. For this material, the mechanical properties in all directions perpendicular to the
longitudinal (fibre) direction are considered the same. Since a transversely isotropic material has more
planes of symmetry, the number of independent elastic constants will be further reduced.
If we consider the fibre direction to coincide with the x1 -direction, making the x2 x3 -plane the plane
of isotropy, it can be shown that:

C22 C23
C22 = C33 , C12 = C13 , C55 = C66 , C44 = (1.59)
2
whereby the Voigt form of the relation between the stresses and strains can be written as:


1 C11 C12 C12 0 0 0
1
C C22 C23 0 0 0 2

2 12







3 C12 C23 C22 0 0 0
3

= C22 C23 (1.60)

23

0 0 0 2 0 0 23

13 0 0 0 0 C66 0 13






12 0 0 0 0 0 C66 12

Thus, only five independent elastic constants are necessary to describe the constitutive response of a
transversely isotropic material.

10
1.8 Constitutive relations for a fibre reinforced lamina
In order to establish the constitutive relations for a unidirectional composite (or a single unidirectional
lamina/ply of a composite laminate), we make the basic assumption that it is in the state of plane
stress. Thus, if a coordinate system x1 x2 x3 is placed such that the x1 -axis points in the direction of
the fibres and the x3 -axis points out of the lamina plane we have 3 = 23 = 13 = 0 and Eq. (1.60) is
reduced to

1 Q11 Q12 0 1
2 = Q12 Q22 0 2 (1.61)
12 0 0 Q66 12

| {z }
[Q]

or equivalently in terms of longitudinal and transverse stresses and strains (1 L , 2 T , 12 LT

etc.

L L
T = [Q] T (1.62)
LT LT

The plane stress assumption is generally a valid assumption since in the majority of the structural
applications, composite laminates are loaded in the plane of the laminate. And even if there are normal
stresses present in the out-of-plane direction (e.g. caused by internal or external pressure) these stresses
are often much smaller than the in-plane stresses.
To arrive at the general relation between stresses and strains, i.e. to establish the expressions for the
components of [Q] we note that it can be considered as the superposition of load cases in which only
one of the in-plane stress components are non-zero. Thus, as a starting point, we consider each of these
states individually.

1.8.1 Longitudinal stress L nonzero

Figure 1.5: Deformation due to nonzero longitudinal stress (from Agarwal et al., Figure 5-3).

The resulting strains when only the longitudinal stress L acts on the lamina are:

L
L = (1.63)
EL
L
T = LT L = LT (1.64)
EL
LT = 0 (1.65)
where EL is the longitudinal stiffness and LT is the so-called major Poissons ratio relating the longitu-
dinal stress to the transverse strain.

11
1.8.2 Transverse stress T nonzero

Figure 1.6: Deformation due to nonzero transverse stress (from Agarwal et al., Figure 5-3).

The resulting strains when only the transverse stress T acts on the lamina are:

T
T = (1.66)
ET
T
L = T L T = T L (1.67)
ET
LT = 0 (1.68)
where ET is the transverse stiffness and T L is the so-called minor Poissons ratio relating the transverse
stress to the longitudinal strain.

1.8.3 Shear stress LT nonzero

Figure 1.7: Deformation due to nonzero shear stress (from Agarwal et al., Figure 5-3).

The resulting strains when only the shear stress LT acts on the lamina are:

L = 0 (1.69)
T = 0 (1.70)
LT
LT = (1.71)
GLT
where GLT is the shear modulus of the lamina.

12
1.8.4 Total constitutive relation
If we superimpose the strains from these three states one obtains:

L T
L = T L (1.72)
EL ET
T T L
T = LT L = LT (1.73)
ET ET EL
LT
LT = (1.74)
GLT
or on Voigt matrix form
1
ETTL

L EL 0 L
LT 1
T = EL

ET 0 T (1.75)
LT 1 LT
0 0

GLT

If we invert this relation, we obtain:

EL T L EL

1 0

L LT T L 1 LT T L
L

T =
LT ET ET
T (1.76)
0
LT 1 LT T L 1 LT T L LT


0 0 GLT
| {z }
[Q]

If we make use of the fact that the lamina/ply is transversely isotropic we can also conclude that the
matrix [Q] is symmetric whereby

T L EL LT ET
= LT ET = T L EL (1.77)
1 LT T L 1 LT T L

1.9 Constitutive relations in terms of global coordinates

In order to make use of the relation (1.77) in analyses of composite laminates, constitutive relations need
to be transformed from fibre oriented coordinate axes LT to global coordinate axes xy. To do so, consider
a lamina with fibres oriented by a rotation angle + around the out-of-plane coordinate axis according
to Figure 1.8. It can be shown that the stresses and strains can be transformed from the fibre oriented
coordinate system to the global coordinate system as:

L x
T = [T1 ] y (1.78)
LT xy


L x
T = [T2 ] y (1.79)
LT xy

where [T1 ] and [T2 ] are the stress transformation matrix and the strain transformation matrix respectively
defined in Eq. (1.24) and Eq. (1.47).
Thus, the relation between stresses and strains in the global coordinates is obtained as

x x
1
y = [T ] [Q][T2 ] y (1.80)
| 1 {z
xy xy
}
[Q]

13
Figure 1.8: UD lamina with fibres oriented + degrees.

14
Chapter 2

Lamina theory

2.1 Volume and weight fractions

A very important factor influencing the properties of a composite is the relative portions of the matrix
and the fibre materials. This can be characterised in (at least) two different ways: by weight fractions,
Wi or volume fractions Vi where i stands for either the matrix (m) or the fibre (f ) material. The
weight fractions are easier to determine from the manufacturing or by subsequent experiments (due to
mass conservation), whereas the volume fractions is used exclusively in the theoretical analysis of the
properties of the composite material.

2.1.1 Weight fraction

If the total weight of the composite material is denoted wc and the total weights of the matrix and fibre
material are denoted wm and wf respectively, we have that

wc = wm + wf (2.1)

and the weight fractions of matrix material (Wm ) and fibre material (Wf ) can be defined as:
wm wf
Wm = , Wf = (2.2)
wc wc

2.1.2 Volume fractions

If the total volumes of the matrix and fibre material are denoted vm and vf respectively, and if one in the
initial stage neglects the existence of pores in the composite material, the total volume of the composite
material may be expressed as
vc = vm + vf . (2.3)

Consequently, the volume fractions of matrix material (vm ) and fibre material (vf ) can be defined as:
vm vf
Vm = , Vf = (2.4)
vc vc

Since the weight fractions are easier to determine from experiments, relations between these weight
fractions and the volume fractions are of importance. Introducing c , m and f for the density of the
composite material, the matrix material and the fibre material, the relation between the volume and
weight fractions are obtained as
wf f vf f
Wf = = = Vf (2.5)
wc c vc c
m
Wm = Vm (2.6)
c

15
 Remains then to determine the density of the composite material (neglecting the existence of pores).
To do this, we first note that the total volume of the composite may be written as
wc wm wf wf m + f wm
vc = = + = . (2.7)
c m f f m

From this expression, we can obtain c as

wc m f 1 1
c = = wf wm = . (2.8)
wf m + f wm wc f + wc m
(Wf /f ) + (Wm /m )

It should be remarked that in the presence of voids, the actual composite density is somewhat lower.
In general, this discrepancy might be rather small (less than 1% for a good composite) but a difference of
up to 5% can be expected for a poorly produced composite. However, the actual effect of a higher amount
of voids may have significant influence on some of the properties such as a lowered fatigue resistance and
strength.

2.2 Analytical and semi-empirical methods for predicting lam-

ina properties of unidirectional composites
The properties of a composite material depend on the properties of its constituents, their concentrations,
distributions and orientations as well as their physical and chemical interaction. The most direct method
to determine the properties of a composite material is by experimental methods. In many cases,
rather simple methods can be used to determine e.g. longitudinal stiffness and (tensile) strength of a
lamina. However, the experimental results obtained are only valid for that particular fibre matrix system
in terms of volume fractions, constituent properties and production method. If any one of these are
changed, new experiments need to be performed in order to establish the properties of the new system.
Such extensive testing may be very time consuming and costly, whereby analytical or semi-empirical
models can be a very valuable tool to predict some of the properties. Some of the available models, to
obtain homogenised properties of the composite, will be discussed below in this section.

2.2.1 Longitudinal properties

2.2.1.1 Longitudinal stiffness prediction based on the constant stress assumption - the
Voigt assumption
To predict the longitudinal stiffness of a unidirectional composite lamina, the fibres are assumed as
uniform in properties and size and parallel throughout the composite, cf. Figure 2.1. Furthermore,
prefect bonding between the matrix and the fibres is assumed. This implies that the longitudinal strains
experienced by the fibres f L , the matrix mL and the composite cL are the same:

f L = mL = cL (2.9)
This is in general homogenisation theory denoted the Voigt assumption which serves as an upper limit
of the possible stiffness of a composite in any direction, given a certain fibre volume fraction.
Based on the assumptions above, the longitudinal load Pc carried by the composite will be shared
between the fibres Pf and the matrix Pm so that

Pc = c A c = Pf + Pm = f A f + m A m (2.10)
where c , f and m are the stresses experienced by the composite, the fibres and the matrix respectively,
and where Ac , Af and Am are the corresponding total cross sectional areas. From Eq. (2.10), we can
get the expression for the composite stress c which, by using the fact that the area fractions equals the
volume fractions of a unidirectional composite, takes the form

Af Am
c = f + m = f Vf + m Vm . (2.11)
Ac Ac

16
Figure 2.1: Model for predicting longitudinal behaviour of unidirectional composite (from Agarwal et al.,
Figure 3-3).

If we now assume that both the fibres and the matrix behaves linear elastic ( = E) we obtain the
expression of the longitudinal elastic stiffness (modulus) of the composite as

EL = Ef Vf + Em Vm (2.12)

which generally is denoted the rule of mixtures (ROM). For a composite of n constituents, the expression
for the longitudinal stiffness can be generalised as
n
X
EL = Ei Vi . (2.13)
i=1

From Eqs. (2.12)-(2.13) and in Figure 2.2 it can be seen that, for the case of significantly stiffer fibres
(compared to the matrix), most of the (longitudinal) load is carried by the fibres already at rather low
fibre volume fractions.
It should be remarked that Eqs. (2.12)-(2.13) are only valid as long as both the matrix and the
fibre material behaves linear elastic. This may however constitute only a small portion of the stress-
strain behaviour and generally, the longitudinal deformation of a unidirectional composite proceed in
four stages:

1. The matrix and the fibres both deform linear elastically

2. The fibres still deform linear elastically whereas the matrix deforms nonlinearly elastic or even
plastically

3. The fibres and the matrix both deform in a nonlinear fashion

4. The fibres break followed by composite failure

Whereas stage 3 can be observed only for ductile fibres, stage 2 may occupy the largest portion
of the composite stress-strain curve (which is no longer linear) and the longitudinal composite stiffness
(modulus) relating an incremental change in strain to the corresponding change in stress (L = EL L )
must be predicted at each composite strain level c as
 
m
EL = Ef Vf + Vm (2.14)
m c
 
m
where is the slope of the matrix stress-strain curve at strain c . However, in practice, the
m c
non-linearity of the stress-strain curve for the matrix material has low effect on the composite stiffness,
especially at significant fibre fraction. Thus, in many cases Eq. (2.12) is a good approximation.

17
 It should be remarked that the rule of mixtures is accurate for longitudinal tensile stiffness. But
when loaded in compression, the response of the composite observed in experiments may
deviate from the analytical predictions (ROM). This because the combined compressive properties
is strongly dependent on the matrix material properties, such as its shear stiffness (cf. buckling), whereas
the tensile response is much more governed by the fibre properties.

Figure 2.2: Graph showing the percentage of the total force carried by the fibres in a unidirectional
composite as function of fibre volume fraction (from Agarwal et al., Figure 3-5).

2.2.1.2 Tensile strength

Failure initiates when the fibres are subjected to their fracture strain, assuming that the fibre failure strain
f is less than the matrix failure strain, which is generally the case. In most practicable applications with
a sufficiently high volume fraction of fibres and where the matrix material behaves linearly elastic up to
fibre failure, the failure stress of the composite cu (cu for composite ultimate) can be expressed by the
rule of mixtures as:

cu A = f u Af + (m ) Am
f
 
Af Am Af Am
cu = f u + (m ) = Vf , = Vm
A f A A A
 
n o Em
cu = f u Vf + (m ) Vm = matrix linear elastic up to f = f u Vf + Vm (2.15)
f Ef

where f u is the fibre fracture stress (or ultimate strength of the fibres) and (m ) is the stress in a
f
matrix subjected to the strain at which the fibres fail (f ), cf Figure 2.3.
It can be seen that Eq. (2.15) predicts a lower composite strength for the composite compared to the
unreinforced matrix material for a certain level of fibre volume fraction (below Vcrit ). Vcrit is defined as
when the failure strength of the composite equals the failure strength of the unreinforced matrix material,
i.e.

mu (m )
f
f u Vcrit + (m ) (1 Vcrit ) = mu Vcrit = . (2.16)
f | {z } f u (m )
f
Vm
However, for most practical applications (Vcrit ) is really small. Consider e.g. the common composite
material consisting of carbon fibres and epoxy matrix (or resin) material. From mean values of the
respective strength for carbon fibres (cf. Table 1-1 in Agarwal et al.) f u 2.3 GPa and epoxy (cf.
Table 2-11 in Agarwal et al.) mu 92.5 MPa and assuming that the epoxy is linear elastic up to the

18
Figure 2.3: Graph showing the composite longitudinal failure stress cu as function of fibre volume
fraction (from Agarwal et al., Figure 3-7).

9
failure strain of the fibres f = 2.3 109 /315 109 = 0.0073 (m ) = |3.425
{z 10} 0.0073 = 25 MPa, we
f
Eepoxy
can obtain Vcrit as:
92.5 106 25 106
Vcrit = 0.03 (2.17)
2.3 109 25 106
That is, the critical volume fraction of fibres Vcrit is about 3% which is considerably lower than in any
practical applications (50-60%).

2.2.1.3 Compressive strength

It is unlikely that the fibres themselves would break due to compressive stresses. But when being subjected
to compressive loads the fibres acts as long columns and either so-called micro-buckling or fibre kinking
can occur. Of course, the buckling load of a fibre embedded into a matrix is significantly larger than
for a free fibre or fibre bundle. But still, fibre (bundle) buckling can occur even when the corresponding
matrix stresses are in the elastic range. However, for practical fibre volume fractions (Vf > 0.4), fibre
buckling of often preceded by other failure phenomena to be discussed further in the chapter covering
failure.

2.2.2 Factors influencing longitudinal stiffness and strength

There are a number of factors that may influence the longitudinal stiffness and strength of a composite
leading, in some cases, to significant deviations from the predictions discussed above.

Misorientation of fibres
Fibre orientation directly influences the properties of the composite. Naturally, the contribution
from the fibres is maximised only when the fibres are aligned with the loading direction. As a
consequence, the stiffness and strength will be reduced when the fibres are not parallel to the
loading direction. However, the discrepancy is small and no corrections have to be made if the
misalignment is small, i.e. limited to a few degrees. It should be remarked that in the general
case, a composite material is often composed of a number of unidirectional lamina with different
orientations stacked on top of each other. For this type of composite structure a laminate there
are several appropriate theories (co-called laminate theories) available to described the structural
response, cf. Chapter 3 and Chapters 6-8 in the course book.
Fibres of nonuniform strength
First of all, it should be remarked that any reduction of the fibre strength directly results in a
lowered strength of the composite material. Consequently, a composite of high strength is obtained
when the fibres are uniform in strength. Reasons for a non-uniform distribution may be due to at
least the following two reasons:

19
 1. A variation of cross section as function of length
2. As a result of initial damage due to handling of the fibres before manufacturing the composite

In addition, the fibre strength decreases as the fibre length increases. This is due to the fact that
statistically, the existence of any strength reducing factor (flaw) increases with length and the fibre
will always first break at its weakest link. In any case, it should be remarked that fibres start
to break at loads lower than the composite strength and that this accumulates up to final failure.
Thus, for detailed predictions of the composite strength, statistical methods need to be incorporated
which, however, is considered as out of scope in the current course.

Discontinuous fibres
When the length of the fibres are in the same order of size as the length over which the load is
transmitted through the matrix material, we speak of discontinuous fibre reinforcements. In these
materials, the end effects (variation of stress along the fibre and stress concentrations) cannot be
neglected and must be dealt with. This is however not included in the current course. Interested
students are referred to Chapter 4 in the course book for reference.

Interfacial conditions

More pronounced importance for discontinuous-reinforced-fibre composites.

Influences the transversal strength.
Good adhesion between fibre and matrix enhances water resistance.
A lower interface strength may lead to a higher composite ductility and, hence, a high fracture
toughness.

Residual stresses
Originates predominantly from the manufacturing process due to e.g.

a variation between thermal expansion coefficients for the fibre and the matrix material,
a significant difference in temperature at manufacturing and temperature at use, cf. also
Computer Assignment 1, and
chemical shrinkage in the manufacturing process

It should be remarked that the existence of residual stresses significantly impact the response
and the properties (strength) of the composite material and should be included in an advanced
assessment analysis of a composite structure. How this is to be incorporated practically is not always
straightforward, but this should definitively be kept in mind! The residual stresses in composites
is much more importance for the component performance than say in the case of traditional metal
components in which the residual stress often decrease with time due to local (and often harmless)
plastic deformations leading to stress relaxation.

2.2.3 Transverse properties

2.2.3.1 Transverse stiffness prediction based on the constant stress assumption - the Reuss
assumption
In order to derive a simple expression for predicting the transverse stiffness of a lamina, let us consider a
lamina with unidirectional (UD) continuous fibre reinforcements loaded by a transversal tensile stress c .
Schematically, we study a Representative Volume Element (RVE) as a generic block consisting of fibre
material bonded to matrix material, as shown in Figure 2.4.
For this case, the elongation of the composite in the loading direction (c ) is obtained by the sum of
the elongation of the matrix material and of the fibre material in the loading direction according to

c = cT L2 = f T Lf + mT Lm (2.18)

20
Figure 2.4: Model for predicting transverse stiffness of unidirectional composite (from Gibson, Principles
of Composite Material Mechanics 2nd ed., 2007).

where cT , f T and mT are the transversal strains in the composite RVE, the fibre materials and the
matrix material respectively. Furthermore, since the RVE do not change along the longitudinal direction,
the length fractions must equal the volume fractions such that

cT = f T Vf + mT Vm . (2.19)

By using Hookes law, (and neglecting any Poisson strains) this can be written as
c f m
= Vf + Vm . (2.20)
ET Ef Em

If we now finally assume that the stress in each of the constituents (fibre and matrix) is the same (the
so-called Reuss assumption), i.e. we have that

c = f = m , (2.21)

we end up with the inverse rule of mixtures for the transverse modulus as
1 Vf Vm
= + (2.22)
ET Ef Em
It should be remarked that by assuming equal stress in both fibres and matrix, the resulting expression
for the (transverse) stiffness can be shown to be the lower bound given a certain fibre volume fraction,
cf. also Figure 2.5.

inverse rule of mixtures

Figure 2.5: Longitudinal and transverse stiffness as function of fibre volume fraction (from Agarwal et
al., Figure 3-9a). Please note that the rule of mixtures (Voigt assumption) serves as an upper limit of the
composite stiffness whereas the inverse rule of mixtures (Reuss assumption) serves as the lower bound

21
2.2.3.2 Halpin-Tsai semi-empirical model for transverse stiffness prediction
The inverse rule of mixtures may lead to poor predictions (and an underestimation of the transverse
stiffness) since, in reality, the stresses are not equal in the matrix and fibre material.
To improve the predictions, Halpin and Tsai developed simple and generalised equations to fit more
advanced micromechanical models for transverse stiffness. The Halpin-Tsai equation for the transverse
composite stiffness (modulus) can be written as

ET 1 + Vf (Ef /Em ) 1
= , = (2.23)
Em 1 Vf (Ef /Em ) +
in which is a measure of the reinforcements and relates to the fibre geometry. Please note the special
cases
1
=0: ET = Vm Vf
(inverse rule of mixtures) (2.24)
Em + Ef
=: ET = Vm Em + Vf Ef (Rule of mixtures) (2.25)

Thus, it is clear that can be viewed as a curve fitting parameter since the real transverse stiffness will
lie between the two limiting cases (closer to the inverse rule of mixture prediction). Halpin and Tsai have
proposed = 2 for circular or quadratic cross sections and = 2 ab for rectangular cross sections
in which a and b are side lengths of the cross section and where a is to be taken as the side length in the
loading direction.
Predictions made by the Halpin-Tsai equations have proven to be adequate in many practical situa-
tions, cf. Figure 2.6 for a comparison made against experimental data for a boron-epoxy lamina reported
by Whitney and Riley (J. IAAA, 4:1537, 1966), and it generally gives a better prediction than the pre-
diction obtained based on the equivalence of stress in the fibre and the matrix material (the inverse rule
of mixtures).

inverse rule of mixtures

Figure 2.6: Transverse stiffness predicted by the inverse rule of mixtures and the Halpin-Tsai equation
compared with experimental data reported by Whitney and Riley for boron-epoxy (from Agarwal et al.,
figure 3-9b).

2.2.3.3 Transverse strength

The transverse strength of a composite is reduced by the existence of fibres. The reason is that due to
their geometry, the fibres cannot carry a large portion of the load which instead is distributed between
the fibre and the matrix material. Instead, the existence of fibres place restrictions on the transverse
deformations, causing strain and stress concentrations in their vicinity, cf. Figure 2.7 which results in an
overall decreased strength in comparison to the unreinforced matrix material.
Thus, the composite transverse failure strength T U can be expressed as
mu
T U = (2.26)
S

22
Figure 2.7: Stress distribution in matrix surrounding a single cylindrical inclusion EF /Em = 10, m =
0.35, f = 0.3 (from Agarwal et al., figure 3-12a).

in which S is a strength reduction factor. This reduction factor can be predicted by several methods.

The strength-of-materials method

Based on the strength-of-materials method, the factor S is assumed e.g. to be equal to the stress
concentration factor SCF which can be predicted by neglecting Poisson effects, cf. Greszczuk Society of
Plastics Industry, 21st Annual Conference, 1966 for details, as:

1 Vf (1 (Em /Ef ))
SCF = 1/2
. (2.27)
1 (4Vf /) (1 (Em /Ef ))

Based on Finite Element Analysis

Another alternative is to analyse the stress distribution in a UD lamina by FEM, cf. the results in
shown in Figure 2.8 obtained by e.g. Chen and Lin (Mater. Res. Stand. MTRSA, 9:29-33, 1969).

Figure 2.8: Principal stress in matrix surrounding multiple fibres obtained by FEM by Chen and Lin
m = 0.35, f = 0.2 (from Agarwal et al., figure 3-12b).

2.2.4 Shear modulus

The derivations of a simple rule for predicting the in-plane shear stiffness of a lamina is similar to those
for the transverse stiffness (EL ). Study Section 3.4 in the course book on your own.

23
2.2.4.1 Major results
Inverse rule of mixtures:
1
GLT = (2.28)
Vf /Gf + Vm /Gm

Improved predictions by means of Halpin-Tsai:

GLT 1 + Vf (Gf /Gm ) 1

= , = , =1 (2.29)
Gm 1 Vf (Gf /Gm ) +

2.2.5 Poissons ratio

Study subsection 3.5 in the course book on your own.

2.2.5.1 Major results

Major Poissons ratio LT
For this case, the underlying assumptions in the derivations are:

The longitudinal strain is the same in both constituents (as for the stiffness prediction)

The Poissons ratios for each constituent is (transversely) isotropic

The total transverse contraction is the sum of the transverse contraction of the fibres and matrix
respectively

Major Poissons ratio LT relating the longitudinal stress to the transverse strain:

LT = f Vf + m Vm (2.30)

Minor Poissons ratio T L

For this case it is used that the stiffness matrix of a linear elastic material is symmetric.

Minor Poissons ratio T L relating the transverse stress to the longitudinal strain:

LT T L ET
= T L = LT (2.31)
EL ET EL

2.3 Expansion coefficients

Study 3.7.1 and 3.7.2 in the course book on your own.

For the thermal case, the underlying assumptions (not explained in the text) in the derivations are:

The fibres are assumed to be isotropic and linear elastic

The matrix is assumed to be isotropic and linear elastic

The strains in the longitudinal direction are the same in the matrix and in the fibres (Voigt as-
sumption)

The stresses in the transverse direction are constant, i.e. the Reuss assumption is used

The homogenised macroscopic stresses in the lamina are zero, i.e. the assumption L = T = 0 is
used.

24
Major results
Thermal expansion coefficients (cf. Scharpery, J. Compos Mater., 2:280-404, 1968 for details):

1
L = (f Ef Vf + m Em Vm ) (2.32)
EL
T = (1 + f ) f Vf + (1 + m ) m Vm L LT (2.33)

Moisture expansion coefficients:

L 0 (given that the fibres are much stiffer than the moisture absorbing matrix) (2.34)
c
T = (1 + m ) m (cf. e.g. Tsai and Hahn, Introduction to Composite Materials, 1980)(2.35)
m

2.4 Hygrothermal degradation of polymer matrix properties in-

fluencing the composite parameters (included for your in-
formation and understanding)
Review of Subsection 5.1-5.2 in Gibson (Principles of Composite Material Mechanics 2nd ed., 2007):

Increased temperature of a polymer causes a gradual degradation of the stiffness up to a certain point,
the so-called glass transition temperature at which the matrix behaviour transforms from being glassy
to being rubbery, cf. also Figure 2.9

Figure 2.9: Variation of stiffness for a typical polymer showing the glass transition temperature Tg and
the effect of absorbed moisture on Tg . Note: Tg0 = dry, Tg0 = wet (from Gibson, figure 5.1.).

In addition, increased moist absorption in the matrix leads to a lowered glass transition temperature
and a decrease in stiffness. Moist saturation 3-4 % leads to a lowered glass temperature of approx 20%, c.f
also Figure 2.10 how the stress-strain curve for 3501-5 epoxy resin is affected by temperature and moist
and Figure 2.11 how this in turn influence the stress-strain of a AS/3501-5 graphite/epoxy composite
(same matrix material!) with fibre volume fraction Vf = 0.63.
The hygrothermal degradation of a composite strength and/or stiffness is often estimated by empirical
models. Chamis and Sinclair (Composite Materials: Testing and Design (Sixth Conference), ASTM STP
787, pp. 498512, 1982) and Chamis (Engineers Guide to Composite Materials, pp. 3-83-24, 1987)
have demonstrated such a procedure. In their approach, the degraded stiffness Em (or strength) is
0
approximated by the nominal stiffness Em multiplied by a degradation factor Fm such that:

0
Em Fm Em (2.36)
and the degraded value Em is then used in e.g. the Halpin-Tsai equation to approximate the temperature
and moisture dependent transverse stiffness of the composite. As for the degradation factor, the following

25
Figure 2.10: Stress-strain curves for 3501-5 epoxy resin at different temperatures and moisture contents
(from Gibson, figure 5.2.)

Figure 2.11: Stress-strain curves for AS/3501-5 graphite/epoxy composite with fibre volume fraction
Vf = 0.63 at different temperatures and moisture contents (from Gibson, figure 5.3.)

expression is proposed:
 1/2
Tgw T
Fm = (2.37)
Tg0 T0
where Tg0 is the dry glass transition temperature of the matrix material (i.e. without moisture), Tgw
is the wet glass transition temperature that depends on the moisture concentration, T0 is the reference
temperature at which the nominal matrix stiffness is measured (E0 ) and T is the current temperature.
Furthermore, Chamis used data for six different epoxy resins to fit the parameters of an empirical curve
in order to get an approximate expression for Tgw for epoxy resin as:


2

32 50Mr2 10Mr
Tgw [C] = 50Mr 10Mr + 1.0 Tg0 [C] + . (2.38)
1.8
where Mr is moist contents in weight percent (4% Mr = 0.04), cf. also Figure 2.12 for tabulated
data for some matrix materials. Please note that the somewhat strange format of the expression for Tgw
stems from the fact that the parameters obtained by Chamis where found for temperatures measured in
Fahrenheit whereas the expression in Eq. (2.38) is valid for temperatures measured in Celsius. Thus, it

26
has been used that
T [F ] = T [C] 1.8 + 32 (2.39)

Figure 2.12: Hygrothermal properties of various polymer matrix materials (from Gibson, Table 5.1.).

27
28
Chapter 3

Laminate theory

Figure 3.1: A four-ply laminate (from Agarwal et al, figure 6-1.)

3.1 Laminate strains

Figure 3.2: Deformation during bending of the laminate in the xz plane (adapted from Agarwal et al,
figure 6-2.)

Consider first the deformation of a laminate in the xz-plane, cf. Figure 3.2. In this case, the displace-
ment u in the x-direction of a point C that is located on the undeformed normal ABCD, at a distance z
from the midplane, is given by
u = u0 z sin u0 z (3.1)
where u0 is the midplane displacement in the x-direction and is the slope of the deformed section
A0 B 0 C 0 D0 or the (negative) rotation of the normal ABCD around the y-axis due to deformation. Sim-
ilarly, the displacement v of a corresponding point when the laminate is deformed in the yz-plane, cf.
Figure 3.3, can be obtained as
v = v0 z (3.2)

29
where v0 is the midplane displacement in the y-direction and is the slope of the deformed section
A0 B 0 C 0 D0 or the (positive) rotation of the normal ABCD around the x-axis due to deformation in the
yz-plane).

Figure 3.3: Deformation during bending of the laminate in the yz plane (adapted from Agarwal et al,
figure 6-2.)

From these kinematical expressions for the displacements, we can derive the in-plane strains as

u u0
x = = z (3.3)
x x x
v v0
y = = z (3.4)
y y y
 
u v u0 v0
xy = + = + z + (3.5)
y x y x y x

The preceding strain-displacement relation can be written as

x 0x kx
y = 0 +z ky (3.6)
0y
xy xy kxy

or
{} = {0 } + z{k} (3.7)

where the midplane strains 0 and plate curvatures k are defined as

u0




0
x
x


v0

0 0
{ } = = (3.8)
0y y
xy


u0 v0


+

y x
and






kx


x




{k} = ky = (3.9)

kxy


y



+

y x

It should be remarked that the definition (and notation) of the curvatures differs in the literature. In these
notes, we follow the notations used in the course book. Please note that we assume that adjacent
lamina across the thickness of the laminate do not slip over each other. Hence, the strain
varies linearly through the thickness even though the properties vary, cf. Figure 3.4. Please
also note that both 0 and k are functions of x and y only.

30
 If we now assume that the normal ABCD remains planar and perpendicular to the midplane also after
deformation (the Kirchhoff assumption), it can be found that

w0
= (3.10)
x
w0
= (3.11)
y

whereby one obtains the plate curvature as

2
w0

x2



kx

2w


0
ky = 2
(3.12)
y
kxy Kirchhof f



2 w0




2


xy

where w0 is the vertical displacement of the midplane (in the z-direction).

Please note that laminate theory based on the Kirchhoff assumption, which also implies that the
transverse shear strains xz and yz are negligible, is often denoted Classical Lamination Theory, CLT.

3.2 Laminate stresses

For a material with identical properties through the thickness, the stresses will, as the strains, vary
linearly. Thus, within each lamina, the stresses are a linear function of z. However, since the properties
in general are different for each lamina the properties are considered constant within each lamina the
stress distribution through the thickness will be discontinuous but piecewise linear in each lamina, cf.
Figure 3.4.

Figure 3.4: Principal variation of stresses and strains in a composite laminate consisting of three unidi-
rectional laminae/plies with different orientation (from Agarwal et al, figure 6-3.)

Given the relation between stresses and strains from Chapter 1:

x Q11 Q12 Q16 x
y = Q12
Q22 Q26 y (3.13)
xy Q16 Q26 Q66 xy

we can now obtain the stresses at any point in lamina k as

0
x Q11 Q12 Q16 x Q11 Q12 Q16 kx
y = Q12 Q22 Q26 0y + z Q12 Q22 Q26 ky (3.14)
0
xy k Q16 Q26 Q66 k xy Q16 Q26 Q66 k kxy

31
3.3 Resultant forces and moments
For convenience, we define the stress resultants in terms of normal forces (per unit length) and moments
(per unit length), obtained by integration of stresses through the thickness of the laminate, as (cf. Fig-
ure 3.5 below for sign conventions):

Figure 3.5: Sign convention for resultant forces and moments.

Resulting force per unit length in x-direction:

Z h/2
Nx = x dz (3.15)
h/2

Resulting force per unit length in y-direction:

Z h/2
Ny = y dz (3.16)
h/2

Resulting in-plane shear force per unit length:

Z h/2
Nxy = Nyx = xy dz (3.17)
h/2

Resulting moment per unit length, acting on the edge with normal in the x-axis direction causing a
(positive) rotation around the y-axis (integration of stress times moment arm):
Z h/2
Mx = x z dz (3.18)
h/2

Resulting moment per unit length, acting on the edge with normal in the y-axis direction causing a
(negative) rotation around the x-axis (integration of stress times moment arm):
Z h/2
My = y z dz (3.19)
h/2

Resulting moment per unit length, acting on the edge with normal in the x-axis direction causing a
(negative) rotation around the x-axis (integration of stress times moment arm):
Z h/2
Mxy = xy z dz (3.20)
h/2

Resulting moment per unit length, acting on the edge with normal in the y-axis direction causing a
(positive) rotation around the y-axis (integration of stress times moment arm):

Myx = Mxy (3.21)

32
 Let us now consider an orthotropic laminate consisting of n laminae (plies), cf. Figure 3.6. Grouping
the normal forces (per unit length) in a vector one obtains:

Nx Z h/2 x Xn Z hk x
{N} = Ny = y dz = y dz (3.22)

Nxy
h/2
k=1 hk1

xy k

xy

where hk is the z-coordinate of the upper (in positive z-direction) part of lamina k.

Figure 3.6: Description of a multilayered laminate geometry (from Agarwal et al, figure 6-5).

Now, combining Eq. (3.14) and Eq. (3.22), N can be written as:
n Z
X hk    0
{N} = Q k + zk dz (3.23)
k=1 hk1

Realising that the midplane strains 0 , the curvatures k and the material properties Qij are constant
within each lamina, we can further rewrite the expression for the normal forces as:

n
!
X Z hk Z hk
0
   
{N} = Q k dz { } + Q k z dz {k}
k=1 hk1 hk1
" n # " n #
X 1 X  (3.24)
Q k (hk hk1 ) {0 } + Q k h2k h2k1 {k}
 

=
2
k=1 k=1
= [A]{0 } + [B]{k}

where

n
" #
X  
[A] = Q k (hk hk1 ) (3.25)
k=1
" n #
1 X  2 2


[B] = Q k hk hk1 (3.26)
2
k=1

(hk hk1 ) is always positive (equal to the thickness of the lamina), but that the
Please note that
term h2k h2k1 is positive for a lamina situated above the midplane and negative for a lamina situated
below the midplane. In this way, B = 0 for laminates where the laminae (plies) are placed symmetrically
around the midplane, thus in that case there is no coupling between curvatures and normal forces
which may occur in the general case.
Performing the same type of analysis for the moments, cf. e.g. Agarwal et al., one obtains:

33

Mx
{M} = My = [B]{0 } + [D]{k} (3.27)
Mxy

where " n #
1 X  3 3


[D] = Q k hk hk1 (3.28)
3
k=1

To summarise, we can write the relation between N-M and 0 -k as:

0
    
N A B
= (3.29)
M B D k

3.4 Laminate orientation code

Each lamina in a laminate is normally defined by its rotation around the positive out-of-plane axis
(positive z-axis) with respect to a global xyz coordinate system. Thus, referring to Figure 3.7 it is clear
that a proper definition of the global coordinate system is required to avoid ambiguity in the definition
of the lay-up sequence of the laminate.

Figure 3.7: Sign convention for orientation of laminae in a laminate. Please note that in the figure, the
laminate is identical in both the left and the right part. The only difference is the orientation of the
coordinate axes (from Agarwal et al, figure A3-1).

With a given coordinate system at hand, the lay-up sequence of a laminate is easily defined according
to the following:

Each lamina is represented by a number representing the angle (in degrees) of rotation around the
z-axis

Individual laminae are separated by a slash sign (/)

The laminae are listed in the sequence they are laid up.

Adjacent laminae with the same orientation are denoted by a numerical subscript.

Adjacent laminae with the same magnitude of orientation, but in the two different directions, are
denoted by a plus-minus () or minus-plus () sign.

For symmetric laminates, only half of the lay-up sequence needs to be defined followed by a subscript
S denoting symmetry.

For symmetric laminates, if the number of plies are uneven, the mid-lamina is denoted by an
overbar.

For a couple of examples of laminate definitions, cf. Figure 3.8.

34
 s

Figure 3.8: Four examples of lay-up sequences and their corresponding laminate orientation code (from
Agarwal et al).

3.5 Read Subsection 6.6 on your own

3.6 Determining relevant stresses and strains in the laminae

It should be remembered that the interesting stresses (and strains) of a composite laminate (in terms of
strength) are (L , T , LT ) (and (L , T , LT )). Given the A, B and D matrices as well as the applied
loads at a certain point (x0 , y0 ), the procedure to determine these stresses (and strains) in each lamina
of the laminate is as follows:

3.6.1 Calculate the strain distribution through the thickness

Invert the relation between 0 , k and N, M to get the strain distribution:

1 
0 (x0 , y0 )
   
A(x0 , y0 ) B(x0 , y0 ) N(x0 , y0 )
= (x0 , y0 , z) = 0 (x0 , y0 ) + zk(x0 , y0 )
k(x0 , y0 ) B(x0 , y0 ) D(x0 , y0 ) M(x0 , y0 )
(3.30)

3.6.2 Transform strains into fibre oriented coordinate system

Transform the strains of each lamina from the global coordinate system to the fibre oriented coordinate
system using the transformation matrix [T2 ]k of each lamina as:

L (x0 , y0 , z) x (x0 , y0 , z)
T (x0 , y0 , z) = [T2 ]k y (x0 , y0 , z) (3.31)
LT (x0 , y0 , z) k xy (x0 , y0 , z) k

where we remember [T2 ]k from lecture 2 as:

35
 cos2 k sin2 k

sin k cos k
[T2 ]k = sin2 k cos2 k sin k cos k (3.32)
2 sin k cos k 2 sin k cos k cos2 k sin2 k

3.6.3 Calculate stresses in global coordinate system for each lamina/ply

To compute the stress distribution through the thickness of each lamina, use the constitutive matrix [Q]k
as
0
x (x0 , y0 , z) x (x0 , y0 ) + zkx (x0 , y0 )
{ k (x0 , y0 , z)}k = y (x0 , y0 , z) = [Q]k 0 (x0 , y0 ) + zky (x0 , y0 ) (3.33)
0y
xy (x0 , y0 , z) k xy (x0 , y0 ) + zkxy (x0 , y0 )

Note that the stress varies linearly with z within the lamina!! Please also note that this is an
unnecessary step if you want to assess the strength of the lamina since you then would need the stresses
expressed in the fibre oriented coordinate system, obtained directly through

L L (x0 , y0 , z)
T = [Q]k T (x0 , y0 , z) (3.34)
LT (x , y , z)

k LT 0 0 k

3.6.4 Transform stresses into fibre oriented coordinate system

Transform the stresses of each lamina from the global coordinate system to the fibre oriented coordinate
system using the transformation matrix [T1 ]k of each lamina as:

L x
T = [T1 ]k y (3.35)
LT

k xy k
where we remember [T1 ] from lecture 2 as:

cos2 k sin2 k

2 sin k cos k
[T1 ]k = sin2 k cos2 k 2 sin k cos k (3.36)
sin k cos k sin k cos k cos2 k sin2 k
Note that the stress varies linearly with z within the lamina!!

3.7 Hygrothermal effects

There are two major sources of change of mechanical behaviour due to temperature (thermal effects) and
moist absorption (hygroscopic effects) of polymer composites:
Hygrothermal stresses and strains
Hygrothermal degradation of properties (matrix dominated properties such as (transverse) stiffness
and strength)

3.7.1 Hygrothermal stresses

Hygrothermal stresses are stresses induced by temperature change or moist absorption/desorption. Please
note that hygrothermal stresses can occur within the different laminae (plies) even though the laminate
itself is free to expand and can at the first glance appear to be stress-free (although the stress resultants
are in fact left unchanged), cf. also Figure 3.9 for an explanatory sketch of a laminate consisting of
isotropic laminae (the principle is the same also for transversely isotropic plies).
Temperature change and moist absorption/desortpion does not, as mentioned above, influence the
stress resultants, but it introduces stresses and strains in each laminae and may cause e.g. warping
for unsymmetric composites. The strains that produce stresses are the mechanical strains:
M = T H . (3.37)

36
Figure 3.9: Concepts of thermal stresses and strains in a three-ply symmetric laminate (from Agarwal et
al, figure 6-19.)

To study the thermal and hygroscopic strains, the hygrothermal strains, we first consider the prop-
erties of an transversely isotropic material, such as an UD lamina. In this case, the coefficients of thermal
and moisture expansion change with direction, cf. e.g. Subsection 3.7.1 and 3.7.3 in Agarwal et al.. Thus,
changes in temperature and/or moist contents produce strains according to:

TL = L T (3.38)
TT = T T (3.39)
H
L = L C = 0 (since L is taken as 0, cf. Subsection 3.7.2) (3.40)
H
T = T C (3.41)
(3.42)

Since the thermal and hygroscopic strains transform in the same way as the total strain (using the
transformation matrix T2 ) it is clear that the coefficients of thermal and hygroscopic expansion transform
in the same way, i.e. we obtain for each lamina:

x L
1
y = [T2 ]k T (3.43)
xy 0

k k
x 0
1
y = [T2 ]k T (3.44)
xy 0

k k

Thus, the thermal and hygroscopic strains may be expressed in the global coordinate system as:

37
 T
x x
T = y T (3.45)
Ty
xy xy

H k
x x
H
y = y C (3.46)
H
xy xy

k

For clarity, due to the similarities between hygroscopic and thermal strains, we will only consider
the thermal strains in the following. The derivations below can easily be extended to include also the
hygroscopic effects, cf. the book, but it involves more terms.
Given the expressions for the thermal strains, and neglecting any moist absorption/desorption, one
obtains the mechanical strains M of each lamina as:
M 0
x x + zkx x
M = 0 + zky y T (3.47)
yM 0y
xy k xy + zkxy xy k

Inserting this in the expression for the stresses in each lamina k, one obtains
0
x   x0 + zkx   x
y = Q k + zky Q k y T (3.48)
0y
xy k xy + zkxy xy k

Now, in order to obtain the relation between the total midplane strains and curvatures and the forces
produced by mechanical and thermal loads, we insert Eq. (3.48) in the expression for the stress resultants
(normal forces and moments due to mechanical loading) one obtains:

{NM } = [A]{0 } + [B]{k} {NT } (3.49)

M 0 T
{M } = [B]{ } + [D]{k} {M } (3.50)

or equivalently
NxT
0 M
x Nx
0 NM NT







0y yT

YM







A B xy Nxy Nxy

= M + (3.51)
B D kx0

M x

MxT

M
k0 My MyT







0y



M

T

kxy Mxy Mxy

where, as stated above, NxM , NyM , ..., Mxy

M
are normal forces and moments due to mechanical loads and
where the thermal loads are obtained as:

NxT Xn
  x
{NT } = NT = T Q k y (hk hk1 ) (3.52)
yT
Nxy xy k

k=1

MxT 1 Xn
  x
{MT } MyT h2k h2k1


= = T Q k y (3.53)
T 2
Mxy xy k

k=1

Remark. Please note that in Agarwal et al., the authors refer to something called thermal stresses
(denoted by a superscript T in e.g. Equation 6.62). This may be somewhat misleading since stresses are
always cased by mechanical straining ( = (M ) and nothing else. Stresses can not be anything but
mechanical. A more proper description would be thermally induced stresses.

38
Chapter 4

Kirchhoff-Love plate theory for

anisotropic laminated plates

4.1 Assumption in the Kirchhoff-Love plate theory

The Kirchhoff-Love plate theory is based the following assumptions (also discussed in Chapter 3 for the
classical laminate theory)
A line originally in the z-direction before deformation remains straight also after deformation.
The deformation of a line in the z-direction follow the rotation of the middle surface during defor-
mation, i.e. they remain perpendicular to the midsurface.
There is no change of thickness of the plate.
.

4.2 Repetition of kinematics for the Kirchhoff-Love plate theory

We recall from Chapter 3 that the strains under the Kirchhoff assumptions become:

u u0 2 w0
x = = z
x x x2
v v0 2 w0
y = = z
y y y 2
u v u0 v0 2 w0
xy = + = + 2z
y x y x x y
It is worth pointing out the two contradiction that
1. the shear strains are assumed to be zero xz = yz = 0 although the corresponding shear stresses
are not (cf. Section 4.3 below), e.g.
u w0 w0 w0
xz = + = + =0
z x x x

2. A simultaneous state of plane stress and constant thickness (after deformation) of the laminate is
assumed.

4.3 Equilibrium equations for a plate section

Consider an infinitesimally small element x-y of a rectangular plate shown in Figure 4.1. Stress
resultants (per unit length) acting on this small plate element are:

39
 out-of-plane shear forces Ryz and Rxz (note that these are resultants from the shear stresses yz
and xz respectively that previously have been neglected but now have to be considered to balance
the applied surface load)
in-plane normal and shear forces Nx , Ny and Nxy = Nyx and
moments Mx , My and Mxy = Myx .
In addition, the outer vertical force acting on the element is the distributed load p(x, y) (force per unit
area and positive in the positive z-direction).
z

y
y

Rxz
Nxy Nx Mx Mxy
Ryz + Ryz
Ryz

Ny Nyx + Nyx p Myx My + My

Nyx Ny + Ny My Myx + Myx

Nxy + Nxy Mx + Mx
Rxz+ Rxz
Nx + Nx
Mxy + Mxy

Figure 4.1: Definition of cross-sectional quantities in a rectangular plate.

As before, the in-plane forces and moments are resultants of the stresses x , y and xy as follows:
Z h/2 Z h/2 Z h/2
Nx = x dz , Ny = y dz , Nxy = xy dz = Nyx
h/2 h/2 h/2
Z h/2 Z h/2 Z h/2
Mx = z x dz , My = z y dz , Mxy = z xy dz = Myx
h/2 h/2 h/2

Correspondingly, the internal shear forces per unit length are the resultants of the out-of-plane shear
stresses Z h/2 Z h/2
Rxz = xz dz, Ryz = yz dz (4.1)
h/2 h/2
We now study equilibrium of the small plate element

4.3.1 Force equilibrium in x-direction

(Nx + Nx )y Nx y + (Nyx + Nyx )x Nyx x = 0
{Nxy = Nyx + divide by xy and let x 0,y 0}
(4.2)
Nx Nxy
+ =0
x y

4.3.2 Force equilibrium in y-direction

(Ny + Ny )x Ny x + (Nxy + Nxy )y Nxy y = 0
{divide by xy and let x 0,y 0}
(4.3)
Nxy Ny
+ =0
x y

40
4.3.3 Force squilibrium in z-direction
(Rxz + Rxz )y Rxz y + (Ryz + Ryz )x Ryz x + pxy = 0
{divide by xy and let x 0,y 0}
(4.4)
Rxz Ryz
+ = p
x y

4.3.4 Moment equilibrium around the x-axis

(My + My )x + My x (Mxy + Mxy )y + Mxy y + (Ryz + Ryz ) xy
y y y
+ (Rxz + Rxz ) y Rxz y + pxy =0
2 2 2
(4.5)
{divide by xy and let x 0,y 0}
Mxy My
+ = Ryz
x y

4.3.5 Moment equilibrium around the y-axis

(Mx + Mx )y Mx y + (Myx + Myx )x Myx x (Rxz + Rxz ) xy
x x x
(Ryz + Ryz ) y + Ryz y + pxy =0
2 2 2
(4.6)
{Mxy = Myx + divide by xy and let x 0,y 0}
Mx Mxy
+ = Rxz
x y

4.3.6 Governing equilibrium equations for the Kirchhoff-Love plate theory

Please note that the above five equations are valid for both plate theories to be considered. However,
since we in the Kirchhoff-Love plate theory only have three unknown fields (displacements of the midplane
u0 , v0 and w0 ), only three equations are needed. This can be accomplished by utilising Eqs. (4.2) and
(4.3) and inserting the results from Eqs. (4.5) and (4.6) into Eq. (4.4) leading to the three governing
equilibrium equations for a Kirchhoff-Love plate as:

Nx Nxy
+ = 0 (4.7)
x y
Nxy Ny
+ = 0 (4.8)
x y
2 2 2
Mx Mxy My
+2 + = p (4.9)
x2 xy y 2
For future purposes, we note that we also can write the equations above in a contracted way as:

TN = 0
(4.10)
T
M = p (4.11)
T
T and defined as:
with

2

0



x x2













2


x

=
0 , =
= with = (4.12)

y

y 2










2


y
2

y x

xy

41
Please note that it is the exact same equilibrium equations as for an isotropic plate. How-
ever, in that case there is no coupling between in-plane loading and bending due to the
material symmetry (or rather isotropy).

4.4 Equilibrium equations in terms of displacements for the

Kirchhoff-Love theory
The equilibrium equations can also be written in terms of displacements using the constitutive relations.
Inserting the following expressions for the normal forces and moments from Chapter 3:

{N} = [A]{0 } + [B]{k} (4.13)

0
{M} = [B]{ } + [D]{k} (4.14)

with 2
u0
w0

x x2










2 w0

v0
=
{0 } = u , {k} = w0 = (4.15)

y

y 2


u0 + v0 2 w0






2

y x

xy
into the equilibrium equations Eqs. (4.7)-(4.9), one obtains the equilibrium equations in terms of the
displacements as:

2 u0 2 u0 2 u0 2 v0 2 v0 2 v0
A11 + 2A 16 + A 66 + A 16 + (A 12 + A 66 ) + A 26
x2 xy y 2 x2 xy y 2
3 3 3 3
(4.16)
w0 w0 w0 w0
B11 3B16 2 (B12 + 2B66 ) B26 =0
x3 x y xy 2 y 3

2 u0 2 u0 2 u0 2 v0 2 v0 2 v0
A16 2
+ (A12 + A66 ) + A26 2
+ A66 2
+ 2A26 + A22 2
x xy y x xy y
3 3 3 3
(4.17)
w0 w0 w0 w0
B16 (B12 + 2B66 ) 2 3B26 B22 =0
x3 x y xy 2 y 3

4 w0 4 w0 4 w0 4 w0 4 w0
D11 + 4D 16 + 2(D12 + 2D 66 ) + 4D 26 + D22
x4 x3 y x2 y 2 xy 3 y 4
3 3 3 3
u0 u0 u0 u0 3 v0
B11 3B 16 (B 12 + 2B 66 ) B 26 B 16 (4.18)
x3 x2 y xy 2 y 3 x3
3 3 3
v0 v0 v0
(B12 + 2B66 ) 2 3B26 B22 3 = p
x y xy 2 y
which are the general equilibrium equations for a laminated plate with any lay-up sequence.

4.4.1 Symmetric and specially orthotropic laminates

In the case of a symmetric laminate, i.e. Bij = 0 the equilibrium equations are simplified to be:

2 u0 2 u0 2 u0 2 v0 2 v0 2 v0
A11 + 2A 16 + A 66 + A 16 + (A 12 + A 66 ) + A 26 =0 (4.19)
x2 xy y 2 x2 xy y 2
2 u0 2 u0 2 u0 2 v0 2 v0 2 v0
A16 + (A 12 + A 66 ) + A 26 + A 66 + 2A 26 + A 22 =0 (4.20)
x2 xy y 2 x2 xy y 2
4 w0 4 w0 4 w0 4 w0 4 w0
D11 + 4D 16 + 2(D 12 + 2D 66 ) + 4D 26 + D22 =p (4.21)
x4 x3 y x2 y 2 xy 3 y 4

42
 Furthermore, specially orthotropic laminates (such as e.g. unidirectional laminates with all fibres
oriented with an angle of 0 or 90 or symmetric cross-ply laminates) the laminate that behave like a
single layer of an orthotropic material satisfy the following additional conditions:

A16 = A26 = 0 (4.22)

D16 = D26 = 0 (4.23)
whereby the equilibrium equations are further simplified to be:

2 u0 2 u0 2 v0
A11 2
+ A66 2
+ (A12 + A66 ) = 0 (4.24)
x y xy
2 2
u0 v0 2 v0
(A12 + A66 ) + A66 2
+ A22 2 = 0 (4.25)
xy x y
4 4
w0 w0 4 w0
D11 4
+ 2(D12 + 2D66 ) 2 2 + D22 = p (4.26)
x x y y 4
It should be remarked that since Bij = 0 the in plane and transverse motions can be treated separately
in this case. In the following, we will study a method for analytically solving the bending part for a simply
supported specially orthotropic plate.

4.4.2 Analytical solutions - Simply supported specially orthotropic plates

Analytical approaches to solve the plate equation (bending part, i.e. Eq (4.26)) for simple boundary
conditions exist. One example is the Naviers solutions for simply supported rectangular plates
dy
y

dx

In this method, the deflection w0 is proposed to be of the following form (Fourier series expansion):
X
X  m x   n y 
w0 (x, y) = wmn sin sin (4.27)
m=1 n=1
a b

where wmn are unknown (Fourier) coefficients. The boundary conditions for a fully simply supported
plate are (
w0 (0, y) = w0 (a, y) = 0 , w0 (x, 0) = w0 (x, b) = 0
(4.28)
Mx (0, y) = Mx (a, y) = 0 , My (x, 0) = My (x, b) = 0
The bending moments can, via (4.14) be expressed as
2 w0 2 w0
Mx = D11 D12 (4.29)
x2 y 2
2 w0 2 w0
My = D12 D22 (4.30)
x2 y 2
and with the assumed expression for w0 :
X X   m 2  n 2   m x   n y 
Mx = wmn D11 D12 sin sin (4.31)
m=1 n=1
a b a b

43
 X
X   m 2  n 2   m x   n y 
My = wmn D12 D22 sin sin (4.32)
m=1 n=1
a b a b

Clearly, the adopted expression for w0 in (4.27) fulfills the boundary conditions since

sin (n ) = 0 , for n = 0, 1, 2, . . .

The loading p must also be expressed in a similar fashion as w0

X
X  m x   n y 
p(x, y) = pmn sin sin
m=1 n=1
a b

where the coefficients pmn can be found from the following expression
Z b Z a
4  m x   n y 
pmn = p(x, y) sin sin dx dy
ab 0 0 a b

To show this we insert p(x, y) into the integral

Z bZ a X X    
i x j y sin m x sin n y dx dy
   
pij sin sin
0 0 i=1 j=1
a b a b

and use the orthogonality of sin functions, i.e.

(
Z a  
i x  m x  a/2 if i = m
sin sin dx =
0 a a 0 else

which proves the expression for pmn .

The next step is now to insert the expressions for w0 and p in the plate equation (4.26) and solve for
the unknown coefficients wmn . The result is
pmn
wmn =   m 4  m 2  n 2  n 4  (4.33)
4
D11 + 2(D12 + 2D66 ) + D22
a a b b

Some examples of loading situations:

Uniform loading p(x, y) = p0 :

4 p0 b a
Z Z  m x   n y 
pmn = sin sin dx dy =
ab 0 0 a b
4 p0 a b h  m x ia h  n y ib
= cos cos
ab m n a 0 b 0
16 p0
= , m, n = 1, 3, . . .
mn 2
Please observe that, despite the fact that the sum in the Fourier series goes to infinity, it is often
for practical reasons sufficient to truncate these series (i.e. use a finite value for the upper limit of
the sum). It is generally the first modes (lowest values of m and n that have the most significant
contributions). As an example, please consider Figure 4.2 where the uniform pressure load is
approximated with truncated Fourier series with different maximum values.

Hydrostatic loading p(x, y) = p0 y/b

4 p0 b a
Z Z  m x   n y 
pmn = sin y sin dx dy =
ab2 0 0 a b
8 p0
= (1)n+1 , m = 1, 3, 5, . . . , n = 1, 2, 3, . . .
2 m n

44
 0.2 0.2 0.2
p(x,y)

p(x,y)

p(x,y)
0.1 0.1 0.1

0 0 0
1 1 1

0.5 0.5 0.5

1 0.8 1 0.8 1
0.6 0.8 0.6 0.6
0.4 0.2 0.4 0.2 0.4
y 0 0.2 y 0 0 y 0 0
0
x x x

a) b) c)

Figure 4.2: Approximation of uniform pressure load p(x, y) = 0.2 using truncated Fourier series with
maximum values according to a) m = n = 3, b) m = n = 7 and c) m = n = 11

Point load Q0 at x0 , y0 . The distributed force can be described by impulse function (Dirac-delta)
p = Q0 (x x0 , y y0 ). With the property
Z bZ a
(x x0 , y y0 ) f (x, y) dx dy = f (x0 , y0 )
0 0

Thereby, pmn become

4 Q0 b a
Z Z  m x   n y 
pmn = (x x0 , y y0 ) sin sin dx dy =
ab 0 0 a b
4 Q0  m x 
0
 n y 
0
= sin sin , m, n = 1, 2, 3, . . .
ab a b

4.5 Buckling of a Kirchhoff-Love plate

a) b)

c) d)

Figure 4.3: Views of the small plate element xy showing the contributions from in-plane forces Nx ,
Ny , Nxy and Nyx to the vertical equilibrium equation. a) Contribution from Nx , b) contribution from
Ny and c)-d) contribution from in Nyx (Nxy can be treated in the same way)

Moderate in-plane loads on a flat symmetric laminate cause in-plane displacements but no out-of-
plane displacements. However, it is known that in-plane compressive loads, when high enough, cause

45
out-of-plane displacements (deflections) that may be excessive and lead to failure. This is called elastic
instability or buckling. In buckling, out-of-plane displacements are caused by in-plane loads and, hence,
the classical Kirchhoff-Love equilibrium equations cannot be used to predict this behaviour because
interaction between in-plane loads and out-of-plane displacements were suppressed in the derivation (it
was assumed that the out-of-plane displacements are so small that the in-plane force resultants Nx , Ny
and Nxy act in their original direction in the xy-plane. Therefore, to study buckling, we need to modify
the governing equations to take into account the out-of-plane displacements and thereby the vertical
projection of the in-plane forces. To do so, we study the in-plane forces acting on the edges of the same
small plate element in Figure 4.1, now accounting for the effects of the out-of-plane deformations as
shown in Figure 4.3.

4.5.1 The vertical projection of Nx

Relate to Figure 4.3a:
w Nx w 2 w 2w Nx w Nx 2 w
Nx y + (Nx + x)( + x)y = N x xy + xy + x2 y
x x x x2 x2 x x x x2
2
w Nx w
Nx 2 xy + xy
x x x
(4.34)
where it was used that
x2 y  xy.

4.5.2 The vertical projection of Ny

Relate to Figure 4.3b:
w Ny w 2 w 2w Ny w Ny 2 w
Ny x + (Ny + y)( + 2
y)x = Ny 2 xy + xy + xy 2
y y y y y y y y y 2
2w Ny w
Ny 2 xy + xy
y y y
(4.35)
where it was used that
xy 2  xy.

4.5.3 The vertical projection of Nyx

Relate to Figure 4.3c)-d):
w Nyx w 2w 2w Nyx w Nyx 2 w
Nyx x + (Nyx + y)( + y)x = Nyx xy + xy + xy 2
x y x xy xy y x y xy
2w Nyx w
Nyx xy + xy
xy y x
(4.36)
where it was used that
xy 2  xy.

4.5.4 The vertical projection of Nxy

In the same way as for Nyx we can obtain the contribution from Nxy as:
w Nxy w 2w 2w Nxy w Nxy 2 w
Nxy y + (Nxy + x)( + x)y = Nxy xy + xy + x2 y
y x y xy xy y y x xy
2w Nxy w
Nxy xy + xy
xy y y
(4.37)

46
where it was used that
x2 y  xy.

4.5.5 Total equation

The total vertical contribution of the in-plane stress resultants becomes (via summation of Eqs. (4.34)-
(4.37)):

2w Nx w 2w Nyx w 2w
Nx xy + xy + N yx xy + xy + N xy xy
x2 x x xy y x xy
Nxy w 2w Ny w
+ xy + Ny 2 xy + xy = {Nyx = Nxy } =
y y y y y
2w 2w 2w (4.38)
=Nx 2 xy + 2Nxy xy + Ny 2 xy
x xy y
   
Nx Nyx w Nxy Ny w
+ + xy + + xy
x y x y y y
| {z } | {z }
=0(in-plane equilib.) =0(in-plane equilib.)

Adding the vertical contributions from in-plane force resultants to Eq. (4.4) yields (when x 0,
y 0):
Rxz Ryz 2w 2w 2w
+ + Nx 2 + 2Nxy + Ny 2 + p = 0 (4.39)
x y x xy y

2 Mx 2 Mxy 2 My 2w 2w 2w
+ 2 + + Nx + 2Nxy + Ny +p=0 (4.40)
x2 xy y 2 x2 xy y 2

4.5.6 Buckling equation for a symmetric laminate

Using that Bij = 0 for a symmetric laminate, the buckling equation in terms of displacements can be
obtained as:
4 w0 4 w0 4 w0 4 w0 4 w0
D11 4
4D16 3 2(D12 + 2D66 ) 2 2 4D26 3
D22
x x y x y xy y 4
(4.41)
2w 2w 2w
+ Nx 2 + 2Nxy + Ny 2 + p = 0
x xy y

4.5.7 Buckling equation for a specially orthotropic laminate

As seen above, specially orthotropic laminates have apart from Bij = 0 also the property that D16 =
D26 = 0 which further reduces the buckling equation to:

4 w0 4 w0 4 w0
D11 4
2(D12 + 2D66 ) 2 2 D22
x x y y 4
(4.42)
2w 2w 2w
+ Nx 2 + 2Nxy + Ny 2 + p = 0
x xy y

By assuming the applied compressive loading according to Nx = Nx0 , Ny = Ny0 , p = 0 and

Nxy = 0 one obtains:

4 w0 4 w0 4 w0 2w 2w
D11 + 2(D12 + 2D66 ) + D 22 = N x0 N y0 (4.43)
x4 x2 y 2 y 4 x2 y 2

For the special case of simply supported specially orthotropic laminates we can again use Naviers
solution. Hence, the deflection w0 is once again given by:

47
 X
X  m x   n y 
w0 (x, y) = wmn sin sin
m=1 n=1
a b
whereby the buckling equation can be rewritten as:
X
X   m 4  mn 2  n 4   m x   n y 
4 wmn D11 + 2(D12 + 2D66 ) + D22 sin sin
m=1 n=1
a ab b a b
  m 2  n 2   m x   n y  (4.44)
XX
= 2 wmn Nx0 + Ny0 sin sin
m=1 n=1
a b a b

Finally, for non-trivial solutions wmn 6= 0, we must have:

  m 4  mn 2  n 4    m 2  n 2 
2
D11 + 2(D12 + 2D66 ) + D22 = Nx0 + Ny0 (4.45)
a ab b a b

4.5.8 Example: Buckling under uniaxial compression

Lets consider the special case of a rectangular plate with side lengths a and b under uniaxial compression
Nx0 = N0 , Ny0 = 0. We then obtain the following expression for the buckling load:
  m 2  n 2  n 4  a 2 
N0 (m, n) = 2 D11 + 2(D12 + 2D66 ) + D22 (4.46)
a b b m
The critical buckling load is then the lowest value of N0 for any combination of m and n which
generally depends on the ratio a/b and the bending stiffness Dij . Please note that m and n characterise
the shape of the buckling modes in terms of half wavelengths in the x and y directions respectively. Thus,
e.g. for a full sine-wave shape in x-direction we have m = 2.
In this particular case, the buckling load is obtained for n = 1. Furthermore, for a laminate with the
properties D11 /D22 = 10 and (D12 + 2D66 )/D22 = 1 this is reduced to:
"  2 #
2 D22 b 1  a 2
N0 (m, 1) = 10m2 +2+ 2 (4.47)
b2 a m b
or written on a non-dimensional form:
"  2 #
N0 (m, 1)b2 b 1  a 2
N= = 10m2 +2+ 2 (4.48)
2 D22 a m b
cf. Figure 4.4. It can be seen that, depending on the geometry and lay-up sequence of the plate, the
most critical load can be obtained also for m > 1. In addition, the corresponding buckling modes for the
first three buckling modes are shown in Figure 4.5.

4.6 Free vibrations of a Kirchhoff-Love plate

4.6.1 Study Subsection 7.3.3 in the course book on your own

48
Figure 4.4: Non-dimensional buckling load N = (N0 [m, 1]b2 )/( 2 D22 ) for a rectangular laminate with
properties D11 /D22 = 10 and (D12 + 2D66 )/D22 = 1 subjected to uniaxial compression.

0.1
w

0
1

0.5
1
0.8
0.6
0.4
y 0 0.2
0
x

a)

0.1

0
w

0.1

0.5
1
0.8
0.6
0.4
y 0 0.2
0
x

b)

0.1
w

0
1

0.5
1
0.8
0.6
0.4
y 0 0.2
0
x

c)

Figure 4.5: Buckling modes for a quadratic laminate (a = b) subjected to uniaxial compression. The first
three modes obtained for n = 1 and a) m = 1, b) m = 2 and c) m = 3.

49
50
Chapter 5

Mindlin - Reissner plate theory for

anisotropic laminated plates

5.1 Extensions of the plate theory to account for out-of-plane

shear: The Mindlin - Reissner theory
Experimental evidence exist in the literature that supports the Kirchhoff-Love plate theory, especially
for thin plates with a high length-to-thickness ratio. Also numerical results where the plate is modelled
as a 3D continuum (higher level of accuracy) agree well for thin plates. However, as the plate thickness
increase in relation to other dimensions, experimental and 3D simulations generally show larger out-
of-plane deflection than what is obtained by using Kirchhoffs plate theory. The discrepancy is largely
dependent on the plate deformations induced by out-of-plane shear (which is neglected in Kirchhoff-
Love plate theory). Thus, for moderately thick plates, the theory needs to take these shear strains into
consideration.

5.2 First-order Mindlin-Reissner shear deformation theory

In order to account for the out-of-plane shear strains, the Kirchhoff-Love theory needs to be modified.
This is done by relieving some of the restrictions imposed by the Kirchhoff-Love assumptions. Thus, we
still assume that lines initially in the z-direction will remain straight after deformation but they do no
longer have to remain perpendicular to the midplane (due to the shear strains).

5.2.1 Kinematics of the Mindlin-Reissner plate theory

To see how the assumptions influence the kinematics of the Mindlin-Reissner theory, we consider first the
deformed plane in Figure 5.1. From this, we can conclude that the displacement in the direction of the
x-axis of a point on the planar cross section can be written as:

u(x, y, z) = u0 (x, y) + zx (x, y) (5.1)

where x is the (positive) rotation around the y-axis 6= w/x. In the same way (if we would consider
deformation in the yz-plane) one obtains the displacement in the y-direction as

v(x, y, z) = v0 (x, y) + zy (x, y) (5.2)

with y being the (negative!) rotation around the x-axis 6= w/y. Finally, the vertical displacement (in
the z-direction) is as before given by
w(x, y, z) = w0 (x, y). (5.3)

By definition, we obtain the strains as:

51
 y x y
y
y y

v
Figure 5.1: Kinematics of Mindlin-Reissner plate deformed in the a) xz plane and b) yz plane (from J.N.
Reddy, Theory and Analysis of Elastic Plates and Shells (2nd ed.), CRC Press, 2007).

u u0 x
x = = +z (5.4)
x x x
v v0 y
y = = +z (5.5)
y y y
 
u v u0 v0 x y
xy = + = + +z + (5.6)
y x y x y x
u w w0
xz = + = x + (5.7)
z x x
v w w0
yz = + = y + (5.8)
z y y
or

{} = {0 } + z{k}M indlin (5.9)

w0
x +


x

{ z } = (5.10)
w0
y +


y
with

u0

x

x x











{0 } = v0 , {k}M indlin = y (5.11)

y

y

u0 v0 x y





+ +






y x y x
As can be observed, the definition of the curvature is extended from the Kirchoff theory to the
Mindlin-Reissner theory. However, it is clear that if we would neglect the shear strains, we would obtain
the relation between the rotation angle and the derivative of the out-of-plane displacement w0 as in the
Kirchhoff-Love theory yielding the same curvature.

5.2.2 Transverse shear forces and stresses

Please note that the shear strains in the Mindlin-Reissner theory are assumed constant through the
thickness which implies that constant shear stresses are predicted within each lamina k through
     
yz Q44 Q45 yz
= (5.12)
xz Q45 Q55 k xz

52
or equivalently
     
xz Q55 Q45 xz
= . (5.13)
yz Q45 Q44 k
yz

Please note that in the last step, the order of the of the shear stresses and strains were reversed for
reasons that will be more clear below (cf. the contracted notation of the equilibium equations, the weak
formulation and in later sections also the FE formulation). Please note that this also means that elements
of the stiffness matrix needs to change place, cf. Eq (5.13).
In order to obtain Q44 , Q45 and Q55 we need to consider the coordinate transformation of the out-of-
plane shear stresses and strains due to a fibre orientation + around the positive z-axis. For this case,
the out-of-plane axis for the global coordinate system (the z-axis) and the out-of-plane axis for the fibre
0
oriented local coordinate system (the so-called T -axis) coincide, whereby the relation between transverse
stresses may be derived from either a standard change of basis or from equilibrium of forces in the vertical
direction as:
    
T T 0 cos sin yz
= . (5.14)
LT 0 sin cos xz
| {z }
[T1 ]

As a consequece, we can obtain Q44 , Q45 and Q55 , knowing that the shear strains z transform in the
same way
    
T T 0 cos sin yz
= , (5.15)
LT 0 sin cos xz
| {z }
[T1 ]

as
   
Q44 Q45 Q44 Q45
= [T1 ]1
k [T1 ]k (5.16)
Q45 Q55 k
Q45 Q55

where     
T T 0 Q44 Q45 T T 0
= (5.17)
LT 0 Q45 Q55 LT 0

with
     
Q44 Q45 GT T 0 0 GT T 0 0 0
= = , GLT 0 = GLT (the T T plane is a plane of isotropy)
Q45 Q55 0 GLT 0 0 GLT 0
k
(5.18)
and where we need either experimental data or a model to predict GT T 0 .
Micromechanical models to predict GT T 0 generally are rather complicated and thereby out of scope
in this course1 . However, interested readers are referred to e.g. the book by Christensen (Mechanics of
Composite Materials, Krieger Publishing company, 1979). It should be known that GT T 0 depends on
fibre volume fraction and the shear properties of the two constituents (fibre and matrix material) and is
generally somewhat lower than the in-plane shear modulus GLT , cf. also Figure 5.2 in which predictions
of GT T 0 (denoted G23 in the figure) from a model developed at Swerea SICOMP denoted CCA (equivalent
to the Christensen model for a pure fibre-matrix material) is compared to the derived upper and lower
bounds by Hashin and the prediction of the in-plane shear modulus GLT using the Halpin-Tsai equations
0
with = 1. Since the T T plane is a plane of isotropy, we know for a fact that:

ET
GT T 0 = (5.19)
2 (1 + T T 0 ) .

As a rule of thumb for carbon fibre/epoxy composites, T T 0 normally lies in the interval 0.4-0.45.

1 As an alternative, FEM can also be used to predict the transverse shear modulus based on an analysis of a representative

volume element

53
Figure 5.2: Predicted values (obtained by Erik Marklund, Swerea SICOMP) for the transverse shear
modulus GT T 0 = G23 of a carbon fibre/epoxy composite based on the self-consistent model developed
at Swerea SICOMP (for this case equivalent to the model by Christensen, Mechanics of Composite
Materials, Krieger Publishing company, 1979) compared to the upper and lower bounds derived by
Hashin and predictions of the in-plane shear modulus GLT = G12 by the Halpin-Tsai equations with
= 1. Parameters used in the predictions are: Em = 3 GPa, m = 0.38, Ef L = 230 GPa, Ef T = 23
GPa, Gf LT = 20 GPa and f LT = f T T 0 = 0.2

5.2.2.1 The shear forces and the shear correction factor

Please note that the representation of shear stress through the thickness in unphysical in the sense that
traction continuity no longer is preserved. This is generally accounted for in the laminate analysis by
introducing a shear correction factor K such that the resulting out-of-plane shear forces are obtained as:
  Z h/2    
Rxz Q55 Q45 xz
=K dz . (5.20)
Ryz h/2 Q45 Q44 yz
| {z }


A55 A45
A45 A44

Please note the slight modification compared to the course book in which the shear correction factor is
included in the Aij components!
By considering a homogeneous specially orthotropic plate and enforcing the condition that the work
done by the external forces Rxz and Ryz

2 2
Ryz
Z Z
1 1 Rxz
Rxz xz + Ryz yz dA = ... = + dA (5.21)
2 Aplate 2K Aplate A55 A44

equals the internal strain energy produced by the shear stresses and strains
h/2 2 2
Ryz
Z Z Z
1 3 Rxz
xz xz + yz yz dA dz = ... = + dA (5.22)
2 h/2 Aplate 5 Aplate A55 A44

the shear correction factor gets the value K = 5/6. It should be noted that even thought the value of 5/6
for K is derived only for a homogeneous orthotropic plate it generally provides good approximations also
for other types of laminates (except for sandwich structures with a thick middle layer for which K 1).

54
5.3 Governing equations for the Mindlin-Reissner theory
The governing plate equations (forces and moments) for the Mindlin-Reissner theory are the same as for
Kirchhoff-Love plate theory. But, since there are now five unknowns (u0 , v0 , w0 , x and y ) the vertical
equilibrium equation cannot be combined with the two moment equations. Thus, the governing equations
for the Mindlin-Reissner (or first-order shear deformation) theory are:

Nx Nxy
+ =0 (5.23)
x y
Nxy Ny
+ =0 (5.24)
x y
Rxz Ryz
+ = p (5.25)
x y
Mx Mxy
+ = Rxz (5.26)
x y
Mxy My
+ = Ryz (5.27)
x y

or on a contracted form:

TN
= 0 (5.28)
T
Rz = p (5.29)
T
M = Rz
(5.30)

with



0
x




Nx Mx

  x



Rxz


N= N y , Rz =
Ryz
, M= My = , = 0 y
(5.31)
Nxy Mxy


y






y x

5.4 Equilibrium equations in terms of displacements and rota-

tions for the Mindlin-Reissner theory
Given the definition of strains and stress resultants, we obtain their inter-relationship as:

+ [B]
[A] 0 + [B] {k} = [A]u

N = (5.32)
{ z } = K[A]
Rz = K[A] ( + w0 ) (5.33)
+ [D]

 0
M = [B] + [D] {k} = [B]u (5.34)

with N, M, A, B and D as before and

       
A = A55 A45 , Rz = Rxz
, u=
u0
, =
x
(5.35)
A45 A44 Ryz v0 y

By insertion of Eqs. (5.32)-(5.35) into the governing equations, Eqs. (5.23)-(5.27), one obtains the equi-
librium equations in terms of the displacement and rotation fields. However, this yields rather lengthy
expressions that can be solved only for simple anisotropic conditions and boundary conditions (in analogy
with Naviers solution for the Kirchhoff-Love plate). Instead, to allow for the solution of more general
problems by FEM, we derive the weak (or variational) formulation of Eqs. (5.23)-(5.27).

55
5.5 Weak form of equilibrium for the Mindlin-Reissner theory
It is first noted that in order to obtain the weak form of the governing equations, each equation is
multiplied by an arbitrary test functions (u = {u v}T , w and = {x y }T ) (often also
denoted virtual displacements) and integrated over the domain . We treat the in-plane equations, the
vertical equation and the moment equations separately to obtain:

Z  T  Z Z  T
N d
uT = uT P d
u N d = 0 (5.36)

Z   Z Z Z
T
w T Rz d = wRn d
(w) Rz d = wp d (5.37)

Z  T  Z
Z  T Z

M d
T = T Mn d
M d = T Rz d (5.38)

T
where P = {Px Py } is the external in-plane force (per unit length) acting on , Rn is the external
T
out-of-plane force (per unit length) and Mn = {Mnx Mny } is the external moment (per unit length).
Please note that P can be expressed in terms of the in-plane normal n and tangent m vectors as
Z h/2 Z h/2
P = Pn n + Pnm m, Pn = n dz, Pnm = nm dz (5.39)
h/2 h/2

where n is the in-plane normal stress acting on the surface and nm is the corresponding in-plane
shear stress (acting on ). In the same way, we have
Z h/2
Rn = nz dz (5.40)
h/2
Z h/2 Z h/2
Mn = Mnn n + Mnm m, Mnn = nn z dz, Mnm = nm z dz (5.41)
h/2 h/2

where nz is the out-of-plane shear stress acting on in the z-direction.

Please note that the weak formulation of the problem does not depend on if the plate is isotropic (also
derived in the course Finite Element Method - Structures) or anisotropic since only stress resultants and
displacements are involved. Please also note that there are a number of steps to be taken in order to
arrive at the weak form of the governing equations (i.e. to derive Eqs. (5.36)-(5.38)) but that the steps
are very similar for the three equations whereby we focus only on the first equation below.

5.5.1 Derivation of Eq. (5.36)

To derive Eq. (5.36), we follow the steps for 2D elasticity derived e.g. in the Chalmers courses MHA021 -
Finit elementmetod (FEM) and VSM167 - Finite element method - basics, both given by the Department
of Applied Mechanics. The procedure is also very nicely described (for 3D elasticity) in the book by
Ottosen and Petersson (Introduction to the Finite Element Method, Prentice Hall, 1992.).
As a first step, Eqs. (5.23) and (5.24) are multiplied by arbitrary test functions u and v respectively
which yields the following two equations:

Z   Z   Z  
Nx Nxy u u
u + d = (uNx ) + (uNxy ) d Nx + Nxy d = 0
x y x y x y
Z   Z   Z  
Nxy Ny v v
v + d = (vNxy ) + (vNy ) d Nxy + Ny d = 0
x y x y x y

Please also note that this corresponds to the multiplication of Eq. (5.28) with u = {u 0}T and
u = {0 v}T respectively.
Utilising that

56
 Z Z

d = nx d (5.42)
x
Z Z

d = ny d (5.43)
y

where n = {nx ny }T is the in-plane normal of the boundary of the expressions may be rewritten
as
Z Z  
u u
u (Nx nx + Nxy ny ) d = Nx + Nxy d (5.44)
x y
Z Z  
v v
v (Nxy nx + Ny ny ) d = Nxy + Ny d (5.45)
x y

Taking the sum of Eqs. (5.44) and (5.45) yields

Z Z  
u u v v
u (Nx nx + Nxy ny ) +v (Nxy nx + Ny ny ) d = Nx + Nxy + Nxy + Ny d

| {z } | {z } x y x y
Px Py

or equivalently Z Z  T
T
u P d =
u N d (5.46)

By taking the sum of Eqs. (5.44) and (5.45), it appears as if we reduced the number of in-plane
equations from two to one. However, as indicated in the beginning of this subsection, u is arbitrary and
by using both u = {u 0}T and u = {0 v}T we retain the two governing equations we started with.
In fact, we can use any two arbitrary combinations of u as long as they satisfy the condition uT1 u2 = 0
(vector multiplication).
In order to see
Px = Nx nx + Nxy ny , Py = Nxy nx + Ny ny
we utilise that P is the projection of the stress in the normal direction of the surface integrated over the
thickness, i.e. in 2D we have

  Z h/2    Z h/2    
Px x xy nx x nx + xy ny Nx nx + Nxy ny
P= = dz = dz =
Py h/2 xy y ny h/2 xy nx + y ny Nxy nx + Ny ny

5.5.2 Boundary conditions

By examining the surface integral (boundary terms) of Eqs. (5.36)-(5.38), we conclude that there are
three pairs of boundary conditions:

5.5.2.1 Boundary conditions for the in-plane force equilibrium

The boundary term is Z
uT P d (5.47)

whereby we can conclude that we can subdivide the boundary into two parts:
D,u on which the in-plane displacements are prescribed: (u = {u0 v0 }T = {u0 v 0 }T = u)
where u0 and v 0 are the prescribed displacement components in x- and y-direction respectively.
Furthermore, for this part of the boundary, the corresponding test functions are set to zero.
N,u on which the external in-plane forces are prescribed: (P = {Px Py }T = {P x P y }T = P)
where P x and P y are the prescribed force components (per unit length) in the x- and y-direction
respectively, acting on N,u .

57
Hence, the boundary term reduces to
Z Z
uT P d = uT P d (5.48)
N,u

5.5.2.2 Boundary conditions for the out-of-plane force equilibrium

The boundary term is

Z
wRn d (5.49)

whereby we can conclude that we can subdivide the boundary into two parts:

D,w on which the out-of-plane displacements are prescribed: (w = w0 = w0 ) where w0 is the

prescribed displacement component in z-direction. In analogy to the in-plane force balance equation,
for this part of the boundary, the corresponding test functions is set to zero.

N,w on which the external out-of-plane force is prescribed: (Rn = Rn ) where Rn is the prescribed
vertical force (per unit length) in the z-direction acting on N,w .

Hence, the boundary term reduces to

Z Z
wRn d = wRn d (5.50)
N,w

5.5.2.3 Boundary conditions for the moment equilibrium

The boundary term is

Z
T Mn d (5.51)

whereby we can conclude that we can subdivide the boundary into two parts:

D, on which the rotations are prescribed: ( = {x y }T = {x y }T = ) where x and y

are the prescribed rotation components. As before, for this part of the boundary, the test functions
are set to zero.

N, on which the external moments are prescribed: (Mn = {Mnx Mny }T = {M nx M ny }T = Mn )

where M nx and M ny are the prescribed moment components (per unit length) acting on N, .

Hence, the boundary term reduces to

Z Z
T Mn d = T Mn d (5.52)
N,

58
5.5.3 Final weak form
Rewriting the boundary terms and using the kinematical relations one obtains the weak form as:

Find u = {u0 v0 }T Vu , w Vw and = {x y }T V where

Vu = {sufficiently regular; u = u on D,u }

Vw = {sufficiently regular; w = w on D,w }
V = {sufficiently regular; = on D, }

such that:

Z  T   Z

u + [B]
[A]u d = uT P d, u Vu
N,u
Z Z Z
T
( + w) d =
K (w) [A] wp d + wRn d, w Vw (5.53)
N,w
Z  T   Z Z


[B]u + [D] d + T
K[A] ( + w) d = T Mn d, V
N,
u = u on D,u
w = w on D,w
= on D,

where

Vu = {sufficiently regular; u = 0 on D,u }

Vw = {sufficiently regular; w = 0 on D,w }
V = {sufficiently regular; = 0 on D, }

which is the starting point for the FE formulation of the problem.

59
60
Chapter 6

FEM formulation of laminated

composite plates

6.1 Preliminaries
The point of departure is the weak (or variational formulation) of the governing equations for the Mindlin
plate as:

Find u = {u0 v0 }T Vu , w Vw and = {x y }T V where

Vu = {sufficiently regular; u = u on D,u }

Vw = {sufficiently regular; w = w on D,w }
V = {sufficiently regular; = on D, }

such that:

Z  T   Z

u
[A]u + [B] d = uT P d, u Vu
N,u
Z Z Z
T
( + w) d =
K (w) [A] wp d + wRn d, w Vw
N,w
Z  T   Z Z

+ [D]
[B]u d + ( + w) d =
T K[A] T Mn d, V
N,
u = u on D,u
w = w on D,w
= on D,

where

Vu = {sufficiently regular; u = 0 on D,u }

Vw = {sufficiently regular; w = 0 on D,w }
V = {sufficiently regular; = 0 on D, }

and where we recall that

0
x
x




= , 0
= (6.1)

y

y
y x

61
6.2 FE-approximations of fields
So far, we have considered only planar problems such that everything is just a function of x and y since
the xy-plane is placed in the plane of the plate. Proceeding in this way, are now in the position of
introducing the FE approximation of the different fields in terms of unknown node variables so-called
degrees-of-freedom (e.g. w
i for the out-of-plane displacement) and shape functions Ni (x, y) as:

  nno  
u0 (x, y) X ui
u(x, y) = = Ni (x, y) i =
ui , where u (6.2)
v0 (x, y) vi
i=1
nno
X
w(x, y) = w0 (x, y) = Ni (x, y)w
i (6.3)
i=1
  nno  
x (x, y) X
i , where
i = x,i
(x, y) = = Ni (x, y) (6.4)
y (x, y) y,i
i=1

or equivalently

u = N
u (6.5)
w = w
N (6.6)
= N (6.7)

with
 
N1 0 N2 0 ... Nnno 0 =
 
N= , N N1 N2 ... Nnno
0 N1 0 N2 ... 0 Nnno
 T  T (6.8)
=
u u
1 v1 u
2 v2 ... u
nno vnno , w
= w
1 w
2 ... w
nno
T
= x,1 y,1 x,2 y,2 ... x,nno y,nno


Pnno
where we dropped (x, y) for brevity. Please note that the sum i=1 is defined over all nodes (nno) in
the domain and that the shape functions have the properties
 nno
1 in node i X
Ni = , Ni = 1 in any point of the domain (6.9)
0 in all other nodes j 6= i
i=1

Furthermore, from the weak form, we know that we also need u, which we, via the
w and
approximation above, obtain as:

u (N
= u) = B
u, B = N (6.10)
w w, = N


w = N =B B (6.11)
 

N
= = B, B = N (6.12)

since all degrees-of-freedom are constant nodal values and not functions of x or y.

6.3 FE-formulation of laminated plates

On order to end up with the FE formulation of the laminated Mindlin plate, we start of by considering
each equation separately.

6.3.1 FE formulation of in-plane equilibrium

Before we insert the FE approximation of the unknown displacement field into the weak form of the
in-plane equilibrium
Z  T   Z

u + [B]
[A]u d = uT P d, u Vu (6.13)
N,u

62
we first note this should hold for all arbitrary test functions u Vu . By using the so-called Galerkins
method, and choosing u1 = {u 0}T and u2 = {0 v}T (where u and v can be arbitrary, given that
they satisfy the necessary restrictions in terms of e.g. boundary conditions) we will obtain two coupled
equations for the in-plane displacements. Furthermore, according to Galerkins method we use the shape
functions Ni to represent also the arbitrary test functions
 
u  T
u1 = = Nc1 , c1 = c1,1 0 c1,2 0 ... c1,nno 0 (6.14)
0
 
0  T
u2 = = Nc2 , c2 = 0 c2,1 0 c2,2 ... 0 c2,nno (6.15)
v

whereby
1
u (Nc1 ) = Bc1 , B = N
= (6.16)
2
u (Nc2 ) = Bc2 , B = N
= (6.17)

If we start by choosing u = u1 (with arbitrary non-zero coefficients in c1 ) and inserting the FE

approximations of u, and u1 we obtain:
Z Z Z !
cT1 BT [A]B d
u+ BT [B]B d NT P d = 0. (6.18)
N,u

In the same way, if we choose u = u2 (with arbitrary non-zero coefficients in c2 ) we obtain:

Z Z Z !
c T T
B [A]B d u+ T
B [B]B d T
N P d = 0 (6.19)
2
N,u

which if combined can be written as:

!
cT1
Z Z Z  
c T T
B [A]B d
u+ T
B [B]B d T
N P d = 0 with cT = (6.20)
N,u cT2

Now, since we know that this should hold for arbitrary test functions, i.e. for any choice of c we must
have Z Z Z
T
B [A]B d u + T
B [B]B d NT P d = 0 (6.21)
N,u
| {z } | {z } | {z }
Kuu Ku f b,u
or
Kuu u = f b,u
+ Ku (6.22)

6.3.2 FE formulation of out-of-plane equilibrium

w = Bc
Following the procedure above, now simply by using w = Nc, we obtain the
and = N
FE-formulation of the out-of-plane equation as:
Z Z Z Z !
c T T
K B [A]N d + T
K B [A]B dw T
N p d T
N Rn d = 0 (6.23)
N,w

Again, since this should hold for arbitrary cT , we obtain the FE-formulation of the out-of-plane equilib-
rium as
Kww w + Kw = f b,w + f l,w (6.24)
with
Z Z Z Z
Kww = T [A]
KB B d, Kw = T [A]N
KB d, f b,w = T Rn d, f l,w =
N T p d (6.25)
N
N,w

63
6.3.3 FE formulation of moment equilibrium
T T
Again, we choose the arbitrary test function in two ways (1 = {x 0} and 2 = {0 y } ) which
results in
Z Z Z  Z !

c T T
B [B]B du+ T
KN [A]B dw + T T
B [D]B + KN [A]N d T
N Mn d = 0
N,
(6.26)
again with
cT1
 
with cT = . (6.27)
cT2
Now, since we know that this should hold for arbitrary test functions, i.e. for any choice of c, we
must have:
+ Kw w
Ku u + K = f b, (6.28)
with:
Z Z
Ku = Ku = BT [B]B d, Kw = KTw = B
KNT [A] d

Z   Z (6.29)
K =
BT [D]B + KNT [A]N d and f b, = NT Mn d
N,

6.3.4 Total formulation

T }T we can write the total system of equations
uT wT
By collecting the unknowns in one vector a = {
as:

Kuu 0 Ku u f b,u

0 Kww Kw
w = f b,w + f l,w (6.30)

Ku KTw K f b,

| {z } | {z } | {z }
K a f
Please note that K is symmetric since Ku = KTu = Ku . Please also note that for symmetric laminates,
can be solved for independently of w and .
Bij = 0 which implies that Ku = 0 leading to that u

6.4 Element-wise approximation

In the derivations so far, we have assumed that the support of the shape functions Ni span over the
entire domain . However, in practice, these shape functions only have local support in the elements
surrounding the node j to which the shape function Nj is associated. This means that the shape functions
are only defined locally in the domain. Often, we speak of element shape functions Nie which is the part
of Ni with support in the element e. Thereby, we can write the approximation of our unknown variables
as a local approximation on each element as:
neno
X
ue
= ei = Ne u
Ni u e = Be u
e u e (6.31)
i=1
neno
X
we = Ni w ew
ie = N ew
we = B (6.32)
i=1
neno
X
e = e = Ne
Ni e e
e = Be (6.33)
i
i=1

where neno is the number of nodes associated with the element and u ei , w e are the associated
ie and i
degrees-of-freedom in node i. Considering the special case of a four node quadrilateral element, u e, w
e,

64
 e , Ne and N
e take on the form shown in Figure 6.1 (one example):


u1
x1

v1 y1








u2 w1 x2









v2 w2 y2

e
=
u e
=
, w e
=
u3 w3 x3


v



w



y3
(6.34)
3 4




u4 x4








v4 y4

 
e N1 0 N2 0 N3 0 N4 0 e = N1
 
N = , N N2 N3 N4
0 N1 0 N2 0 N3 0 N4

and we have
e, B
Be = N e = N
e (6.35)

x4 x3

v4 v3
y4 y3

w4 u4 w3 u3

x1 x2

v1 v2
y1 y2

w1 u1 w2 u2

Figure 6.1: Ordering of the degrees of freedom for the prototype element according to Eq. (6.36)

This means that we can obtain the FE equations on one element as:
e
Keu e f eb,u

Kuu 0 u
0 Keww Kew e
w = f eb,w + f el,w (6.36)
Keu Kew T Ke e
f eb,

| {z }| {z } | {z }
Ke ae fe

where Keuu is obtained by replacing B by Be and so on in the equations above. Finally, the global
stiffness matrix and force vector are obtained by assembling the element contributions into the global
system, i.e. by adding the contributions of each element stiffness matrix and element force vector to the
corresponding positions in the global matrices, cf. any basic course in FEM.

6.4.1 Prototype example: The four node isoparametric quadrilateral element

6.4.1.1 Preliminaries
The isoparametric quadrilateral (bilinear) element fulfil both requirements for convergence:

Compatibility: The approximation of all fields is continuous across the common boundary of two
neighbouring elements

Completeness: The approximation over the element is able to represent both an arbitrary constant
gradient as well as an arbitrary constant value of any of the fields.

6.4.1.2 Parent and global domain

If we consider a four-noded finite element in the xy-plane as part of the discretisation (or mesh) of a
plate, cf. Figure 6.2 (right), the basics of the isoparametric mapping is that this element in the so-called

65
global domain (xy-plane) corresponds to an element in the so-called parent domain (-plane) with edges
parallel with the coordinate axes = 1 and = 1, cf. Figure 6.2 (right). Furthermore, there exist
a one-to-one mapping from the parent domain to the global domain such that any point (x, y) can be
found as:
x = x(, ), y = y(, ) (6.37)
Furthermore, if we differentiate x and y we obtain:

x x y y
dx = d + d, dy = d + d (6.38)

or
x x

   
dx d
= (6.39)
dy y y d

| {z }
J
In order for the mapping to be one-to-one, we need to be able to invert the Jacobian J in Eq. (6.39).
Thus we have the requirement that
det J > 0

Figure 6.2: Mapping of four-node isoparametric quadrilateral element from the parent domain to the
global (FE) domain. (from Ottosen and Petersson, Introduction to the Finite Element Method, Prentice
Hall, 1992)

6.4.1.3 Element shape functions and their derivatives

The basics of the isoparametric mapping is that the mapping of a point (x, y) is represented via the
element shape functions Nie such that:

x(, ) = N1e (, )x1 + N2e (, )x2 + N3e (, )x3 + N4e (, )x4 (6.40)
y(, ) = N1e (, )y1 + N2e (, )y2 + N3e (, )y3 + N4e (, )y4 (6.41)

or
e (, )xe , y(, ) = N
x(, ) = N e (, )ye (6.42)
with
x1 y1
e = N1 (, ) N2 (, ) N3 (, ) N4 (, ) , xe =
  x 2
e y2

N x3 , y = y3

(6.43)
x4 y4

66
where xi , yi are the coordinates of node i numbered counter-clockwise and where the shape functions are
(considering the numbering according to Figure 6.2):

1
N1e = ( 1) ( 1) (6.44)
4
1
N2e = ( + 1) ( 1) (6.45)
4
1
N3e = ( + 1) ( + 1) (6.46)
4
1
N4e = ( 1) ( + 1) (6.47)
4
In the FE-formulation above, it is clear that we need the derivatives of the shape functions with
respect to x and y. By noting that

e e x N e y
e
N N x y N

x +
y x

= = (6.48)


N e N e x N e y x y N e
+
x y y
| {z }
JT

the derivatives with respect to x and y may be obtained as

e


N e N e N2e N3e N4e

N
1
1
x x x x x
 
= T
e = J (6.49)

N e N1e N2e N3e N4

Ne


y y y y y

Finally, JT can be found as

e

N N1e N2e N3e N4e

x1 y1
 x2 y2

JT = xe ye = (6.50)

N e N e N2e N3e N4e x3 y3
1
x4 y4

6.4.1.4 Numerical integration

We start by considering that the integral over an element e can be written as:
Z e Z Z
d = (x, y) dxdy (6.51)

It can be shown that by using the isoparametric mapping (x = x(, ) and y = y(, )) the integral can
be rewritten as: Z Z Z 1 Z 1
(x, y) dxdy = (, ) det Jdd (6.52)
1 1

Furthermore, in any textbook containing numerical integration, it is shown that the integral in one
dimension Z 1
f () d (6.53)
1

may be approximated, with close accuracy as:

Z 1 n
X
f () d f (i )Hi (6.54)
1 i=1

67
where a finite number (n) of so-called integration points i (in which the function to be integrated is
evaluated) is used, and where Hi is the integration weight associated with that point. This can be
expanded to two dimensions as:
Z 1 Z 1 n X
X m
f (, ) dd f (i , j )Hi Hj (6.55)
1 1 i=1 j=1

where each integration point ij now has the coordinates (i , j ).

There exists several ways of choosing the positions of the integration points and the corresponding
weights. A common (and the most accurate) scheme is the so-called Gauss integration scheme for which it
can be proven that in one dimension, the Gauss integration using n integration points exactly integrates
a polynomial of order 2n 1.
In our case, if we want to integrate the stiffness matrix contributions Kuu , Ku , ... and the force vector
contribution f l,w as accurately as possible, we start by noting that the highest order of accuracy necessary
is for a polynomial of order 3 (since that is approximately what get by multiplying N (bilinear) with B
(linear)). Thus, if the integration should be performed in one dimension we would need two integration
points, Now, since it is an area integral performed in two dimensions, we need 2x2 integration points, cf.
Figure 6.3.

Figure 6.3: Location of integration points of the Gauss scheme for 1 x 1, 2 x 2 and 3 x 3 point integration
in the parent domain (, ). (from Ottosen and Petersson, Introduction to the Finite Element Method,
Prentice Hall, 1992)

For the Gauss integration scheme, it can be shown that for a 2D integration with n x n integration
points, we have the positioning of the integration points and the corresponding weights as:

For n = 1
Point 1 (n = 1, m = 1) : 1 = 0, 1 = 0, H1 = 2, H1 = 2

For n = 2  
1 1
Point 1 (n = 1, m = 1) : (, ) = , , H1 = 1, H1 = 1
3 3
 
1 1
Point 2 (n = 1, m = 2) : (, ) = , , H1 = 1, H2 = 1
3 3
 
1 1
Point 3 (n = 2, m = 1) : (, ) = , , H1 = 1, H2 = 1
3 3
 
1 1
Point 4 (n = 2, m = 2) : (, ) = , , H2 = 1, H2 = 1
3 3

6.4.1.5 Reordering the degrees-of-freedom before assembly

In order to derive the expression for the stiffness matrix and force vector for the laminate element, it was
convenient to arrange the degrees-of-freedom in the element degree-of-freedom vector (cf. also Figure 6.1)

68
ae as

u1
x1

v1 y1









u w1

2 x2
e

u







v2 w2

e
a =
w e e
=
with u e
=
, w =
e y2
e u 3 w3 x3








v w4 y3

3






u

4 x4






v4 y4

which yields the format of the element stiffness matrix Ke and element force vector f e as in Eq. (6.36).
However, for computational reasons, it is more convenient to number the degrees of freedom as for the
element in Figure 6.4. Please note that we keep the y degrees-of-freedom defined as positive
around the negative x-axis.

a 20 a 15

a 17 a 12
a 19 a 14

a 16 a 11
a 18 a 13

a5 a 10

a2 a7
a4 a9

a3 a1 a8 a6

Figure 6.4: Ordering of the degrees of freedom for the prototype element according to the global num-
bering.

Therefore, first determine the element contributions according to Eq. (6.36) and then reorder the rows
of the force vector and rows and columns of the stiffness matrix to meet the numbering in Figure 6.4.
This final step can in MATLAB be accomplished e.g. by:

Ke = Ke([1 2 9 14 13 3 4 10 16 15 5 6 11 18 17 7 8 12 20 19],...
[1 2 9 14 13 3 4 10 16 15 5 6 11 18 17 7 8 12 20 19])

fe = fe([1 2 9 14 13 3 4 10 16 15 5 6 11 18 17 7 8 12 20 19])

6.4.2 Reduced integration of out-of-plane shear parts

It can be shown that the current formulation with linear (bi-linear) approximations of u, w and
experience locking (i.e. the element become too stiff) when fully integrated. This can be seen by
considering a plate in pure bending, cf. Figure 6.5 for a view in the xz-plane, for which the out-of-
plane shear stresses and strains should be zero. However, if we use a linear (or bi-linear) approximation
of the out-of-plane displacement, as indicated by red piecewise linear curves (one for each element) in
Figure 6.5a, and the same approximation for the rotation angles, zero shear strains xz = x + w/x
can only be accomplished over the whole element in the case of both x and w/x being constant over
the element. This means that the plate can not have any curvature. Hence, locking occurs for this case,
cf. the sketch in Figure 6.5b.
The locking effect can be avoided by so-called selective integration where parts of the stiffness matrix
are integrated by reduced integration (1 x 1 integration point in the plane of a plate). This means that
if the terms involving e.g. xz are integrated using just one integration point in each element, cf. the
cross in Figure 6.5b, it is possible to achieve xz = 0 (in this point) for these elements. The rest of the
terms (not involving out-of-plane shear strains) can still be fully integrated (2 x 2 integration points in
the plane).

69
 -x,2
z
x

-x,1
xz x
x,2
w2
x,1
w1
-w/x = -(w2 - w1 )/x
x

a) b)

Figure 6.5: Sketch explaining the cause of locking for bilinear quadrilateral Mindlin plate elements with
full integration. In a), the red curves indicate the linear approximation of the out-of-plane displacement
w for the plate which is sketched in black. Also the resulting rotations in the nodes are indicated by x,1
and x,2 . In b), the resulting x variation (black), w/x variation (red) and the resulting shear strain
xz are indicated together with the point in which xz = 0 can be obtained (black cross).

To generalise this discussion, it can be shown that in order to avoid locking, all parts of K involving
need to be integrated with a reduced order. Thus, if we subdivide
the out-of-plane shear stiffness ([A])
K as: Z Z
K = T
B [D]B d +
KNT [A]N d = K,1 + K,2 (6.56)

the order of integration should be like:

1 x 1 integration of:
K,2 , Kw , Kw , Kww

2 x 2 integration of:
K,1 , Kuu , Ku , Ku

Please note that selective (under) integration of the stiffness matrix may lead to spurious zero energy
modes under certain conditions. Therefore, one should always carefully investigate the results obtained.
The treatment (or stabilisation) of such behaviour is however considered as outside the scope of this
course.

70
Chapter 7

Failure in composites

The definition of failure depend on the application. This means that a component can be considered as
failing long before it breaks. As an example, if a laminated composite components main purpose is to
ensure a certain stiffness, it has failed already when the first ply breaks since this will have a direct impact
on the stiffness, cf. the first kink in the load-displacement curve in Figure 7.1. However, the component
may still be able to carry an increase in load, meaning that the function is lost whereas collapse or
catastrophic failure still can be avoided. On the other hand, some components, such as crash boxes or
other components intended for energy absorption during a car crash, are designed to fail in a reliable
way. In this case, failure is considered first when the component looses all load carrying capacity.

loss of load carrying capacity

loss of stiffness

Figure 7.1: Load-displacement behaviour of a hypothetical laminate (from Agarwal et al., Figure 6-13).

In addition, in the case of composite materials, internal material failure occurs much before any change
in macroscopic behaviour can be observed. Examples of types of internal failure are:
Breakage of fibres (but not causing total ply-failure)
Microcracking of the matrix
Separation of fibres from the matrix (debonding)
Local separation of plies (local delamination)
The effects of such internal material failure is however only observable if it occurs to a large extent.

7.1 Intralaminar failure - failure within a lamina

7.1.1 Lamina failure modes
Most unidirectional materials (including UD fibre reinforced plies) exhibit a linear elastic behaviour up
to failure. Thus, for a single lamina the decrease (loss) of stiffness and loss of all load carrying capacity

71
 Loading type Possible failure modes
LU Longitudinal tension Brittle (concentrated fibre breakage)
Brittle + fibre pullout (stress concentration at fibre ends)
Brittle + interface shear failure or debonding
0
LU Longitudinal compression Micro-buckling
Shear failure
Fibre crushing
Splitting
T U Transverse tension Matrix tensile failure
Constituent debonding and/or fibre splitting
0
T U Transverse compression Shear failure of the matrix
-In combination with constituent debonding
-In combination with fibre crushing
LT U In-plane shear Matrix shear failure
Matrix shear failure + constituent debonding
Pure constituent debonding

Table 7.1: Possible lamina failure modes

occurs simultaneously. In the presentation material for this lecture, a number of different failure modes
(depending on loading type) is presented, cf. also a summary of these failure modes in Table 7.1.

7.1.2 Multi-axial lamina failure criteria

There exist a number of failure criteria to predict failure of a composite ply subjected to multi-axial
loading, some of them being presented below. All these criteria require a number of material parameters
to characterise failure under certain one-dimensional loading conditions, cf. the failure modes of a ply.
These material strength parameters, summarised in Table 7.1, can either be predicted by micromechanical
models (such as the rule of mixtures for predicting the longitudinal tensile failure stress, cf. the lecture
notes for lecture 3) or determined experimentally (interested readers are referred to Chapter 10 in the
course book).

7.1.2.1 The maximum-stress criterion

This criterion states that fracture will occur if the in-plane shear stress or any of the normal stresses
0
in the direction of the principal material axes (the L, T and T axes of a ply) exceed the corresponding
allowable stress. Thus, in order to avoid failure, the following inequalities need to be satisfied:
0
LU < L < LU (7.1)
T0 U < T < T U (7.2)
abs(LT ) < LT U (7.3)

Please note that, of course, the sign of the shear stress does not have an influence. According to this
theory, failure occurs when any of the inequalities are violated. Thereby, it is not really one fracture
criterion but five subcriteria not taking into consideration the effects of coinciding loads. Thereby, it
becomes non-conservative under certain conditions, e.g. combined tensile and shear loading.

7.1.2.2 The maximum-strain criterion

This criterion states that fracture will occur if the in-plane shear strain or any of the normal strains in
0
the direction of the principal material axes (the L, T and T axes of a ply) exceed the corresponding
allowable strains. Thus, in order to avoid failure, the following inequalities need to be satisfied:

0LU < L < LU (7.4)

0T U < T < T U (7.5)
abs(LT ) < LT U (7.6)

72
 According to this theory, failure occurs when any of the inequalities are violated. Thereby, it is neither
in this case really one fracture criterion but five subcriteria without consideration of interaction. This
criterion is indeed very similar to the maximum stress criterion, especially if we assume linear elastic
behaviour of the composite up to failure. Thereby, only minor differences can be noted that are attribute
to Poisson effects, cf Figure 7.2.

Figure 7.2: Off-axis strength predicted my maximum stress and maximum-strain theories of failure (from
Agarwal et al., Figure 5-13).

7.1.2.3 The maximum-work criterion (Tsai-Hill)

The maximum-work criterion, also denoted the Tsai-Hill criterion, states that a lamina in the state of
0
plain stress (T 0 = LT 0 = T T 0 = 0 where T is the coordinate axis pointing out of the lamina plane)
fails when the following criterion is violated:
 2  2  2   
L T LT L T
+ + <1 (7.7)
LU T U LT U LU LU

The criterion is derived on the basis of the theory for anisotropic plasticity according to Hill:

2 2
2 2
H (L T ) + F (T T 0 ) + G (T 0 L ) + LLT + M T2 T 0 + N LT
2
0 = 1 (7.8)
where F N are material constants. By adding the restrictions that

Failure in pure shear occurs when the shear stress reaches the maximum strength value (for all
2 2 2
three directions) L = (LT U ) , M = (T T 0 U ) , N = (LT 0 U )

Failure under uniaxial loading occurs when the normal stress reaches the maximum strength value
2 2
(for all three directions) H + G = (LU ) , H + F = (T U )

The ply is transversely isotropic H = G

The out of plane normal stress T 0 is zero

the parameters F N can be expressed in the material strength parameters as

 2  2  2  2  2   
L T LT T T 0 LT 0 L T
+ + + + <1 (7.9)
LU T U LT U T T 0 U LT 0 U LU LU
which is further reduced to Eq. (7.7) if we assume that also the transverse shear stresses are zero. For a
visualisation of the Tsai-Hill criterion, refer to Figure 7.3 where the failure envelope is plotted under the
condition LT = 0.

73
Figure 7.3: Visualisation of the Tsai-Hill criterion for zero shear stresses (from Zenkert and Battley,
0
Foundations of Fibre composites, Figure 4-10). 1c = LU , 2c = T0 U ,
1t = LU , 2t = T U

It should be realised that this criterion in its original form does not distinguish between the tensile
and compressive strength in longitudinal and transverse direction. Therefore, it cannot in this form
adequately model materials with different strength in tension and compression (as is the case for many
polymer composite materials). However, to improve the performance of the criterion, one should carefully
take into consideration the sign of L and T meaning that if any of these are negative, the
corresponding compressive strength should be used.
As an example, if the stress state is such that L > 0 and T < 0, the Tsai-Hill criterion should take
the form:
 2  2  2   
L T LT L T
+ 0 + <1 (7.10)
LU T U LT U LU LU
It should be remarked that the Tsai-Hill criterion indeed takes the interaction between different stress
components into account, but really without any micromechanical motivation to support the actual
interaction. Thereby, this criterion can be considered more in the form of a curve-fitting criterion which
matches the single-mode conditions. Generally, this criterion is more conservative than the maximum-
stress criterion (important in industrial applications) except in small regions in the third and fourth
quadrant (when T is negative (depending on the strength values).

7.1.2.4 The Tsai-Wu failure criterion

Another common criterion used for evaluating possible failure of composite plies is the Tsai-Wu criterion.
It states that a lamina in the state of plain stress (T 0 = LT 0 = T T 0 = 0) fails when the following
criterion is violated:
2
L T2 2
LT L L T T
0 + + + 2F12 L T + 0 + 0 <1 (7.11)
LU LU T U T0 U 2
LT U LU LU TU TU

where F12 is another material parameter that needs to be determined experimentally in addition to the
separate strength parameters. In order for the criterion to represent a closed ellipse (and not parallel
lines or a hyperbola) we need the restriction
q
1 < F12 LU LU0 0
T U T U < 1

It is very difficult to perform specific experiments to determine F12 and it is thereby often necessary
to combine several experiments to calibrate the most suitable value. It should however be remarked that
the parameter F12 influence both the slenderness ratio and the inclination of the major axis of the failure
ellipse described by the Tsai-Wu criterion, cf. Figure 7.4 for a visualisation of the Tsai-Wu criterion when
LT = 0.

74
Figure 7.4: Visualisation of the Tsai-Wu criterion for zero shear stresses (from Zenkert and Battely,
0
Foundations of Fibre composites, Figure 4-11).
1c = LU , 2c = T0 U ,
1t = LU , 2t = T U

As for the Tsai-Hill criterion, the Tsai-Wu criterion also takes multi-axial loading effects into account,
but also in this case without any real motivation for the interaction. Instead, it is a mathematically
easy-to-use criterion adapted for a limited number of measure points. This criterion is widely used, it has
however some significant limitations. Most important, if the transverse compressive strength is decreased,
the predicted strength in compressive-compressive loading is actually increased (unphysical!).

7.1.2.5 The Hashin criterion

In an attempt to account for load interaction and at the same time having a physical basis (failure mode
basis) for the criterion, Hashin (Journal of Applied Mechanics, 47:329334, 1980) developed a set of
equations to predict failure based on four failure modes (details are omitted here but can be found in the
reference):
Fibre direction tensile failure (combined effect from shear and tensile loading):
 2  2
L LT
+ < 1 when L > 0 (7.12)
LU LT U
Fibre direction compressive failure (not considering the failure mode of matrix shear failure):
 2
L
0 < 1 when L < 0 (7.13)
LU
Transverse direction matrix tensile failure:
 2  2
T LT
+ <1 (7.14)
T U LT U
Transverse direction matrix compressive failure
2 2 ! 2
T0 U
  
T T LT
+ 1 + <1 (7.15)
2LT 0 U 2LT 0 U T0 U LT U

where LT 0 U in the last criterion is the out-of-plane shear strength and where is an additional material
parameter (representing the effect of shear on the longitudinal tensile failure) that needs to be calibrated
against experimental data. Often, however, is chosen as 1. I should also be pointed out that the two
matrix failure modes also comprise the loading case of pure shear, for which we set T = 0 such that
failure in pure shear occurs when the in-plane shear stress reaches the critical shear stress value LT U .

7.1.2.6 The LaRC set of criteria

To further consider the physical mechanisms behind failure initiation in a laminate, a set of failure criteria
denoted the LaRC criteria have been proposed in the literature, cf. Pinho et al. (Pinho et al., Failure

75
Models and Criteria for FRP Under In-Plane or Three-Dimensional Stress States Including Shear Non-
Linearity. NASA/TM-2005-213530. NASA Langley Research Center. Hampton, VA 23681, 2005). The
criteria are based on physical models for each failure mode and also take into consideration non-linear
matrix shear behaviour. Furthermore, the model for matrix compressive failure is based on the so-called
Mohr-Coulomb criterion and it predicts the fracture angle. Fiber kinking is triggered by an initial fiber
misalignment angle and by the rotation of the fibers during compressive loading. The plane of fiber
kinking is also predicted by the model. LaRC consists of 6 expressions that can be used directly for
design purposes. This set of criteria are however quite advanced and considered out-of-scope in the
current course. Interested readers are referred to the reference mentioned above (also available on the
course homepage).

7.2 Progressive failure of a laminate - The analysis of laminates

after initial intralaminar failure
7.2.1 Preliminaries
Given the failure criteria presented in Section 7.1.2, it is rather straightforward to predict whether any
of the plies will fail under certain loading conditions. This means that the load when the first ply will
fail can be computed. However, since the strength of a ply is directly related to its orientation, not all
plies will fail at the same time. Instead, one will see progressive ply failure as indicated in Figure 7.1.
Failure will most likely first occur in plies oriented perpendicular to the loading (due to the lower
strength). This may in fact occur at rather low loads meaning that the laminate at this point is in no
real danger of complete failure. Sometimes even, the effect of the failure of the first couple of plies may
be difficult to observe from the macroscopic viewpoint. However, as the number of failed plies increase,
the macroscopic behaviour will change (e.g. manifested as a decrease in stiffness). Finally, when enough
plies have failed, a critical point will be reached after which no further increase in the load can be made
without total failure and fracture of the component. To analyse the ultimate strength of a laminate and
to be able to follow the entire load-displacement history, a procedure where progressive failure is taken
into account is thereby necessary.

7.2.2 Simplified approach to handle progressive failure

One approach to invoke the progressive failure in the modelling of laminates is to update the lamina
properties whenever failure is detected. For instance, if longitudinal tensile failure occurs, it is very
likely that the stiffness in this direction is reduced (or vanishes completely). This can be modelled by
modification of the [Q]-matrix of the current ply in such a way that EL is set to zero (or rather a
very low number to avoid numerical difficulties) whereas the other properties, such as the transverse
modulus ET and the in-plane shear modulus GLT can be considered as unaffected by the failure mode.
In a general sense, the modification of particular material parameters should be made considering the
associated failure mode. One such procedure, could be as described in Table 7.2 (from Zenkert and
Battley (Zenkert and Battley, Foundations of Fibre Composites, Paper 96-10, KTH, 2003)).
It should be remarked that in order for such a failure assessment procedure to work well, it is
essential to involve a physically based intralaminar failure criterion (or set of criteria) such that the
associated failure mode in each ply is indicated. If instead a more pragmatic curve-fitting failure criterion
(such as Tsai-Hill and Tsai-Wu) is used, the only option is to reduce all stiffness properties (since the
mode of failure is NOT indicated) following the so-called ply discount method. This ply discount method
however, provides a more conservative estimation of the total load carrying capacity of the laminate.
Please note that whenever any of the material parameters for a ply is updated, new [A], [B], and
[D] matrices need to be generated, cf. the work flow described in Figure 7.5. Please also note that, due
to failure, an initially symmetric ply will most probably become un-symmetric whereby a coupling is
introduced between membrane and bending action (i.e. [B] 6= 0)
It should also be remarked that by applying such a stiffness modification approach in a finite
element analysis, based on a stress (or strain) based failure criterion, the results may be mesh dependent.
Therefore, the mesh convergence of the results should always be carefully considered.

76
 Material state Material constants to be reduced (or set to zero)
No failure
Transverse matrix failure ET
Fibre failure (tension) EL
Fibre/matrix shear failure LT , GLT
Shear failure GLT
Matrix failure + fibre/matrix shear failure ET , LT , GLT
Fibre/matrix shear + shear failure ET , LT , GLT
Total failure EL , ET , LT , GLT

Table 7.2: Stiffness reduction scheme for progressive failure analysis (Reproduced from Zenkert and
Battley (Foundations of Fibre Composites, Paper 96-10, KTH, 2003).

Figure 7.5: Work flow for progressive failure analysis for laminated fibre composites (From Zenkert and
Battley (Foundations of Fibre Composites, Paper 96-10, KTH, 2003).

7.2.3 Application to in-plane loading of cross-ply laminates

Consider a cross-ply laminate made up of n equal plies of which l have the fibres oriented in the direction
of the load, here equal to the x-direction, and m have the fibres oriented perpendicular to the load
(l + m = n). If we set the longitudinal and transverse stiffnesses (moduli) to EL and ET respectively, we
have by the rule of mixtures that the composite modulus E (in the loading direction x) is

l m
E= EL + ET (7.16)
n n
which sometimes is denoted the primary modulus. Adopting the maximum-strain criterion, failure in the
plies will occur when the composite strain is equal to T U at which stage the composite stress C , defined
as the stress resultant force Nx divided by the laminate thickness h as

Nx
C = , (7.17)
h
becomes
C = A = ET U . (7.18)
Let us now first assume that there is no stress relaxation in the failing plies. Then, when the load
increase further, the transverse plies will fail and the remaining secondary modulus ES becomes (by
setting ET = 0)
l
ES = EL . (7.19)
n

77
Final fracture will then occur when the load is such that the composite fracture strain is reached ( = LU )
resulting in the composite failure stress F as:

F = A + ES (LU T U ) (7.20)

The total response can be illustrated by studying Figure 7.6. Based on the stress-strain curve in Figure 7.6,
we can also write the composite strain as:

C

E C A
= (7.21)
A + C A
A < C < F
E ES

or in the latter stage using an effective stiffness Ee as

C
C A
E

= (7.22)
C
A < C < F
Ee

where
E
Ee (C ) = (7.23)
1 + [(E/ES ) 1] [1 (C /A )]

Figure 7.6: Stress-strain diagram for a cross-ply laminate up to failure (from Agarwal et al, Figure 6-14).

It should be remarked that, due to the assumption made above that there is no stress relaxation
in the failing plies, the resulting composite stress-strain (or load-displacement) curve is continuous, cf.
Figure 7.6. However, in many relevant cases, (partial) stress relaxation will occur in the failing plies which
will yield a sudden change in the magnitude of the stress in the case of displacement (or strain) controlled
loading (path C) or the magnitude of the strain for load controlled loading (path B), cf. Figure 7.7. Please
also note that the intermediate path (path D) is possible.
The level of discontinuity in the laminate response will depend how pronounced the stress relaxation
in the failed plies is. If the plies were to be considered more or less independent of each other, full stress
relaxation would be expected. But in most real applications, when the interlaminar bond is strong the
adjoining 0 plies will restrain the failure in the 90 plies. As a result, the failure of a 90 ply is localised,
and only partial stress relaxation is to be expected.
Read the remaining part of Chapter 6.8 on your own! Pay special attention to the
experiments by Hahn and Tsai that better explains the true stress relaxation in failing
plies

78
Figure 7.7: Stress-strain diagram for a cross-ply laminate with (partial) stress relaxation in the plies that
have failed. (from Agarwal et al, Figure 6-15).

7.3 Fracture mechanics of fibre composites

7.3.1 Preliminaries
Up until now we have treated the composites as more or less free of defects, flaws, microcracks etc.
whereby a continuum approach has been valid. However, in order to account for both defects and
stress concentrations, e.g. holes, we need to resort to the theory of fracture mechanics with emphasis
on anisotropic material. This theory can be used to determine the criticality of an existing crack and
thereby also maximum allowable crack size within a component (normally to be related to the size of
detectable cracks during inspection).
Generally, failure within composites (just as in metals) will emanate from existing flaws or microcracks
such as broken fibres, matrix cracks and debonded interfaces between fibres and matrix. From there,
cracks will propagate very much like in the schematic in Figure 7.8. Hence, at some distance in front of
the crack, the fibres remain intact. Closer to the vicinity of the crack tip, fibres will break due to the high
stresses, although not necessarily in the plane of the crack. As a consequence, fibre pullout will occur
behind the crack tip. In the case of strong bonding between the matrix and the fibres, it is more likely
that the fibres break in the crack plane in front to the crack tip, leaving a pure matrix dominated zone
which will break in a ductile (plastic deformation) manner.

Figure 7.8: A model of the fracture-process zone (from Agarwal et al, Figure 8-11).

It should be remarked that for laminated composites, the direction of crack propagation will be
strongly affected by the fibre orientation (as opposed to metals where the crack propagation direction
is more or less entirely governed by the stress state in the vicinity of the crack tip). In fact, cracks
may very well propagate in different directions in different plies with alternating fibre

79
orientation. Another conclusion that can be made is that practically all laminates are highly notch
sensitive and that several of the following parameters may affect this notch sensitivity:
Fibre volume fraction
Fibre properties
Manufacturing conditions
Environmental conditions
There exists several methods to assess the notch sensitivity and strength of a notched composite of
which the approach proposed by Whitney and Nuismer is one of the more commonly used, partly because
it is simple to operate.

7.3.2 Whitney-Nuismer criteria for holes

The basic motivation for the Whitney-Nuismer failure criteria for notched laminates (with holes) is the so-
called hole-size effect which means that larger holes causes a greater strength reduction of the laminate.
This can be explained by the fact that even though the stress concentration factor (and thereby the
maximum stress) is independent of the size of the hole, a larger hole will produce large stresses over a
larger area with higher chance for the existence of flaws that will initiate failure, cf. Figure 7.9 for a
sketch of the hole-size effect for isotropic plates.

Figure 7.9: Sketch of the motivation behind the so-called hole-size effect, i.e. that larger holes causes a
greater strength reduction of the laminate (from Zenkert and Battley, Foundations of Fibre composites,
Figure 6-8).

As a starting point for the derivation of these criteria, let us consider an infinite orthotropic plate
with a hole with radius R subjected to a uniform load far away from the hole, as shown in Figure 7.10.
If the axes x and y are assumed to be normal to the planes of elastic symmetry (the laminate is specially
orthotropic with respect to in-plane loading, cf. Chapter 6.6.2 in the course book), the normal stress y
acting along the line y = 0 can be approximated by
(  2  4 "    8 #)
6
R R R R
y (x, 0) = 2+ +3 (kT 3) 5 7 (7.24)
2 x x x x

where kT is the orthotropic stress concentration factor

s
A11 A22 A212
p 
2
kT = 1 + A11 A22 A12 + (7.25)
A11 2A66

where index 2 relates to the loading (y) direction opposed to the way it is written in the book
where, in that case, 1 denotes the direction of the applied load. Please note that for a uni-
directional ply loaded in the longitudinal direction (in this case L = 2), the expression for the stress

80
concentration factor takes the form
v !
u r
u E L EL
kT = 1 + t2 LT + (Note the error in the course book!) (7.26)
ET GLT

Figure 7.10: An orthotropic plate with a circular hole of radius R (from Agarwal et al, Figure 8-17).

Now, Whitney and Nuismer propose two ways to predict failure:

1. The point-stress criterion saying that failure will occur when the stress at a distance d0 from the
notch (hole edge or crack tip) reach the unnotched strength 0 of the material

2. The average stress criterion saying that failure will occur when the average value of the stress over
some fixed distance a0 from the notch (hole edge or crack tip) reach the unnotched strength 0 of
the material

For both cases, the underlying assumption is that d0 and/or a0 are laminate strength parameters which
are constant and applicable to any type of notch for a given material system, volume fraction of fibres
and lay-up sequence. Please note that, thereby, the in-plane strength of the laminate 0 will depend on
the material system as well as the orientation of the fibres.

7.3.2.1 Point-stress criterion for a hole

According to the statement above, the point-stress criterion for holes states that failure occurs when the
stress at a distance d0 from the hole edge reach the unnotched strength 0 such that

y (R + d0 , 0) = 0 , (7.27)

cf. also Figure 7.11 for an illustration.

Inserting Eq. (7.27) in Eq. (7.24) yields the relation between the notched (N ) and unnotched (0 )
strength (corresponding to the critical stress applied far away from the hole) as

N 2
= (7.28)
0 2 + p1 + 3p1 (kT 3)(5p61 7p81 )
2 4

where
R
p1 = . (7.29)
R + d0
Interesting to note is that, for very large holes p1 1 and consequently N /0 (1/kT ) whereas for
vanishingly small holes p1 0 and consequently N /0 1.

81
Figure 7.11: Illustration of the Whitney-Nuismer point-stress criterion for laminates with a hole (from
Zenkert and Battley, Foundations of Fibre composites, Figure 6-28).

7.3.2.2 Average-stress criterion for a hole

The average-stress criterion for holes states that failure occurs when the averaged stress over a distance
a0 from the hole edge reach the unnotched strength 0 such that
Z R+a0
1
y (x, 0) dx = 0 , (7.30)
a0 R

cf. also Figure 7.12 for an illustration (where y is the averaged y -stress over the distance a0 and

1t = 0 ).

Figure 7.12: Illustration of the Whitney-Nuismer average-stress criterion for laminates with a hole (from
Zenkert and Battley, Foundations of Fibre composites, Figure 6-28).

Inserting this condition in Eq. (7.24) yields the relation between the notched (N ) and unnotched
(0 ) strength as
N 2(1 p2 )
= (7.31)
0 2 p22 p42 + (kT 3)(p62 p82 )
where
R
p2 = . (7.32)
R + a0
Interesting to note is that also in this case, for very large holes p2 1 and consequently N /0 (1/kT )
whereas for vanishingly small holes p2 0 and consequently N /0 1.

82
 Figure 7.13: An orthotropic plate with a sharp crack of length 2c (from Agarwal et al, Figure 8-18).

7.3.3 Whitney-Nuismer criteria for cracks

In order to derive similar criteria for the case of a cracked laminate, let us now instead consider an infinite
orthotropic plate with a sharp crack of length 2c subjected to a uniform load y far away from the crack,
as shown in Figure 7.13. If the axes x and y are assumed to be normal to the planes of elastic symmetry
(of the laminate), the normal stress y acting along the line y = 0 can be approximated by

x k1 x
y (x, 0) = =p x>c (7.33)
2
x c2 c (x2 c2 )

where, in the latter stage, the orthotropic stress intensity factor k1 was introduced as

k1 = c. (7.34)

By applying the two failure criteria to the case of a sharp crack, we get the following:

7.3.3.1 Point-stress criterion for a crack

The point-stress criterion for cracks in a laminate which is specially orthotropic with respect to in-plane
loading states that failure occurs when the stress at a distance d0 from the crack tip reach the unnotched
strength 0 such that
y (c + d0 , 0) = 0 (7.35)
which yields the relation between the notched (N ) and unnotched strength (0 ) as

N
q
= 1 p23 (7.36)
0
where
c
p3 = . (7.37)
c + d0
p
Interesting to note is that, for very large cracks (N /0 ) 2d0 /c whereas for vanishingly small cracks
p3 0 and consequently N /0 1.

7.3.3.2 Average-stress criterion for a crack

The average-stress criterion for cracks in a laminate which is specially orthotropic with respect to in-plane
loading states that failure occurs when the averaged stress over a distance a0 from the hole edge reach
the unnotched strength 0 such that
Z c+a0
1
y (x, 0) dx = 0 (7.38)
a0 c

83
which yields the relation between the notched and unnotched strength (N ) as
r
N 1 p4
= (7.39)
0 1 + p4
where
c
p4 = (7.40)
c + a0
p
Interesting to note is that, for very large cracks (N /0 ) 0.5a0 /c whereas for vanishingly small
cracks p4 0 and consequently N /0 1.

7.3.4 Concluding remarks regarding the Whitney-Nuismer criteria

It should be remarked that the benefit and usefulness of the Whitney-Nuismer approach to both holes and
cracks is that the characteristic distances remains constant for all hole or crack sizes in at least a
particular laminate of a particular material system. In this case, only two test are necessary (one
for 0 and one for a0 /d0 ) to characterise the notch strength. This does not of course hold in the general
case but there are some experimental evidence that it holds for fibre dominated laminates in glass-epoxy,
boron-epoxy and graphite-epoxy, cf. e.g. Figure 7.14 and Figure 7.15. Consequently, these criteria
appears to a suitable rough estimate of the strength of notched laminates. For further reading on the
validity of the Whitney-Nusimer criteria, interested readers can consult e.g. the paper by Eriksson and
Aronsson (Eriksson and Aronsson, Journal of Composite Materials, Vol 24, pp 456482, 1990). However,
in real applications the appearance of e.g. cracks that penetrate the entire laminate, as in the idealised
cases above, is in fact limited. Most often, crack initiate and propagate differently in the different plies
depending on the fibre orientation.

Figure 7.14: Comparison of experimentally measured and theoretically predicted strengths of [0/
45/90]2S glass-epoxy laminates containing holes (from Agarwal et al, Figure 8-19).

Figure 7.15: Comparison of experimentally measured and theoretically predicted strengths of [0/
45/90]2S glass-epoxy laminates containing holes (from Agarwal et al, Figure 8-20).

84
Chapter 8

Viscoelasticity and damping

8.1 Introduction to viscoelasticity

In general, polymers behave viscoelastically. This will have a direct impact on fibre composites with the
matrix material being a polymer. Viscoelastic materials are capable of storing elastic energy, while at the
same time dissipating a certain amount of energy, especially when loaded at high rates or during a long
time interval. Being elastic means that, opposed to energy dissipating plastic materials, the viscoelastic
material will return to the original form (after some time) when all load is relieved. Furthermore,
viscoelastic materials have a fading memory which means that the current stress state depends on the
total strain history, but with a larger emphasis on the most recent part. Thus, what happened long time
ago fades away.
Generally, viscoelasticity has (at least) the following four characteristics:

When applying a constant stress, the strain will vary over time. This is denoted creep and a
representative creep curve can be seen in Figure 8.1a. Of interest is that an initial elastic strain
is obtained at the time the stress is applied. Thereafter, the strain will gradually increase up to a
limiting value at infinite time.

When applying a constant strain, the stress will vary over time. This is denoted relaxation and
a representative relaxation curve can be seen in Figure 8.1b. Of interest is that an initial elastic
stress response is obtained at the time the constant strain is applied. Thereafter, the stress will
gradually decrease (relax) down to a limiting value at infinite time.

A viscoelastic material will dissipate energy when subjected to cyclic loading, characterised by the
hysteresis loop as shown in Figure 8.1c.

A viscoelastic material is strain-rate dependent such that a stiffer (and sometimes more brittle)
response is obtained when the strain-rate is increased, cf Figure 8.1d.

8.1.1 Creep compliance

In the case of creep where the stress is constant = 0 , the strain will vary over time = (t).
Consequently, we can define the relation between the stress and strain by the time-dependent creep
compliance S(t):

(t)
(t) = S(t)0 S(t) = (8.1)
0

The shape of the creep compliance as function of time will be the same as for the strain history shown
in Figure 8.1a (the scale factor 0 relates the two).

85
Figure 8.1: Characteristic behaviour of a viscoelastic material (from Gibson, Principles of Composite
Material Mechanics (2nd ed.), CRC Press, 2007, Figure 8.1).

8.1.2 Relaxation stiffness

In a similar fashion, we define the relaxation stiffness (or modulus) as the relation between the constant
strain = 0 and time-varying stress = (t) as

(t)
(t) = C(t)0 C(t) = (8.2)
0
Please note that the shape of the creep compliance as function of time normally takes the same form
as the stress history curve shown in Figure 8.1b (the scale factor 0 relates the two).

8.1.3 Boltzman superposition principle

It can be shown that, by knowing the creep compliance S(t) (and the relaxation stiffness C(t)) and a
given stress (strain) history, the corresponding time-varying strain (stress) can be computed as
Z t
d( )
(t) = S(t ) d (8.3)
d

Z t
d( )
(t) = C(t ) d (8.4)
d

The above equations are denoted as the Boltzman superposition principle and can easily be generalised
to
Z t
dj ( )
i (t) = Sij (t ) d (8.5)
d

Z t
dj ( )
i (t) = Cij (t ) d (8.6)
d

where Sij (t) and Cij (t) are the creep compliance tensor and the relaxation stiffness tensor respectively.

86
8.1.4 Basic models for viscoelasticity
A normal way to visualise the models for viscoelasticity is by a system of springs (to represent the elastic
response) and dashpots (for the rate-dependence causing energy dissipation). Below, three of the most
basic models for viscoelasticity will be presented and discussed.

8.1.4.1 The Maxwell model

The simple model by Maxwell can be visualised by one spring with stiffness k [N/m2 ] and one dashpot
with viscosity [N s/m2 ] in series, as shown in Figure 8.2a. In order to find a way to describe the
relation between stresses, strains and their mutual rates, we start of by noting that the total strain over
the system will be split into one part over the spring 1 and one part over the dashpot 2 as:

= 1 + 2 (8.7)

Taking the rate of this equation:

d d 1 d 2
= + (8.8)
dt dt dt
and realising that we, by equilibrium, have the same stress in both components such that:
d 1 1 d
= (8.9)
dt k dt
d 2
= (8.10)
dt
we end up with the following differential equation to describe the stress-strain relation:
d 1 d
= + (8.11)
dt k dt

Figure 8.2: The Maxwell model for viscoelasticity (from Gibson, Principles of Composite Material Me-
chanics (2nd ed.), CRC Press, 2007, Figure 8.8).

Creep compliance of the Maxwell model

To study the response of the Maxwell model in creep, we set = 0 (constant) and = (t). Thereby,
the differential equation reduces to:
d 0
= (8.12)
dt
which integrated gives the following form for the strain
0
(t) = t + C1 (8.13)

In order to determine the constant C1 we make use of the initial condition (0) = k0 (the stress is applied
so quickly that the dashpot has no chance to react). Thus, we obtain the strain as
0 0 0
= C1 (t) = t+ (8.14)
k k

87
and, consequently, the creep compliance as

(t) t 1
S(t) = = + (8.15)
0 k

which has the shape as shown in Figure 8.2b. Note that this does not match what is generally observed
in experiments, cf. Figure 8.1a, which leads to the conclusion that the Maxwell model is not adequate to
describe creep in a realistic way.

Relaxation stiffness of the Maxwell model

To study the response of the Maxwell model in relaxation, we set = 0 (constant) and = (t).
Thereby, the differential equation reduces to:

1 d
0= + . (8.16)
k dt

Integration and utilisation of the initial value (0) = k 0 (cf. motivation above) we end up with the
relaxation stiffness C(t) as:
(t)
C(t) = = k et/ (8.17)
0
where = /k is the so-called relaxation time. The resulting shape of C(t) can be seen in Figure 8.2c
and it is clear that it corresponds to what is generally observed in experiments, cf Figure 8.1b. Thus, it
can be concluded that the Maxwell model can represent relaxation in a physically realistic way.

8.1.4.2 The Kelvin-Voight model

The simple model by Kelvin and Voight can be visualised by one spring with stiffness k [N/m2 ] and one
dashpot with viscosity [N s/m2 ] in parallel, as shown in Figure 8.3a. In order to find a way to describe
the relation between stresses, strains and their mutual rates, we start of by noting that the strain over
the individual components will be the same:

= 1 = 2 (8.18)

and that the total stress will be the sum of the contributions from the spring and the dashpot. Thereby,
we directly obtain the differential equation to describe the Kelvin-Voight model as:

d
= k+ (8.19)
dt

Figure 8.3: The Kelvin-Voight model for viscoelasticity (from Gibson, Principles of Composite Material
Mechanics (2nd ed.), CRC Press, 2007, Figure 8.9).

Creep compliance of the Kelvin-Voight model

It is rather straightforward to prove the expression of the creep compliance to be:

1 
S(t) = 1 et/ (8.20)
k

88
where = /k is the retardation time. Show this on your own! (use the same argumentation for the
dashpot as in the description of the Maxwell model). It can from Figure 8.3b be concluded that the shape
of the creep compliance is reasonable, cf. Figure 8.1a, except for the lack of an initial elastic response.

Relaxation stiffness of the Kelvin-Voight model

It is straightforward to prove the expression of the relaxation stiffness to be:

C(t) = k. (8.21)

It can from Figure 8.3c be concluded that the shape of the relaxation stiffness is unphysical (should
decrease and reach an asymptotic value), cf. Figure 8.1b. Thus, it can be concluded that the Kelvin-
Voight cannot represent creep in a physically realistic way.

8.1.4.3 The Zener model (or three parameter model)

As can be realised from the analysis of the two most basic models for viscoelasticity above, there is a
need for a somewhat more advanced model in order to have a model that can represent both limiting
cases (creep and relaxation) in a physically acceptable way. One such model is the Zener single relaxation
model (sometimes denoted the three-parameter model) as shown in Figure 8.4a. Thus, this model can
be described by the parallel coupling between one spring with stiffness k0 and one spring (k1 ) and one
dashpot (1 ) in series. Furthermore, it can be shown that the differential equation for this model can be
written as:
1 d 1 d
+ = k0 + (k0 + k1 ) (8.22)
k1 d t k1 dt

Figure 8.4: The Zener model for viscoelasticity (from Gibson, Principles of Composite Material Mechanics
(2nd ed.), CRC Press, 2007, Figure 8.10).

Creep compliance of the Zener model

It can be shown (cf. the tutorial session) that the creep compliance for the Zener model can be written
as
 
1 k1
S(t) = 1 et/1 (8.23)
k0 k0 + k1

where the retardation time is given by 1 = k0 k1 1 (k0 + k1 ). From the shape of this curve, cf. Figure 8.4b,
it can be concluded that the Zener model can represent creep in a physically realistic way.

Relaxation stiffness of the Zener model

It can be shown (cf. the tutorial session) that the relaxation stiffness for the Zener model can be
written as  
C(t) = k0 + k1 et/1 (8.24)

where the relaxation time is 1 = 1 /k1 . From the shape of this curve, cf. Figure 8.4c, it can be conclude
that the Zener model can represent also relaxation in a physically realistic way.

89
8.1.5 Viscoelastic material subjected to sinusoidal loading
In many practical structural applications, the response under oscillatory loading is of major importance.
Therefore, we will in the following section investigate the response of a viscoelastic material under sinu-
soidal loading
(t) = 0 sin (t) (in the following we use
(t) and (t) to denote sinusoidal variations in
time). We will use the Zener model as a representative of the viscoelastic response since it is the simplest
possible model with physically acceptable behaviour.
We start by noting that we can, by making use of complex numbers, write the time-varying stress as1 :

(t) = Imag 0 eit

 
(8.25)

where Imag [] represents the imaginary part of and where i = 1. However, as will be shown below,
if we drop the Imag-part, more valuable information concerning the energy dissipation and damping can
be extracted, whereby we in fact consider the applied load as

(t) = 0 eit
(8.26)

in the following. Furthermore, we realise that the resulting strain will oscillate with the same frequency,
but that it due to the rate-dependence may follow the stress with a certain delay, or phase-lag, such
that
(t) = 0 eit (t) = 0 ei(t)
(8.27)
If we insert the expression for
(t) and (t) into the differential equation of the Zener model, Eq. (8.22),
we obtain    
1 1

(t) 1 + (i) = k0 + (k0 + k1 ) (i) (t) (8.28)
k1 k1
1
which, if we multiply by 1 k1 (i), can be reformulated to obtain the following relation between the
stress and strain:
k0 + 21 (k0 + k1 ) 2
 
1 k1

(t) = +i (t) (8.29)
1 + 21 2 1 + 21 2
or
(t) = E ()
(t) (8.30)
where E () is the complex Youngs modulus

E () = E 0 () + iE 00 () (8.31)

with
k0 + 21 (k0 + k1 ) 2
E 0 () = (8.32)
1 + 21 2
1 k1
E 00 () = . (8.33)
1 + 21 2

Here, E 0 () is the so-called storage modulus and E 00 () is the so-called loss modulus. By introducing
the frequency-dependent loss-factor (), which is a measure of the relative damping properties of the
material, Eq. (8.31) can be reformulated as

E () = E 0 () (1 + i()) (8.34)

where, consequently, () is defined as

E 00 ()
() = . (8.35)
E 0 ()
In Figure 8.5, the generic dependence of E 0 () and E 00 () on is shown. It can be seen that for low
frequencies ( 0) the storage modulus approach the long-term stiffness (E 0 () k0 ) whereas the loss
modulus approach zero (E 00 () 0). Thus, as expected the rate-independent response is obtained for low
loading rates (=low frequencies). On the other hand, for large frequencies ( ) the storage modulus
approach the dynamic stiffness E 0 () k0 + k1 meaning that the response is physically acceptable also

90
Figure 8.5: Variation of storage modulus, E 0 (), and loss modulus, E 00 (), with frequency for the Zener
model (from Gibson, Principles of Composite Material Mechanics (2nd ed.), CRC Press, 2007, Figure
8.33).

at high loading-rates. In addition, it can be seen that the loss modulus is maximised for an angular
frequency of 1/1 .
It should also be noted that the complex Youngs modulus can be written as:

E () = |E ()| ei (8.36)

where
p
|E ()| = E 0 2 + E 00 2 (8.37)
 00 
E ()
= arctan (8.38)
E 0 ()

(t) = E ()
which, inserted into (t), yields

(t) = 0 eit = |E ()| ei 0 ei(t) (8.39)
 00 
E ()
= = arctan (8.40)
E 0 ()
0
0 = (8.41)
|E ()|

8.1.5.1 Generalisation
It can be shown that Eq. (8.30) can be generalised to the following form


i (t) = Cij ()
j (t) (8.42)

where

Cij () = iCij () (8.43)
in which Cij () is the Fourier transform of the relaxation stiffness tensor Cij (t) as
Z
Cij () = Cij ( ) ei d (8.44)

Thus, in the 1D case, E () can be directly obtained e.g. for the Zener model by taking the Fourier
transform of Eq. (8.24) as Z  

E () = i k0 + k1 e /1 ei d (8.45)

1 Since 0 eit = 0 (cos(t) + i sin(t))

91
8.2 Application to fibre composites - damping
Making use of the derivations above, we can now address the dynamic response and damping properties
of a fibre composite laminate in which the fibres can be seen as more or less purely elastic whereas the
matrix material (if being a polymer) behaves viscoelastically.
First, we note that in order to enable a simple procedure for the combination of the material properties
(homogenisation of properties), the following conditions must be satisfied (cf. also Figure 8.6):

The fibre dimension (diameter) d must be much smaller than the characteristic length of the
laminate L, i.e. d << L. This is in fact a necessary condition also for the preceding discussion
regarding homogenisation of the elastic properties.

The dimension (diameter) d must be much smaller than the wave length of the oscillations, i.e.
d << , meaning that the applicability of the approach presented below is limited to sufficiently low
frequencies. What the actual limit is can only be determined by comparisons between the predicted
response and experimental measurements of the same.

Figure 8.6: Critical dimensions for the validity of the homogenisation (from Gibson, Principles of Com-
posite Material Mechanics (2nd ed.), CRC Press, 2007, Figure 8.6).

Assuming that the requirements specified above are met, we can continue to characterise the dynamic
response and damping properties of the laminate. First, we note that, when subjected to sinusoidal
loading, we have the following situation

if (t) = Cij fj (t) = [Q

ij ]f f (t)
j for the fibres (8.46)

im (t) = Cm,ij
()m
j (t) for the matrix (8.47)
c
cj (t)
i (t) = Cc,ij () for the composite (8.48)

in which [Q ij ]f is the normal tensor relating fibre strains to fibre stresses, C () is the complex
m,ij

stiffness for the matrix and Cc,ij () is the homogenised complex stiffness for the composite, which can be
obtained by the ordinary homogenisation rules used for the elastic properties. What this means is that

Em is replaced by Em in the models to predict the different moduli etc. Thus, the general procedure is
as follows:
0 00
1. Determine Em () = Em () + iEm ()

2. utilise Em , m , Ef , f , Vf to determine [Q ](EL , ET , GLT , LT

, T L )
0 00
3. For each ply k, transform [Q ] to obtain [Q ]k = [T1 ]1
k [Q ][T2 ]k = [Q () + iQ ()]k

92

4. Utilise [Q ]k to compute [A ] = [A0 () + iA00 ()], [B ] = [B0 () + iB00 ()] and
[D ] = [D0 ()+ iD00 ()] matrices
  0 toobtain the relation
A B

N(t) (t)

=

M(t) B D
k(t)

5. Applied to a specific kinematical representation of e.g. a plate (e.g. Mindlin theory), also a complex
stiffness matrix can be computed on the structural level as

K = K0 () + iK00 ().

In this case, K00 () is nothing but a frequency dependent structural damping matrix Cmat () due
to internal (material) damping., cf. structural dynamics. Is should however be remarked that this
material damping due to the viscoelastic material response is only one part of the total structural
damping.

It should be remarked that it is not straightforward to realise or prove that the approach above holds,
but several experimental comparisons prove a very good agreement between predicted and experimentally
measured behaviour, cf. e.g. Figure 8.7.

a) b)

Figure 8.7: Predicted and measured a) off-axis storage modulus ratio, Ex0 /Em 0
and b) loss factor
of graphite/epoxy for various fibre orientations (from Suarez, Gibson, Sun, Chaturvedi. Experimental
Mechanics, 26(2):175184, 1986).

8.2.1 Example: In-plane sinusoidal loading of a symmetric laminate

T
Considering a symmetric laminate subjected to an in-plane load N(t) = {1, 0, 0} eit N/m and zero

moment (M(t) = 0), we note that the time-varying mid-plane strains can be obtained as

0 (t) = [A ]1 N(t)
(8.49)

Using the following notation

A1 A1 A1

11 12 16
[A ]1 = A1
12 A1
22 A1
26

A1
16 A1
26
1
A66
we obtain the different strain components as:

0x (t) = A1
11 e
it
= |A1
11 | e
i(t11 )
(8.50)
0y (t) = A1
12 e
it
= |A1
12 | e
i(t12 )
(8.51)
0
xy (t) = A1
16 e
it
= |A1
16 | e
i(t16 )
(8.52)

As can be seen, the different strain component will oscillate with the same frequency but with different
amplitudes (e.g. |A1 0x (t)) and different phase-lags (e.g. 11 for 0x (t)) where the latter are defined
11 | for

93
as
Imag [A1
 
11 ]
11 = arctan (8.53)
Real [A1
11 ]
Imag [A1
 
12 ]
12 = arctan 1 (8.54)
Real [A12 ]
Imag [A1
 
16 ]
16 = arctan (8.55)
Real [A1
16 ]

In the present case, since the amplitude of the force is unity, the amplitudes of the strain components
correspond to the so-called amplification factors (|A1
ij |), which in general relate the amplitude in force
to the amplitude in strain.

To investigate the frequency dependency of the amplification factors, we consider the special case of
a [90/0/45]S -laminate of thickness h = 12 mm with volume fraction of fibres Vf = 0.6 with the following
properties:

Fibres:
Ef = 350 GPa, f = 0.2

Matrix (represented by the Zener model):

k0 = 3.5 GPa, k1 = 3.5 GPa, 1 = 1 MPa s, m = 0.35

It can be seen from Figure 8.8a that damping results in that the amplification factors decrease with
angular frequency giving the largest amplitudes for really slow loading. Furthermore, from Figure 8.8b,
the phase-lag of 0x (t) is evident (there is a time-shift between the oscillations in force and in strain) in
the plot showing the normalised force and normalised strain.

10
x 10 strain amplification factors 1
10
0.8 normalised strain
normalised force
0.6
8 H()
normalised force and strain

x 0.4
H()
y
0.2
H()[m/N]

H()xy
6 0

0.2

4 0.4

0.6

2 0.8

1
0 200 400 600 800 1000 0 0.002 0.004 0.006 0.008 0.01
time [s]
[rad/s]

a) b)

Figure 8.8: a) Predicted amplification factors (H()x = |A1 1

11 |, H()y = |A12 | and H()xy = |A16 |)
1
0 0
and b) predicted normalised x-component of strain x (t)/|x | plotted together with the normalised force
Nx (t)/|Nx |.

8.2.2 Example: Sinusoidal bending of a symmetric laminate

In order to study the effect of laminate lay-up on the damping properties, we consider two symmetric
T
laminates [90/0]S and [0/90]S subjected to a bending moment load M(t) = {1, 0, 0} eit Nm/m and

zero in-plane load (N(t) = 0). For this case, we note that the time-varying curvatures can be obtained
as
= [D ]1 M(t)
k(t) (8.56)

94
 Using the following notation
1 1 1

D11 D12 D16
1 1 1 1
[D ] = D12 D22 D26
1 1 1
D16 D26 D66

we obtain the different curvature components as:

kx (t) = D11
1 it 1 i(t11 )
e = |D11 |e (8.57)
1 it 1 i(t12 )
ky (t) = D12 e = |D12 | e (8.58)
1 it
kxy (t) = D16 1 i(t16 )
e = |D16 |e (8.59)

As can be seen, the different curvature component will oscillate with the same frequency but with
1
different amplitudes (e.g. |D11 | for kx (t)) and different phase-lags (e.g. 11 for kx (t)) where the latter
are (in the bending case) defined as

Imag [D1
 
11 ]
11 = arctan (8.60)
Real [D1
11 ]
Imag [D1
 
12 ]
12 = arctan 1 (8.61)
Real [D12 ]
Imag [D1
 
16 ]
16 = arctan (8.62)
Real [D1
16 ]

To investigate the lay-up dependency on the amplification factors (and thereby the damping proper-
ties), we consider two laminates with different lay-up but with the same total thickness of 8 mm and the
constituent properties:

Fibres:
Ef = 350 GPa, f = 0.2

Matrix (represented by the Zener model):

k0 = 3.5 GPa, k1 = 3.5 GPa, 1 = 10 MPa s, m = 0.35

4
x 10 curvature amplification factors
9

8 H()x [90/0]S
7 H()x [0/90]S

6
H()[m/N]

1
0 200 400 600 800 1000
[rad/s]

1
Figure 8.9: Predicted amplification factor (Hx () = |D11 |) for the [0/90]S and [90/0]S laminates.

What is interesting to note in Figure 8.9 is the expected weaker response when the upper and lower
plies has fibres in the 90 -direction ([90/0]S ), but also that a pronounced frequency dependency is obtained
1
whereas in the reversed case ([0/90]S ), the amplification factor H()x = |D11 | is more or less constant.

95
This means that more damping is obtained whenever plies with fibres oriented 90 are placed at (or
close to) the top and bottom of the laminate. Thereby, it can be concluded that, not only the stiffness
properties but also the damping properties may be tailored by changing the order of the plies.

96