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

Glass fibre-reinforced plastic (GFRP) composite laminates are used in many industries due to their excellent mechanical and thermal properties. However, these materials are prone to the initiation and propagation of delamination crack growth between different plies forming the laminate. The crack propagation may ultimately result in the failure of GFRP laminates as structural parts. In this research, a comprehensive mathematical model is presented to study the delamination crack growth in GFRP composite laminates under fatigue loading. A classical static damage model proposed by Allix and Ladevèze is modified as a fatigue damage model. Subsequently, the model is implemented in commercial finite element software via UMAT subroutine. The results obtained by the finite element simulations verify the experimental findings of Kenane and Benzeggagh for the fatigue crack growth in GFRP composite laminates.


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 [1].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 [4], while the damage mechanics deals with the propagation of the cracks as well as the simulation of the crack initiation process [5][6][7][8].
Many authors have also proposed different mathematical formulations to analyze the delamination crack growth in composite laminates under static loadings [9][10][11][12].Martin and Murri [13], Paas et al. [14], Asp et al. [15], and Blanco et al. [16] performed the fatigue-driven delamination crack growth experiments.Robinson et al. [17], Tumino and Cappello [18], Turon et al. [19], and Ijaz et al. [20] 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 [21].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 Ladevèze 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 1 , d 2 , and d 3 ), which correspond to the three orthotropic directions.Here, d 1 represents mode I delamination crack growth, while d 2 and d 3 present mode II and mode III crack propagations, respectively.The damage variable d i is divided into the static (d iS ) and fatigue components (d iF ).Hence, the total damage (d i ) can be calculated by taking the sum of the two aforementioned damage variables ( 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][23][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.

Static Damage Model
The relative displacement of the laminate layers corresponds to the three mutually perpendicular vectors (N 1 , N 2 , and N 3 ) in an orthotropic reference frame and can be expressed in the following form: The interfacial stresses corresponding to the three damage variables may be expressed as In the above equation, k 0 1 , k 0 2 , and k 0 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

3
where σ 33 + 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 [6]: Here, γ 1 and γ 1 are the coupling parameters and α 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]: where the damage evolution material function ω Y is defined as [6] Here, Y O and Y C are the threshold and critical damage energy release rates, and n is termed as the characteristic function of material.Y R is defined as the damage energy associated with rupture and can be calculated by Identification of different parameters Y C , γ 1 , and γ 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]:

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. [14] and Peerlings et al. [21].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 [17].
In the present work, the damage evolution is governed by a single equivalent damage energy release rate Y t 2 International Journal of Aerospace Engineering (see ( 4)); hence, the same is used to predict the fatigue damage, as follows: Here, f is a loading function and is defined as f = Y − Y * .Y * 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 t and the total damage d.This parameter is expressed as Here, C, λ, and β 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 Δt [20]:

International Journal of Aerospace Engineering
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 ΔN.Hence, (10) takes the form 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: International Journal of Aerospace Engineering Now, the total damage due to the fatigue loading can be calculated by summing (11), ( 12) and ( 13) as

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) [25].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 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 [26].
Different parameters like Y O , Y C , γ 1 , n, k 0 1 , and k 0 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 O is taken as zero in all simulations.If the values of mode I and mode II critical energy release rates (G IC and G IIC ) are already known from LEFM experiments, then Y C , γ 1 , and γ 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 [20]: In (15), σ 33 and σ 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 [29]: 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 [29] G II = 3 4 The proposed fatigue damage model is able to produce the linear crack growth rate as obtained by the classical Paris law [23,24]: In the above equation, B and m are the material parameters and are determined from different fatigue tests (ΔG = G max − G min ).Here, G max and G min correspond to the maximum and minimum energy release rates during the cyclic variation of load.G 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 [22] 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 0 = 25 and width is b = 20.Similarly, for the ELS specimen, dimensions are L = 65, h = 3 0, a 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 2 and 2.9 kJ/m 2 , respectively [22].
The modulus, E11, used in the fibre direction is 36.2GPa, and the major Poisson ratio is 0.33 for the GFRP composite laminate [22].The normal and in-plane shear stiffnesses, k 0 3 and k 0 1 , are found to be 9000 MPa/mm and 1500 MPa/mm, respectively, by using (15) for the maximum stress value of  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 (λ, β, 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 max = 0 35 k J/m 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 max = 0 07 kJ/m 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 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 max = 0 3 kJ/m 2 (Figure 8).On the other hand, for G max = 0 1 kJ/m 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 max = 0 3 kJ/m 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.

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

Figure 3 :
Figure 3: Evolution of stress with displacement (a) normal σ 33 and displacement U 3 and (b) shear σ 13 and displacement U 1 .

Table 1 :
Fracture toughness and associated fatigue parameters.The evolution of the normal and shear stresses with the interfacial displacement of identified damage parameters is shown in Figure3.

Figure 8 :
Figure 8: Damage variable evolution due to mode I fatigue loading: maximum fatigue load corresponds to 0.1 kJ/m 2 .
0E + 00 1.0E + 04 2.0E + 04 3.0E + 04 4.0E + 04 5.0E + 04 6.0E + 04 7.0E + 04 8.0E + 04Figure 6: Crack growth rate versus number of fatigue loading cycles (mode I).Figure 7: Damage variable evolution due to mode I fatigue loading: maximum fatigue load corresponds to 0.3 kJ/m 2 .6InternationalJournal of Aerospace Engineering 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.