A Fatigue Model to Predict Interlaminar Damage of FRP Composite Laminates Subjected to Mode I Load

In fiber-reinforced polymer (FRP) composite laminate structures operating under fluctuating stresses, interface delamination is seen as one of the significant damage mechanisms. The constant degradation of their relatively low interlaminar strength and stiffness are the primary reasons for delamination. This study develops an interlaminar fatigue damage model to quantify the mechanics of the damage process and address the reliability of composite structures. The model considers the failure process in two stages: (1) damage due to degradation of interlaminar elastic properties, and (2) damage due to dissipation of fracture energy through the damage evolution process. The model is examined for a case study of mode I fatigue loading of a carbon-fiber-reinforced polymer (CFRP) composite laminate. The results show that the interlaminar normal stress is confined to the crack front region, with tensile stress peaks at 70% of the interlaminar strength. Furthermore, a stable interface crack growth is predicted initially, followed by a sudden crack “jump” at 14,000 cycles. The simulation results are compared with the experimental data, with very good agreement, showing a successful validation of the fatigue model.


Introduction
Fiber-reinforced polymer (FRP) composite laminates have attracted attention for many industrial applications in recent decades. Their high strength-to-weight ratio makes them desirable for aerospace applications. Typical applications include the fuselage and spoilers of aircraft, the skin of wind turbine blades, and vehicle components. These composite structures undergo static and fatigue loading throughout their service life. Under such loading, various types of inter-and intralaminar fatigue damage are predominantly caused by cyclic loading, especially due to the low strength and toughness of the different constituents of laminated composites, affecting the reliability of FRP composite structures [1,2]. As a result, the development of fatigue damage models is valuable as a tool for complete insight into the mechanics of deformation and failure of the material, as well as for the reliability assessment of composite structures.
Composite structures as a part of vehicles' bodies are often designed to bear complex loads that appear in the form of lateral loads, bending, twisting, etc. [3,4]. Such loading causes various models of failure in laminated composites, including matrix yielding and cracking, fiber/matrix interface debonding, fiber pullout, fiber buckling, fracture, and interface delamination [5][6][7]. Multiple interface delamination is among the most frequently seen failure modes, often occurring under mixed-mode loading conditions. The constitutive damage models of interface failure are derived to consider damage modes I and II, along with their interaction; therefore, the knowledge of the mechanics and mechanisms of proposed elastic-plastic cohesion laws for polycrystalline low-cycle fatigue, taking the cohesive-frictional behavior into account [44].
In this study, a new cyclic cohesive zone model (CCZM) was developed and employed for the reliability prediction of FRP composite laminates under mode I crack loading conditions. The theoretical damage model of the interface was developed by combining two models, i.e., (1) a property degradation model, and (2) a fatigue life model. The models encompass the interface's typical fatigue life and recognize the fatigue degradation of interlaminar strength and stiffness properties during the damage progression to the onset of crack nucleation. The dissipation of fracture energy governs the separation of the material point in the interface. The model was embedded in a finite element analysis (FEA) package using the subroutine UMAT and used to investigate the failure and damage during mode I fatigue loading of CFRP composite laminates. A comparison of the model's FEpredicted data and experimental results (from the literature) showed a successful validation of the model.

Cohesive Zone Model for Static Loading
The interfaces between two laminas in an FRP composite laminate are assumed to have cohesive behavior. In addition, the interface is relatively thin, such that only the out-of-plane normal and in-plane shear stress components act on the mathematically "zerothickness" layer. Bilinear traction-displacement softening law is assumed for the cohesive interface under the quasi-static loading, as schematically shown in Figure 1. This cohesive zone model (CZM), describing the criterion for the onset of interlaminar damage, is given as follows [18]: where σ 33 is the tensile normal stress and (τ 31 , τ 32 ) are the components of shear stress on the interface constituent. The interlaminar normal tensile and shear strength are represented by T 0 and S 0 , respectively. The normal stress as a result of compressive load, is not affecting the interlaminar damage process. In addition, the shear stress terms diminish under the mode I crack loading. The subsequent damage evolution process causing the separation of the interface material point is governed by the dissipation of interlaminar strain energy. The fracture energy under individual mode I is given by the area below the traction-relative displacement curve, as follows: where T 0 is the interface normal strength and δ n f is the relative displacement of the interface at fracture in mode I. Similarly, for isolated mode II, the fracture energy is given as follows: where S 0 is the interface shear strength and δ S f is the relative displacement at fracture in mode II. The strain energy release rate G T under the mixed-mode (modes I and II) crack loading condition is described as follows [45]: where G IC and G IIC are critical strain energy release rates in mode I and mode II loading, respectively. The exponent η represents the degree of interaction of each loading mode, taken as 1.45. G II and G I represent the energy release rates, and their ratio in Equation (4) indicates the participation of individual modes. In the isolated mode I loading condition, G I I = 0 and the equation is reduced to G T = G IC . In the isolated mode II crack loading conditions, G I = 0 and the equation is reduced to G T = G I IC . The critical material point on the interface separates when the total fracture energy is dissipated.
where GIC and GIIC are critical strain energy release rates in mode I and mode II loading, respectively. The exponent η represents the degree of interaction of each loading mode, taken as 1. 45. GII and GI represent the energy release rates, and their ratio in Equation (4) indicates the participation of individual modes. In the isolated mode I loading condition, = 0 and the equation is reduced to = . In the isolated mode II crack loading conditions, = 0 and the equation is reduced to = . The critical material point on the interface separates when the total fracture energy is dissipated. The CZM properties' parameters were determined through the isolated mode testing of CFRP composite laminate specimens. At the damage initiation point, the local maximum stress attained is referred to as the interface strength. A validated hybrid FEexperimental approach was employed to extract the CZM parameters. The hybrid FEexperimental approach is detailed in the published literature [23]. Table 1 lists the resulting interface properties.  The CZM properties' parameters were determined through the isolated mode testing of CFRP composite laminate specimens. At the damage initiation point, the local maximum stress attained is referred to as the interface strength. A validated hybrid FE-experimental approach was employed to extract the CZM parameters. The hybrid FE-experimental approach is detailed in the published literature [23]. Table 1 lists the resulting interface properties.

Interlaminar Fatigue Damage Model
Interlaminar strength and stiffness degradation, the fatigue life of the interface, and characteristic damage evolution under fatigue loading conditions were employed to quantify the process of the mechanics of interlaminar fatigue damage. Figure 2 illustrates the mechanics of the interlaminar fatigue damage model through the evolution of normal traction-relative displacement curves at the critical interface material point experiencing fluctuating mode I crack loading. The interface was loaded with the maximum traction of σ max due to the application of constant-amplitude external loading with a load ratio (κ > 0). The load ratio κ is the ratio of minimum to maximum stress. The interlaminar fatigue damage process-from pristine condition to full separation-of any material point of interface is composed of two stages: The first stage of fatigue damage is caused by the degradation of the interface properties due to fatigue loading. Curves a-b-c-d indicate this fatigue damage stage as the fatigue load cycles advance. The end of this curve, point d, represents the onset of the crack nucleation event. The stage of fatigue damage evolution follows the first step of damage. This second stage of damage is governed by fracture energy dissipation (∆0-d-g) that starts at point d (nucleation onset) and ends at point g (separation of the material point). This implies that the governing damage mechanisms for both stages of damage are different and, therefore, should be treated separately. As soon as the material point on the interface separates, the local stress reduces to zero and the load is redistributed to nearby elements to maintain the forces' equilibrium. All of the subsequent material points on the interface experience a similar fatigue damage process. These separated material points collectively create a fatigue crack that advances with the increasing number of load cycles until a final fracture of the whole interface.

Interlaminar Fatigue Damage Model
Interlaminar strength and stiffness degradation, the fatigue life of the interface, and characteristic damage evolution under fatigue loading conditions were employed to quantify the process of the mechanics of interlaminar fatigue damage. Figure 2 illustrates the mechanics of the interlaminar fatigue damage model through the evolution of normal traction-relative displacement curves at the critical interface material point experiencing fluctuating mode I crack loading. The interface was loaded with the maximum traction of due to the application of constant-amplitude external loading with a load ratio ( > 0). The load ratio is the ratio of minimum to maximum stress. The interlaminar fatigue damage process-from pristine condition to full separation-of any material point of interface is composed of two stages: The first stage of fatigue damage is caused by the degradation of the interface properties due to fatigue loading. Curves a-b-c-d indicate this fatigue damage stage as the fatigue load cycles advance. The end of this curve, point d, represents the onset of the crack nucleation event. The stage of fatigue damage evolution follows the first step of damage. This second stage of damage is governed by fracture energy dissipation (Δ0-d-g) that starts at point d (nucleation onset) and ends at point g (separation of the material point). This implies that the governing damage mechanisms for both stages of damage are different and, therefore, should be treated separately. As soon as the material point on the interface separates, the local stress reduces to zero and the load is redistributed to nearby elements to maintain the forces' equilibrium. All of the subsequent material points on the interface experience a similar fatigue damage process. These separated material points collectively create a fatigue crack that advances with the increasing number of load cycles until a final fracture of the whole interface.

Cyclic Cohesive Zone Model for Mode I Crack Loading
The hypothesis of the fatigue damage model is based on the gradual deterioration of the properties of the interface (i.e., normal strength 0 , fracture energy 0 , and penalty stiffness 0 ). Curve 0ae is the reference bilinear softening response at the beginning of the load cycles or, alternatively, a quasi-static response. The interlaminar tensile strength 0 , penalty stiffness 0 , and fracture energy 0 degrade with the elapsed of the load cycles (n1), to new values of ( 1 ), ( 1 ), and ( 1 ) signifying the accu-

Cyclic Cohesive Zone Model for Mode I Crack Loading
The hypothesis of the fatigue damage model is based on the gradual deterioration of the properties of the interface (i.e., normal strength T 0 , fracture energy G 0 IC , and penalty stiffness k n 0 ). Curve 0ae is the reference bilinear softening response at the beginning of the load cycles or, alternatively, a quasi-static response. The interlaminar tensile strength T 0 , penalty stiffness k n 0 , and fracture energy G 0 IC degrade with the elapsed of the load cycles (n 1 ), to new values of T(n 1 ), k n (n 1 ), and G IC (n 1 ) signifying the accumulation of fatigue damage. The initial critical strain energy release rate G 0 IC represented by the area under the ∆0ae is reduced to the area beneath the ∆0bf. Further fatigue load cycles (n 2 ) would continue the degradation of these properties to T(n 2 ), k n (n 2 ), and G IC (n 2 ), while the apex of the traction-relative displacement curve follows the path a-b-c-d. The interlaminar tensile strength could only degrade to the maximum stress level σ max of the load cycles (point d), where the interface crack nucleation process begins. The evolution of interlaminar fatigue Polymers 2023, 15, 527 6 of 17 damage due to property degradation D f , which follows the path a-b-c-d from the pristine condition (D f = 0) to the onset of the crack nucleation (D f = 1.0), is given as follows: The variables T(n) and S(n) are the residual tensile and shear strength, respectively, corresponding to n elapsed load cycles. At this stage (point d), the damage due to fracture energy dissipation D e = 0. Following the onset of the interface crack nucleation process, the fatigue damage due to fracture energy dissipation D e under mode I crack loading is governed by the dissipation of the residual interlaminar strain energy release rate G I (n).
Separation of the critical interface material point occurs when G I (n) reduces to zero or D e = 1.0. This signifies the nucleation of the interlaminar fatigue crack. Adjacent separated material points collectively form the structural interface crack. This enables the interface crack propagation process in the cohesive interface plane to be simulated.

Interlaminar Property Degradation Model
The interrupted fatigue tests under different applied stress conditions (κ, σ max ) and the number of accumulated stress cycles n was used to establish the degradation of the interface properties. These properties include normal stiffness k n 0 , fracture energy G 0 IC , and tensile strength T 0 for the interface. The tests were carried out on CFRP composite laminates with DCB ENF specimen geometry under specific fatigue loading conditions (κ, σ max , n). These pre-fatigued specimens were then loaded in the tensile testing machine for a quasi-static loading until fracture, and the load-displacement curves were documented. A validated hybrid FE-experimental technique was then employed for each test specimen to establish the normal penalty stiffnesses and interface strength. Experimental load and displacement were used for the calculation of fracture energy for the isolated loading mode. A hybrid FE-experimental technique has already been published [23] that guides the extraction of interface properties of FRP composite laminates under damage due to cyclic loads. It should be noted that throughout the interface's fatigue life, the properties degrade in a similar manner. As a result, the degraded strength property data can be stated in their normalized form as follows [48]: By having static strength T 0 , maximum applied stress σ max , the number of cycles to failure N f (related to the specific loading state), and curve fitting parameters α and β, the residual interface strength at any number of cycles T(n) can be calculated. By executing a similar process, the normalized residual penalty stiffness can be calculated as follows [48]: where k n (n) is residual interface stiffness, k n 0 is interface stiffness for quasi-static cases, δ n 0 is the relative displacement of the interface at damage onset, and γ and λ are curve-fitting parameters.
Similarly, the normalized residual fracture energy can be calculated as follows [48]: where δ n f is the relative displacement at fracture, µ, and φ are curve fitting parameters. The outcome of the interface normalized property model for CFRP composite laminates is displayed in Figure 3a-c. Different curve fitting parameters for the above-mentioned models are listed in Table 2.

Interlaminar Fatigue Life Model
The fatigue life of an FRP composite laminate's interface is affected by the mean stress of the loading cycle and is considered in this model through a series of fatigue life tests on FRP specimens. Double-cantilever beam (DCB) specimens were used until fracture for mode I loading [21]. Load-control-based fatigue testing at various stress levels (κ, σ max ) was conducted on end-notched flexure (ENF) specimens until fracture. For mode I, identical test data for fatigue life was acquired from the published literature [34]. The life parameter χ is used to include the effect of mean stress in the model, as follows [3,49]: where A and B are curve-fitting constants. The constant c is optimized such that all of the data points are aligned with the best fit in a straight line when presented in the χ − log 10 N f plot, as shown in Figure 4 with c = 1.939. Such a linear fit to the test data is valuable for computational life prediction exercises. By using Figure 4, the number of cycles to failure N f can be calculated for any stress ratio. The fatigue life of an FRP composite laminate's interface is affected by the mean stress of the loading cycle and is considered in this model through a series of fatigue life tests on FRP specimens. Double-cantilever beam (DCB) specimens were used until fracture for mode I loading [21]. Load-control-based fatigue testing at various stress levels ( , ) was conducted on end-notched flexure (ENF) specimens until fracture. For mode I, identical test data for fatigue life was acquired from the published literature [34]. The life parameter χ is used to include the effect of mean stress in the model, as follows [3,49]: where A and B are curve-fitting constants. The constant c is optimized such that all of the data points are aligned with the best fit in a straight line when presented in the χ − log10 Nf plot, as shown in Figure 4 with c = 1.939. Such a linear fit to the test data is valuable for computational life prediction exercises. By using Figure 4, the number of cycles to failure Nf can be calculated for any stress ratio. The model was coded in Fortran language and integrated with the FE package [24] through a user-written subroutine (UMAT). The flowchart of the subroutine is shown in Figure 5. The model was coded in Fortran language and integrated with the FE package [24] through a user-written subroutine (UMAT). The flowchart of the subroutine is shown in Figure 5.

Fatigue Failure Process of Mode I Interface Loading
The interface fatigue cracking process of the CFRP composite laminates under mode I loading was demonstrated through an FE simulation case study. Necessary conditions for simulating the response of a CFRP composite laminate beam with the double-cantilever beam (DCB) are discussed, including model geometry, boundary conditions, applied load cycles, and mesh convergence analysis. The results are presented here and discussed with respect to the damage evolution of the interface due to fatigue loading, the stress distribution of the interface, and the number of cycles to failure N f . FE Simulation Case Study The DCB specimen was simulated under mode I fatigue loading conditions. The geometry and test setup of the DCB specimen are shown in Figure 6. The specimen's nominal dimensions (in mm) were length L = 100, width B = 21, thickness 2h = 9.6, and a o = 45.

Fatigue Failure Process of Mode I Interface Loading
The interface fatigue cracking process of the CFRP composite laminates unde I loading was demonstrated through an FE simulation case study. Necessary con for simulating the response of a CFRP composite laminate beam with th ble-cantilever beam (DCB) are discussed, including model geometry, boundary tions, applied load cycles, and mesh convergence analysis. The results are present and discussed with respect to the damage evolution of the interface due to fatigu ing, the stress distribution of the interface, and the number of cycles to failure Nf.

FE Simulation Case Study
The DCB specimen was simulated under mode I fatigue loading conditio geometry and test setup of the DCB specimen are shown in Figure 6. The spe nominal dimensions (in mm) were length L = 100, width B = 21, thickness 2h = 9.6, 45. An initial interface crack was introduced in the middle plane by overlying without applying a bond between them to create a traction-free crack surface at ginning of the simulation. The reference static material properties of the CFRP com laminas and the critical interface of the DCB specimen employed in this simula shown in Table 1. The pair of loading blocks attached to the specimen were assu be rigid bodies. Figure 7 shows the specimen's mesh, boundary conditions, and geometry discretized into 3D solid elements. The lower loading block was con non-moveable (Ux = Uy = Uz = 0). The load cycles were applied to the top loadin thus, it was moveable in the z-direction (Ux = Uy = 0). Both of the loading block rotate about the y-axis ( = = 0). The load cycle blocks simulated the app tigue loading with Pmax = 148 N and load ratio = 0.11. An initial loading step fined to bring the load to the minimum level (Pmin = 28 N) at the beginning of the simulation. The DCB specimen was modeled with 32 layers of laminas above the plane and 32 layers below the middle plane interface containing the initial crack cohesive layer (a total of 64 layers). The lamina was discretized into a total o eight-node continuum shell elements (Abaqus SC8R element). This implies that 3 of the half-laminate are defined in 4 layers of SC8R elements. A layer of 8000 eig cohesive elements (Abaqus COH3D8 element) with matching nodes to the a laminar surfaces is therefore prescribed. An initial interface crack was introduced in the middle plane by overlying nodes without applying a bond between them to create a traction-free crack surface at the beginning of the simulation. The reference static material properties of the CFRP composite laminas and the critical interface of the DCB specimen employed in this simulation are shown in Table 1. The pair of loading blocks attached to the specimen were assumed to be rigid bodies. Figure 7 shows the specimen's mesh, boundary conditions, and model geometry discretized into 3D solid elements. The lower loading block was considered non-moveable (Ux = Uy = Uz = 0). The load cycles were applied to the top loading block; thus, it was moveable in the z-direction (Ux = Uy = 0). Both of the loading blocks could rotate about the y-axis (UR x = UR z = 0). The load cycle blocks simulated the applied fatigue loading with P max = 148 N and load ratio κ = 0.11. An initial loading step was defined to bring the load to the minimum level (P min = 28 N) at the beginning of the fatigue simulation. The DCB specimen was modeled with 32 layers of laminas above the middle plane and 32 layers below the middle plane interface containing the initial crack and the cohesive layer (a total of 64 layers). The lamina was discretized into a total of 31,552 eight-node continuum shell elements (Abaqus SC8R element). This implies that 32 layers of the halflaminate are defined in 4 layers of SC8R elements. A layer of 8000 eight-node cohesive elements (Abaqus COH3D8 element) with matching nodes to the adjacent laminar surfaces is therefore prescribed.
The mesh convergence analysis is shown in Figure 8. This was performed for the interface crack front region to ensure that the largest element employed had an insignificant influence on the FE-calculated variables. The interlaminar tensile stress σ 33 was used as the monitoring variable in the mesh convergence analysis. In this respect, the FE model of the interface region ahead of the crack front was discretized into cohesive elements, each with an edge length of 0.1 mm along the crack growth direction. The mesh convergence analysis is shown in Figure 8. This was performed for the interface crack front region to ensure that the largest element employed had an insignificant influence on the FE-calculated variables. The interlaminar tensile stress σ33 was used as the monitoring variable in the mesh convergence analysis. In this respect, the FE model of the interface region ahead of the crack front was discretized into cohesive elements, each with an edge length of 0.1 mm along the crack growth direction.

Results and Discussion
The focus of this study was to develop a model to capture the fatigue failure process in FRP composite laminates. The results of the case study are presented and discussed in the following sections.

Stress Distribution of the Interface
Fatigue damage accumulates due to interlaminar property degradation before the onset of interface crack nucleation. However, the stress experienced by the degrading interface material points remains unchanged. A typical interlaminar stress field corresponding to the peak applied load cycle is shown in Figure 9. The loading of the DCB  The mesh convergence analysis is shown in Figure 8. This was performed for the interface crack front region to ensure that the largest element employed had an insignificant influence on the FE-calculated variables. The interlaminar tensile stress σ33 was used as the monitoring variable in the mesh convergence analysis. In this respect, the FE model of the interface region ahead of the crack front was discretized into cohesive elements, each with an edge length of 0.1 mm along the crack growth direction.

Results and Discussion
The focus of this study was to develop a model to capture the fatigue failure process in FRP composite laminates. The results of the case study are presented and discussed in the following sections.

Stress Distribution of the Interface
Fatigue damage accumulates due to interlaminar property degradation before the onset of interface crack nucleation. However, the stress experienced by the degrading interface material points remains unchanged. A typical interlaminar stress field corresponding to the peak applied load cycle is shown in Figure 9. The loading of the DCB

Results and Discussion
The focus of this study was to develop a model to capture the fatigue failure process in FRP composite laminates. The results of the case study are presented and discussed in the following sections.

Stress Distribution of the Interface
Fatigue damage accumulates due to interlaminar property degradation before the onset of interface crack nucleation. However, the stress experienced by the degrading interface material points remains unchanged. A typical interlaminar stress field corresponding to the peak applied load cycle is shown in Figure 9. The loading of the DCB beam specimen induced normal stress with a peak at 49.9 MPa. The interlaminar stress field was highly concentrated in the vicinity of the crack tip, as expected. The corresponding normal strain magnitude at the peak stress cycle was 5.12 × 10 −5 mm. Similar stress contours but at lower magnitudes were predicted following the start of the crack nucleation process, as governed by the strain energy release rate of the interface. beam specimen induced normal stress with a peak at 49.9 MPa. The interlaminar stress field was highly concentrated in the vicinity of the crack tip, as expected. The corresponding normal strain magnitude at the peak stress cycle was 5.12 × 10 −5 mm. Similar stress contours but at lower magnitudes were predicted following the start of the crack nucleation process, as governed by the strain energy release rate of the interface. Figure 9. The interlaminar normal stress field in the vicinity of the starter crack tip corresponds to the peak applied load cycle at the start of the fatigue loading.
The stress field at the crack front was analyzed further by plotting the variation in the normal stress of the interface against the true distance along the length of the interface. This effect is shown in Figure 10. The vertical axis is the stress magnitude normalized by the normal static strength of the interface. It can thus be concluded that the normal stress on the interface is highly concentrated at the crack front while sharply declining in magnitude after the length of a few elements.  The stress field at the crack front was analyzed further by plotting the variation in the normal stress of the interface against the true distance along the length of the interface. This effect is shown in Figure 10. The vertical axis is the stress magnitude normalized by the normal static strength of the interface. It can thus be concluded that the normal stress on the interface is highly concentrated at the crack front while sharply declining in magnitude after the length of a few elements.
beam specimen induced normal stress with a peak at 49.9 MPa. The interlaminar stress field was highly concentrated in the vicinity of the crack tip, as expected. The corresponding normal strain magnitude at the peak stress cycle was 5.12 × 10 −5 mm. Similar stress contours but at lower magnitudes were predicted following the start of the crack nucleation process, as governed by the strain energy release rate of the interface. The stress field at the crack front was analyzed further by plotting the variation in the normal stress of the interface against the true distance along the length of the interface. This effect is shown in Figure 10. The vertical axis is the stress magnitude normalized by the normal static strength of the interface. It can thus be concluded that the normal stress on the interface is highly concentrated at the crack front while sharply declining in magnitude after the length of a few elements.

Evolution of the Interlaminar Fatigue Damage
In this study, the interlaminar fatigue damage process was divided into two stages: A damage variable D f was used to calculate the damage until the onset of the crack nucleation in the first stage. From the onset of crack nucleation to the final separation of the material point, another damage variable D e was computed. The distribution of these damage variables is shown in Figure 11 for the different numbers of cycles. The magnitude of both variables exceeding 0.99 is assumed to satisfy the failure criteria. Figure 11a shows the first event of damage due to property degradation D f for the first row of elements at the crack front at almost 8200 cycles. This marks the crack nucleation onset. Figure 11b shows the separation of these damaged elements at almost 10,000 cycles, forming the first crack increment where the local stress was diminished to zero. Additional fatigue cycles forced the neighboring row of cohesive elements to separate (D e = 1). Beyond almost 14,100 cycles, the process of separation of elements became faster, with two or more elements separating in each load block increment. Figure 11c shows the contour of damage at almost 14,100 cycles. The combination of all of these separating elements forms the fatigue crack growth stage. After this many load cycles, the crack advances at a higher rate and forms a structural interface crack. Figure 11d shows the larger interface crack increment after enduring almost 2000 more fatigue load cycles.

Conclusions
In this study, a fatigue damage model was developed to predict and analyze the interface damage initiation and propagation. The fatigue model predicts the damage due to the degradation of interface properties up to fatigue crack nucleation onset. Fracture energy is the governing mechanism in this model for the calculation of the subsequent separation of the interface material point at the interface. As a case study, this model was used to predict the fatigue failure behavior of CFRP laminated composite subjected to mode I fatigue load and boundary conditions. A 3D simulation of the model was carried out based on an experimentally tested CFRP DCB sample, and the FE results were compared with experimentally measured data, in which a good agreement is observed, indicating the successful validation of the model and the simulation process. The interlaminar fatigue damage model was used to illustrate the evolution of stress and quantify the fatigue crack growth in the DCB specimen. The selective results of the FE model for describing the fatigue damage evolution are as follows: • The model successfully confined the high-stress gradient at the interface crack front region, with a normal tensile stress level of 70% of the respective interlaminar strength.

•
The first event of onset of nucleation at the material point in the crack front occurred at almost 8200 cycles. • The first increase in crack growth (i.e., the first row of elements separation) occurred at almost 10,000 cycles.

•
After almost 14,000 load cycles, the crack advanced at a higher rate and formed a structural interface crack.

•
A similar damage evolution process as predicted by the FE model was observed in the experimental case.

•
As indicated by the results, it can be concluded that interface properties' degradation and the dissipation of fracture energy are appropriate physical properties to be employed for interlaminar fatigue damage modeling of FRP composite laminates. Future research should aim at using a similar model for the fatigue analysis of FRP composites under mode II and mixed-mode I/II loading conditions. Another important future research direction will be to investigate the full fatigue life of FRP composite structures where fatigue damage occurs in both lamina and interface constituents, where the current interlaminar model may be integrated with a laminar fatigue damage model to predict the life of the composite structure.

Conclusions
In this study, a fatigue damage model was developed to predict and analyze the interface damage initiation and propagation. The fatigue model predicts the damage due to the degradation of interface properties up to the onset of fatigue crack nucleation. Fracture energy is the governing mechanism in this model for the calculation of the subsequent fatigue fracture of the interface material point. As a case study, this model was used to predict the fatigue failure behavior of CFRP laminated composite subjected to mode I fatigue load and boundary conditions. A 3D simulation of the model was carried out based on an experimentally tested CFRP DCB sample, and the FE results were compared with experimentally measured data, in which a good agreement is observed, indicating the successful validation of the model and the simulation process. The interlaminar fatigue damage model was used to illustrate the evolution of stress and quantify the fatigue crack growth in the DCB specimen. The selective results of the FE model for describing the fatigue damage evolution are as follows:

•
The model successfully confined the high-stress gradient at the interface crack front region, with a normal tensile stress level of 70% of the respective interlaminar strength.

•
The first event of onset of nucleation at the material point in the crack front occurred at almost 8200 cycles. • The first increase in crack growth (i.e., the first row of elements separation) occurred at almost 10,000 cycles. • After almost 14,000 load cycles, the crack advanced at a higher rate and formed a structural interface crack.

•
A similar damage evolution process as predicted by the FE model was observed in the experimental case. • As indicated by the results, it can be concluded that interface properties' degradation and the dissipation of fracture energy are appropriate physical properties to be employed for interlaminar fatigue damage modeling of FRP composite laminates. Future research should aim at using a similar model for the fatigue analysis of FRP composites under mixed-mode I/II loading conditions. Another important future research direction will be to investigate the full fatigue life of FRP composite structures where fatigue damage occurs in both lamina and interface constituents, where the