Review and Assessment on Fatigue Delamination Damage of Laminated Composite Structures

1. Introduction

Composite materials have many advantages such as high specific strength and strong designability. With the development of preparation and detection technology of composite material, the application of composite materials on civil aircraft is widely used. The composite materials which are applied in B757 and B767 account for 4%, which is 11% in B777 and even 50% in B787. Furthermore, the percent of composite materials in A350 has amazingly reached at 52%. In pursuit of lightweight for aircraft, a trend using more and more advanced composite materials is formed. Thus, the application of composite materials has become one of the important indicators to measure the performance of aircrafts. Although composite materials have many advantages, there are also many problems which are required to be solved in the servicing process.

Delamination is one of the most important damage modes of laminated composite structures. Delamination is very dangerous because it is visually undetectable within the laminated structures. Once delamination occurs, the stiffness and strength of the structure will greatly degrade, especially in compression mode. While the composite laminated structures of aircraft are subjected to the complex fatigue loading conditions, the damage propagation process of fatigue delamination is very difficult to be detected and monitored. Moreover, the damage of fatigue delamination is irreversible, which often leads to a serious decline in the performance of components and even flight accidents. Therefore, the researches on initiation and propagation law [1,2], the failure mechanism [3,4] and the numerical simulation [5,6] of delamination under fatigue loads are the significant issues to the damage tolerance and durability design of composite laminated structures.

As we all know, Fiber Reinforced Plastic (FRP) is often stacked into laminates by different directions laminas which has a single fiber direction, so the mechanical properties of laminates can be flexibly designed. However, the failure mechanism and failure modes of laminated composite structures largely differ from those of the metal structures. When dealing with angle-ply laminates, the mismatch of elastic modulus between the adjacent layers can result in the interlaminar stresses at the free edge of laminates. Furthermore, the interlaminar performance of FRP composite laminates is relatively weak, so that the larger interlaminar stress may lead to the failure mode of delamination. Even on the low stress level, it will cause fatigue delamination damage (FDD) as the performance of the interlaminar material degrades. In addition, the problem of FDD is more serious and critical in engineering structures with openings. Therefore, it is desirable to find a universal model of fatigue delamination damage for different composite laminated structures, which can simply and efficiently evaluate span fatigue life of delamination.

By reviewing a large number of literatures, the fatigue delamination damage model (FDDM) of the composite laminated structures can be divided into the strength model based on stress/strain field analysis [7,8], the fracture mechanics model [9,10], the damage mechanics model [11,12,13] and the hybrid model [14,15,16,17] according to the theory framework. In general, the performance of the FDDM can be verified according to the results of experiment including the double cantilever beam (DCB), the single lap joint (SLJ), the free-edge delamination and so on. At present, the Paris law [18] in metals based on fracture mechanics is commonly used to characterize the fatigue delamination growth (FDG) in the FDDM of composite materials. Furthermore, the phenomenological model is also established by utilizing the data from experiments. Considering about the complexity of composite laminated structures and loading conditions, various ways have been used to create and express a phenomenological model. Based on the influence factors of fatigue delamination failure, researchers have proposed a series of empirical and semi-empirical correction formulas [19,20,21], which provided methods of damage evolution for the simulation of FDDM. In recent years, the physical model, cohesive element based on damage mechanics [14,22,23] has received extensive attention and research. The cohesive element possesses the simple principle and high efficiency in numerical calculation, which greatly reduces the difficulty in simulation of FDG. In addition, some models in conjunction with the strength criterion with the fatigue progressive damage theory [24] have also been applied to estimate damage evolution of fatigue delamination, which shows a good agreement with experiments. Furthermore, investigators have already studied different hybrid frameworks of theories. A lot of valuable explorations have been made in the FDDM. The performances of these models have been verified through experiments and numerical simulations, and some practical conclusions have been also obtained. However, there are a number of problems without breakthrough.

Therefore, this article is to be divided into 5 sections. Section 2 introduces the characterization methods and significant influence factors of fatigue delamination. In Section 3, the general formula and procedure of the typical phenomenological models are summarized, including delamination onset and delamination propagation under fatigue loadings. Based on the framework of theories, Section 4 classifies the existing numerical simulations of fatigue delamination. Furthermore, the procedure, application and limitations of various FDDM have been compared in details. Finally, Section 5 sums up the up-to-date development of the FDDM, and gives suggestions on the research trends of the fatigue delamination and technical challenges.

2. Experiment and Mechanism of Fatigue Delamination

In order to investigate the delamination of composite laminated structures under fatigue loadings, many researchers have carried out numerous experimental studies based on the quasi-static delamination. In order to quantitatively characterize the damage of fatigue delamination, the relationship between the crack growth rate da/dN and the fatigue loading is established by combining the fatigue crack growth of metal with the experimental phenomenon of fatigue delamination. However, due to various of parameters and expressions of fatigue loadings, it is necessary to make a profound investigation of the influence factors on the crack propagation of fatigue delamination. In conclusion, researches on the experimental phenomenon and delamination mechanism of composite laminated structures under fatigue loadings have laid the foundation of the subsequent phenomenological model based on the experiment.

2.1. Delamination Test

According to the deformation characteristics of crack tip, the delamination crack growth of composite laminated structures can be divided into the open mode (mode I), the sliding mode (mode II) and the tearing mode (mode III) under the framework of fracture mechanics, as are shown in Figure 1.

2.1.1. Pure Mode Delamination Test

In order to obtain the fracture toughness of mode I [25] delamination propagation, ASTM specially formulated DCB test [26], with respect to unidirectional FRP composite laminates, as shown in Figure 2 (a). Because of the simple procedure, intuitive methods and fewer interference, DCB test is one of the most commonly used tests for researchers to verify the proposed FDDM.

Although ASTM has put forward a three-point bend end notched flexure (3ENF) test [27] for mode II, as shown in Figure 2 (b). Nevertheless, delamination crack growth is very unstable and the initial strain energy release rate (SERR) can only be obtained. Furthermore, the complete resistance curve (R-curve) cannot be obtained, which has significant impact on input parameters of the FDDM. Therefore, the test methods of mode II are still incomplete. In addition to 3ENF, other test methods such as the center notched flexure (CNF) [28], the over notched flexure (ONF) [29], the four-point bend end notched flexure (4ENF) [30] and the end load split (ELS) [31] also have different limitations. When verifying and analyzing the performance of the FDDM, it is required to select an appropriate test method.

There is no recognized standard experiment for mode III delamination test up to now. The experimental results obtained by edge crack torsion (ECT) [32] and split cantilever beam (SCB) [33] are sensitive to the geometric parameters of structures. In addition, the results are affected by mixed-mode of delamination, which result in incorrect measurement. Because the delamination test of mode III is still an open investigation subject, it is conservatively assumed that the delamination mechanism of mode III is similar to mode II. However, the fracture toughness of mode III is slightly higher than mode II.

2.1.2. Mixed-mode Delamination Test

In general, composite laminated structures are in a mixed-mode loading. Thus, the test methods of I/II mixed-mode delamination have been developed, including single leg bending (SLB) [34], fixed ratio mixed mode (FRMM) [35], cracked lap shear (CLS) [36], asymmetric double cantilever beam (ADCB) [37] and mixed-mode bending (MMB) [38]. It requires to be noted that the specimen form of MMB is the same as DCB and 3ENF. Moreover, MMB has developed a specific analytical solution, Therefore, it is included in ASTM [39], as shown in Figure 3.

2.1.3. Other Delamination Test

In engineering, the stacking sequence and ply orientation of composite laminated structures are complex. Thus, researchers often carry out fatigue delamination tests on the following specific composite laminated structures in order to verify the universality of FDDM. Firstly, fatigue delamination tests are studied by using the single-lap joints (SLJ) [40], the double-lapped joint (DLJ) [41] and the inclined and wedge joint. Second, fatigue delamination test of open-hole composite laminated structures [42] with different diameters and locations is also explored. Third, fatigue delamination test for angle-ply composite laminates [43] is further investigated. Figure 4 shows the finite element model of these three types of laminated structures, which are finished independently by authors.

2.2. Characterization for Fatigue Delamination

Based on the above tests, the delamination mechanism of composite laminated structures is gradually revealed. By summarizing the experimental phenomena, the researchers have obtained the characterization methods for fatigue delamination of FRP, which is similar to fatigue crack growth of metal. However, the preparation process of composite laminated structures is complex and different from metal. Furthermore, various initial damage modes such as voids, wrinkles and delamination are unavoidable. Composite laminated structures have not the stage of nucleation of cracks which is different from that of metal cracks. Therefore, the review article focuses on the description of experimental phenomena. Moreover, it is significant to characterize fatigue delamination onset and propagation of composite laminated structures under fatigue loads.

2.2.1. Fatigue Delamination Onset

First of all, considering the mode I fatigue delamination onset of FRP composites [44]. The onset of crack growth is determined by monitoring crack tip or flexibility of DCB specimens. The purpose of the test is to relate the maximum SERR G_Imax to the cycle of fatigue delamination onset N_onset . What’s more, G-N curve which is similar to S-N curve commonly used in fatigue of metal can be established. The critical fracture toughness G_IC is equivalent to the static strength σ_s , and the fatigue SERR threshold G_th is equivalent to the fatigue strength limit σ_f . A typical fatigue delamination growth onset curve [44] can be obtained as shown in Figure 5.

2.2.2. Fatigue Delamination Propagation

During the test of FDG, the length of delamination a and cycle of loading N are under a whole-process monitor. The data of experiment is utilized to fitting the relationship between the crack growth rate and the fatigue load range. Furthermore, the expression of fatigue delamination propagation adopts the general formula of classical Paris law like Eq.(1). Therefore, the typical curve of normalized fatigue delamination crack growth (da/dN-∆K/K_c) [18] is shown in Figure 6.

$$da/dN=C\times {[f(G \mathrm{or} K)]}^m$$

(1)

Same as fatigue delamination onset, the delamination of mode I is selected as an example. As shown in Figure 6, the fatigue crack growth rate develops from the high-rate region III to the low-rate region I, which is following the direction of arrow. As the cracks propagate along the fiber direction, the loading range of crack tip and the fatigue crack growth rate decrease gradually. It is noted that the loading range can be described by different parameters, such as the stress amplitude ∆σ, stress intensity factor (SIF) range ∆K and SERR range ∆G.

In medium-rate region II (Paris law region), the crack growth rate is linearly logarithmic with the loading range. Furthermore, propagation of crack stops when reaching the lower limit of fatigue delamination load range ∆K_th. The fatigue crack growth rate of FRP composites is measured according to ASTM standard of metal crack growth [45]. It is worth noting that the slope of fatigue crack growth rate curve of most FRP composites is much larger than that of most metals [46]. Therefore, it is indicated that little disturbance of load range can produce large errors in fatigue crack growth rate. In other words, the variation of fatigue crack growth rate is more sensitive than metals.

2.3. Influencing Factors of Fatigue Delamination

The relationship between fatigue crack growth rate and loading range is usually used to characterize the propagation of fatigue delamination. Thus, the parameters involved in loading range [47,48,49,50] and different forms of load range [51] may directly influence the FDG. These main influence factors include the mixed-mode ratio, the fiber bridging, the load ratio, the temperature, the humidity, the load sequence and the expression of f (K or G) in Eq.(1). What’s more, there are several secondary factors such as loading frequency and wave form.

2.3.1. Mixed-Mode Ratio

According to the definition of the mixed-mode ratio in terms of energy, the mixed-mode ratio is ϕ=G_II/G (G=G_I+G_II). In general, G_IIC is obviously greater than G_IC in FRP composites. As shown in Figure 7 (a) and (c), ϕ ranges from 0 to 1 with the gradual transition from mode I to mode II. In addition, the fracture toughness of mixed-mode evolves from G_IC to G_IIC. It can be found that upper limits of G_max-N curves and log(da/dN)-log(G_max) curves changed with different ϕ. However, lower limits of SERR remains unchanged. In other words, mixed-mode ratio directly affects the delamination onset and the fatigue crack growth rate [52]. Nevertheless, the influence of mixed-mode ratio can be ignored under low fatigue crack growth rate. These rules have been confirmed by investigation on experiment of Obrien [53].

In order to more accurately describe the relationship between delamination onset /fatigue crack growth rate and mixed-mode ratio, researchers [54,55] proposed a general curve fitting model based on a series of fatigue delamination tests of FRP composites. It is the most important that the fitting coefficient is a non-monotonic function of ϕ, which makes the tracking process of delamination onset and propagation more precise. As shown in Figure 7 (b) and (d), the change of maximum SERR G_max and fatigue crack growth rate da/dN is divided into two stages. In the first stage, G_max and da/dN monotonically change with the increase of mixed-mode ratio. In the second stage, G_max and da/dN non-monotonically change along the arrow direction with the increase of mixed-mode ratio.

2.3.2. Fiber Bridging effect

In existing test standards, most test objects of fracture toughness and fatigue crack growth are unidirectional FRP composite laminates. The influence of fiber bridging [56] on the critical SERR G_C is artificially introduced, which can obtain a complete R-curve [57,58]. From the initial stage of fiber bridging region [59], the critical SERR of delamination gradually increases and eventually tends to a stable value. Research of Zhang et al. [50] shows that the fiber bridging effect will lead to more complex characterization of delamination onset and propagation. It means that the lower limit G_th of curves will no longer be independent of the mixed-mode ratio.

As shown in Figure 8 (a) and (b), Yao et al. [60] carried out the test of mode I fatigue delamination and fitted curves of G_IC and G_th with length of delamination. It can be found that fracture toughness of mode I G_IC and lower limit SERR G_th can be constructed as piecewise functions of fatigue crack length a-a₀. G_IC and G_th gradually develop to a stable stage, which have a significant promotion when compared with the initial values. If the influence of fiber bridging is neglected, and the stable value of interlaminar parameters is taken for numerical simulation, it will lead to a conservative estimation of the fatigue delamination life of FRP composites.

2.3.3. Loading Ratio

When studying fatigue of metals, load ratio R has been always the focus of studies, and load ratio R has influence on the fatigue delamination of composite materials. If two of three loading parameters G_max, G_min and R are known, the form of fatigue loading can be determined. Therefore, while G_max is fixed, the effect of load ratio on delamination onset and propagation is investigated. As shown in Figure 9 (a) and (c), the experimental study of Matsubara et al. [61] indicated that the higher load ratio is, the higher the lower limit of the delamination onset curve is, which increases the number of cycles of delamination onset. Meanwhile, the larger the average load is, the greater the slope of the delamination propagation curve is, which is increasing the fatigue crack growth rate’s sensitivity to disturbance of G_max.

In order to simultaneously consider the effects of mixed-mode ratio and load ratio except fiber bridging, Allegri et al. [62] attempted to establish a relationship of coupling influence curve in form of Figure 9 (b) and (d). At the beginning of fatigue load, the higher the mixed-mode ratio, the greater the maximum SERR. With the evolution of fatigue delamination, the higher the loading ratio, the greater the maximum SERR under the same fatigue crack growth rate. Therefore, when considering the coupling effects of mixed- mode ratio and loading ratio, the contribution of both to fatigue delamination onset and propagation must be reasonably allocated. It is significant to obtain a phenomenological model matching through a series experiment.

2.3.4. Other Major Factors

Due to long-term service in hygrothermal conditions, the delamination resistance of composite laminates is weakened. The tests of Springer [63], Jin et al. [21] and Ramirez et al. [64] showed that the polymer matrix will be modified under high temperature and humidity especially for FRP composites. Therefore, matrix and inter-laminar properties of composite laminates will be directly affected. Hence, it is necessary to take temperature and humidity as parameters which influences the fatigue crack growth rate in special environments. Based on experiments, the quantitative expression of fatigue delamination behavior is obtained, with corresponding material system and loading cases.

There is no consensus on mechanism of the effect of fatigue load sequence on delamination propagation of composite laminates. On the basic of micromechanics, Li tried to study the interfacial fatigue debonding of unidirectional ceramic matrix composites under different loading sequences. In addition, Van Paepem and Degrieck [49] reviewed analytical methods of high-low fatigue load and low-high fatigue load of FRP composites. Thus, some issues and ideas about studying the influences of load sequence on fatigue delamination have been proposed.

When Paris Law is extended to composite materials, the problem of stress singularity at the crack tip is inevitable. Therefore, in the framework of linear elastic fracture mechanics (LEFM), the SERR G is established from the perspective of energy. In addition, the relationship between the SERR G and the SIF K is derived, which realize the substitution of load range. After fitting the data obtained from the test of fatigue crack growth rate, there are various function forms such as ∆G, G_max, $\sqrt{G_{\max }}$,${\left(\sqrt{G_{\max }}-\sqrt{G_{\max }}\right)}^2$and so on [66,67,68], which have completely different intrinsic physical meanings.

2.3.4. Other Secondary Factors

In general, loading frequency [49,69] and loading rate [70] are secondary factors on strain energy release rate. As long as the frequency below 100Hz, the limitation of frequency ensures that high temperature will not locally occur in composite laminates, and fatigue crack growth rate will be not affected. The influence of load waveform is much smaller [69], so researchers rarely investigated on it.

3. Phenomenological Model of Fatigue Delamination

By Referring to section 2, the analytical method of fatigue delamination is established by qualitative analysis of observation phenomenon in experiments. Furthermore, the important factors affecting fatigue delamination onset and propagation are clarified. Therefore, calculation of SERR is briefly summarized in this section. Then, a number of phenomenological models are introduced from the different perspectives of influencing factor and data processing.

3.1. Calculation of SERR

Two methods are usually employed to solve SERR, which are analytical solutions and finite element method (FEM). For analytical solutions, please refer to the related standards of ASTM [26,27,39], which is including modified beam theory (MBT), modified compliance calibration (MCC) and so on. The area method can also be used to estimate SERR [71,72]. The FEM includes the virtual crack closure technique (VCCT), the J-integral method, the virtual crack extension and so on. Among them, VCCT and J-integral method are most widely used.

3.1.1. VCCT

VCCT has been rapidly developed since Rybicki and Kanninen [73] proposed it in 1977. The basic assumption of VCCT [74] is that the opening displacement of the virtual crack tip is approximately equal to that of the initial crack tip. In two-dimension (2D) finite element model shown in Figure 10, SERR of pure mode I and mode II of delamination are calculated by Eq.(2).

$$\begin{array}{l}{G}_I=\frac{1}{2B\Delta a}{F}_y\Delta v\hfill \\ G_{II}=\frac{1}{2B\Delta a}{F}_x\Delta u\hfill \end{array}$$

(2)

where, B is the crack thickness, ∆a is the crack length, F_x and F_y denote the force of node e and node f along the x and y directions respectively, ∆u and ∆v denote nodal displacement difference between the node c and node d along the x and y directions respectively.

However, the path of damage evolution of VCCT needs to be set up in advance. Unless the prefabricated laminated specimens which are used for testing, the path of delamination evolution cannot be determined in structures, much less to complex geometric characteristics and load conditions. Moreover, VCCT can only track the propagation of delamination but cannot predict onset of delamination.

3.1.2. J-Integral

In 1968, Rice and Rosengren [75] proposed the J-integral for 2D problems. The stress and displacement far-field solved by line integral, which are used to describe the mechanical properties of the crack tip. Since the phenomenological models involved in this review are all based on theoretical framework of LEFM, J-integral and SERR are completely equivalent for problems of plane stress or plane strain. Therefore, SERR can be easily calculated by line integral, which is independent with the path from the lower surface to the upper surface of the crack.

3.2. Delamination Onset Phenomenological Model

The simplest functional relationship like Eq.(3), which is based on Basquin formula to fit data from delamination onset test of FRP composites.

$$G_{max}=A{N}_{\mathrm{onset}}^{-m}$$

(3)

where m and A are fitting parameters, N_onset is the cycles of fatigue delamination onset at a certain loading level. However, the model does not consider the main influencing parameters such as mixed-mode ratio, loading ratio and temperature. Any change of structure, environment and load condition will result in different curve fitting results. Thus, easy formula like Eq.(3) has poor universality in practice.

Jagannahan et al. [76] introduced the concept of constant life diagram (CLD) which is commonly used in fatigue of metals. Furthermore, the constant onset life diagram (COLD) of mode II fatigue delamination of carbon fiber composite was established. After fitting the test data of three loading ratios R = 0, R = 0.5 and R = 1 (that is pure mode II) by Basquin formula, the curve between amplitude of SERR and the average SERR (G_a -G_m curve) was achieved. The COLD is obtained by connecting points of constant life. Then, cycles of fatigue delamination onset under the corresponding loading level can be obtained by linear interpolation.

However, the latest research on mode I and mode II fatigue delamination onset of carbon fiber reinforced polymer (CFRP) materials in hygrothermal conditions, Ramirez et al. [64] still adopted the model same as Eq.(3). It is obvious that there is no general model which can comprehensively consider influencing factors of fatigue delamination onset. There are few phenomenological models for fatigue delamination onset in composite materials, and there is a lack of experimental verification.

3.3. Delamination Propagation Phenomenological Model

Up to now, the phenomenological models of fatigue delamination propagation in composite materials are mostly based on classical Paris Law. According to the discussion on influencing factors of fatigue crack growth rate in Section 2.3, the existing modified Paris Law model is unified and summarized. The original formula of the modified Paris Law model is as follows Eq.(4).

$$\frac{da}{dN}=C f{\left(K or G\right)}^m$$

(4)

3.3.1. Mixed-Mode Ratio Modified Model

For the Paris Law considering the mixed-mode ratio, SERR is generally used to characterize load level of composite materials. Thus, the similar parameters f (K or G) in Eq.(4) is simplified to f (G). When mode III is assumed equal to mode II, four issues need to be considered while establishing the mixed-mode ratio modified model.

Firstly, it needs to present the methods defining the contribution of different modes to fatigue crack growth rate in mixed-mode loading. Therefore, the uncoupled models [77,78] and coupled models [79] have been proposed. Secondly, it requires to determine the approach normalizing the load level [63,80]. Thirdly, coefficients of Paris Law may use the material constants or the functions of mixed-mode ratio. If the functions of mixed-mode ratio are adopted, the monotone [81,82] or the non-monotonic [19] function relationships have been presented. Fourthly, an approach introducing upper and lower limits [83] of fatigue crack growth rate to Paris Law model is also required.

From the simplest uncoupled models of mode I and mode II, modified Paris Law which considers uncoupled mixed-mode can be proposed by linear superposition. The general formula is given by

$$\frac{da}{dN}={\displaystyle \sum _i{C}_i}f{\left(G_i\right)}^{m_i}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} i=\mathrm{I},II$$

(5)

where, C_i and m_i are material constants of mode i, f (G_i) are functions of loading levels. Different models can be obtained by changing forms of function. For example, the models of Gustafson and Hojo [78] can be obtained when f (G_i) = G_imax - G_imin. Nevertheless, the models of Ramkumar and Whitcomb [77] can be obtained when f (G_i) = G_i /G_ic, which is normalized by critical SEER of mode i. Therefore, the fatigue crack growth rates of different FRP composites can be compared.

However, laminated composite structures under fatigue loading are usually in mixed mode. What’s more, the interaction between mode I and mode II/III in process of delamination propagation can’t be ignored. The simplest general interaction model between mode I and mode II is as follows.

$$\frac{da}{dN}=C(\varphi )f{(G)}^{m(\varphi )}$$

(6)

When considering the general linear interaction,

$$\begin{array}{l}C(\varphi )={\displaystyle \sum _i{C}_i}{\varphi }_i,\hspace{1em}{\displaystyle \sum _i{\varphi }_i}=1\hfill \\ m(\varphi )={\displaystyle \sum _i{m}_i}{\varphi }_i,\hspace{1em}{\displaystyle \sum _i{\varphi }_i}=1\hfill \end{array}$$

(7)

It is observed that the essence of considering coupling effect of mode is to regard coefficients of Paris Law as a function of the mixed-mode ratio. However, C(ϕ) and m(ϕ) in Eq.(6) which is following the linear coupling relationship like Eq.(7) are monotone functions. For the epoxy matrix composite, it has been proved [80] that the coefficients of Paris Law are non-monotone functions of the mixed-mode ratio. Therefore, a non-monotone function can be constructed which needs to satisfy two conditions. Firstly, values of endpoint are equal to Paris Law coefficients of mode I and mode II respectively. Secondly, the stationary point of non-monotone function is within the range of 0 and 1. The quadratic function of ϕ can satisfy requirements by controlling boundary conditions.

$$\begin{array}{l}m(\varphi )=d+e\varphi +f{\varphi }^2,\hspace{1em}m(0)=m_I,\hspace{1em}m(1)=m_{II}\hfill \\ C(\varphi )=a+b\varphi +c{\varphi }^2,\hspace{1em}C(0)=C_I,\hspace{1em}C(1)=C_{II}\hfill \end{array}$$

(8)

The most general non-monotonic modified model including the effect of mixed-mode ratio is proposed by Blanco et al. [19]. The model is based on experimental observation in Section 2 and procedure of construction above.

As we all known, the classical Paris Law only describes middle-rate region of fatigue delamination. In order to simultaneously describe the thresholds of fatigue crack growth rate curve, it is possible to multiply Eq.(6) by functions which is representing lower and upper horizontal asymptote. The general formula of this kind of model is given by

$$\frac{da}{dN}=C(\varphi )f{(G)}^{m(\varphi )}{g}_{th}(G)h_c(G)$$

(9)

where, the loading level G can be expressed in various forms and normalized or not. The multiplier g_th(G) characterizes the lower horizontal asymptote of the fatigue crack growth rate curve. The value of function needs to satisfy two conditions. When load level is infinitely close to the fatigue crack growth threshold G_th, the value of function tends to be 0. When load level is slightly larger than the fatigue crack growth threshold G_th, the value of function tends to be 1. In a similar way, the multiplier h_c(G) characterizes the upper horizontal asymptote of the fatigue crack growth rate curve. The value of function needs to satisfy two conditions. When load level is infinitely close to the critical fatigue crack growth G_c, the value of function tends to be positive infinite. When load level is slightly less than the critical fatigue crack growth G_c, the value of function tends to be 1.

Finally, the requirements of g_th(G) and h_c(G) are satisfied simultaneously. Therefore, the modified Paris Law which can describe the complete fatigue crack growth rate curve is obtained. The classical form of the function g_th(G) and h_c(G) can refer to the phenomenological model of Shivakumar et al. [83].

3.3.2. Fiber Bridge Modified Model

Based on the extensive application of unidirectional composite laminates in practice, the influence of the fiber bridge on fatigue crack growth rate can’t be neglected. By observing the phenomenological models which are considering the fiber bridge effect based on Paris Law, four issues need to be solved when establishing the model.

Firstly, it needs to select proper parameter which is similar to f (G). Secondly, the way considering the effect of fiber bridging by normalization is required. Thirdly, it requires to be determined that the coefficient of Paris Law is a material constant or a function of fatigue crack length. Fourthly, it needs to present a method which introduces the upper and lower threshold to fatigue crack growth curve. Based on the above considerations, the following general formula is given by

$$\frac{da}{dN}=C(a)f{\left(G(a),G_{CR}(a)\right)}^{m(a)}$$

(10)

where c(a) and m(a) are coefficients of Paris Law which depend on the fatigue crack length. G_CR(a) denotes the R-curve obtained from the fracture toughness experiment.

In the model of Shivakumar et al. [83], the definition of fatigue crack growth rate is simple. The function f = G_max(a) / G_CR(a), C and m denote the material constants. However, Gong et al. [68] have deduced the reasonable form of SERR by superposition principle from the concept of SERR. ${\left(\sqrt{G_{max}}-\sqrt{G_{min}}\right)}^2$ has a good agreement with the similar parameters ∆K in classical Paris Law. Yao et al. [66,67] elaborate on how to consider the effect of fiber bridge on mode I fatigue crack growth rate in a modified Paris Law. On the basis of above similar parameters, the Hartman-Schijve equation [60,84] and the idea of equivalent SERR at the crack tip [85] are used. Therefore, similar parameters which can consider the upper and lower threshold of fatigue crack growth curve and fiber bridge are respectively constructed by Eq.(11). These two phenomenological models based on modified Paris Law have been verified by various fatigue delamination tests. In addition, it is also concluded that the Paris Law coefficient m of mode I fatigue delamination is independent of fatigue crack length in FRP composite materials.

$$\begin{array}{c}f\left(G(a),G_{CR}(a)\right)=\frac{\sqrt{\Delta G(a)}-\sqrt{\Delta G_{th}(a)}}{\sqrt{[1-\sqrt{G_{\max }}(a)/\sqrt{G_{CR}(a)}}]}\\ f\left(G(a),G_{CR}(a)\right)=\frac{G_0}{G_{\mathrm{CR}}(a)}\Delta G(a)\end{array}$$

(11)

where ∆G(a), ∆G_th(a) and G_max(a) are functions of SERR range, the fatigue crack growth threshold range and the maximum SERR which varies with the fatigue crack length a respectively. G₀ denotes the initial SERR without fiber bridge effect. This model can well capture the fatigue crack growth rate at different fatigue crack lengths with high accuracy and simple fitting.

3.3.3. Loading Ratio Modified Model

The influence of loading ratio on fatigue delamination has always been a focus. From the analysis in Section 2.3, it is well known that as long as the combination of maximum and minimum loading levels such as G_mean, ∆G and G_min/G_max are involved in modified Paris Law, which is equal to consider the influence of loading ratio on the fatigue crack growth rate. Therefore, five issues need to be concentrated on when establishing a general model.

Firstly, it needs to select similar parameters in classical Paris Law. Secondly, the method of normalization is required to be determined. Thirdly, it is required to build the expression of loading level. Fourthly, the coefficient of Paris Law may use a material constant or the function of loading ratio. Finally, an approach introducing the upper and lower threshold of modified Paris Law is also required. Therefore, the general formula for developing such a model is as follows

$$\frac{da}{dN}=C(R)f{\left(G_{\max },G_{\min }\right)}^{m(R)}$$

(12)

where C(R) and m(R) are coefficients of Paris Law which depends on loading ratio. In simplified model, C and m can also represent material constants. f (G_max, G_min) denotes all combined expressions of maximum and minimum loading levels including normalized models. The simplest expression is ∆G = G_max - G_min.

Considering loading ratio in fatigue crack growth rate, the model proposed by Hojo et al. [47] is the most commonly used. The model represents the relative contribution of the SIF range ∆K and maximum SIF K_max by parameter γ which ranges from 0 to 1. The construction of function f (G_max, G_min) is shown in Eq.(13).

$$f\left(G_{\max },G_{\min }\right)=K_{\max }^{\gamma }\Delta K^{1-\gamma }$$

(13)

In fact, after taking the logarithm on both sides of Eq.(13)，

$$\log f=\gamma \log K_{\max }+(1-\gamma )\log \Delta K$$

(14)

The kernel of constructing the function f (G_max, G_min) is the logarithmic weighted average of ∆K and K_max. Actually, the concept is reflected in many classical modified Paris Law considering loading ratio. For instance, the expression of Eq.(15) was specially proposed for glass fiber reinforced laminated composites by Atodaria et al. [86].

$$f={(\sqrt{G})}_{\mathrm{mean}}^{\gamma }{(\Delta \sqrt{G})}^{1-\gamma }$$

(15)

However, the model mentioned can’t consider the effect of high and low rate regions of the fatigue crack growth curve. By directly incorporating the horizontal asymptote function into f (G_max, G_min) refer to the Eq.(9), the models of Shintarou et al. [87] considering the low rate regions and Blanco et al. [19] considering both the high and low rate regions are obtained. The specific functions are given as follows.

$$\begin{array}{c}f\left(G_{\max },G_{\min }\right)=\Delta K-\Delta K_{th}\\ f\left(G_{\max },G_{\min }\right)=\frac{2K_{\mathrm{mean}}\left(\Delta K-\Delta K_{\mathrm{th}}\right)}{K_{\mathrm{C}}^2-K_{{\max }^2}}\end{array}$$

(16)

In fact, the model which is considering fiber bridging proposed by [60,85,88] also take the effect of loading ratio into account. In addition to replacing similar parameters with${\left(\sqrt{G_{max}}-\sqrt{G_{min}}\right)}^2$, the coefficient of Paris Law C (a, R) is set as a function of the fatigue crack length a and the loading ratio R. However, m(R) only depends on loading ratio R. By this way, the fatigue crack growth curves under different loading ratios are unified under the influence of fiber bridge.

3.3.4. Temperature Modified Model

As for the influence of hygrothermal conditions on fatigue crack growth rate, the influence of temperature is simply discussed here and the research of humidity can refer to temperature. Miura et al. [89] have studied I/III mixed mode fatigue crack growth curves of woven glass/epoxy composite laminates at 295K (room temperature), 77K and 4K. The results indicate that fatigue crack growth rates under three different temperatures satisfy da/dN_295K > da/dN_4K > da/dN_77K. However, the modified Paris Law is primitive just like Eq.(17).

$$\frac{da}{dN}=C{\left(\Delta G_T\right)}^m$$

(17)

Above all, there is issue in similarity of ∆G_T and ∆K. Secondly, the coefficient of Paris Law is not considered as a function of temperature. Therefore, not only accuracy and application of the model can’t be guaranteed, but also the influence of environmental temperature under same material system can’t be unified in a single fatigue crack growth curve. The model of Zhong et al. [90] also has similar defects.

In latest investigations, Jin et al. [21] measured the tensile-tensile fatigue crack growth rate of unidirectional and orthogonal Ti/CFRP Laminates at 25℃, 80℃, 120℃ and 150℃, respectively. The fatigue crack length nonlinearly increases with the increase of temperature. What’s more, acceleration of fatigue crack growth is obvious during the temperature difference from 120℃ to 150℃. Based on the concept of SIF, the general formula of modified Paris Law which is considering the effect of temperature is given by

$$\frac{da}{dN}=C(T){(\Delta K)}^{m(T)}$$

(18)

According the experimental phenomena, the fitting functions which are satisfying requirement of nonlinear monotone behavior of C(T) and m(T) are as follows

$$\begin{array}{l}C(T)=C_0+C_1{e}^{D_{0}T}\hfill \\ m(T)=m_0+m_1{e}^{n_{0}T}\hfill \end{array}$$

(19)

where C₀, C₁, D₀, m₀, m₁ and n₀ are fitting parameters of experimental data. When material system and environment condition change, this kind of model with coefficient of Paris Law as the temperature function is very effective. The most important is that the reasonable fitting function should be adapted through observation of experiment.

3.3.5. Other Modified Models

With respect to influence of loading sequence on fatigue crack growth rate, investigators have different opinions because of the lack of experimental data about different material systems and loading conditions. Specific progress in this field can refer to references [49,65,69,91].

In terms of secondary factors, the FDM which considering the influences of loading frequency and loading waveform on fatigue crack growth rate can be reconstructed by Eq.(12)~(14). The term which is characterizing the effect of frequency and waveform can be directly added into expression of the model. Furthermore, the steps in details can be referred to relevant articles [92].

4. Finite Element Simulation of Delamination Fatigue

In this section, FDDMs are classified into four categories based on finite element model of fatigue delamination onset and propagation under fatigue loadings. In which, the framework of theories including strength model, fracture mechanics model, damage mechanics model and hybrid model. In order to comprehensively compare existing FDDMs, establishment process, application scope, prediction ability and limitation of various models are all taking into account.

4.1. Strength Model

Generally, the fatigue model of FRP composites can be classified into three types [93]: the model based on S-N curve, residual stiffness/strength model and progressive damage model. The similarity of these models is that they are based on classical strength theory and different forms of strength failure criterion. A similar method can be used to establish FDDMs of laminated composite structures.

In early stage, the methods of delamination simulation are mainly based on strength theory. After obtaining stresses field of laminated composite structures, the delamination damage is calculated by strength failure criterion expressed by stress components. Among macroscopic strength failure criterion of composites, criterions such as maximum stress, maximum strain, Tsai-Hill, Tsai-Wu and Hoffman are commonly used in engineering. For example, Bernasconi et al. [7] directly verified the correlation between the maximum shear stress and the number of failure cycles. Thus, τ-N curve which was obtained through SLJ test provided a simple and practical criterion for bonding structures.

However, these criteria are not able to consider failure mode of delamination alone. Hence, investigators have proposed microscopic strength failure criterion for delamination, such as Hashin criterion [94], Chang-change criterion [95], Hou criterion [96], Zou criterion [97] and LaRC criterion [98], which are shown in Table 1. These criteria are usually expressed by quadratic polynomials. Furthermore, the essence of criteria is transforming complex stress state into equivalent stress state. In many cases, in order to reduce the effect of stress singularity at crack tip, stress strength criteria are also rewritten as strain strength criteria.

While considering the influence of fatigue delamination in engineering, the maximum stress due to cyclic loading is usually less than 50% of ultimate strength of material. It is unreasonable that the structures without degradation of material properties never suffer from strength failure. Thus, material properties must be reduced in simulation of fatigue delamination. The residual stiffness/strength model and fatigue progressive damage model (FPDM) [99] are essentially models within the framework of damage mechanics for predicting fatigue delamination onset using strength theory. Furthermore, it will be elaborated in Section 4.3.

4.2. Fracture Mechanics Model

It should be emphasized that the fracture mechanics models are commonly based on framework of LEFM [100,101,102]. LEFM assumes that the complexity of the crack tip is not considered and the stress singularity of the crack tip is avoided. Instead, stress states of outer region of the crack tip are used to characterize the fracture behavior.

Compared with strength criterion, fracture mechanic model has developed a series of standard tests and formulas for solving the stress field at the crack tip under various crack forms. However, the unit of SIF K which is necessary for developing fracture criterion is [N][L]^-2/3. The quantity lacks of explicit physical meaning. Thus, SERR G defined by Irwin [103] represents the energy consumption per unit area of crack growth and the unit is [N][L]^-1. Through the strict proof of equivalence between G and K, fracture criterion can be automatically transformed into energy criterion. For instance, the relationship between SERR G_I and SIF K_I of mode I delamination is shown in Eq.(20). The loading levels can be characterized by SERR with more specific physical meaning in the process of fatigue crack growth. The main numerical computation method can refer to the brief introduction and references in Section 3.1.

$$G_I=\frac{K_I^2}{E_1}$$

(20)

In order to establish the relationship between fatigue crack growth rate and load level, almost all phenomenological models of fatigue delamination are based on the modified Paris Law. Based on the experiment, the influence of parameters and the methods of construction can refer to sections 2.2~2.3 and 3.2~3.3. The schematic diagram of FDM [104,105] based on Paris Law of fracture mechanics is shown in Figure 11. The steps in details are as follows.

• Artificial pre-cracks are added to laminated composite structures. Set the growth path of crack propagation and adapt displacement loading mode.
• When structure is applied cyclic loading, it is assumed that f (K or G) is constant during the increasement of fatigue crack length ∆a. The general Paris Law is converted into incremental form as shown in Eq.(21).

$$\Delta N=\frac{\Delta a}{Cf{\left(K\hspace{0.17em}or\hspace{0.17em}G\right)}^m}$$

(21)

3.

By means of VCCT or J-integral method in FEM, f (K or G)_n near the crack tip of delamination in current increment step is calculated.

4.

(∆N)_n is acquired by substituting known coefficients of Paris Law C and m, ∆a and f (K or G)_n which are obtained in previous step into Eq.(21). It should be noticed that C and m can be material constants or functions of parameters such as fatigue crack length a.

5.

Update cycles of fatigue loading N_n+1 by Eq(22). The stiffness of structure declines because of delamination damage. Go back to step (2) and calculate delamination propagation of next incremental step until it stops when f (K or G) ≤ f (K or G)_th. Record final fatigue crack length a and cycles of fatigue loading N. f (K or G)_th is inherent properties of material which is obtained in test of fatigue crack growth rate.

$$N_{n+1}=N_n+{(\Delta N)}_n$$

(22)

Paris Law is directly applied to fatigue delamination in finite element software. The life prediction of laminated composite structures can be completed as long as the general formulas summarized in Section 3.3 are used to construct the modified Paris Law. In addition, various influence parameters mentioned before can be considered. Then, the data of material system are fitted by experiments.

However, the fracture mechanic model also has limitations [106]: (1) The essence of finite element is discretizing structure, and increment of fatigue crack length ∆a is usually equal to characteristic length of element L_c. The higher the mesh density is, the higher the accuracy of fatigue crack growth curve da/dN-∆G is. Therefore, the model is strongly sensitive to mesh size. (2) Coefficients of Paris Law are parameters of material properties. Thus, it is difficult to accurately simulate fatigue delamination at the interface of mixed materials. (3) The FEM needs the growth path of pre-crack in advance. What’s more, it is not able to predict onset of delamination and simulate the spatial cracking because there are no parameters in traditional Paris Law to determine the direction of delamination.

4.3. Damage Mechanics Model

The damage mechanics model [107] can directly predict delamination onset and propagation. Thus, compared with fracture mechanics model, it has been paid more attention in recent years. At present, numerical methods of fatigue delamination based on damage mechanics can be classified into three types: FDDMs based on degradation of intra-laminar element, FDDMs based on degradation of inter-laminar element and FDDMs based on degradation of hybrid intra-laminar element and inter-laminar element. Taking delamination of mode I as an example, the schematic diagrams are shown in Figure 12.

4.3.1. Intra-laminar Element Model

Based on degradation of intra-laminar element in FDDMs, stresses of intra-laminar element around the crack tip are substituted into strength criterion to estimate the delamination onset. Furthermore, delamination propagation is simulated by damage evolution. The concept and procedure of residual stiffness/strength model and FPDM are similar. Therefore, the specific analysis process of residual stiffness/strength model is given in Figure 13.

As the number of cycles increases, properties of the intra-laminar element gradually degrade [99]. Stiffness degradation leads to redistribution of stress fields. Furthermore, laminated composite structure achieves a higher stress state. Strength degradation changes parameters of strength in failure criterion. Therefore, local static strength failure occurs first in structure, and it finally leads to overall structural failure.

In addition to strength failure criterion of delamination, many investigators have studied degradation model of material parameter, damage variable [108,109,110] and damage stiffness matrix [111] on intra-laminar elements. Although intra-laminar element model can implement the simulation of delamination onset and propagation, the failure mechanism is essentially different from the interface of actual structures. Moreover, stiffness/strength degradation model depends on a large number of experiment data, which increases cost. Besides, the empirical assumptions of damage variable and constitutive model reduce the prediction accuracy of fatigue delamination.

4.3.2. Inter-laminar Element Model

The cohesive zone model (CZM) is the main method of FDDMs based on degradation of inter-laminar element. The history of CZM can be traced back to the D-B model [112,113] which is established by Dugdale and Barenblatt. CZM is essentially a phenomenological model. The traction-separation law of cohesive element is a phenomenological description of fracture process near the crack tip. CZM has simple constitutive equation and clear physical meaning. What’s more, it is easy to be realized in FEM software. Therefore, simulation of fatigue delamination based on cohesive elements has been rapidly developed. The basic procedures of inter-laminar element model under fatigue loading are as follows.

(1)

Derive the stiffness matrix and the constitutive equation (that is cohesive law) of cohesive element.

Camanho et al. proposed an interface element between solid elements, which can simulate onset and propagation of mixed-mode delamination. This element has been incorporated into ABAQUS software. The cohesive element established by Jiang et al. [114] considered the change of mixed-mode ratio during delamination propagation. What’s more, the numerical results are in good agreement with the experimental results. In terms of cohesive law, Needleman [115] proposed the exponential law, Tvergaard and Hutchinson [116] proposed the trapezoidal law, Mi et al. [117] proposed the bilinear law, Heidari-Rarani et al. [118] proposed the trilinear law. These constitutive relations are shown in Figure 14. However, relevant studies [119] have indicated that different forms of constitutive relations have no obvious effect on the load-displacement curve of laminated composite structures. The bilinear law has good computational efficiency, accuracy and convergence. Therefore, most of FDDMs based on cohesive elements generally use the simplest bilinear softening constitutive relations.

(2)

Determine the damage initiation and evolution criterion of mixed-mode fatigue delamination.

At present, the damage initiation criterions of delamination which are integrated in ABAQUS include the maximum principal stress/strain criterion, the maximum nominal stress/strain criterion and the quadratic nominal stress/strain criterion. The damage evolution criteria of delamination include Power law, B-K law [120] and Reeder law [121], in which Reeder law can automatically degenerate to Power law when G_IIC=G_IIIC [6]. Turon et al. [122] proposed a delamination initiation criterion which takes consideration about the effect of crack closure and satisfies the condition of thermodynamic consistency. By controlling the relationship between interface strength and stiffness, a method for accurately predicting the propagation of mixed-mode delamination has been developed [102]. The common criteria introduced above are summarized in Table 2.

Figure 15 shows the quadratic nominal stress criterion and B-K law of mixed-mode delamination which are the most commonly used. The black solid line derived from the arrow in Figure 15 is the envelope of damage initiation and evolution of cohesive element. Therefore, different damage initiation and evolution criteria of delamination should be constructed according to material systems. It is the guarantee of accurately tracking the damage initiation and critical fracture toughness under mixed-mode delamination.

(3)

Choose the strategy of FDDMs based on loading-unloading hysteresis [123,124,125] or envelope load (cycle jump approach) [126,127].

The patterns of inter-laminar properties degradation under fatigue loading are shown in Figure 16. It illustrates the difference between two kinds of model. It is necessary to calculate degraded parameters of each cycle in loading-unloading hysteresis damage model. This method is not desirable under high-cycle fatigue delamination because of extreme high computational cost. However, only the maximum value of load block is considered in envelope load damage model. It is assumed that the damage is continuous and average in each loading block [128], which is greatly improving the efficiency of simulation.

At present, most damage models which are based on CZM adopt the method of envelope load. Zhang et al. [129] proposed an explicit FEM for predicting the fatigue delamination of composite laminates by using twin cohesive models and a combined static and fatigue cohesive formulation. Compared with the traditional envelope load method, twin cohesive model is subjected to peak and valley envelope load respectively which means that the effect of loading ratio is considered in the model. At the same time, it is significant that the novel model improves the computational efficiency of fatigue delamination.

(4)

Establish the differential relationship between damage variables and the number of loading cycles. Furthermore, the reasonable formula of fatigue delamination damage rate can be obtained.

In situation of applying the envelope load method to simulate fatigue delamination, the total damage rate of delamination dD/dN is the sum of quasi-static damage rate dD_s/dN and fatigue damage rate dD_f /dN. Assuming that the current cycles of loading is N, the damage variable is updated by Eq.(23) after increment of loading cycles ∆N.

$$D(N+\Delta N)=D(N)+{\displaystyle {\int }_N^{N+\Delta N}\frac{dD}{dN}}dN$$

(23)

Generally, left rectangle formula [14,130,131] and trapezoid formula [132] are used to numerically integrate Eq.(23) to obtain the damage variable under current cycles of loading.

Robinson et al. [133] proposed a fatigue damage rate model with three parameters. The curves of fatigue crack growth rate obtained by simulation are consistent with the fitted Paris curves obtained from experiment data of mode I, mode II and mixed-mode. Tumino and Cappello [132] developed a fatigue damage rate model with two parameters, which also reproduced phenomenon of fatigue delamination tests. The difference is that one of parameters is a function of mixed-mode ratio. In three parameters model established by Khoramishad et al. [130], both the crack growth threshold of fatigue delamination and the mixed-mode ratio were taken into account simultaneously. Thus, it is more applicable than the former models. It should be noted that the parameters of FDDMs purely based on damage mechanics and can’t be directly obtained by existing standard tests. Therefore, fitting parameters of the model can only be adjusted by FE model calibration method [134,135,136] to reproduce the crack growth curve of fatigue delamination. At present, the correlation between this model and geometry and properties of structures has not been clearly clarified. Hence, it hinders application of directly applying the model to engineering problem for life predicting of fatigue delamination.

(5)

Select the inter-laminar parameters of fatigue delamination damage that need to be degraded. Then, establish the functional relationship between degradation parameters and damage variables.

In CZM, there are only three inter-laminar parameters. Penalty stiffness K, interfacial strength σ^o and critical fracture toughness G_c are required to describe bilinear cohesive relationship. The essence of degradation of inter-laminar cohesive elements is to establish a functional relationship between K, σ^o, G_c and damage variable D(N). Furthermore, the function is applied for softening cohesive law. The simplest softening curve of bilinear constitutive model [136] is shown in Figure 17.

The process of bilinear softening law is simple. It inherits the concept of residual strength/stiffness of metal fatigue and defines the variables of residual penalty stiffness/ interface strength/fracture toughness, which are naming after $K_{{\mathrm{N}}_{\mathrm{i}}}$,${\sigma }_{{\mathrm{N}}_{\mathrm{i}}}^{\mathrm{o}}$,$G_{\mathrm{C}}^{{\mathrm{N}}_{\mathrm{i}}}$ in Figure 17. The model assumes that the initial bilinear constitutive relationship of cohesive element is always maintained in process of fatigue loading. With the cycles of fatigue loading increasing, the corresponding parameters such as penalty stiffness, interfacial strength and fracture toughness continuously degrade. In other words, these residual parameters defined based on FDDMs are related to fatigue life N. Therefore, as long as the relationship between damage parameters of inter-laminar element and fatigue life is given, a complete FDDM can be established.

Residual and initial inter-laminar parameters are set to be R(N) and R(N)⁰respectively. What’s more, the most commonly used degradation function of inter-laminar performance is as follows.

$$R(N)=(1-D)R{(N)}^o$$

(24)

The model established by Turon et al. [102] and Chen et al. [137] respectively degraded the inter-laminar parameters of interface strength and penalty stiffness according to Eq.(24). The model established by Mazaheri and Hosseini et al. [11] simultaneously degraded six inter-laminar properties corresponding to mode I and mode II delamination. The performance of models are consistent with those of Roe and Siegmund [138]. Since the function of degradation is empirical, there is no specific correlation between the degradation of inter-laminar performance parameters and physical mechanism of fatigue delamination.

(6)

Realize loading/unloading cycles of FDDMs by FEM software and formulate the standard of integral structure failure. Finally, fatigue life of delamination under specific loading can be calculated.

Taking fatigue delamination of DCB as an example, a modified FDDM is established by referring to the model proposed by Koloor et al. [136]. As shown in Figure 18, the results of FEM are obtained by UMAT subroutine in ABAQUS. The model is loaded with prefabricated crack within the first second. Afterwards, the zigzag wave with frequency of 5Hz is applied. The loading mode is displacement control. Figure 18 (a) shows the local stiffness degradation at the fatigue crack front of the DCB specimen. Fig 18 (b) presents Mises stress of 10 cohesive elements at the crack front. The smaller the element number is, the closer it is to the crack front. It is obvious that first five elements have already failed in stage of pre-cracking. The last five elements gradually fail under the fatigue loading of delamination. Corresponds to Figure 18 (b), Figure 18 (c) clearly demonstrates the evolution process of damage variable, which is representing the damage degree of the element at the crack tip. Because of displacement control, as shown in Figure 18(d), the descending speed of load-displacement curve gradually slows down and finally tends to stop. This is consistent with the delamination crack growth rate which is observed in Figure 18 (e).

In order to further realize the implement of fatigue delamination crack growth with cohesive element model, the sensitivity analysis of parameters is carried out on the basis of quasi-static delamination damage. Gustafson and Waas et al. [139] studied the influence of penalty stiffness, inter-laminar strength, fracture toughness and constitutive shape on FEM results of DCB and ENF. The results indicated that the output of DCB model mainly depends on one parameter, while the ENF model depends on multiple parameters. The same is that the constitutive shape has little influence on results of the model. Turon et al. [140] found that inappropriate parameters might lead to large error between simulation and experiment. Therefore, a penalty stiffness model with variation of mixed-mode ratio was developed to adjust the model parameters, which is dynamically matching the result of experiment. In order to better illustrate the influence of parameters on numerical examples, the constitutive relationship of cohesive element has been realized in UMAT subroutine. The accuracy of the program is verified by comparing with the models in ABAQUS help documentation. Figure 19 gives the sensitivity of mesh size in numerical simulation. In addition, it visually presents the influence of the mesh size on load-displacement curve of MMB under different mixed-mode ratios. It can be seen that the influence of mesh size is significant and non-negligible. The convergence of critical load and displacement is sensitive to the element size.

4.3.3. Mixed Intra-laminar/Inter-laminar Element Model

Detailed numerical procedures of fatigue delamination based on mixed intra-laminar/inter-laminar element model are not repeated any more. The degradation of intra-laminar and inter-laminar elements are considered simultaneously in finite element model. In addition, the fatigue damage of two different kinds of element affects each other.

Du et al. [42] combined intra-laminar element of Linde et al. [111] with the exponential cohesive law by referring to an example of tensile failure of fiber-metal laminates (FML) in ABAQUS help documentation. The failure analysis of Ti/Cf/PEEK laminates with delamination was carried out. Tarfaoui et al. [141] studied the progressive delamination damage failure of CFRP composite laminates with open hole by different failure criteria of intra-laminar elements and cohesive element. Aoki et al. [24] classified FPDMs of laminated composite structures into two types: the intra-laminar damage model and the inter-laminar damage model. The fatigue tests of DCB, ENF and curved laminates verified the accuracy of the model in respect to life prediction of CFRP composite laminates in different levels of loading.

The mixed intra-laminar/inter-laminar element model takes the fatigue damage of intra-laminar element into account in while simulating damage evolution of inter-laminar element. In particular, with regard to laminated composite structures which have complex stacking sequence and geometries, coupling effects always accelerate the process of fatigue delamination crack growth and make fatigue life prediction closer to actual conditions. However, this model also has the similar issues of inter-laminar element model and intra-laminar element model. It means that the model has limitations in obtaining parameters of numerical simulation and empirical degradation assumptions, which is restricting development of the model to some extent.

4.4. Hybrid Model

In study of FDDMs, models which are based on damage mechanics are not able to essentially reflect the properties degradation of laminated composite structures under fatigue loading. What’s more, it does not establish the relationship between fatigue delamination growth rate in numerical simulation and standard experiment. Then, investigators start to combine CZMs based on damage mechanics with Paris Law models based on fracture mechanics. They have explored a general method of establishing FDDMs based on fatigue delamination growth test of laminated composite structures.

The procedure based on hybrid fracture and damage mechanics model can refer to damage mechanics models in Section 4.3. In step 4, fatigue damage rate dD/dN is related to fatigue crack growth rate da/dN of Paris Law by chain rule [142]. The final general formula of fatigue damage rate is as Eq.(25).

$$\frac{dD}{dN}=\frac{dD}{da}\frac{da}{dN}=\frac{dD}{da}Cf{(G)}^m$$

(25)

Turon et al. [14] proposed a damage model to simulate delamination propagation under high-cycle fatigue loading. The relationship between damage variable and crack growth rate in Paris law is established by cohesive element. The verification test of element is carried out by UMAT subroutine in ABAQUS. Furthermore, accuracy of the model is verified by comparing with the experiment of fatigue delamination. At the same time, it is also the basis procedure for many researchers in this field to deeply study the fatigue delamination damage of interface elements. The performances of cohesive element reproduced by UMAT subroutine in ABAQUS are given in Figure 20. In which, Figure 20 (a) shows the cohesive element model and Traction-Logarithmic Cycles curve. The traction of cohesive element experiences a sudden drop and then slow drop. Figure 20 (b) shows the Traction-Displacement curve of cohesive element, which fails after about 28K cycles of fatigue loading.

Harper and Hallett [143] have developed a degradation law of cohesive elements for analysis of delamination propagation under fatigue loading. The model directly related the SERR of cohesive zone to the Paris curve which was obtained from the experiment. Thus, it reproduces the phenomenon of delamination propagation in DCB, ENF and MMB tests. The results are in good agreement with experimental data and analytical solutions in the literature. Nixon-Pearson et al. [144] established a cohesive element based on Paris law to simulate relationship of delamination and matrix cracking in composite laminates with open hole. The S-N curve is consistent with the experimental results. Because of the sensitivity of hybrid fracture and damage mechanics model to element size, Tao et al. [145] proposed a method of using virtual fatigue damage variable to locate the element at the crack tip. By combining the modified Paris law and extended cohesive element, a hybrid model with weak dependent on mesh size was obtained. The numerical simulation is able to reproduce the results of DCB and 4ENF under fatigue loading. Ebadi-Rajoli et al. [146] developed a new fatigue damage accumulation (FDA) model by correlating cohesive element based on fracture mechanics and damage mechanics, and adopting degradation of interface stiffness.

Compared with damage mechanics models, the hybrid fracture and damage mechanics models which are based on Paris Law does not need to be calibrated by parameters of numerical simulation to match fatigue delamination growth rate. On account of Paris Law based on LEFM framework is an inherent property of material and it does not depend on geometry of specimen. Therefore, fatigue delamination growth rate curve of the standard specimen, that is input parameters of hybrid model, reflects the phenomenological model summarized from the experiment. In addition, the general modified Paris Law which considers the influence factors of fatigue delamination propagation in Section 3.3 can arbitrarily form hybrid fracture and damage mechanics model with CZM. The high generality and operability can satisfy requirements of research and engineering.

The damage of angle-ply laminates with open hole under fatigue loading is not only the delamination propagation of cohesive elements, but also including coupling effect of matrix cracking and fiber-matrix shear failure of delamination interface elements. Therefore, it is natural to combine intra-laminar element based on strength criterion and inter-laminar cohesive element based on modified Paris Law. Zhou et al. [8] developed such a hybrid model and carried out the fatigue life prediction of GFRP and CFRP composite joints to validate accuracy of the model. This kind of model which synthesizes strength criterion, fracture mechanics and damage mechanics has strong structural universality and engineering practicability in theory. However, more standard experiments are needed to support input parameters of model, which has the limitations of low efficiency and high cost.

5. Summary and Prospect

Based on delamination phenomena of laminated composite structures under fatigue loading, this review discusses the general methods of establishing phenomenological model for fatigue delamination. It focuses on the influence factors and the correction techniques of Paris law. Furthermore, classification, procedure, performance and limitations of existing FDDMs are also included.

As a phenomenological model based on framework of LEFM, Paris law has been widely studied in considering the effects of mixed-mode ratio, loading ratio, threshold of fatigue delamination growth rate, fiber bridging, temperature and humidity, and loading sequence on fatigue delamination life. The cohesive element model based on damage mechanics has clear physical meaning and simple constitutive relation. Therefore, the FDDMs based on modified Paris Law within framework of hybrid fracture and damage mechanics is the most effective and stable method so far to simulate fatigue delamination onset and propagation.

However, the phenomenological model is not able to explicitly explain the fatigue delamination mechanism of composite materials. The establishment of FDDMs mostly depend on accumulation of experiment data and assumption of semi-empirical formula. Thus, it is difficult to obtain accurate results when studying complex structures and loading conditions, which greatly limits its application in practical engineering. Furthermore, the following five directions of FDDM are suggested according to aforementioned technical problems to be solved: (1) Multi-scale CZM fatigue delamination model simultaneously considering macro and micro model. (2) Use extended finite element method (XFEM) to analyze spatial fatigue delamination onset and propagation. (3) Use NURBS to study fatigue delamination life of irregular geometric model. (4) Construct higher-order composite beam element or plate element to predict fatigue delamination onset and critical loading of laminate structure. (5) Integrate the artificial intelligence algorithm into the data processing of the phenomenological model.