Three-dimensional nonlinear micro/meso-mechanical

response of the fibre-reinforced polymer composites

Z. Ullaha,1 , L. Kaczmarczyka, C. J. Pearcea

arXiv:1610.04610v1 [[Link]] 14 Oct 2016

a School of Engineering, Rankine Building, The University of Glasgow, Glasgow, G12 8LT,
UK

Abstract

A three-dimensional multi-scale computational homogenisation framework is

developed for the prediction of nonlinear micro/meso-mechanical response of
the fibre-reinforced polymer (FRP) composites. Two dominant damage mecha-
nisms, i.e. matrix elasto-plastic response and fibre-matrix decohesion are consid-
ered and modelled using a non-associative pressure dependent paraboloidal yield
criterion and cohesive interface elements respectively. A linear-elastic trans-
versely isotropic material model is used to model yarns/fibres within the repre-
sentative volume element (RVE). A unified approach is used to impose the RVE
boundary conditions, which allows convenient switching between linear displace-
ment, uniform traction and periodic boundary conditions. The computational
model is implemented within the framework of the hierarchic finite element,
which permits the use of arbitrary orders of approximation. Furthermore, the
computational framework is designed to take advantage of distributed memory
high-performance computing. The accuracy and performance of the compu-
tational framework are demonstrated with a variety of numerical examples,
including unidirectional FRP composite, a composite comprising a multi-fibre
and multi-layer RVE, with randomly generated fibres, and a single layered plain
weave textile composite. Results are validated against the reference experimen-
tal/numerical results from the literature. The computational framework is also
used to study the effect of matrix and fibre-matrix interfaces properties on the

1 Correspondence to: Z. Ullah, E-mail: [Link]@[Link]

Preprint submitted to Composite Structures October 18, 2016

homogenised stress-strain responses.
Keywords: Finite element analysis, Fibre reinforced polymer, Multi-scale
computational homogenisation, Elasto-plasticity, Cohesive interface elements,
Transverse isotropy

1. Introduction

Compared to conventional materials, fibre-reinforced polymer (FRP) com-

posites can offer exceptional physical and chemical properties (including high
strength, low specific weight, fatigue and corrosion resistance, low thermal ex-
pansion and high dimension stability), making them ideal for a variety of en-
gineering applications, including aerospace, marine, automotive industry, civil
structures and prosthetics [1, 2, 3]. Phenomenological or macro-level models
cannot accurately describe the complex behaviour of FRP composites due to
their underlying complicated and heterogeneous microstructure. Furthermore,
nonlinearities associated with the matrix elasto-plasticity and fibre-matrix de-
cohesion make the computational modelling even more challenging. Multi-scale
computational homogenisation (CH) provides an accurate modelling framework
to simulate the behaviour of FRP composites and determine the macro-scale
homogenised (or effective) response, based on the physics of an underlying, mi-
croscopically heterogeneous, representative volume element (RVE) [4, 5, 6, 7,
8, 9, 3]. The homogenised properties calculated from the multi-scale CH are
subsequently used in the numerical analysis of the macro-level structures.
A variety of numerical techniques have been developed to model the nonlin-
ear micro-mechanical response of unidirectional (UD) FRP composites, mostly
based on finite element analysis. For UD glass/carbon (G/C) FRP composites,
a computational model was developed in [10, 11] within the framework of fi-
nite deformation. Both in-plane shear and compressive loading scenarios were
considered. The Mohr-Coulomb yield criterion and cohesive interface elements
were used to model the response of epoxy matrix and fibre-matrix interfacial
decohesion respectively. Fibres were generated randomly within the RVEs using

2
the algorithm presented in [12] and were modelled as a linear-elastic and iso-
topic material. A parametric study, including the effect of matrix and interface
properties on the stress-strain response, was also conducted. The idea of [10]
was extended further in [13] by incorporating thermal residual stresses (due to
cooling of FRP composites after curing process, caused by the mismatch in ther-
mal expansion coefficients of matrix and fibres) in the simulation, in addition
to transverse tensile and cyclic loading for the CFRP composites. The nearest
neighbour algorithm (NNA) [14] developed by the same authors, was used to
randomly generate the fibres within the RVEs. Using the same constitutive
models for matrix, fibres and fibre-matrix decohesion as in [10], a multi-layer
multi-fibre (M2 ) RVE was used in [15] for laminates. Each lamina was modelled
as a cube with randomly distributed but axially aligned fibres, generated using
a fibres randomisation algorithm in DIGIMAT FE [16]. Both cross [0/90]ns
and angle [±45]ns (where the subscript ns represents n layers with the same se-
quence and symmetric about the mid plane) GFRP composites were considered
with in-plane shear loading and results of stress-strain behaviour were validated
against the experimental results. A combined transverse compression and axial
tension loading scenario was considered in [17] for UD GFRP composite. In ad-
dition to matrix plasticity and fibre-matrix decohesion, fibre breakage was also
included in the FE simulation. The pressure dependent, Drucker-Prager yield
criterion was used to model matrix plasticity and both fibre breakage and fibre-
matrix interfacial decohesion were modelled with cohesive interface elements. A
simple periodic, hexagonal fibre arrangement was assumed. In [18], a modified
von Mises yield criterion was used to model the behaviour of the matrix mate-
rial, while a maximum tensile stress criterion was used to model fibre breakage.
Fibre-matrix decohesion was also included in the simulation and was modelled
with cohesive interface elements. The random distribution of the fibres was also
included within the RVE based on the optical microscopy of real composites. A
variety of loading conditions was used subsequently to study the response of the
UD FRP composite. The limitations of different plasticity models for modelling
matrix materials, including Mohr-Coulomb and Drucker-Prager were argued in

3
[19, 20], especially in complex loading scenarios. Instead of the conventional
plasticity models, a pressure dependent thermodynamically consistent plastic-
ity model [21] was used. A statistically proven random distribution algorithm
proposed by the same authors in [22] was used to randomly generate UD fibres
within the RVEs. Similar to previous studies, fibres were modelled as linear-
elastic and isotropic material and fibre-matrix decohesion was modelled with
cohesive interface elements. A variety of RVE loading scenarios was considered
including transverse tension and compression, transverse and longitudinal shear
and combined transverse compression and transverse shear.
A number of numerical modelling approaches have been used to simulate
the behaviour of textile composites subjected to different loading scenarios. A
comprehensive review of these methods can be found in [23]. Continuum dam-
age mechanics (CDM) was used in [24] to model both matrix and yarns for
glass and carbon plain weave textile composites. Dissipated energy density was
used as damage parameter and both material and geometric nonlinearities were
included in the simulation. Further use of CDM in the simulation of textile com-
posites can also be found in [25, 26, 27]. Moreover, a three-dimensional CDM
based approach was used to simulate the progressive damage in laminated FRP
composites in [28, 29]. A variety of failure mode, including matrix tensile and
compressive cracking, fibre tensile and compressive failure, fibre-matrix shearing
and delimitation between the layers were included in the simulations. For a twill
weave textile CFRP subjected to in-plane loading, a meso-mechanical analysis
was performed in [30]. The matrix was modelled as elasto-plastic material with
the same plasticity model as in [19, 20, 21], while yarns were modelled as linear-
elastic and transversely isotropic material. Results of the RVE strain fields
and homogenised stress-strain response were validated against the experimental
results and found in a good agreement.
These numerical simulations, described above, of FRP composite behaviour
are limited to specific RVE type (2D or 3D, UD or woven/textile) or load-
ing scenarios (normal or shear). In contrast, this paper develops a fully gen-
eralised three-dimensional micro/meso-mechanical framework, which is subse-

4
quently implemented in the authors’ open source FE software, MOFEM [31].
The dominant damage mechanisms (observed experimentally [10]), i.e. ma-
trix elasto-plasticity and fibre-matrix decohesion, are included in the compu-
tational framework. Matrix material is modelled using a pressure dependent
paraboloidal yield criterion [19, 20, 30, 21] with an exponential hardening law.
Fibre-matrix decohesion is modelled with zero thickness cohesive interface ele-
ments. Yarns are modelled as linear-elastic and transversely isotropic materi-
als. Rather than simplified fibre arrangements for UD FRP composites, e.g. in
[32, 33, 17, 34], which are not the actual representation of the real FRP compos-
ites and can lead to erroneous results. this study adopts a statistically proven
random distribution algorithm proposed in [22] to generate fibre arrangements
within the RVE. The RVE boundary conditions are imposed in a unified manner
which allows convenient switching between displacement, traction and periodic
boundary conditions [35]. Hierarchic finite elements are adopted, which per-
mits the use of arbitrary order of approximation, leading to accurate results for
relatively coarse meshes. The computational framework is designed to take ad-
vantage of distributed memory high-performance computing. Moreover, CUBIT
[36] and ParaView [37] are used as pre- and post-processor respectively.
This paper is organised as follows. The computational framework is fully
described in §2. The material models are given in §2.1, consisting of material
model for matrix §2.1.1, yarns/fibres 2.1.2 and fibre-matrix interfaces §2.1.3.
The nonlinear multi-scale CH with corresponding RVE boundary conditions
are explained in §2.2. Calibration and validation of the matrix plasticity model
is given in §3. Three numerical examples are given in §4, including UD GFRP
composites §4.1, M2 RVE §4.2 and plain weave textile composites §4.3. Finally,
the concluding remarks are given in §5.

2. Computational framework

The computational framework developed for FRP composites consists of a

set of constitutive models for individual components including matrix, yarns/fibres

5
and fibre-matrix interface and implemented within the formulation of first-order
multi-scale CH.

2.1. Material constitutive models

Typical RVEs in the case of UD FRP and textile composites are shown in
Figure 1(a) and (b) respectively, consisting of yarns/fibres embedded within a
polymer matrix. The constitutive model for FRP composites is a combination
of constitutive models for these individual components, together with fibre-
matrix interface decohesion. In the following, each of these constitutive model
is explained in detail.

Full RVE Full RVE

Matrix Matrix

Fibres Yarns

(a) Unidirectional RVE (b) Textile RVE

Figure 1: Typical RVEs for unidirectional FRP and textile composites

6
2.1.1. Matrix
The polymer matrix is modelled as an elasto-plastic material using a non-
associative pressure dependent paraboloidal yield criterion [21, 19, 20, 30, 9].
This plasticity model can incorporate different yield strengths in tension and
compression and is shown Figure 2(a) in the principal stress space. The yield
function is written as

f (σ, σc , σt ) = 6J2 + 2I1 (σc − σt ) − 2σc σt , (1)

where σ is Cauchy stress tensor, I1 = tr(σ) is the first invariant of Cauchy stress
tensor, J2 = 12 η : η is the second invariant of deviatoric stress η = σ − 31 I1 and
σt and σc are yield strengths in tension and compression respectively. A non-
associative flow rule is used, for which the plastic potential function is written
as
1 − 2νplas
g (σ, σc , σt ) = 6J2 + 2αI1 (σc − σt ) − 2σc σt , α= , (2)
1 + νplas

where νplas is a material parameter and is known as plastic Poisson’s ratio.

Furthermore, the Helmholtz free energy in the case of linear isotropic hardening
is written as

1 1 1
ψ= λtr[ε]2 + µε : ε + σt0 α0 + Ht α20 + σc0 α1 + Hc α21 , (3)
2 2 2

where λ and µ are the Lame parameters, ε is the strain tensor, σt0 and σc0
are the initial yield strengths in tension and compression respectively, α0 and
α1 are internal kinematic variables and Ht and Hc are hardening parameters
in case of tension and compression respectively. Following Equation (3), yield
strengths in tension and compression are written as

∂ψ ∂ψ
σt = = σto + α0 Ht , σc = = σco + α1 Hc (4)
∂α0 ∂α1

7
A more realistic, exponential hardening law is presented in this paper, due to
which Equation (3) is rewritten as

1 Ht −nt α0 Hc −nc α1
ψ= λtr[ε]2 + µε : ε + (σt0 + Ht ) α0 + e + (σc0 + Hc ) α1 + e ,
2 nt nc
(5)
where Ht , Hc are the difference between the ultimate and yield strengths, nt
and nc are material parameters and determine the rate of convergence between
yield and ultimate strengths. Following Equation (5), yield strengths in tension
and compression are written as

∂ψ ∂ψ
= σto +Ht 1 − e−nt α0 , = σco +Hc 1 − e−nc α1 . (6)



σt = σc =
∂α0 ∂α1

A similar hardening law as a function of equivalent plastic strain was also used
in [30].

Traction

2
Fracture energy

1
3

Separation

(a) Paraboloidal yield surface (b) Transversely isotropic material (c) Traction-separation law

Figure 2: Material models for matrix and fibre-matrix decohesion

8
2.1.2. Yarns/fibres
In textile composites, as shown in Figure 1(b), yarns are interwoven together
to form a textile structure. On the micro-level, yarns are the same as UD FRP
composites consisting of bundles of glass/carbon fibres within the polymer ma-
trix. On the meso-level yarns are modelled as homogenous and transversely
isotropic material with homogenised or effective properties obtained from the
multi-scale CH of the UD FRP composites. Five material parameters are re-
quired for transversely isotropic materials, which are Ep , νp , Ez , νpz and Gzp ,
where z and p are fibres and transverse directions respectively as shown in Fig-
ure 2(b). Ep and νp are Young’s modulus and Poisson’s ratio in the transverse
direction respectively, while Ez , νpz and Gzp are the Young’s modulus, Pois-
son’s ratio and shear modulus in the fibre directions respectively. In order to
re-orient the known stiffness matrix for the transversely isotropic material from
the local coordinates to global coordinates, the yarns directions at each Gauss
point need to be determined. To do this it is possible to simply use the cubic
splines that were used to construct the yarns. However, this can lead to inac-
curacies in the case of yarns with non-uniform cross-sections along their length.
An alternative approach is used in this paper, in which the yarns directions are
determined by solving the potential flow along these yarns. A detailed descrip-
tion of this approach and how to transfer the stiffness matrix from the local to
global coordinate axes is given in [6, 3, 7, 8, 9].
On the micro-level in the case of UD FRP composites, fibres are modelled
as linear-elastic and isotropic material (a special case of transversely isotropic
material model), for which only two material parameters, i.e. Young’s modulus
and Poisson’s ratio (Ef , and νf ) are required. In this case, the following material
parameters are used:

Ef
Ep = Ez = Ef , νp = νpz = νf , Gp = . (7)
2 (1 + νf )

9
2.1.3. Fibres/matrix interfaces
Fibre-matrix decohesion is modelled using standard cohesive interface ele-
ments with a straightforward material model, i.e. linear traction-separation law
(shown in Figure 2(c)). Only three material parameters are required for the
material model, including cohesive strength ft , fracture energy Gf (shaded area
in Figure 2(c)) and material parameter β, which assign different weight to open-
ing and shear displacements. Mathematically the material model for cohesive
interface elements is written as

 E δ if δ < δ0 ,
 0


t= (1 − ω) E0 δ if δ0 ≤ δ < δmax , (8)


 0

if δ < δmax ,

p
where E0 is the initial stiffness, δ = 2 + δ 2 ) is the displacement jump
δn2 + β(δs1 s2

with δn and δs as its normal and shear components and ω is damage parameter.
κ is a history parameter and is equal to the highest value of displacement jump
δ. Furthermore, δ0 and δmax are respectively the displacement jumps at the
onset of damage (ω = 0) and when the interface is fully damaged (ω = 1). E0
and δ0 are written as
Em ft
E0 = , δ0 = , (9)
h E0

where Em is the Young’s modulus of matrix material and h is the interface

thickness. Furthermore, from Figure 2(c), the damage parameter ω is written
as
2Gf E0 + ft2 κ


ω= . (10)
2Gf (ft + κE0 )

The constitutive matrices for cohesive interface elements in the local and global
coordinate systems are written as

Dloc = (1 − ω) IE0 , Dglob = RT Dloc R, (11)

10
where I is the unit matrix and R is the transformation matrix [38]. Equations
for stiffness matrix and corresponding internal forces are written as
Z Z
el T glob
K = Φ D ΦdA, Fel
int = ΦT RT tloc dA, (12)
A A

where Φ is a matrix of shape functions and tloc is the traction vector in local
coordinates, details of which are given in [38].

2.2. Multi-scale computational homogenisation

In multi-scale CH, a heterogeneous RVE is associated with each Gauss point

of the macro-homogeneous structure, the boundary conditions for which are im-
plemented using the generalised procedure proposed in [3, 35]. Small displace-
ment and small strain formulations are used within the framework of first order
multi-scale CH, the basic concept of which is shown in Figure 3, where Ω ⊂ R3
and Ωµ ⊂ R3 are macro and micro domains respectively. Macro-strain ε =
h iT
ε11 ε22 ε33 2ε12 2ε23 2ε31 is first calculated at each Gauss point
h iT
x = x1 x2 x3 of the macro-structure, which is then used to formulate
the boundary value problem on the micro-level. After solution of the micro-level
h iT
boundary value problem, homogenised stress σ = σ 11 σ 22 σ 33 σ 12 σ 23 σ 31
and stiffness matrix C are calculated.
For a global step n + 1, the discretised system of equations in case of an
iteration i of the Newton-Raphson algorithm is written as
  
Kin+1 C T  △ui
n+1

  = Fin+1 , (13)
C 0  △λi 
n+1

where K and u are the standard FE tangent stiffness matrix and displacement
vector respectively and λ is the unknown vector of Lagrange multipliers re-
quired to impose the RVE boundary conditions. Matrix C in Equations (13)
are calculated over the boundary Γ of the RVE and are constant throughout

11
 Integration point

Homogenised stress Macro strain

Stiffness matrix

Macro-structure

Traction

Solve boundary
value problem 2

1
3

Separation

Micro-structure (RVE)

Figure 3: Multi-scale computational homogenisation

the calculations [3, 35] and are given as

Z
C= HNT Nd∂Ωµ . (14)
∂Ωµ

In Equation (14), N is a matrix of shape functions and H is a matrix that is

specific to the type of boundary conditions used, each row of which represents
an admissible distribution of nodal traction forces on the RVE boundary [35].
The specific choice of H in the case of linear displacement, periodic and uniform
traction boundary conditions can be found in [35, 3] and is not repeated here.
Matrix K comprises contributions of the matrix, yarns and yarn-matrix in-
terface elements. Fin+1 is a vector of residuals and is written as
 
 CT λi − Fint i  Z
n+1 n+1
Fin+1 = , D= HNT Xd∂Ωµ , (15)
 C ui − Dε  ∂Ωµ
n+1 n+1

where X is a matrix of spatial coordinates, evaluated at Gauss points during

12
numerical integration of the surface integrals and is given as
 
2y1 0 0 y2 y3 0
1 
X=  0 2y2 0 y1 0 y3  . (16)
 
2 
0 0 2y3 0 y1 y2

At Newton-Raphson iteration i, variable ξ = u, λ is calculated using ξin+1 =

Pi
ξn + m=1 ξ m int i
n+1 . In Equation (15), Fn+1 is a vector of internal forces. Further-

more, Cuin+1 and CT λin+1 are associated with the RVE boundary conditions
and are written as
Z Z
Cuin+1 = HNT uhn+1
i
dΓ, CT λin+1 = HNT λhn+1
i
dΓ, (17)
Γ Γ

where uh and λh are displacements and Lagrange multipliers calculated at a

Gauss point, i.e. ξ h = uh , λh = Nξen+1
i
, where ξ e is a matrix of displacements
or Lagrange multipliers associated with element e. Finally, the homogenised
stress for global increment n + 1 is written as:

1 T
σ n+1 = D λn+1 , (18)
V

To compute the homogenised stiffness matrix C at the end of global increment

n + 1, the converged matrix K is subjected to six different macro-strain pertur-
bations of unit vector leading to six linear system of equations. This will give a
set of homogenised stresses, i.e.

h i
C= σ1 σ2 σ3 σ4 σ5 σ6 , (19)

where for example:

h iT
σ 1 : for ε = 1 0 0 0 0 0
h iT . (20)
σ 4 : for ε = 0 0 0 1 0 0

13
In each of the six cases, only the right-hand side of the system of Equations (13)
changes, which is solved very efficiently as the left-hand side matrix is factorised
only once.

3. Calibration and validation of plasticity model

Following [30], the plasticity model is first calibrated against the experimen-
tal results from [39, 19, 30] for epoxy resin subjected to tensile and compressive
loading. A list of material parameters used in this case is shown in Table 1,
where E, ν, σto , σto and νplas are given in [30, 19]. Moreover, the hardening pa-
rameters, i.e. Ht , Hc , nt and nc are determined from the numerical simulation
based on the experimental stress-strain curves. The estimated parameters Ht ,
Hc are the same as given in [30] but in contrast, due to the use of hardening
law as a function of internal kinematic variables leads to different nt and nc in
our case.
The geometry considered in this case is a cube of dimension 1mm, which
is discretised with 1191 tetrahedral elements and 299 nodes. The cube is fixed
at the bottom face and subjected to tension, compression and shear loading on
the top face as shown in Figures 4(a), (b) and (c) respectively. A comparison
between numerical and experimental stress-strain responses for all the three
loading scenarios are shown in Figure 4(d). As expected the numerical and
experimental responses in tension and compression are in good agreement, as a
result of parameter fitting. The response in shear is not fitted and also shows
fairly good agreement. All of the three responses are highly non-linear and it
is clear that the plasticity model can capture them well. It must be noted that
the plasticity model requires input data for tension and compression and can
be used subsequently to simulate more generalised loading scenarios.

4. Numerical Examples

Three numerical examples are now given to demonstrate the correct imple-
mentation and performance of the developed computational framework.

14
 Parameter Value
Young’s modulus (E) 3.76 GPa
Poisson’s ratio (ν) 0.39
Plastic Poisson’s ratio (νplas ) 0.3
Initial yield strength in tension (σto ) 29 MPa
Initial yield strength in compression (σto ) 67 MPa
Ht 67 MPa
Hc 58 MPa
nt 170
nc 150

Table 1: Material parameters for epoxy resin

(a) Tension

(b) Compression (c) Shear (d) Stress-strain response

Figure 4: Comparison of numerical and experimental [39] stress-strain response for an epoxy
resin

4.1. Unidirectional GFRP composites

The first numerical example consists of polymer composites reinforced with

unidirectional glass fibres subjected to transverse tension. A similar numerical
example is also considered in [20]. Four different RVE sizes are considered in
this case, consisting of periodic, randomly distributed but axially aligned fibres
and are shown in Figure 5. The algorithm proposed in [22] is used to randomly
generate the fibres within the RVEs with diameter of 5µm and volume fraction

15
of 60%. The four RVEs are discretised with 1,499, 3,239, 13,023 and 21,140
tetrahedral elements. The detailed geometry for RVE-4, showing the individual
matrix, fibres and cohesive interface elements, is presented in Figure 6(a).
A list of material properties used for the elasto-plastic matrix materials are
the same as given in Table 1, while for the linear-elastic and isotropic glass fibres
Young’s modulus and Poisson’s ratio used are 74 GPa and 0.2 respectively. For
the cohesive interface elements, interface strength and fracture energy are 50
MPa and 2 J/m2 respectively. The macro-strain (applied to the RVEs using
periodic boundary conditions) versus homogenised stress response for all four
RVEs are compared to reference numerical results from [20] and are shown in
Figure 6(b). The RVEs are subjected to a transverse strain εxx of 1 percent. It is
clear from Figure 6(b) that the developed computational framework accurately
predicts the stress-strain behaviour in the pre-peak region (up to εxx = 0.65%)
for all of RVEs. The size effect can be clearly seen in the post-peak region
(beyond εxx = 0.65%); increasing the size of the RVE leads to a more brittle
response. A similar behaviour was also reported in [20, 40]. Issues related with
the existence and size of the RVE and pre- and post-peak region behaviour are
described in detail in [40, 41, 42, 43], where the ill-posedness of the macro-level
BVP and its non-objectivity with respect to the size of the RVE is discussed. A
detailed description of dealing with these limitation of the classical CH schemes
is given in [43]. These specialised treatments are not considered in this paper.
The final damaged RVEs with clear localisation zones/debonding are also shown
in Figure 6(b). The damaged zones consist of fully damaged cohesive interface
elements that are perpendicular to the direction of the applied strain. It is clear
that fibre-matrix decohesion interface leads to a reduction in load transfer from
the matrix to the fibres, which results in the overall stiffness. Furthermore,
strain localisation associated with the damaged zones subsequently leads to
severe plastic deformation of matrix material.
A parametric study is also conducted to investigate the effect of different
parameters on the macro-strain versus homogenised stress response. The effect
of fracture energy on the stress-strain response is shown in Figure 7(a), where

16
fracture energies of 2, 3, 4 and 100 J/m2 are used but all other parameters are
kept constant. It is clear from Figure 7(a) that the stress-strain response are
the same for all the four cases up to εxx = 0.6%. Furthermore, lower fracture
energies leads to clear damaged zone with high strain localisation (severe plastic
deformation). The effect of interface strength on the stress-strain response is
shown in Figure 7(b), where interface strength of 20, 35 and 50 MPa are used and
all other parameters are kept constant. Furthermore, the effect of considering
unlimited interface strength, i.e. ft = ∞ in shown in Figure 7(b). A clear
localisation can be seen for both ft = 35 and ft = 50 cases as compared to the
ft = 20 for the applied strain of 1 %.
The effect of a linear elastic material, as opposed to an elasto-plastic ma-
terial, on the strain-stress response is shown in Figure 7(c). In addition, three
cases with different interface strengths of cohesive interface elements, i.e. 20,
35 and 50 MPa are considered. It is clear from Figure 7(c) that in the pre-
peak regions, the use of either linear-elastic or elasto-plastic matrix material
leads to almost similar stress-strain response while in the post-peak region the
use of linear-elastic matrix material leads to relatively stiff response. The fi-
nal damaged RVEs for both ft = 20 and ft = 35 with both linear-elastic and
elasto-plastic matrix materials are also shown in Figure 7(c). The high strain
localisation in the damaged zones leads to severe plastic deformations leading
to a more brittle stress-strain response.

4.2. Multi-fibres multi-layer RVE

A multi-fibre multi-layer RVE subjected to in-plane shear is considered in

the second example. A similar example is also analysed experimentally and
numerically in [15], the stress-strain response from which is used here as a
reference. The UD FRP composite used in this case, consist of E-glass (ER-
459L) and epoxy matrix (EPOFINE-556) with FINEHARD- 951 hardeners. The
RVE geometry is shown in Figure 8(a), consisting of two cubes of dimension
1mm with randomly distributed fibres (generated using the algorithm in [22])
of 24 µm and volume fraction of 28 % and are placed on the top of each other

17
 21.4 m

16.2 m

11.4 m

5.7 m

(a) RVE-1 (b) RVE-2 (c) RVE-3 (d) RVE-4

Figure 5: Different RVE sizes for the UD GFRP example

with an angle of 90o . In Figures 8(b) and (c) individual matrix and fibres are
shown respectively. The RVE is discretised with 32,818 tetrahedral elements
and is shown in Figure 8(d), while fibre-matrix interfaces are discretised with
3,056 cohesive interface elements and are shown in Figure 8(e). Moreover, a
perfect bond is assumed between laminae.
For the linear-elastic and isotropic glass fibres, Young’s modulus and Pois-
son’s ratio are 73 GPa and 0.23 respectively. For the matrix, most of the material
parameters used are the same as given in Table 1 with the only exception of
Young’s modulus and Poisson’s ratio, which are 4.7 GPa and 0.3 respectively.
For cohesive interface elements, interface strength and fracture energy used are
30 MPa and 100 J/m2 respectively [15]. The RVE in this example is subjected
to shear strain γzx = 4%, as shown in Figure 9. The shear stress versus shear

18
Full RVE

Matrix

Fibres

Reference
Cohesive
elements

(a) RVE (b) Stress - strain curves

Figure 6: RVE components and comparison of numerical and reference stress-strain curves
for different RVE sizes for the UD GFRP example

strain (γzx versus τzx ) response is compared with the experimental and nu-
merical results from [15] and is shown in Figure 9, which are in a very good
agreement. Stress-strain response is almost linear up to γzx = 1.5%, beyond
which the response is nonlinear due to the decohesion between fibres and ma-
trix. The difference between the numerical and experimental results (especially
between γzx = 1% and 2.5%), might be due to the assumption of perfect bond-
ing between the 0o and 90o laminae. At the end of the simulation, contours of
γzx over the deformed RVE are also shown in Figure 9. Strain is very small
in the glass fibres as compared to matrix material due to the associated high
stiffness. Decohesion between matrix and fibres can also be seen Figure 9.

19
 (a) Effect of fracture energy (b) Effect of i
 
 

(c) P tic versus elasti c
i materials

Figure 7: Parametric study for the UD GFRP example

4.3. Plain weave textile composites

Finally, a plain weave textile composite subjected to a variety of normal and

shear loading conditions is considered, consisting of E-glass fibres and epoxy
matrix. A similar numerical example is also considered in [24]. An RVE, con-
sisting of similar yarns in warp and weft directions is used in this example, for

20
 0 m
m

(a) Full RVE (b) Matrix (c) Fibres (d) M (e) Cohesive elements

Figure 8: Geometry and mesh for the M2 RVE example

Simulation

Figure 9: Comparison of numerical and reference stress-strain response for the M2 RVE ex-
ample

21
which the geometry with all of the required dimensions are shown in Figure
10(a). Elliptical cross-sections and cubic splines are used respectively to model
the cross sections and paths of the yarns. The volume fraction of fibres within
the yarns is 65% while the total volume fraction of fibres within the RVE is 35%.
The RVE is discretised with 11,516 tetrahedral elements and is shown in Figure
10(b). For the elasto-plastic matrix material, the same properties are used as
given in Table 1, while for the linear-elastic and transversely isotropic yarns,
material properties are given in Table 2 [24]. Furthermore, a perfect bond is
assumed between yarns and matrix.
The yarns direction, calculated from the potential flow analysis, are shown in
Figure 11(a). Four loading conditions, including two normal (εxx and εyy ) and
two shear (γyz and γzx ) are considered, where the RVE is subjected to a macro-
strain of 3%. For the shear case, the stress-strain responses are shown in Figure
11(a). The nonlinear response beyond a strain of 1.5% is due to matrix failure.
The stress-strain response in the case of γyz is also compared with numerical
results from [24], which are in a very good agreement, especially in the linear
region (up to γyz =1.5%). Beyond γyz =1.5%, our simulation result is relatively
stiffer, which might be due to the use of linear-elastic material for the yarns. A
high strain gradient can also be seen in the matrix, especially in regions of the
thin matrix layer. The response in the case of γzx involves shearing of yarns
leading to stiffer behaviour as compared to γyz . Furthermore, response in the
case of εxx and εyy are shown in Figure 11(b). For the given range of applied
strains, stress-strain responses for both εxx and εyy are linear. Moreover, εxx
involves direct tensile load on the yarns and behave stiffer as compared to εyy
case, where strain is applied directly on the matrix.

Parameter Value
Ep 18.06 GPa
νp 0.34
Ez 48.47 GPa
νpz 0.25
Gzp 5.58 GPa

Table 2: Yarns material parameters for the plain weave textile composites example

22
 Yarns

Matrix

3

Full RVE


! 5 
1.7

( eometry (b) Mesh

Figure 10: Geometry and mesh for the plain weave textile composites example (dimensions
in mm)

5. Concluding remarks

A three-dimensional, nonlinear micro/meso-mechanical multi-scale CH frame-

work is developed for FRP composites. The matrix material is modelled as
elasto-plastic, using a paraboloidal yield surface. Decohsion of the fibre-matrix
interface are modelled using cohesive interface elements. The yarns/fibres are
modelled as linear-elastic and transversely isotropic material. It is shown that
the two dominant damage mechanisms, i.e. matrix plasticity and fibre-matrix
interfacial decohesion control the strength of FRP composites. Experimental
stress-strain results for epoxy resin for both tension and compression load cases
are used to calibrate the plasticity model and is validated subsequently for the
shear loading. Three numerical examples with a variety of RVEs and load-

23
 y&

Reference

(a) Shear "#$ %ormal

Figure 11: Stress-strain responses for the plain weave textile composites subjected to different
loading conditions

ing conditions are considered to demonstrate the performance of the proposed

computational framework. In both UD FRP composite and M2 RVE examples,
fibres are randomly generated within the RVEs using a statistically proven ran-
dom distribution algorithm. In the UD FRP numerical example, the developed
computational framework can accurately predict the stress-strain behaviour in
the pre-peak region, while in the post-peak region size dependent response is
observed, which is natural in the case of first-order computational homogenisa-
tion. A parametric study is also conducted for the UD FRP numerical example,
i.e. the effect of different matrix and fibre-matrix interface parameters on the
stress-strain behaviour and it is shown that failure starts at fibre-matrix inter-
face followed by formation of a shear band or matrix plasticity. Furthermore,
from the M2 RVE and plain weave textile composite numerical examples, it
is shown that the computational framework can accurately predict the stress-
strain behaviour of RVEs with complicated geometries subjected to different
loading scenarios. The developed computational framework is implemented in
the authors’ open-source FE software MOFEM; this has additional capabilities,
including generalised RVE boundary conditions, hierarchic finite elements and

24
optimisation for high-performance computing. The developed computational
framework provides the nonlinear micro/meso-mechanical response at lamina
level, which will be used subsequently to simulate FRP composites at both
laminate and structure level.

Acknowledgements

• The authors gratefully acknowledge the support of the UK Engineering

and Physical Sciences Research Council through the Providing Confidence
in Durable Composites (DURACOMP) project (Grant Ref.: EP/K026925/1).

• The first author also acknowledges A. R. Melro (University of Bristol, UK)

for providing the MATLAB code for generation of randomly distributed
fibres within the RVE for the UD FRP composites.

References

References

[1] L. Tong, A. P. Mouritz, M. K. Bannister, 3D Fibre Reinforced Polymer

Composites, Elsevier Science Ltd, The Boulevard, Langford Lane, Kidling-
ton Oxford, UK, 2002.

[2] A. Mouritz, M. Bannister, P. Falzon, K. Leong, Review of applications for

advanced three-dimensional fibre textile composites, Composites Part A:
Applied Science and Manufacturing 30 (12) (1999) 1445 – 1461.

[3] Z. Ullah, L. Kaczmarczyk, S. A. Grammatikos, M. C. Evernden, C. J.

Pearce., Multi-scale computational homogenisation to predict the long-
term durability of composite structures, Computers and Structures (Under
Review).

[4] S. Nemat-Nasser, M. Hori, Micromechanics: Overall Properties of Hetero-

geneous Materials, Elsevier Science Publishers, North Holland, Amster-
dam, 1993.

25
 [5] M. Geers, V. Kouznetsova, W. Brekelmans, Multi-scale computational ho-
mogenization: Trends and challenges, Journal of Computational and Ap-
plied Mathematics 234 (7) (2010) 2175 – 2182.

[6] X.-Y. Zhou, P. D. Gosling, C. J. Pearce, Z. Ullah, L. Kaczmarczyk,

Perturbation-based stochastic multi-scale computational homogenization
method for woven textile composites, International Journal of Solids and
Structures 80 (2016) 368 – 380.

[7] Z. Ullah, L. Kaczmarczyk, M. Cortis, C. J. Pearce, Multiscale modelling

of the textile composite materials, in: 22nd UK Conference of the Associ-
ation for Computational Mechanics in Engineering (ACME), University of
Exeter, Exeter, UK, 2014, pp. 214–217.

[8] Z. Ullah, L. Kaczmarczyk, C. J. Pearce, Multiscale computational ho-

mogenisation to predict the long-term durability of composite structures,
in: 23rd UK Conference of the Association for Computational Mechanics
in Engineering (ACME), University of Swansea, Swansea, UK, 2015, pp.
95–98.

[9] Z. Ullah, L. Kaczmarczyk, C. J. Pearce, Nonlinear micro-mechanical re-

sponse of the fibre-reinforced polymer composites including matrix damage
and fibre-matrix decohesion, in: 24th UK Conference of the Association for
Computational Mechanics in Engineering (ACME), University of Cardiff,
Cardiff, UK, 2016, pp. 256–259.

[10] C. González, J. LLorca, Mechanical behavior of unidirectional fiber-

reinforced polymers under transverse compression: Microscopic mecha-
nisms and modeling, Composites Science and Technology 67 (13) (2007)
2795–2806.

[11] F. Naya, C. Lopes, C. González, J. LLorca, Computational micromechanics

strategies for the analysis of failure in unidirectional composites, in: P. P.
Camanho, S. R. Hallett (Eds.), Numerical Modelling of Failure in Advanced

26
 Composite Materials, Woodhead Publishing Series in Composites Science
and Engineering, Woodhead Publishing, 2015, pp. 411 – 433.

[12] J. Brockenbrough, S. Suresh, H. Wienecke, Deformation of metal-matrix

composites with continuous fibers: geometrical effects of fiber distribution
and shape, Acta Metallurgica et Materialia 39 (5) (1991) 735 – 752.

[13] T. Vaughan, C. McCarthy, Micromechanical modelling of the transverse

damage behaviour in fibre reinforced composites, Composites Science and
Technology 71 (3) (2011) 388 – 396.

[14] T. Vaughan, C. McCarthy, A combined experimental–numerical approach

for generating statistically equivalent fibre distributions for high strength
laminated composite materials, Composites Science and Technology 70 (2)
(2010) 291 – 297.

[15] G. Soni, R. Singh, M. Mitra, B. G. Falzon, Modelling matrix damage and

fibre-matrix interfacial decohesion in composite laminates via a multi-fibre
multi-layer representative volume element (M2 RVE), International Journal
of Solids and Structures 51 (2) (2014) 449 – 461.

[16] Digimat Users’ Manual, Version 4.2.1 (2011).

[17] M. Romanowicz, A numerical approach for predicting the failure locus of

fiber reinforced composites under combined transverse compression and
axial tension, Computational Materials Science 51 (1) (2012) 7 – 12.

[18] N. K. Parambil, S. Gururaja, Micromechanical damage analysis in lami-

nated composites with randomly distributed fibers, Journal of Composite
Materials - (-) (2015) 1–14.

[19] A. R. Melro, P. P. Camanho, F. M. A. Pires, S. T. Pinho, Micromechanical

analysis of polymer composites reinforced by unidirectional fibres: Part
I – constitutive modelling, International Journal of Solids and Structures
50 (11–12) (2013) 1897 – 1905.

27
[20] A. R. Melro, P. P. Camanho, F. M. A. Pires, S. T. Pinho, Micromechanical
analysis of polymer composites reinforced by unidirectional fibres: Part II
– micromechanical analyses, International Journal of Solids and Structures
50 (11–12) (2013) 1906 – 1915.

[21] N. Tschoegl, Failure surfaces in principal stress space, Journal of polymer

science Part C: Polymer symposia 32 (1) (1971) 239–267.

[22] A. R. Melro, P. P. Camanho, S. T. Pinho, Generation of random distribu-

tion of fibres in long-fibre reinforced composites, Composites Science and
Technology 68 (9) (2008) 2092 – 2102.

[23] M. Matveev, A. Long, Numerical modelling for predicting failure in textile

composites, in: P. P. Camanho, S. R. Hallett (Eds.), Numerical Modelling
of Failure in Advanced Composite Materials, Woodhead Publishing Series
in Composites Science and Engineering, Woodhead Publishing, 2015, pp.
435 – 455.

[24] M. Kollegal, S. Chatterjee, G. Flanagan, Progressive failure analysis of

plain weaves using damage mechanics based constitutive laws, International
journal of damage mechanics 10 (4) (2001) 301–323.

[25] M. Zako, Y. Uetsuji, T. Kurashiki, Finite element analysis of damaged wo-

ven fabric composite materials, Composites Science and Technology 63 (3–
4) (2003) 507 – 516.

[26] D. Blackketter, D. Walrath, A. Hansen, Modeling damage in a plain weave

fabric-reinforced composite material, Journal of Composites, Technology
and Research 15 (2) (1993) 136–142.

[27] D. S. Ivanov, F. Baudry, B. V. D. Broucke, S. V. Lomov, H. Xie, I. Ver-

poest, Failure analysis of triaxial braided composite, Composites Science
and Technology 69 (9) (2009) 1372 – 1380.

28
[28] B. G. Falzon, P. Apruzzese, Numerical analysis of intralaminar failure
mechanisms in composite structures. Part I: FE implementation, Compos-
ite Structures 93 (2) (2011) 1039–1046.

[29] B. G. Falzon, P. Apruzzese, Numerical analysis of intralaminar failure

mechanisms in composite structures. Part II: Applications, Composite
Structures 93 (2) (2011) 1047–1053.

[30] B. Stier, J. W. Simon, S. Reese, Comparing experimental results to a nu-

merical meso-scale approach for woven fiber reinforced plastics, Composite
Structures 122 (2015) 553 – 560.

[31] Mesh Oriented Finite Element Method (MoFEM) Version 0.3.8, Univer-
sity of Glasgow, Glasgow, UK, [Link]
(2016).

[32] J. Whitcomb, X. Tang, Micromechanics of moisture diffusion in composites

with impermeable fibers, Journal of Composite Materials 36 (9) (2002)
1093–1101.

[33] T. Hobbiebrunken, M. Hojo, T. Adachi, C. De Jong, B. Fiedler, Evaluation

of interfacial strength in CF/epoxies using FEM and in-situ experiments,
Composites Part A: Applied Science and Manufacturing 37 (12) (2006)
2248–2256.

[34] A. R. Maligno, N. A. Warrior, A. C. Long, Finite element investigations

on the microstructure of fibre-reinforced composites, eXPRESS Polymer
Letters 2 (9) (2008) 665–676.

[35] L. Kaczmarczyk, C. J. Pearce, N. Bićanić, Scale transition and enforcement

of RVE boundary conditions in second-order computational homogeniza-
tion, International Journal for Numerical Methods in Engineering 74 (3)
(2008) 506–522.

[36] CUBIT, The Geometry and Mesh Generation Toolkit, Version 13.0, Sandia
National Laboratories, USA, [Link]

29
[37] ParaView , Open-source, multi-platform data analysis and visualization
application, Version 4.4.0, [Link]

[38] J. Segurado, J. LLorca, A new three-dimensional interface finite element to

simulate fracture in composites, International Journal of Solids and Struc-
tures 41 (11–12) (2004) 2977 – 2993.

[39] B. Fiedler, M. Hojo, S. Ochiai, K. Schulte, M. Ando, Failure behavior of an

epoxy matrix under different kinds of static loading, Composites Science
and Technology 61 (11) (2001) 1615 – 1624.

[40] V. P. Nguyen, O. Lloberas-Valls, M. Stroeven, L. J. Sluys, On the ex-

istence of representative volumes for softening quasi-brittle materials – a
failure zone averaging scheme, Computer Methods in Applied Mechanics
and Engineering 199 (45–48) (2010) 3028 – 3038.

[41] I. Gitman, H. Askes, L. Sluys, Representative volume: Existence and size

determination, Engineering Fracture Mechanics 74 (16) (2007) 2518 – 2534.

[42] Z. P. Bazant, Can multiscale-multiphysics methods predict softening dam-

age and structural failure?, International Journal for Multiscale Computa-
tional Engineering 8 (1) (2010) 61–67.

[43] V. P. Nguyen, M. Stroeven, L. J. Sluys, Multiscale continuous and discon-

tinuous modeling of heterogeneous materials: A review on recent develop-
ments, Journal of Multiscale Modelling 03 (04) (2011) 229–270.

30