Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

57, 1, pp. 17-26, Warsaw 2019
DOI: 10.15632/jtam-pl.57.1.17

MATHEMATICAL MODELLING AND SIMULATION OF DELAMINATION
CRACK GROWTH IN GLASS FIBER REINFORCED PLASTIC (GFRP)

COMPOSITE LAMINATES

Hassan Ijaz

University of Jeddah, Mechanical Engineering Department, Jeddah, Saudi Arabia

e-mail: hassan605@yahoo.com

Delamination crack growth is a major source of failure in composite laminates under static
and fatigue loading conditions. In the present study, damagemechanics based failuremodels
for both static and fatigue loadings are evaluated via UMAT subroutine to study the dela-
mination crack growth phenomenon in Glass Fiber Reinforced Plastic (GFRP) composite
laminates. A static local damage model proposed by Allix and Ladevèze is modified to an
non-local damagemodel in order to simulate the crackgrowthbehaviordue to static loading.
Next, the same classical damage model is modified to simulate fatigue delamination crack
growth.The finite element analysis results obtained by the proposedmodels are successfully
comparedwith the available experimental data on the delamination crack growth forGFRP
composite laminates.

Keywords: finite element analysis, GFRP, damage mechanics, non-local, fatigue, delami-
nation

1. Introduction

Composite laminates are frequentlyused inmodernstructuralmaterials dueto thehigh strength-
-to-weight ratio. Moreover, by adjusting the orientation of fibers one can also get desired me-
chanical properties in desired loading directions (Herakovich, 1997). Carbon and glass fibers are
commonly used to manufacture composite laminates. Carbon fibers have better strength and
less density than glass fibers, but they are not cost effective. Glass Fiber Reinforced Plastic
(GFRP) composite laminates are used in avionic, automobile, ship andwind turbine industries.
In the present study, delamination crack growth simulations for theGFRP composite laminates
are performed under static and fatigue loadings. Delamination may be defined as a crack like
an entity between the composite laminates. The cracks can grow within laminates under static
and fatigue loadings andmay result in failure of structural parts (Davies et al., 1989; Allix and
Ladevèze, 1992). Normally, damage or fracture mechanics based approaches are used to study
the cark growth behavior in different structural elements. Fracture mechanics deals with the
propagation of already existing crack (Meng andWang, 2014) while, on the other hand, damage
mechanics can not only simulate the propagation of cracks but also deals with initiation of the
crack (Allix et al., 1995, 1998; Allix and Ladevèze, 1996; Ijaz et al., 2016).
Damagemechanics based formulations have been used to simulate the crack growth behavior

in composite laminatesmostly forCFRP(Corigliano, 1993;Corigliano andAllix, 2000;Chaboche
et al., 1997; Alfano and Crisfield, 2001). In the present study, delamination crack growth in the
GFRP composite laminates is focused using the damage mechanics based formulation.
Classical static damagemodels proposedby earlier authorsweremostly local in nature (Allix

and Ladevèze, 1992; Corigliano, 1993; Chaboche et al., 1997; Alfano and Crisfield, 2001). Loca-
lization means that damage tends to localize in a narrow zone in front of the crack tip rather
than a uniform distribution over a certain region (Jirasek, 1998). Bažant and Pijaudier-Cabot



18 H. Ijaz

(1988, 1989) proposed an integral type non-local damage model for brittle concrete materials.
Similarly, a rate dependent damage model is also proposed to avoid the localization issues
in CFRP composite laminates by introducing a time delay in the damage evolution formula-
tion (Allix et al., 2000; Marguet et al., 2007). To counter the localization problem, Peerlings
introduced a gradient enhanced damage evolution model (Peerlings et al., 2001). Borino ga-
ve the idea of using the integral type non-local damage model for the interface damage mo-
dels for composite laminates (Borino et al., 2007). Ijaz used the idea of an integral type non-
-local interface damage model for the study of delamination crack growth in CFRP composite
laminates (Ijaz et al., 2014). GFRP composite laminates also show a considerable amount of
fiber bridging during crack growth (Davidson and Waas, 2012). In the present study, an inte-
gral type non-local damage is used to accommodate the spurious localization and fiber bridging
issues duringdelamination crack growth inGFRPcomposite laminates under static loading con-
ditions. The classical damage model proposed by Allix and Ladevèze (1992, 1996) is modified
to a non-local one.

This article is organized as follows: inSection 2, basics of the classical interface damagemodel
are recalled. The proposed non-local static interface damage model is discussed in Section 3.
Finite element simulation results and their comparison with the experimental data are detailed
in Section 4. Finally, some concluding remarks are given in Section 5.

2. Introduction to the classical local interface damage model

Simulation of delamination crack growth in composite laminates is performed by coupled in-
terface damage modelling. The interface is a crack like entity that exists between two adjacent
lamina layers. The relative displacement of the two adjacent layers with respect to each other
can be described as

U=U+−U− =U1N1+U2N2+U3N3 (2.1)

whereN1,N2 andN3 aremutual perpendicular vectors in an orthotropic reference frame for the
interface. The failure or deterioration of the interface is taken into account by the introduction
of three damage variables, d1, d2 and d3 correspond to orthotropic direction vectors. Here, d3
corresponds to the out-of-plane opening mode (Mode I), whereas d1 and d2 correspond to the
in-plane shearing and tearing failure modes (Mode II and Mode III). The damage variable is
divided into two parts, i.e. static damage variable diS and fatigue damage variable diF . Hence,
the total damage di can be calculated by taking the sum of the two aforementioned damage
variables di = diS +diF , i=1,2,3.

Ifσ13,σ23 andσ33 are interfacial stress components inN1,N2 andN3 directions, respectively,
then the damage variables are related to the interfacial displacements as







σ13
σ23
σ33






=







k01(1−d1) 0 0
0 k02(1−d2) 0
0 0 k03(1−d3)













U1
U2
U3






(2.2)

here k01, k
0
2 and k

0
3 are defined as interface rigidities corresponding to three failure modes. The

damage model is built by considering thermodynamic forces combined with damage variables
and are associated with three modes of delamination as follows (Allix et al., 1995; Allix and
Ladevèze, 1996)

Yd3 =
1

2

〈σ33〉
2
+

k03(1−d3)
2

Yd1 =
1

2

σ213
k01(1−d1)

2
Yd2 =

1

2

σ232
k02(1−d2)

2
(2.3)



Mathematical modelling and simulation of delamination crack growth... 19

where 〈σ33〉+ represents the positive value of σ33, i.e. damage will not grow during compression
loading when a normal loading is applied. Now the three damage variables are assumed to be
strongly coupled and are governed by a single equivalent damage energy release rate of the
following form (Allix and Ladevèze, 1996)

Y (t)=max
r¬t

(

(Yd3)
α+(γ1Y d1)

α+(γ2Y d2)
α
)

1

α

(2.4)

where γ1 and γ2 are coupling parameters, and α is a material parameter which governs the
damage evolution under mixed mode loading conditions. Now, the damage evolution law is
defined as an isotropic material function of the following form

if [(d3S < 1) and (Y <YR)] then d1S = d2S = d3S =ω(Y )

else d1S = d2S = d3S =1
(2.5)

where the damage evolution material function ω(Y ) is defined as (Allix and Ladevèze, 1996)

ω(Y )=
( n

n+1

〈Y −YO〉+
YC −YO

)n

(2.6)

where YO and YC are threshold and critical damage energy release rates. n is termed as a
characteristic function of thematerial. Higher values ofn correspond to a brittle interface. YR is
defined as damage energy associated to rupture and can be calculated using following formula

YR =YO+
n+1

n
d
1

n

c (YC −YO)

Now the identification of parametersYC,γ1 andγ2 can bedoneby comparing the energy dissipa-
tion yielded by the damagemechanics approach andLEFM(LinearElastic FractureMechanics).
For a puremode case energy release rateGiC (i= I,II,III) obtained from fracturemechanics,
the experiments considering LEFM can be compared to the damage mechanics approach using
the following relation (Allix et al., 1995; Allix and Ladevèze, 1996)

GIC =YC GIIC =
YC

γ1
GIIIC =

YC

γ2
(2.7)

For a mixed-mode loading case, a standard LEFM model can be recovered as (Allix and La-
devèze, 1996)

( GI

GIC

)α

+
( GII

GIIC

)α

+
( GIII

GIIIC

)α

=1 (2.8)

The equivalence between damage mechanics and LEFM also leads to the following relation
during the complete debonding process (DP)

GIC =

∫

DP

σ33 dU3 GIIC =

∫

DP

σ13 dU1 GIIIC =

∫

DP

σ23 dU2 (2.9)

Equation (2.9) states that for any puremode debonding case, the area under the curve obtained
from the damage mechanics approach is equal to the experimentally obtained critical energy
release rateGiC.



20 H. Ijaz

3. Non-local interface damage model for static loading

The classical interface damagemodel shows the strain softening phenomenon during the degra-
dation process (Borino et al., 2007; Ijaz et al., 2014). Due to this softening behaviour, the stress
tends to localize in a narrow region in front of the crack tip. This localization phenomenon is
more obvious for 3D delamination simulations over 2D analysis. Moreover, fibre bridging also
occurs during delamination crack growth in GFRP composite laminates (Yao, 2015; Davidson
andWaas, 2012). Themathematical non-local interface damagemodel presented in this Section
will also inherently accommodate the fibre bridging process.

In the proposedmethodology, the damage variable ismade non-local by taking spatial avera-
ging over a certain domain using theGaussian distributionmethodology. The averaging domain
is dictated by the characteristic lengthparameter l.Ahigher value of lmeans thatmore elements
are used for the averaging. This domain dependent characteristic length l will also simulate the
fibre bridging behaviour since bridging occurs over a certain domain during delamination crack
growth in composite laminates.

Now one can write the average of damage variable d over the surrounding domain using the
weight function α0(r) (Ijaz et al., 2014)

d(x)=
1

Vr(x)

∫

α0(‖x− ζ‖)d(ζ) dζ (3.1)

where

Vr(x)=

∫

α0(‖x− ζ‖) dζ (3.2)

The damage variable calculated using Eq. (2.5) is made non-local using Eq. (2.10) over a pre-
scribed selected domain. The isotropic weight function α0(r) is calculated using the Gaussian
distribution function of exponential form as follows

α0(r)= exp
(

−
r2

2l2

)

(3.3)

From the above equation it is clear that the weight function α0(r) depends on the distance
between two points r= ‖x− ζ‖ and the characteristic length parameter l. A smaller value of l
corresponds to a less number of elements available for averaging. Higher values mean thatmore
elements will take part in the for averaging.

Equation (2.10) describes the damage variable d which is made non-local by taking spatial
averaging. Similarly, likeEq. (2.10), twoother variables like the equivalent damage energy release
rate Y and interfacial displacementU are alsomade non-local, and then the damage variable is
calculated using these as

Y (x)=
1

Vr(x)

∫

α0(‖x− ζ‖)Y (ζ) dζ

U(x)=
1

Vr(x)

∫

α0(‖x− ζ‖)U(ζ) dζ

(3.4)

Once the values of Y (x) and U(x) are known then the non-local damage variable can be cal-
culated as a function of Y (x) and U(x), respectively. Once the value of the averaged non-local
variable d is calculated either using Eqs. (2.10) and (2.13)1 or (2.13)2 then the updated value
of dwill be used in Eq. (2.2) to calculate the interfacial stresses.



Mathematical modelling and simulation of delamination crack growth... 21

4. Finite element analysis and results

In this Section, details and results of finite element simulations are presented and discussed
separately for both static non-local and fatigue interface damage models for GFRP composite
laminates.DoubleCantileverBeam(DCB) specimen forMode I is shown inFig. 1.The specimen
has total length L, initial crack length a0 and total height 2h, as shown in Fig. 1.

Fig. 1. DCB specimen for Mode I delamination crack growth

All simulations are performed in finite element software Cast3M (CEA) (Verpeaux et al.,
2000). The geometry of the beam is modelled with 2D plane strain quadrangles. The interface
between the specimen arms is modelled with 2D interface element JOI2 to simulate the debon-
ding process (Beer, 1985). Different parameters like YO, YC, γ1, n, k

0
1 and k

0
3 are needed to be

identified for the finite element analysis ofMode I andMode II delamination crack growth. The
value of threshold damage energy YO is taken zero, i.e. YO = 0 for all the finite element simu-
lations. The identification of YC and γ1 can be done by using Eq. (2.7) provided that Mode I
and Mode II critical energy release rates GIC and GIIC are already determined from LEFM
experiments. The identification of value of n depends on the brittleness of the interface. The
value ofn varies between 0-1.0, and a good value can be identified bymatching the experimental
and numerical results. The values of interfacial rigidities can be calculated using the following
relation (Ijaz et al., 2011)

k03 =
(2n+1)

2n+1

n

8n(n+1)Yc
σ233 k

0
i
=
γi(2n+1)

2n+1

n

8n(n+1)Yc
σ23i i=1,2 (4.1)

In Eq. (4.1), σ33 and σ3i (i = 1,2) are the maximum interfacial normal and in-plane shear
stresses. The energy release rate calculated using fracture mechanics theory for pureMode I is
given as (Willams, 1988)

GI =
M2

bEI
(4.2)

whereM is the applied moment, b is width of the specimen,E is flexure modulus and I is the
secondmoment of area of the bear arm.

5. FE analysis of delamination crack growth under static loading

In this Section, finite element analysis ofMode I delamination crack growth forGFRPcomposite
laminates is performed using the non-local interface damage model. The experimental work of
Davidson andWaas (2012) onMode I delamination crack growth forGFRPcomposite laminates
is used for finite element analysis. Nominal dimensions forDCB specimen are:L=130, h=2.5,
a0 =50andwidth is b=25.4.All thedimensionsmentioned above forDCBspecimenare inmm.
theMode I critical energy release rate GIC value is 1.45KJ/m

2. The modulusE11, in the fibre
direction is 11.5GPa and themajor Poisson’s ratio is 0.3 forGFRP (Davidson andWaas, 2012).
In the interface damagemodel, different values are identified as discussed above. The interfacial
rigidity k03 is found to be 9000MPa/mm for maximum normal stress value of 50MPa. Figure 2



22 H. Ijaz

shows the evolution of normal stress with respect to the interfacial displacement for different
identified damage parameters. From Fig. 2, one can also observe that the area under the curve
is always equal to the critical energy release rate for the complete debonding process, Eq. (2.9).

Fig. 2. Evolution of normal stress σ33 with interfacial displacementU3

Three different formulations on averaging of the variables d, Y and U are discussed earlier
for the non-local interface damage model. These three different non-local interface damage for-
mulations are implemented in finite element software Cast3M via procedure PERSO1 and user
material subroutine UMAT.

Figure 3 shows the normalized value of evolution of the reaction force with crack opening
displacement forMode I delamination crack growthusingdbasednon-local damage formulation.
Reaction force values are normalized to 120N for ease of presentation and comparison with the
experimental results. Figure 3 presents the curves for the non-local model for three different
characteristic lengths l with values of 0.2, 0.4 and 0.6.

Fig. 3. (a) Normalized reaction force and (b) crack growth vs crack opening displacement for d based
non-local formulation

From Fig. 3a, one can notice that as the characteristic length value reduces, the results
coincide with the classical local damage model proposed by Allix and Ladevèze. The reason
is that as the characteristic length value reduces, a fewer number of elements is available for
the averaging and the results come close to those predicted by the local damage model. The
finite element simulation results are also in good agreement with the experimental results of
(Davidson andWaas, 2012). Peng andXu (2013) proposed a damagemodel that accommodates



Mathematical modelling and simulation of delamination crack growth... 23

thebridging effect bydividing thedamage variable into staticdS andbridgingdb parts.However,
in the present work, the non-local model accommodates the bridging effect by controlling the
characteristic length l.ThemodelofPengandXu(2013) alsodemonstrated the similarbehaviour
of the increasing reaction forcewith an increase in the bridging force, whereas in the present case
a larger value of the characteristic length causes an increase in the reaction force, see Fig. 3a.
From Fig. 3a, one can also note that although the peak reaction force is different for various
characteristic length values, but they all tend to converge close to each other once a stable crack
growth is established.

Figure 3b presents the evolution of delamination crack growth as a function of the crack ope-
ning displacement. Figure shows the evolution of the crack for the non-local model using three
different characteristic length values, i.e. 0.2, 0.4 and 0.6. In Fig. 3, the finite element analysis
results are also compared with the experimental results. The start of crack growth correspon-
ding to different l values is in accordance with the behaviour depicted in Fig. 3. The start of
crack growth for the non-local damagemodel with a small value of l (0.2) and the classical local
damagemodel is almost the same. Similarly, crack growth starts late for a larger l (0.6) value in
comparison to smaller l (0.2,0.4) values, and this behaviour is also in accordance with the one
predicted in Fig. 3. The late start of crack growth is due to a relatively large value of l, which
means thatmore elements are available for the averaging in the non-local zone and, thus, impar-
ting extra resistance to crack growth resulting in an increase in the reaction force and, hence,
the late start of crack growth. The actual fibre bridging behaviour during delamination crack
growth in GFRP composite laminates is simulated by introducing the characteristic length l
into the classical local damagemodel. The larger value of l indicates awider fibre bridging zone,
whichmeans extra resistance to crack growth and a higher value of the reaction force.Moreover,
there is a good agreement between numerical and experimental results for delamination crack
growth with crack opening displacement, Fig. 3b.

Figure 4 shows the evolution of the reaction force and crack growth with crack opening
displacement using theU based non-local damage model. Trend of the reaction force and crack
growth for the U based non-local damage model is very similar to that predicted by d based
formulation. In the case of d based non-local formulation, the reaction force and crack growth
converge very close to each otherwhile forU based formulationbothalmost converge to the same
path for different values of l once the stable crack growth rate establishes. The results presented
in Fig. 4 also show a similar trend as predicted by Peng andXu (2013) for the reaction force in
the GFRP composite laminate delamination crack growth for different bridging forces.

Fig. 4. (a) Normalized reaction force and (b) crack growth vs. crack opening displacement forU based
non-local formulation



24 H. Ijaz

Similarly,Fig. 5 shows theevolutionof the reaction force andcrackgrowthwith crackopening
displacement using theY based non-local damagemodel.The results predicted by this non-local
damage model show erratic behaviour. The results start to deviate from the experimental ones
as the characteristic length l value increases for both the reaction force and crack growth, see
Fig. 5.

Fig. 5. (a) Normalized reaction force and (b) crack growth vs crack opening displacement for Y based
non-local formulation

From the above discussions, it can be inferred that d andU based non-local damage formu-
lations give reasonable results for delamination crack growths in GFRP composite laminates.
The proposed methodology will not only avoid the localization issue but will also compensate
the fibre bridging effect.

Fig. 6. Evolution of the damage variable with crack opening displacement: (a) local damagemodel,
(b) non-local damagemodel, l=0.6

Figure 6 shows the evolution of the damage variablewith crack opening displacement for the
first ten elements fromthe crack tip for local andU basednon-local damagemodels. In thefigure
one can observe that for the first five elements, from the crack tip, the damage evolution starts
as soon the load is applied in the local damagemodel, whereas in the non-local damagemodel,
the damage evolution for the first seven elements starts as soon the load is applied. Similarly, for
the eighth element, the damage evolution starts when the crack opening displacement is 17mm
in the local damage model, and it starts around 8mm in the non-local damage model. The
characteristic length l helps one to involve more elements in the averaging process and makes
more elements be a part of damage initiation in the non-local model in comparison to the local
damagemodel. Hence, one can say that for larger values of l, more elements will be available in



Mathematical modelling and simulation of delamination crack growth... 25

the averaging process andwill take part in initial damage growth as soon as the load is applied.
Since in the local damage model no averaging of damage variable is done over a certain area,
therefore, fewer elements will take part in the initial damage growth in comparison with the
non-local damage model, see Fig. 6.

6. Conclusion

Delamination crack growth in GFRP composite laminates under static loading is simulated
in the present study. The original classical local damage proposed by Allix and Ladevèze is
modified to a non-local static damagemodel and a fatigue damagemodel. The proposedmodels
are implemented via UMAT subroutine and PERSO1 procedure in Cast3M FE software. The
non-local damage model is not only capable of avoiding the accumulation of damage in front of
the crack tip butalso compensates thefibrebridgingphenomenonat the interface by introducing
characteristic length l. The reaction force and crack growth are plotted against the crack opening
displacement for different values of l in pureMode I static loading. Three different formulations
based on averaging of variables d, Y and U are proposed. FE simulation results show that d
and U based non-local models predict good results. The results determined by finite element
analysis for the static loading is successfully compared with the available experimental data of
GFRP composite laminates.

References

1. AlfanoG., CrisfieldM.A., 2001, Finite element interfacemodels for the delamination analysis
of laminated composites:mechanical and computational issue, International Journal for Numerical
Methods in Engineering, 50, 1701-1736

2. AllixO., LadevèzeP., 1992, Interlaminar interfacemodelling for thepredictionofdelamination,
Composite Structures, 22, 235-242

3. AllixO.,LadevèzeP., 1996,Damagemechanics of interfacialmedia:basic aspects, identification
and application to delamination, Studies in Applied Mechanics, 44, 167-188

4. Allix O., Ladevèze P., Corigliano A., 1995, Damage analysis of interlaminar fracture speci-
mens,Composite Structures, 31, 61-74

5. AllixO.,LadevèzeP.,DeuJ.F.,LévêqueD., 2000,Amesomodel for localizationanddamage
computation in laminates,ComputerMethods in AppliedMechanics andEngineering,183, 105-122

6. Allix O., Ladevèze P., Gornet L., Perret L., 1998, A computational damage mechanics
approach for laminates: identification and comparisonwith experimental result, Studies in Applied
Mechanics, 46, 481-500

7. Bažant Z.P., Pijaudier-Cabot G., 1988, Nonlocal continuum damage, localization instability
and convergence, Journal of Applied Mechanics, 55, 287-293

8. Bažant Z.P., Pijaudier-Cabot G., 1989, Measurement of characteristic length of nonlocal
continuum, Journal of Engineering Mechanics, 115, 755-767

9. Beer G., 1985, An isoparametric joint/interface element for finite element analysis, International
Journal for Numerical Methods in Engineering, 21, 585-600

10. BorinoG.,FaillaB.ParrinelloF., 2007,Nonlocal elasticdamage interfacemechanicalmodel,
International Journal for Multiscale Computational Engineering, 5, 153-165

11. Chaboche J.L., Girard R., Levasseur P., 1997, On the interface debondingmodels, Interna-
tional Journal of Damage Mechanics, 6, 220-257



26 H. Ijaz

12. Corigliano A., 1993, Formulation, identification and use of interface models in the numerical
analysis of composite delamination, International Journal Solids and Structures, 30, 2779-2811

13. Corigliano A., Allix O., 2000, Some aspects of interlaminar degradation in composites,Com-
puter Methods in Applied Mechanics and Engineering, 185, 203-224

14. Davidson P., Waas A.M., 2012, Non-smoothMode I fracture of fibre-reinforced composites: an
experimental, numerical and analytical study, Philosophical Transactions of the Royal Society of
London, Series A, 370, 1942-1965

15. DaviesP.,CantwellW.,MoulinC.,Kausch,H.H., 1989,A studyofdelamination resistance
of IM6/PEEK composites,Composites Science and Technology, 36, 153-166

16. HerakovichC.T., 1997,Mechanics of Fibrous Composites, 1st ed., JohnWiley&SonsLtd., New
York, U.S.

17. Ijaz H., Asad M., Gornet L., Alam S.Y., 2014, Prediction of delamination crack growth
in carbon/fiber epoxy composite laminates using non-local interface damage model, Mechanics &
Industry, 15, 293-300.

18. Ijaz H., Khan M.A., Saleem W., Chaudry S.R., 2011, Numerical modeling and simulation
of delamination crack growth in cf/epoxy composite laminates under cyclic loading using cohesive
zonemodel,Advanced Materials Research, 326, 37-52

19. Ijaz H., Zain-ul-AbdeinM., Saleem W., AsadM., Mabrouki T., 2016,A numerical appro-
achonparametric sensitivityanalysis foranaeronauticaluminiumalloyturningprocess,Mechanics,
2, 149-155

20. Jirasek M., 1998, Nonlocal models for damage and fracture: comparison of approaches, Interna-
tional Journal of Solids and Structures, 35, 4133-4145

21. Marguet S., Rozycki P., Gornet L., 2007, A rate dependent constitutive model for carbon-
-fiber reinforcedplasticwoven fabric,Mechanics of AdvancedMaterials and Structures,14, 619-631

22. Meng Q., Wang Z., 2014, Extended finite element method for power-law creep crack growth,
Engineering Fracture Mechanics, 127, 148-160

23. Peerlings R.H.J., Geers M.G.D., De Borst R., Brekelmans W.A.M., 2001, A critical
comparison of nonlocal and gradient enhanced softening continua, International Journal of Solids
Structures, 38, 7723-7746

24. Peng L., Xu J., 2013, Fatigue delamination growth of composite laminates with fiber bridging:
Theory and simulation, [In:]Proceeding of 13th International Conference on Fracture (ICF13), Yu
S. and Feng X.-Q. (Edit.), Beijing, China, 16-21

25. Verpeaux p., CharrasT., MillardA., 1998,Castem 2000: UneApprocheModerne du Calcul
des Structures, Fouet J.M., Ladevèze P., OhayonR. (Edit.), 227-261

26. Williams J.G., 1988, On the calculation of energy release rates for cracked laminates, Interna-
tional Journal of Fracture, 36, 101-119

27. Yao L., 2015, Mode I fatigue delamination growth in composite laminates with fibre bridging,
PhDDissertation, Technische Universiteit Delft, Germany

Manuscript received April 13, 2017; accepted for print May 4, 2018