# Fatigue Delamination Crack Growth in GFRP Composite Laminates: Mathematical Modelling and FE Simulation.

1. Introduction

Composite laminates are frequently used in the modern automobile, aerospace, and manufacturing industries due to their excellent mechanical properties and lightweight characteristics. Moreover, the composite properties can be tailored in the desired directions by adjusting the stacking sequences and the directions of different plies . The list of composite is quite diversified, but carbon, Kevlar, and glass fibres are most common in manufacturing of different structural parts. Among these, the glass fibre composite laminates are the most economical for structural parts, like the windmills, UAV bodies, and the filament tubes/cylinders.

Delamination is considered as a crack like entity between any two plies that can initiate and propagate in the composite laminates under different loading conditions [2, 3], and the situation may become severe since many structural parts may fail in the real-life applications under the cyclic fatigue loading. Many authors proposed different mathematical models based on the damage and the fracture mechanics theories to simulate the delamination crack growth in the glass and the carbon fibre composite laminates. Fracture mechanics can predict the propagation of a crack that already exists in the structural parts , while the damage mechanics deals with the propagation of the cracks as well as the simulation of the crack initiation process [5-8].

Many authors have also proposed different mathematical formulations to analyze the delamination crack growth in composite laminates under static loadings [9-12]. Martin and Murri , Paas et al. , Asp et al. , and Blanco et al.  performed the fatigue-driven delamination crack growth experiments. Robinson et al. , Tumino and Cappello , Turon et al. , and Ijaz et al.  developed the models to simulate the delamination crack growth using finite element (FE) analysis. Peerlings et al. proposed a fatigue damage model for the analysis of crack growth in metals . Most of the proposed simulation work in literature is related with the carbon fibre-reinforced plastic (CFRP) composites. In the present work, however, the glass fibre-reinforced plastic (GFRP) composite laminates under the fatigue loading conditions are examined. A static damage model proposed by Allix and Ladeveze is modified to accommodate the fatigue-driven delamination in the present study.

A delamination crack is, generally, represented by an interfacial entity between two plies of the laminates. Failure of this interface is dictated by the three damage variables ([d.sub.1], [d.sub.2], and [d.sub.3]), which correspond to the three orthotropic directions. Here, [d.sub.1] represents mode I delamination crack growth, while [d.sub.2] and [d.sub.3] present mode II and mode III crack propagations, respectively. The damage variable [d.sub.i] is divided into the static ([d.sub.iS]) and fatigue components ([d.sub.iF]). Hence, the total damage ([d.sub.i]) can be calculated by taking the sum of the two aforementioned damage variables ([d.sub.i] = [d.sub.iS] + [d.sub.iF], i = 1, 2, 3).

In this study, the FE simulation results for mode I and mode II fatigue delamination crack growth are compared with the experimental findings of Kenane [22-24].

This article is organized as follows: the classical static damage model is presented in Section 2, the proposed fatigue damage model is detailed in Section 3, details of finite element analysis and its comparison with the experimental results are discussed in Section 4, and finally, the concluding remarks are given in Section 5.

2. Static Damage Model

The relative displacement of the laminate layers corresponds to the three mutually perpendicular vectors ([N.sub.1], [N.sub.2], and [N.sub.3]) in an orthotropic reference frame and can be expressed in the following form:

U = [U] = [U.sup.+] - [U.sup.-] = [[U.sub.1]][N.sub.1] + [[U.sub.2]][N.sub.2] + [[U.sub.3]][N.sub.3] (1)

The interfacial stresses corresponding to the three damage variables may be expressed as

[mathematical expression not reproducible]. (2)

In the above equation, [k.sup.0.sub.1], [k.sup.0.sub.2], and [k.sub.0.sup.3] are the corresponding interfacial stiffnesses. The damage model is derived by considering the thermodynamic forces. These are combined with the three damage variables [5, 6] as

[mathematical expression not reproducible]. (3)

where [<[[sigma].sub.33]>.sub.+] represents that damage will grow only in the tensile loading (opening mode I) and will not grow during the compression. The above three damage variables are combined to form a single equivalent damage energy release rate for the mixed mode loading :

[mathematical expression not reproducible] (4)

Here, [[gamma].sub.1] and [[gamma].sub.1] are the coupling parameters and [alpha] is the material parameter which governs the damage evolution under the mixed mode loading conditions. For static loads, the damage evolution law is defined in the following form [5, 6]:

[mathematical expression not reproducible], (5)

where the damage evolution material function [omega]([Y.bar]) is defined as 

[omega]([Y.bar]) = [[n/n+1 [<[Y.bar] - [Y.sub.O>.sub.+]/[Y.sub.C] - [Y.sub.O]].sup.n]. (6)

Here, [Y.sub.O] and [Y.sub.C] are the threshold and critical damage energy release rates, and n is termed as the characteristic function of material. [Y.sub.R] is defined as the damage energy associated with rupture and can be calculated by [Y.sub.R] = [Y.sub.O] + (n + 1/n) [d.sup.1/n.sub.c] ([Y.sub.C]-[Y.sub.O]).

Identification of different parameters [Y.sub.C], [[gamma].sub.1], and [[gamma].sub.1] is carried out by comparing the energy dissipation yielded by the linear elastic fracture mechanics (LEFM) and the damage mechanics approaches for different failure modes [5, 6]:

[mathematical expression not reproducible] (7)

3. Fatigue Damage Model

In this section, the mathematical formulations and assumptions used in the fatigue damage model are discussed. The actual cyclic fatigue load varies between a maximum and minimum value as shown in Figure 1. But, for a high cyclic fatigue, the actual numerical load applied in FE analysis corresponds to the maximum value of the cyclic loading; see Figure 1.

The original idea of the fatigue damage is adopted from Paas et al.  and Peerlings et al. . Peerlings proposed a strain-based damage model for cyclic loading for the crack propagation in metals. Robinson adopted the idea of Peerlings for CFRP composite laminates and used for the normalized interfacial displacement .

In the present work, the damage evolution is governed by a single equivalent damage energy release rate [Y.bar](t) (see (4)); hence, the same is used to predict the fatigue damage, as follows:

[mathematical expression not reproducible]. (8)

Here, f is a loading function and is defined as f = [Y.bar]-[Y.sup.*]. [Y.sup.*] is a threshold damage energy release rate. Here, g is a dimensionless parameter that governs the fatigue damage. This parameter depends on the equivalent damage energy [Y.bar](t) and the total damage d. This parameter is expressed as

g(d, [Y.bar]/[Y.sub.C]) = [Ce.sup.[lambda]d] [([Y.bar]/[Y.sub.C]).sup.[beta]]. (9)

Here, C, [lambda], and [beta] are the fatigue damage material parameters and will be identified by comparing the experimental results with the FE results. Equation (8) is in a rate form; therefore, it is integrated over time to get the increase in damage with respect to the time increment [DELTA]t :

[mathematical expression not reproducible]. (10)

In the above equation, t corresponds to the number of cycles N and [DELTA]t corresponds to the number of increment cycles [DELTA]N. P(d, [Y.bar]) is the evolution of damage due to the fatigue per cycle:

P(d, [Y.bar]) = C/1 + [beta] [[e.sup.[lambda]d] [[Y.bar]/[Y.sub.C]].sup.1+[beta]] (11)

The sum over the number of cycles presented in (10) is done by means of the trapezoidal rule for definite integrals. The trapezoidal rule estimates the desired value by taking the average of the initial and final values. At the end of the increment, it is multiplied by [DELTA]N. Hence, (10) takes the form

[mathematical expression not reproducible]. (12)

As stated earlier, for the high-cycle fatigue, the total damage is a sum of the fatigue and static damage components. The following relation expresses the evolution of the static damage over a certain number of cycles:

[mathematical expression not reproducible]. (13)

Now, the total damage due to the fatigue loading can be calculated by summing (11), (12) and (13) as

d(N + [DELTA]N = [d.sub.F] (N + [DELTA]N) + [d.sub.S](N + [DELTA]N). (14)

4. Finite Element Analysis and Results

Mode I and mode II delamination crack growth simulations of GFRP composite laminates are performed in FE software CAST3M (CEA) . The proposed fatigue damage model described in Section 3 is implemented via UMAT subroutine and used for FE analysis. Double cantilever beam (DCB) and end load split (ELS) specimens are modelled for mode I and mode II delamination crack growths, respectively. This is shown in Figure 2. The specimen has a total length L, an initial crack length [a.sub.0], and the total height 2h.

The geometry of beam is modelled with 2D plane strain quadrangles. To simulate the debonding process, the interface between the specimen arms is modelled with an interface element JOI2 .

Different parameters like [Y.sub.O], [Y.sub.C], [[gamma].sub.1], n, [k.sup.0.sub.1], and [k.sup.0.sub.3] need to be identified for the FE analysis of mode I and mode II crack growths [27, 28]. Note that the value of threshold damage energy [Y.sub.O] is taken as zero in all simulations. If the values of mode I and mode II critical energy release rates ([G.sub.IC] and [G.sub.IIC]) are already known from LEFM experiments, then [Y.sub.C], [[gamma].sub.1], and [[gamma].sub.2] can be identified using (7). The value of n varies between 0 and 1.0, and a good value can be estimated by comparing the experimental and numerical results. The values of interfacial stiffnesses can be calculated by using the following relation :

[mathematical expression not reproducible]. (15)

In (15), [[sigma].sub.33] and [[sigma].sub.3i] are the maximum interfacial normal and in-plane shear stresses. The energy release rate is calculated by using the fracture mechanics theory for mode I :

[G.sub.I] = [M.sub.2]/bEI, (16)

where M is the applied moment, b is the width of the specimen, E is the flexure modulus, and I is the second moment of area of bear arm. Similarly, the energy release rate for pure mode II can be expressed as 

GII = 3/4 [M.sup.2]/bEI. (17)

The proposed fatigue damage model is able to produce the linear crack growth rate as obtained by the classical Paris law [23, 24]:

da/dN = B[[DELTA]G/[G.sub.C].sup.m]. (18)

In the above equation, B and m are the material parameters and are determined from different fatigue tests ([DELTA]G = [G.sub.max]-[G.sub.min]). Here, [G.sub.max] and [G.sub.min] correspond to the maximum and minimum energy release rates during the cyclic variation of load. [G.sub.C] is the fracture toughness of the material, which is determined through the static delamination crack growth experiments. In this study, the results of the experimental work of Kenane and Benzeggagh  on the fatigue delamination growth of GFRP composite laminates is used for comparison with the FE simulation results. Nominal dimensions for DCB specimen are L = 150, h = 3.0, and [a.sub.0] = 25 and width is b = 20. Similarly, for the ELS specimen, dimensions are L = 65, h = 3.0, [a.sub.0] = 25, and b = 20. All the dimensions mentioned above are in millimeters for both types of specimens. Mode I and mode II critical energy release rate values are 0.429 kJ/[m.sup.2] and 2.9 kJ/[m.sup.2], respectively .

The modulus, E11, used in the fibre direction is 36.2 GPa, and the major Poisson ratio is 0.33 for the GFRP composite laminate . The normal and in-plane shear stiffnesses, [k.sup.0.sub.3] and [k.sup.0.sub.1], are found to be 9000 MPa/mm and 1500 MPa/mm, respectively, by using (15) for the maximum stress value of 40 MPa. The evolution of the normal and shear stresses with the interfacial displacement of identified damage parameters is shown in Figure 3.

Figure 4 shows the linear Paris plot behaviour for mode I delamination crack growth due to the fatigue loading. The FE results obtained, using the proposed fatigue model, are found to be in good agreement with the experimental results [23, 24]. The suitable values of fatigue parameters ([lambda], [beta], and C) of the proposed model are selected to give the best fit of the experimental results. Similarly, Figure 5 presents the Paris plot behaviour for mode II delamination crack growth under the cyclic loading. The linear line of Paris plot behaviour is in good agreement with experimental results. The identified fatigue parameters of GFRP composite laminates for a cyclic loading are given in Table 1.

Figure 6 shows the crack growth rates with a number of cycles for different energy release rates, which correspond to the maximum amplitude values of the fatigue loading. The methodology of applied loading is explained earlier; see Figure 2. From Figure 6, it can be observed that when applied fatigue loading amplitude is very high, that is, [G.sub.max] = 0.35 k J/[m.sup.2], the crack growth rate is also high. Similarly, when the fatigue crack growth rate is slow, the low amplitude values are applied, that is, [G.sub.max] = 0.07 kJ/[m.sup.2]. The same phenomenon of the large crack growth rates with the larger loading amplitudes is also depicted in Paris plot; see Figures 4 and 5.

Figures 7 and 8 present the damage evolution with respect to the number of cycles for the first ten elements from the crack tip for maximum loading amplitudes corresponding to 0.3 and 0.1 kJ/[m.sup.2] energy release rates for mode I fatigue delamination crack growth. From the figures, one can observe that damage initiates and finishes rapidly for high amplitude loading, that is, [G.sub.max] = 0.3 kJ/[m.sup.2] (Figure 8). On the other hand, for [G.sub.max] = 0.1 kJ/[m.sup.2] loading, damage for the first element completes at a point that corresponds to 1.0E5 cycles, and this point corresponds to the complete damage of 10th element for [G.sub.max] = 0.3 kJ/[m.sup.2] loading. Figures 7 and 8 depict the reason of high delamination crack growth rates due to rapid initiation and completeness of the damage variables for high amplitude fatigue loadings.

5. Conclusion

In this study, delamination crack growth simulations for the GFRP composite laminates under the fatigue loading are presented. Details of the proposed mathematical model are also explained. The model is implemented in CAST3M software via UMAT subroutine. The crack growth rates obtained from FE analysis are plotted against the energy release rates to obtain the linear Paris plot behaviour for mode I and mode II load cases. The simulation results are compared with the available experimental data of GFRP composite laminate and were found to be in good approximation. FE analysis results predicted the large crack growth rates at high amplitude for the applied loading values.

https://doi.org/10.1155/2018/2081785

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

References

 C. T. Herakovich, Mechanics of Fibrous Composites, John Wiley and sons Ltd, New York, NY, USA, 1st edition, 1997.

 P. Davies, W. Cantwell, C. Moulin, and H. H. Kausch, "A study of the delamination resistance of IM6/PEEK composites," Composites Science and Technology, vol. 36, no. 2, pp. 153-166, 1989.

 O. Allix and P. Ladeveze, "Interlaminar interface modelling for the prediction of delamination," Composite Structures, vol. 22, no. 4, pp. 235-242, 1992.

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

 O. Allix, P. Ladeveze, and A. Corigliano, "Damage analysis of interlaminar fracture specimens," Composite Structures, vol. 31, no. 1, pp. 61-74, 1995.

 O. Allix and P. Ladeveze, "Damage mechanics of interfacial media: basic aspects, identification and application to delamination," Studies in Applied Mechanics, vol. 44, pp. 167-188, 1996.

 P. Ladeveze, O. Allix, L. Gornet, D. Leveque, and L. Perret, "A computational damage mechanics approach for laminates: identification and comparison with experimental results," Studies in Applied Mechanics, vol. 46, pp. 481-500, 1998.

 H. Ijaz, M. Zain-ul-abdein, W. Saleem, M. Asad, and T. Mabrouki, "A numerical approach on parametric sensitivity analysis for an aeronautic aluminium alloy turning process," Mechanics, vol. 22, no. 2, pp. 149-155, 2016.

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

 A. Corigliano and O. Allix, "Some aspects of interlaminar degradation in composites," Computer Methods in Applied Mechanics and Engineering, vol. 185, no. 2-4, pp. 203-224, 2000.

 J. L. Chaboche, R. Girard, and P. Levasseur, "On the interface debonding models," International Journal of Damage Mechanics, vol. 6, no. 3, pp. 220-257, 1997.

 G. Alfano and M. A. Crisfield, "Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues," International Journal for Numerical Methods in Engineering, vol. 50, no. 7, pp. 1701-1736, 2001.

 R. H. Martin and G. B. Murri, "Characterization of mode I and mode II delamination growth and thresholds in AS4/PEEK composites," in Composite Materials: Testing and Design ASTM STP, Volume 1059, pp. 251-270, Sparks, NV, USA, January 1990.

 M. H. J. W. Paas, P. J. G. Schreurs, and W. A. M. Brekelmans, "A continuum approach to brittle and fatigue damage: theory and numerical procedures," International Journal of Solids and Structures, vol. 30, no. 4, pp. 579-599, 1993.

 L. E. Asp, A. Sjogren, and E. S. Greenhalgh, "Delamination growth and thresholds in a carbon/epoxy composite under fatigue loading," Journal of Composites Technology and Research, vol. 23, pp. 55-68, 2001.

 N. Blanco, E. K. Gamstedt, L. E. Asp, and J. Costa, "Mixed-mode delamination growth in carbon-fibre composite laminates under cyclic loading," International Journal of Solids and Structures, vol. 41, no. 15, pp. 4219-4235, 2004.

 P. Robinson, U. Galvanetto, D. Tumino, G. Bellucci, and D. Violeau, "Numerical simulation of fatigue-driven delamination using interface elements," International Journal for Numerical Methods in Engineering, vol. 63, no. 13, pp. 1824-1848, 2005.

 D. Tumino and F. Cappello, "Simulation of fatigue delamination growth in composites with different mode mixtures," Journal of Composite Materials, vol. 41, no. 20, pp. 2415-2441, 2007.

 A. Turon, J. Costa, P. P. Camanho, and C. G. Davila, "Simulation of delamination in composites under high-cycle fatigue," Composites Part A: Applied Science and Manufacturing, vol. 38, no. 11, pp. 2270-2282, 2007.

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

 R. H. J. Peerlings, W. A. M. Brekelmans, R. de Borst, and M. G. D. Geers, "Gradient-enhanced damage modelling of high cycle fatigue," International Journal for Numerical Methods in Engineering, vol. 49, no. 12, pp. 1547-1569, 2000.

 M. Kenane and M. L. Benzeggagh, "Mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites under fatigue loading," Composites Science and Technology, vol. 57, no. 5, pp. 597-605, 1997.

 M. Kenane, "Delamination growth in unidirectional glass/epoxy composite under static and fatigue loads," Physics Procedia, vol. 2, no. 3, pp. 1195-1203, 2009.

 M. Kenane, S. Benmedakhene, and Z. Azari, "Fracture and fatigue study of unidirectional glass/epoxy laminate under different mode of loading," Fatigue & Fracture of Engineering Materials & Structures, vol. 33, pp. 284-293, 2010.

 P. Verpeaux, T. Charras, and A. Millard, Castem 2000: Une approche moderne du calcul des structures, J. M. Fouet, P. Ladeveze, and R. Ohayon, Eds., pp. 227-261, 1998, http://www-cast3m.cea.fr.

 G. Bfer, "An isoparametric joint/interface element for finite element analysis," International Journal for Numerical Methods in Engineering, vol. 21, no. 4, pp. 585-600, 1985.

 A. Corigliano and S. Mariani, "Parameter identification of a time-dependent elastic-damage interface model for the simulation of debonding in composites," Composites Science and Technology, vol. 61, no. 2, pp. 191-203, 2001.

 A. Corigliano, S. Mariani, and A. Pandolfi, "Numerical modeling of rate-dependent debonding processes in composites," Composite Structures, vol. 61, no. 1-2, pp. 39-50, 2003.

 J. G. Williams, "On the calculation of energy release rates for cracked laminates," International Journal of Fracture, vol. 36, no. 2, pp. 101-119, 1988.

Hassan Ijaz (iD), (1) Waqas Saleem (iD), (1) Muhammad Zain-ul-abdein, (1) Aqeel Ahmad Taimoor, (2) and Abdullah Salmeen Bin Mahfouz (3)

(1) Department of Mechanical Engineering, University of Jeddah, Jeddah, Saudi Arabia

(2) Department of Chemical and Materials Engineering, King Abdulaziz University, Jeddah, Saudi Arabia

(3) Department of Chemical and Materials Engineering, University of Jeddah, Jeddah, Saudi Arabia

Correspondence should be addressed to Hassan Ijaz; hassan605@yahoo.com

Received 29 September 2017; Revised 9 January 2018; Accepted 22 January 2018; Published 28 March 2018

Caption: Figure 1: Envelope of the applied cyclic load.

Caption: Figure 2: Specimen nomenclature and boundary conditions.

Caption: Figure 3: Evolution of stress with displacement (a) normal [[sigma].sub.33] and displacement [U.sub.3] and (b) shear [[sigma].sub.13] and displacement [U.sub.1].

Caption: Figure 4: Fatigue delamination crack growth (mode I).

Caption: Figure 5: Fatigue delamination crack growth (mode II).

Caption: Figure 6: Crack growth rate versus number of fatigue loading cycles (mode I).

Caption: Figure 7: Damage variable evolution due to mode I fatigue loading: maximum fatigue load corresponds to 0.3 kJ/[m.sup.2].

Caption: Figure 8: Damage variable evolution due to mode I fatigue loading: maximum fatigue load corresponds to 0.1 kJ/[m.sup.2].
```Table 1: Fracture toughness and associated fatigue parameters.

Test        [G.sub.C]             Interface          Damage fatigue
method    (kJ/[m.sup.2])                               parameters
[23, 24]

Mode I        0.429         n = 0.5, [Y.sub.O] =     [lambda] = 3.0
0 kJ/[m.sup.2]

[Y.sub.c] = 0.429      [beta] = 1.0 x
kJ/[m.sup.2]           [10.sup.-4]

[k.sup.3.sub.0] = 9.3       C = 4.0 x
x [10.sup.3] MPa/mm       [10.sup.-5]

Mode II        2.9          n = 0.5, [Y.sub.O] =     [lambda] = 3.0
0 kJ/[m.sup.2]

[Y.sub.c] = 2.9        [beta] = 5.0
kJ/[m.sup.2]

[k.sup.1.sub.0] = 1.5       C = 75.0
x [10.sup.3] MPa/mm

[[gamma].sub.1] = 0.148
```
Title Annotation: Printer friendly Cite/link Email Feedback Research Article; Glass fibre-reinforced plastic; finite element Ijaz, Hassan; Saleem, Waqas; Zain-ul-abdein, Muhammad; Taimoor, Aqeel Ahmad; Mahfouz, Abdullah Salme International Journal of Aerospace Engineering Report 1USA Jan 1, 2018 3845 Effective Mechanical Property Estimation of Composite Solid Propellants Based on VCFEM. Analysis of Global Sensitivity of Landing Variables on Landing Loads and Extreme Values of the Loads in Carrier-Based Aircrafts. Cracking (Materials) Fatigue (Materials) Fatigue testing machines Finite element method Laminated materials Laminated plastics Laminates Materials Mathematical models Plastic laminates