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

Authors: Hassan Ijaz; Waqas Saleem; Muhammad Zain-ul-abdein; Aqeel Ahmad Taimoor; Abdullah Salmeen Bin Mahfouz
Journal: International Journal of Aerospace Engineering
 (2018)

---

## Abstract

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.


---

## Body

## 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 [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–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 [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 (d1, d2, and d3), which correspond to the three orthotropic directions. Here, d1 represents mode I delamination crack growth, while d2 and d3 present mode II and mode III crack propagations, respectively. The damage variable di is divided into the static (diS) and fatigue components (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).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 Section2, 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 (N1, N2, and N3) in an orthotropic reference frame and can be expressed in the following form:
(1)U=U=U+−U−=U1N1+U2N2+U3N3.The interfacial stresses corresponding to the three damage variables may be expressed as(2)σ13σ23σ33=k101−d1000k201−d2000k301−d3U1U2U3.In the above equation,k10, k20, and k30 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)Yd3=12σ33+2k301−d32,Yd1=12σ132k101−d12,Yd2=12σ322k201−d22,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]:
(4)Y¯t=Maxτ≤tYd3α+γ1Yd1α+γ2Yd2α1/α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]:
(5)ifd3S<1andY¯<YRthend1S=d2S=d3S=ωY¯elsed1S=d2S=d3S=1,where the damage evolution material function ωY¯ is defined as [6]
(6)ωY¯=nn+1Y¯−YO+YC−YOn.Here,YO and YC are the threshold and critical damage energy release rates, and n is termed as the characteristic function of material. YR is defined as the damage energy associated with rupture and can be calculated by YR=YO+n+1/ndc1/nYC−YO.Identification of different parametersYC, γ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]:
(7)GIC=YC,GIIC=YCγ1,GIIIC=YCγ2.


## 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 Figure1. 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.Figure 1
Envelope of the applied cyclic load.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 rateY¯t (see (4)); hence, the same is used to predict the fatigue damage, as follows:
(8)d•F=∂dF∂t=gd,Y¯YCY•¯tYCifY•¯≥0andf≥00ifY•¯<0orf<0.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
(9)gd,Y¯YC=CeλdY¯YCβ.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]:
(10)dFt+Δt=dFt+∫tt+ΔtdF•dt=dFt+∑n=NN+ΔNΡd,Y¯.In the above equation,t corresponds to the number of cycles N and Δt corresponds to the number of increment cycles ΔN. Ρd,Y¯ is the evolution of damage due to the fatigue per cycle:
(11)Ρd,Y¯=C1+βeλdY¯YC1+β.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
(12)dFN+ΔN=dFN+12PdN+ΔN,Y¯N+ΔN+PdN,Y¯NΔN.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:(13)dSN+ΔN=dSN+nn+11YC−YOnY¯N+ΔN−YO+n−Y¯N−YO+nifY¯N+ΔN≥Y¯N.Now, the total damage due to the fatigue loading can be calculated by summing (11), (12) and (13) as
(14)dN+ΔN=dFN+ΔN+dSN+ΔN.


## 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) [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 a0, and the total height 2h.Figure 2
Specimen nomenclature and boundary conditions.
(a)
GFRP composite specimen nomenclature
(b)
DCB specimen (mode I)
(c)
ELS specimen (mode II)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 likeYO, YC, γ1, n, k10, and k30 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 YO is taken as zero in all simulations. If the values of mode I and mode II critical energy release rates (GIC and GIIC) are already known from LEFM experiments, then YC, γ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]:
(15)k30=2n+12n+1/n8nn+1Ycσ332,ki0=γi2n+12n+1/n8nn+1Ycσ3i2,withi∈1,2.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]:
(16)GI=M2bEI,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]
(17)GII=34M2bEI.The proposed fatigue damage model is able to produce the linear crack growth rate as obtained by the classical Paris law [23, 24]:
(18)dadN=BΔGGCm.In the above equation,B and m are the material parameters and are determined from different fatigue tests (ΔG=Gmax−Gmin). Here, Gmax and Gmin correspond to the maximum and minimum energy release rates during the cyclic variation of load. GC 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 a0=25 and width is b=20. Similarly, for the ELS specimen, dimensions are L=65, h=3.0, a0=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/m2 and 2.9 kJ/m2, respectively [22].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 [22]. The normal and in-plane shear stiffnesses, k30 and k10, 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 3
Evolution of stress with displacement (a) normalσ33 and displacement U3 and (b) shear σ13 and displacement U1.
(a)
(b)Figure4 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 4
Fatigue delamination crack growth (mode I).Figure 5
Fatigue delamination crack growth (mode II).Table 1
Fracture toughness and associated fatigue parameters.
Test method
GC (kJ/m2) [23, 24]
Interface
Damage fatigue parameters
Mode I
0.429
n=0.5,YO=0kJ/m2
λ=3.0
Yc=0.429kJ/m2
β=1.0×10−4
k03=9.3×103MPa/mm
C=4.0×10−5
Mode II
2.9
n=0.5,YO=0kJ/m2
λ=3.0
Yc=2.9kJ/m2
β=5.0
k01=1.5×103MPa/mm
C=75.0
γ1=0.148Figure6 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, Gmax=0.35kJ/m2, the crack growth rate is also high. Similarly, when the fatigue crack growth rate is slow, the low amplitude values are applied, that is, Gmax=0.07kJ/m2. 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.Figure 6
Crack growth rate versus number of fatigue loading cycles (mode I).Figures7 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/m2 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, Gmax=0.3kJ/m2 (Figure 8). On the other hand, for Gmax=0.1kJ/m2 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 Gmax=0.3kJ/m2 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.Figure 7
Damage variable evolution due to mode I fatigue loading: maximum fatigue load corresponds to 0.3 kJ/m2.Figure 8
Damage variable evolution due to mode I fatigue loading: maximum fatigue load corresponds to 0.1 kJ/m2.


## 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.



---

*Source: 2081785-2018-03-28.xml*