Progressive Fatigue Modelling of Open-Hole Glass-Fibre Epoxy Laminates

Abstract

The failure of composite laminates under cyclic fatigue loads is complex, as multiple failure mechanisms are in play at different scales and interact with each other. Predicting the remaining fatigue life as well as the residual capacities of a composite laminate or component is crucial, particularly for safety-critical applications. A progressive fatigue model is proposed to describe the catastrophic failure of open-hole laminates under tensile cyclic fatigue. To represent both intra-laminar and inter-laminar damage, a combination of a continuum damage mechanics model (CDM) and a discrete cohesive zone model (CZM) is implemented in the finite element (FE) software Abaqus. The CDM combines fibre- and matrix-dominated S-N curves with the Palmgren–Miner accumulation rule and Hashin’s residual strength to form a fatigue failure criterion differentiating between fibre failure (FF) and matrix failure (MF). The CZM implemented in this work is the CF20 model proposed by NASA. Fatigue cycling is simulated using an external cycle-jump scheme, where the stiffness degradation is conducted between the FE simulations outside of the implicit solver [90/0] s. Glass fibre reinforced polymer (GFRP) open-hole specimens were tested in tensile cyclic fatigue at a load ratio of 0.1. The experiments were reproduced numerically and the results compared. After calibration of a set of parameters based on one load level, the model was able to reproduce the experimental S-N curve very well, predicting a slope of −0.10, while the experimental value was −0.11. The failure sequence of the laminate was also successfully reproduced. The growth of the split from the hole, and its interaction with inter-laminar delamination, was successfully captured. The proposed approach was able to describe the fatigue failure of an open-hole laminate with a minimal set of material inputs using a simplified fatigue damage model while avoiding convergence issues.

1. Introduction

The failure of composite laminates under fatigue loads is a complex phenomenon, not so far well described. Multiple failure mechanisms are in play, at different scales, growing at different rates and interacting with each other. But determining the remaining fatigue life and the residual properties of a laminate is a crucial question, especially for safety-critical applications with severe cyclic loads, like wind turbine blades, offshore structures, or more recently, offshore composite pipelines. The industry needs reliable predictive tools to calculate the fatigue life and to evaluate the residual capacity of the composite components with complex geometries. The end of the useful life may be determined by a certain change in stiffness or a reduction in strength, depending on the application. This work is related to pressure-containing equipment where strength is critical and stiffness changes can be tolerated.

The different approaches to describe the fatigue of composites can be classified in three categories [1]: fatigue-life models, phenomenological models, and progressive fatigue models (PFMs). In the present paper, the term static is used to refer to quasi-static loading conditions, and fatigue refers to cyclic fatigue loading conditions. The term lifetime refers to the cycling stage, before the failure of the laminate. The term fatigue life refers to the duration of the specimen’s life, i.e., the number of cycles at failure. Fatigue-life models use S-N curves and constant life diagrams to predict the failure of a laminate under the relevant load conditions. No information is given on the evolution of the properties along the lifetime. The damage mechanisms and the change in the mechanical response are not described either. While providing limited insights, these models are simple, reliable, and are commonly used as design rules [2]. Typical fatigue-life models are based on S-N curves, the Goodman diagram, and the Palmgren–Miner accumulation rule [3,4,5,6]. Recently, some commercial programmes have integrated composite fatigue analysis tools based on fatigue-life models. For example, FEMFAT Laminate [7,8,9] offers fast fatigue analyses for composite structures. It is based on a static stress-field calculation and is limited to the prediction of the fatigue life.

Phenomenological models describe the evolution of macroscopic mechanical properties along the lifetime, for example, the residual strength and the residual stiffness. A damage variable is generally used to describe the current residual properties at one stage of the lifetime. The final failure can be linked to this damage variable, which can also be correlated with the material S-N curve. Hashin proposed an analytical residual strength model, using the Palmgren–Miner sum as a damage variable and describing the sudden-death behaviour of the material [10]. Other models are constructed by fitting experimental data on the residual strength or the residual stiffness with mathematical expressions in order to describe the damage growth rate [11,12,13,14]. These models describe the evolution of macroscopic properties quite well and are easily implemented. However, they usually require a lot of testing to account for the different layups and loading conditions. Also, they do not describe the actual physical failure mechanisms at play in the laminate.

Progressive failure models aim to describe damage and failure arising from the physical failure mechanisms in the material. Fibre failure, matrix failure, and delamination can, for example, be modelled. Typically, a failure criterion on the stress state will govern the onset of a failure mechanism and its evolution can be modelled by a softening law until final failure. Continuum damage mechanics models (CDMs) are widely used to describe intra-laminar failure mechanisms in a continuous and homogenised way [15,16,17]. Cohesive zone models (CZMs) have proven powerful and convenient to model inter-laminar failure, such as delamination or macroscopic cracks [18,19]. CZMs are implemented in discrete predefined locations where cohesive fractures are expected or observed. The XFEM approach can be used to model failure and crack growth within a material without a predefined path [20,21]. Progressive models can be applied both to static and fatigue problems. In the latter case, they are called progressive fatigue models (PFMs). For fatigue load cases, PFMs capture the redistribution of stresses in laminates along the lifetime. This is particularly important for laminates with geometrical singularities (for example, a hole) where the stress field is not homogeneous. In such cases, the failure mechanisms are not spatially even and the redistribution of stresses is significant. This affects the future growth of the damage and the component’s fatigue life. PFMs implemented in a finite element (FE) scheme have the potential to describe the lifetime of such laminates with a high level of fidelity. It is, thus, the approach chosen in this work.

Extensive work has been undertaken developing progressive failure models to describe or predict the failure of laminates under static loads, using CDMs [22,23,24] and CZMs [19,25]. Recently, Nguyen and Waas proposed an innovative semi-discrete framework to accurately model both intra-laminar and inter-laminar failure under static loads [26]. Fibre-failure, including non-linearity, is described by bulk elements. Splitting-bands with a novel mixed-mode law for matrix damage are inserted parallel to the fibres to model intra-laminar failure. Cohesive elements are inserted between plies to model inter-laminar failure. A statistical description of the material properties was integrated to account for the scatter in the material properties. CZMs were used successfully to model the interactions between intra-laminar transverse matrix cracks, inter-laminar delamination, and intra-laminar longitudinal matrix splits in [0/90] open-hole laminates under static load [27]. Static CDMs and CZMs can be combined to model the different failure mechanisms in play in the laminate and their interactions [23,28,29,30].

When it comes to fatigue, many CDMs are already available, usually in the form of residual strength models corresponding to each failure mode [13,14,31]. Recently, several studies have proposed approaches to implement them in FEA, illustrating the growing interest in this kind of model.

Koch et al. used a PFM to describe the catastrophic failure of a carbon-fibre component under cyclic loads [32]. The CDM was based on damage growth rates measured on test coupons for the three failure mechanisms modelled (fibre failure, and transverse and shear matrix failures). From this growth rate, the damage can be predicted along the lifetime. A cycle-jump scheme similar to the one originally presented by Van Paepegem [33] was used. The damage growth rate is calculated based on the peak stresses seen during a loading step. Then, the damage is extrapolated over a jump of cycles and a new loading step is conducted on the newly damaged structure. To ensure a good description of the damage growth, the maximum acceptable jump size is controlled by a criterion.

Samareh-Mousavi et al. proposed another PFM based on a CDM [34]. The emphasis was put on the non-linear stress–strain material response and its influence on the fatigue failure criterion. Hashin’s static failure criterion was used [35]. To minimise the number of parameters needed, the fatigue damage was described in terms of the residual stiffness for the transverse plies and in terms of the residual strength for the longitudinal ones. This allowed a relevant description of the fatigue damage at the laminate scale, while optimizing the number of experiments and inputs needed. The authors successfully applied this model to a pin-loaded laminate in a 2D case, describing both the residual stiffness evolution and the final failure with good agreement.

Gerendt et al. presented an advanced PFM using Puck’s 3D failure criterion and a fracture energy equivalence proposed by Pfanner [36] to evaluate the damage evolution in fatigue [37,38,39]. Coupling this hypothesis with damage evolution functions determined experimentally, they were able to describe very well both the static and the progressive fatigue failure of pin-loaded laminates, using the plies’ residual stiffness functions to quantify local progression of the fibre and matrix damage.

Similar approaches were also proposed in [40,41,42,43,44]. While their descriptive and predictive capacities are satisfying, the approaches typically rely on large sets of material parameters for the residual stiffness or strength models within the CDMs. Also, they usually do not feature any CZM and are limited in terms of inter-laminar failure.

Many models are available to describe intra- and inter-laminar matrix damage with CZM models. Bi-linear softening laws were proposed to model mode I, mode II, and mixed-mode delamination [25], while non-linear empirical laws were developed to account for non-linear effects, for example, in the presence of fibre-bridging or for toughened laminates [45,46,47]. Davila recently proposed two models for cohesive fatigue damage growth based on Turon’s cohesive law, which were able to reproduce Paris Law curves in the case of double cantilever beam (DCB) and mixed-mode bending (MMB) tests with low mode mixity [48,49]. The fatigue damage increment was calculated from an S-N curve of the material in the transverse direction. Harper and Hallett developed a similar model where the fatigue damage was calculated directly from three Paris Law curves (one in each mode) [50]. The model was in good agreement with crack growth experiments on test coupons.

These fatigue models can be combined within the same PFM framework to simultaneously model intra- and inter-laminar cohesive failure, along with intralaminar continuous failure mechanisms, with flexibility and accuracy. For example, the CompDam tool from NASA [23] features several CDMs for the failure of the matrix and the fibres, and the CF20 CZM from Davila is integrated to model inter- and intra-laminar matrix failure in fatigue [48]. A limitation of CompDam is the lack of a model describing the fibre failure under cyclic fatigue loads. Thus, it focuses mainly on stiffness reduction prediction rather than catastrophic failure.

Two other proposed models involve a combination of CDM and CZM under fatigue loading. Aoki et al. implemented a fatigue CDM from Ladevèze [51,52] in combination with a CZM from Yashiro [27,53]. The CZM was limited to inter-laminar delamination. The model was applied to a curved geometry and the residual stiffness and fatigue life was describe with good agreement. Hugaas developed a similar approach to describe the fatigue failure of split disks with a hole [54]. In his work, a discrete CZM layer was also inserted to model an intra-laminar split in the load-bearing hoop plies.

All the PFMs presented here rely on empirical damage evolution laws (e.g., residual strength, residual stiffness, or general fatigue damage growth models, usually at the ply level) [55], which need to be obtained experimentally. They focus on matrix-dominated failure and stiffness loss along the fatigue lifetime. This makes it easier to implement the PFM in implicit FE solvers, where the convergence issues are limited thanks to the smooth stiffness evolution resulting from the damage evolution laws. Such approaches are limited when attempting to predict the fatigue life in terms of the residual strength, and especially catastrophic fibre-dominated failure. This catastrophic failure implies a more abrupt stiffness loss and more load redistribution than matrix-dominated failure, leading to more convergence issues in classical implicit schemes.

This study attempts to reproduce the failure of a laminate with a hole, focusing on fibre-dominated catastrophic failure and fatigue life prediction, with a minimal amount of input parameters and a reasonable calculation time. In contrast to other models, a precise description of the stiffness evolution at the early stages of the fatigue life was not sought after here. This deliberate trade-off aims at keeping the method simple, accessible, and able to be implemented in use-cases where the residual strength and the fatigue life are more critical than the residual stiffness. Few proposed models combine a fatigue CDM and a discrete description of macroscopic matrix failures. But such macroscopic failures may drive the stress redistribution and impact the fatigue life. Thus, a coupled FE progressive fatigue model is proposed, where microscopic matrix failure and fibre failure are modelled using a fatigue CDM approach, while selected macroscopic matrix failures are accounted for with a discrete CZM at predefined locations. The CZM has the advantage of describing efficiently cohesive failures when their propagation path is known in advance [30]. As the addition of cohesive interfaces is computationally costly, the present study proposes the integration of selected cohesive layers where major delamination and splits are expected.

The CDM model uses S-N curves and the material strength to compute a fatigue failure criterion based on the Palmgren–Miner accumulation rule [5,6] and Hashin’s residual strength [10]. The Palmgren–Miner damage accumulation rule features linear fatigue damage growth, while Hashin’s residual strength exhibits a sudden-death behaviour. The use of these two simple approaches is justified by the focus of the study, which is the prediction of catastrophic failure and the fatigue life. They avoid the necessity of phenomenological fatigue damage growth characterisation, as classically performed in other studies. The CZM model implemented here is the CF20 fatigue crack growth model [48]. It relies on the inter-laminar strengths and energy release rates, as well as a matrix-dominated S-N curve, without requiring any Paris Law parameters. The model is implemented in an implicit solver and the fatigue analysis is conducted through a cycle-jump scheme. A novelty of the current implementation is the external cycle-jump scheme proposed to mitigate the convergence issues linked to fibre failure in the model. For the CDM, the fatigue damage calculation and the stiffness degradation are computed by an external Python script outside of the FE solver. The model was compared to tensile cyclic fatigue experiments on glass-fibre/epoxy specimens with a hole. The comparison focused on the description of the ultimate failure. The geometric singularity at the hole allowed evaluation of how well the FEA could cope with an uneven stress field.

2. Model Description

To describe the failure of a composite laminate under tensile cyclic fatigue loads, a combination of a CDM and a CZM is proposed.

The main failure mechanisms to take into account are:

• Microscopic matrix cracking
• Microscopic fibre failure
• Intra-laminar matrix cracking
• Inter-laminar matrix cracking
• Macroscopic fibre failure

Matrix cracking usually onsets at a microscopic level, between fibres. In a ply, the damage coalesces into intra-laminar cracks, which are macroscopic failures parallel to the fibres, usually through the entire thickness of a ply. In transverse plies, they are perpendicular to the main load direction. Thus, the term “transverse matrix crack” is used. Their main impact is a reduction of the shear and transverse load-bearing capacity in these plies. These cracks govern mainly the stiffness reduction of the laminate. In this approach, the focus is on the fibre-dominated catastrophic failure, so transverse matrix cracks are modelled in a homogenised way within a CDM. This might not be as accurate as a discrete or semi-discrete modelling, but the authors believe that it is a good compromise between fidelity and computational efficiency. In longitudinal plies, intra-laminar matrix cracks occur parallel to the main loading direction. Thus, the term “matrix split” is used. Their creation is shear-stress dominated in a mode II or mixed mode I–mode II crack opening. They split the longitudinal plies into separate sections. A split reduces the shear-based stress transfer and generates a significant stress redistribution within the longitudinal plies. Given the influence of the split on the stress field and the stress redistribution, it is important to model it with high fidelity. Thus, a discrete CZM model is used here.

Microscopic fibre failures occur early on in the lifetime, usually starting from fibre microscopic damage. They will also induce some microscopic matrix cracking. However, they are not the driving mechanism behind catastrophic failure. The stress redistribution in the plies due to all the matrix-related damage will generate stress redistribution and some stress concentrations at the scale of a couple of fibres. These stress concentrations coupled with variation in the fibre strengths will induce the first failure of the fibre clusters. These failures then coalesce, propagating to neighbouring clusters, finally leading to the catastrophic failure of the load bearing plies and failure of the laminate. Modelling this chain-mechanism precisely would be ideal but is difficult to achieve without a micro-mechanical model coupled with a statistical description of the fibre properties. In a simpler and faster way, a smeared CDM approach can give meaningful results by applying a failure criterion at the meso-scale without differentiating between microscopic and macroscopic failure mechanisms.

The model presented here is a coupled approach, where microscopic matrix cracks, macroscopic matrix cracks (in transverse plies), and fibre failures are described using a CDM, while matrix split (in the longitudinal plies) and inter-laminar delamination are modelled discretely using the CZM in the cohesive elements. The model principle is presented in Figure 1. The discrete CZM location is, thus, predefined. This limits the predictive capacity of the model, but represents an effective simplification when the location of the split and the delamination are known. For the laminate with a hole studied here, the expected locations of matrix splitting and delamination are known. This approach has already been applied successfully by Yashiro et al. [27]. For more complex cases, or if more fidelity is required, advanced meshing techniques, including numerous splits and transverse cracks in the mesh (similar to the approach proposed by Nguyen and Waas [26]), could be used.

2.1. Continuum Damage Mechanics

The CDM accounts for fibre and matrix failure in the bulk material in a continuous and homogenised way. The fatigue failure criterion is obtained by combining a maximum stress criterion with the Palmgren–Miner fatigue damage accumulation rule and the residual strength formulated by Hashin [10]. The maximum stress failure criterion is quite simple and is usually not used when a precise description of matrix-dominated failure is studied. Failure criteria such as those of Puck [56] or LaRC [57] may be preferred. However, they require additional material parameters not available for the current study. The adaptation of the maximum stress criterion in fatigue using Hashin’s residual strength is straightforward. This approach is chosen here.

The CDM model differentiates between fibre failure (FF) in the longitudinal direction $\left\{11\right\}$ and matrix failure (MF) in the other directions $\{22,33,12,13,23\}$. The maximum stress criterion in static is given by Equations (1) and (2):

• Fibre failure (FF):

  $$-X_C\le {\sigma }_{11}\le X_T$$

  (1)
• Matrix failure (MF):

  $$max(-\frac{{\sigma }_{22}}{Y_C},\frac{{\sigma }_{22}}{Y_T},-\frac{{\sigma }_{33}}{Z_C},\frac{{\sigma }_{33}}{Z_T},\frac{|{\sigma }_{12}|}{S_{12}},\frac{|{\sigma }_{13}|}{S_{13}},\frac{|{\sigma }_{23}|}{S_{23}})\le 1$$

  (2)

The fatigue failure criterion is obtained by substituting the static strength with Hashin’s fatigue residual strength [10], as described below. It is calculated using Palmgren–Miner’s accumulation damage variable M [5,6], based on the material’s S-N curve.

The fatigue life $N_{ij}$ of the material exposed to a given maximum stress ${\sigma }_{ij}$ in the $ij$-direction at a given load ratio R is given by Basquin’s S-N curve equation [3]:

$${\sigma }_{ij}={\sigma }_{ij}^0{N_{ij}}^{-a_{ij}}$$

(3)

where $-a_{ij}$ and ${\sigma }_{ij}^0$ are respectively the slope and intercept of the S-N curve in a log-log plot. $N_{ij}$ is the number of cycles to failure for the applied stress ${\sigma }_{ij}$.

This S-N curve formulation usually exceeds the material static strength in the low-cycle fatigue domain. As this is not physically possible, a cut-off is applied when the S-N curve overshoots the static strength. So, when the material is submitted to a maximum stress ${\sigma }_{ij}$, the number of cycles to failure is given by:

$$N_{ij}\left({\sigma }_{ij}\right)=\left\{\begin{array}{cc}{\frac{|{\sigma }_{ij}|}{{\sigma }_{ij}^0}}^{-\frac{1}{a_{ij}}}\hfill & \mathrm{if}|{\sigma }_{ij}|<X\hfill \\ 1\hfill & \mathrm{if}|{\sigma }_{ij}|\ge X\hfill \end{array}\right.$$

(4)

with X the material strength in the considered loading direction:

$$X\in \{X_T,X_C,Y_T,Y_C,Z_T,Z_C,S_{12},S_{13},S_{23}\}$$

(5)

During the lifetime, damage develops and the stresses in the material are redistributed. The stress redistribution leads to changes in the local maximum stresses seen at any point in the material. To account for this variation in the fatigue stress, the Palmgren–Miner damage accumulation rule is used [5,6]:

$$M_{ij}=\sum _k\frac{\Delta N_k}{N_{ij}\left({\sigma }_{ij,k}\right)}$$

(6)

where k indexes the different load blocks, $\Delta N_k$ is the number of cycles during which the material is exposed to a stress ${\sigma }_{ij,k}$, and $N_{ij}$ is the fatigue life predicted by Equation (4) at this stress level.

In each direction, the value of $M_{ij}$ quantifies the fatigue damage in the material direction $ij$. Hashin’s residual strength provides the updated material strength given this fatigue damage:

$$X_{res}=X{[1-M_{ij}]}^{a_{ij}}$$

(7)

Here, X and $X_{res}$ refer to the initial and residual material strength in any material direction and loading mode (traction or compression). This residual strength model exhibits a sudden-death behaviour–the strength remains quite close to the static value until the last fraction of the lifetime, where it drops rapidly. Under fatigue loading, the maximum stress failure criterion becomes:

Fibre failure in fatigue (FF):

$$-X_{res,C}\le {\sigma }_1\le X_{res,T}$$

(8)

Matrix failure in fatigue (MF):

$$max(-\frac{{\sigma }_{22}}{Y_{res,C}},\frac{{\sigma }_{22}}{Y_{res,T}},-\frac{{\sigma }_{33}}{Z_{res,C}},\frac{{\sigma }_{33}}{Z_{res,T}},\frac{|{\sigma }_{12}|}{S_{res,12}},\frac{|{\sigma }_{13}|}{S_{res,13}},\frac{|{\sigma }_{23}|}{S_{res,23}})\le 1$$

(9)

The exposure factors corresponding to these failure criteria are defined as:

$$f_{FF}=\left\{\begin{array}{cc}\frac{{\sigma }_{11}}{X_{res,T}}\hfill & \mathrm{if}{\sigma }_{11}\ge 0\hfill \\ -\frac{{\sigma }_{11}}{X_{res,C}}\hfill & \mathrm{if}{\sigma }_{11}<0\hfill \end{array}\right.$$

(10)

$$f_{MF}=max(-\frac{{\sigma }_{22}}{Y_{res,C}},\frac{{\sigma }_{22}}{Y_{res,T}},-\frac{{\sigma }_{33}}{Z_{res,C}},\frac{{\sigma }_{33}}{Z_{res,T}},\frac{|{\sigma }_{12}|}{S_{res,12}},\frac{|{\sigma }_{13}|}{S_{res,13}},\frac{|{\sigma }_{23}|}{S_{res,23}})$$

(11)

A classical way to model failure numerically is to reduce the material stiffnesses at failure. Here, this is performed by introducing stiffness reduction factors—$R_{FF}$ for the fibre-dominated properties, and $R_{MF}$ for the matrix-dominated properties:

$$\begin{array}{c}{E}_{11}=E_{11}^0·R_{FF}\\ E_{22}=E_{22}^0·R_{MF}\\ E_{33}=E_{33}^0·R_{MF}\\ G_{12}=G_{12}^0·R_{MF}\\ G_{13}=G_{13}^0·R_{MF}\\ G_{23}=G_{23}^0·R_{MF}\\ {\nu }_{12}={\nu }_{12}^0·R_{MF}\\ {\nu }_{13}={\nu }_{13}^0·R_{MF}\\ {\nu }_{23}={\nu }_{23}^0·R_{MF}\end{array}$$

(12)

where $E_{ii}^0$, $G_{ij}^0$ and ${\nu }_{ij}^0$ are the initial properties. $R_{FF}$ and $R_{MF}$ are defined by empirical softening laws based on the values of $f_{FF}$ and $f_{MF}$, as implemented in [54]. $f_{FF}=1$ corresponds to fully failed fibres and $f_{MF}=1$ corresponds to a fully failed matrix. In the laws used here, the failures induced from $f_{FF}$ and $f_{MF}=0.8$ and the respective mechanical properties are reduced linearly until their minimal values when $f_{FF}$ or $f_{MF}=1$. To prevent any numerical instability in the FE analysis, the stiffnesses are reduced to a fraction of their initial value rather than 0. Here, a value of 10% is used. In reality, matrix and fibre failures interact and affect all the material properties (not only the respective ones in the fibre- or matrix-dominated directions). Typically, microscopic fibre failure will generally initiate microscopic matrix cracks. Hence, $f_{FF}$ will also influence the material properties in directions $\{22,33,12,13,23\}$. Conversely, microscopic matrix cracks will reduce the load transfer along and between fibres, which can isolate some fibres, induce their failure, and lower the local longitudinal stiffness. So an increase in $f_{MF}$ will affect the stiffness in the $\left\{11\right\}$ direction. These couplings are accounted for in an empirical manner. The values of $R_{FF}$ and $R_{MF}$, depending on $f_{FF}$ and $f_{MF}$ and the coupling between fibre and matrix damage, are detailed in Equations (13)–(17). These softening laws are also represented in Figure 2.

• No failure: In the case of exposure factors below 0.8, no reduction in the material properties is applied:

  $$f_{FF},f_{MF}<0.8=>R_{FF}=R_{MF}=1$$

  (13)
• Fibre failure only: In the case of fibre failure, both the fibre and matrix properties are degraded equally.
  In the case of partial failure:

  $$0.8\le f_{FF}<1=>R_{FF}=R_{MF}=1-(f_{FF}-0.8)·4.5$$

  (14)
  In the case of full failure:

  $$f_{FF}\ge 1=>R_{FF}=R_{MF}=0.1$$

  (15)
• Matrix failure only: In the case of matrix failure, the matrix-dominated properties are degraded down to 10% of their initial values, but the fibre-dominated properties are only degraded to an empirical value of 60%:
  In the case of partial failure:

  $$0.8\le f_{MF}<1=>\left\{\begin{array}{c}{R}_{MF}=1-(f_{MF}-0.8)·4.5\hfill \\ R_{FF}=1-(f_{MF}-0.8)·2\hfill \end{array}\right.$$

  (16)
  In the case of full failure:

  $$f_{MF}\ge 1=>\left\{\begin{array}{c}{R}_{MF}=0.1\hfill \\ R_{FF}=0.6\hfill \end{array}\right.$$

  (17)
• Mixed failure: In the case where both the fibres and the matrix are damaged, $R_{FF}$ and $R_{MF}$ are calculated independently, as previously, and the minimum material properties resulting from the fibre and the matrix failures are taken.

2.2. Cohesive Failure Model

Inter-laminar delamination and intra-laminar matrix splits are modelled discretely with the CF20 CZM proposed by Davila [48]. This has proven to reproduce crack growth experiments quite well and requires a limited set of material parameters, which are already available for this study.

The crack opening is described as a combination of the three following modes, as presented in Figure 3. The static failure is modelled by a bi-linear traction-separation law. The onset of the damage is defined by a maximum traction criterion using the interface strengths in the three modes. The softening is defined based on the energy release rate, which is the energy per unit length released during the opening of the crack. The mixed-mode configuration is managed by the Benzeggaggh–Kenane [58,59] mode interaction law.

Once the static state is defined, the fatigue damage is calculated from the opening displacement of the crack. It degrades the stiffness until the element is fully damaged. The details of the calculation can be found in [48]. An overview is given here.

The internal damage variable D is defined. Geometrically, it represents the proportion of the softening domain which has been reached:

$$D=\frac{\lambda -{\Delta }_c}{{\Delta }_f-{\Delta }_c}$$

(18)

D is equal to 0 at damage onset and to 1 at full failure. The stiffness is reduced according to a derived damage variable d:

$$K=K_i\ast \left(1-d\right)$$

(19)

where d and D are related geometrically by the bi-linear traction-separation law, as illustrated in Figure 4. This relation is written:

$$1-d=\frac{(1-D){\Delta }_c}{D{\Delta }_f-(1-D){\Delta }_c}$$

(20)

During a loading step, the peak opening displacement $\lambda$ is recorded. The ratio between this peak displacement and the maximum displacement allowed by the static failure envelope is computed. It is called the relative displacement jump $\frac{\lambda }{\lambda *}$. The damage increment is then calculated according to Equation (21). The full procedure is illustrated in Figure 5.

$$\Delta D=\frac{1}{\gamma }\frac{{(1-D)}^{\beta -p}}{E^{\beta }(p+1)}{\left(\frac{\lambda }{{\lambda }^{*}}\right)}^{\beta }\Delta N$$

(21)

where E and $\beta$ are model variables related to the material fatigue life under the current mixed-mode opening. $\gamma$ and p are model parameters to calibrate the fatigue crack growth rate from experiments.

2.3. Cycle-Jump Numerical Procedure

Simulating the fatigue damage of a laminate raises the challenge of computation time. Progressive fatigue modelling needs to model the evolution of gradual and interacting damage mechanisms all along the lifetime. Simulating each and every cycle would be too computationally costly. To overcome this, Van Paepegem and Degrieck proposed a so-called cycle-jump approach [1]. In this method, the stress state is calculated at given cycles and the damage is then extrapolated over a jump of cycles. Numerous examples of successful implementation of this approach can be found in the literature [34,38,40,42,43,53].

This model implements a cycle-jump approach with a rather low number of cycles, typically ranging from 10–100 steps. This aims at keeping the computational time as low as possible, with a view to applying this model at the component scale. As a reminder, the purpose of this attempt is to assess the descriptive capabilities of the model in terms of catastrophic failure with minimal computation time.

A cycle-jump scheme can be implemented in two different ways: internally in the FEA software, or externally, as presented by Leidermark and Simonsson [60]. The internal implementation is integrated within the FE software and is possibly quicker but needs to handle the material degradation within the simulation. This means that element softening must be handled in the analysis. The inherent negative tangential stiffness resulting from material degradation leads to convergence difficulty in conventional implicit numerical solvers. To overcome this issue, an external implementation of the cycle-jump scheme is chosen, where the stresses are read from the FE simulation and processed outside of the solver to compute the fatigue damage before starting a new simulation with updated material properties. This is further described in the following section.

First, a loading step is conducted from the initial state up to the maximum fatigue load value. During this quasi-static loading step, the bulk material behaves according to its linear-elastic constitutive law without damage. As mentioned earlier, this considerably eases the convergence in the implicit scheme. The cohesive model implemented here for macroscopic matrix splits and delamination integrates a quasi-static damage in the traction-separation law, which is active during the loading step. Thus, cohesive damage can onset and propagate in a static mode during this step.

Once the peak of the fatigue cycle is reached, the loading step stops. The stress state at the peak is recorded in the whole model. The maximum stresses are recorded from all the integration points in the bulk elements, and the maximum opening displacements are recorded from the integration points in the cohesive elements.

From these maximum stresses and displacements, the fatigue damage for both the CDM and the CZM are computed for the current cycle-jump $\Delta N$. The material properties are then degraded accordingly and a new loading step is started. An illustration of the principle is proposed in Figure 6.

In this implementation, the size of the cycle-jump is user-defined prior to the simulation using a regular interval on a logarithmic scale.

3. Experiments and FEA Model Assessment

Tensile cyclic fatigue tests of glass-fibre epoxy laminates with a hole were conducted. The tests were reproduced numerically with the FE model. The flat laminate allows for a predictable plane stress-field in the laminate, while the hole introduces a singularity and a stress concentration. This causes the onset of the fatigue damage in a predetermined area. A layup of [90/0]${}_s$ was chosen, as it exhibits the growth of an intra-laminar matrix split in the longitudinal plies and of inter-laminar delamination between the $0^{°}$ and ${90}^{°}$ plies. The specimen’s residual stiffness, matrix damage, and fatigue life were observed and compared with the numerical results.

3.1. Experimental Setup

Tensile cyclic fatigue tests were performed on R-glass epoxy open-hole laminates produced by vacuum-assisted resin transfer (VART), according to the test standard ASTM D7615 [61]. The layup was [90/0]${}_s$. The reinforcement used is a Saertex unidirectional (UD) stitch-bonded layer (a small amount of fibres, less than 0.1% of the total weight, is used in the transverse direction to keep the fibre tows together) made from 3B’s HiPer-tex W2020 R-glass fibre with a weight of 1150 g/m${}^2$. EPIKOTE™ MGS™ RIMR 135 epoxy resin was used with the hardener EPIKURE™ MGS™ RIMH 137 in a mix ratio of 100:30. The laminate was cured for 24 h at room temperature and post-cured at 80 ${}^{°}$C for 15 h. Specimens were cut with a water-jet to the dimensions L = 300 mm, w = 40 mm. The geometry of the open-hole test specimens was chosen to challenge the FE analysis creating a significant stress concentration at the hole and an uneven stress field between the hole and the edge. Thus, the geometry was different than the one recommended by the test standard ASTM D5766 [62]. The ply thickness was measured to 0.95 mm, giving a total thickness of 3.8 mm for the 4-plies laminate. The open hole at the centre had a diameter of 10 mm. The geometry of the specimens is presented in Figure 7. Defects induced during the machining of the hole can affect the load-bearing capacity of the laminate [63]. In the case studied here, small delaminations could be observed after machining. The impact of such delamination under tensile loading condition was minimal and the size of the defects was negligible compared to the extent of damage developing in the early stage of the fatigue lifetime.

The maximum static tensile load was estimated for one specimen. A quasi-static tensile test was conducted in displacement control mode on an MTS Series 809,100 kN servo-hydraulic test machine, with a speed of 2 mm·min${}^{-1}$. Failure occurred after a displacement of 7.53 mm at a maximum load of 56.84 kN.

Cyclic tensile fatigue tests were conducted in load control mode at a load ratio of R = 0.1. At the start of the fatigue tests, a quasi-static ramp was done up to the mean load. The load rate of this ramp was set to 1 kN·s${}^{-1}$. Then the load cycles were conducted, starting upwards from the mean to the maximum load. The load levels were defined to obtain failures in each decade from ${10}^4$ to ${10}^6$. To minimise the influence of visco-elastic effects, the frequency was scaled proportionally to the load so that the average strain-rate stayed as constant as possible between the different specimens. This point is slightly different than the test method recommended from the test standard ASTM D7615 [61]. Three specimens were tested at the same load level (30 kN) to observe the scatter in the fatigue life.

Table 1. Summary of the tensile cyclic fatigue tests. The frequency was scaled proportionally to the inverse of the load to keep the strain rate as constant as possible between the specimens.

Specimen	Load [kN]	Frequency [Hz]	Fatigue Life [Cycles]
1	22	1.85	Runout 2M
2	26	1.54	1,060,822
3	28	1.43	441,030
4	30	1.33	374,169
5	30	1.33	205,428
6	30	1.33	332,648
7	32	1.25	114,051
8	32	1.25	169,506
9	34	1.18	28,627

summarizes the specimen test loads and frequencies.

The residual stiffness and the number of cycles at failure were recorded by the load frame. The residual stiffness of the specimen was defined as the change in the macroscopic stiffness of the specimen. It is, thus, equal to the change in the crosshead displacement at the peak of the cycle. While it is usually not advised to use the load frame measurement to quantify the specimen’s elongation during the test, the precision of the measurement was sufficient for our study here. Due to the significant stiffness difference between the fixture chain (grips, load-cell, and crosshead) and the specimen, the deviation in the stiffness measured from the crosshead displacement and the true specimen stiffness was estimated to be well below 0.01%.

For two specimens, a camera was installed to capture pictures of the specimen along the lifetime. The matrix damage development was observed qualitatively and the split length was measured along the lifetime. This was performed on two specimens tested at the same load level of 32 kN. The camera was synced to the load frame, and the picture acquisition was performed precisely at the peak of the cycle at a predefined interval. A speckle pattern was applied to the specimen to post-process the images with a digital image correlation (DIC) software. The goal was to obtain a second measure of the stiffness evolution of the specimen. However, due to extensive matrix cracking of the top transverse ${90}^{°}$ ply, the DIC measurement was not reliable enough to provide a correct measure of the residual stiffness. The stiffness measured from the maximum crosshead displacement was, thus, conserved.

3.2. FEA Implementation

The flat specimen with a hole was modelled in the commercial FEA software Abaqus 2019. The geometry was reduced to 1/8th of the specimens thanks to the three symmetry planes of the test setup, as described in Figure 8. Three planar-symmetry boundary conditions (X-, Y-, and Z-planes) were applied on the corresponding faces of the specimen. The specimen was modelled up to the beginning of the machine grips. The load was applied to the tip of the specimen using a reference point, which was connected to the tip surface by a beam multi-point constraint (MPC). On this point, only the translation along the x-axis was allowed, and the axial force was applied along this axis. The PFM requires each ply to be modelled discretely with solid elements.Thus, each ply is modelled as an individual part, using a solid composite layup in Abaqus, and several elements through the thickness. In this way, the through-thickness stresses are accounted for and the CDM is applied on each ply individually. The inter-laminar CZM was inserted between the plies. To model the intra-laminar matrix split in the $0^{°}$ ply, this ply was modelled as two separate parts linked by a vertical band of cohesive elements. The individual ply and cohesive parts in the model are connected using surface-based tie constraints. In order to ease the meshing of the geometry, the split was modelled slightly aside of the hole. This creates a so-called bridge area in the $0^{°}$ ply, as shown in Figure 8. By doing so, a regular mesh without wedged elements can be created. With this geometry, the right part of the $0^{°}$ ply is not totally free-hanging as the elements in the bridge area are constrained in the (y,z) plane by the symmetry boundary condition (BC). To assess the influence of this bridge, simulations were run where these elements were free from this BC. No significant changes in the predicted residual stiffness, fatigue life, or split growth were observed.

3.2.1. Mesh Description

To allow more flexibility in the mesh pattern and refinement, the meshes between the different parts are not compatible. The parts are tied together using Abaqus surface-based tie constraints. The top ${90}^{°}$ ply and the core $0^{°}$ ply are modelled using C3D8R solid elements with reduced integration. This element is a standard continuum solid hexahedral element with eight nodes and a reduced integration scheme, which is more computationally efficient. The element deformations are controlled using the built-in hourglass-enhanced method and distortion controls.

The delamination interface and the split band use fully integrated COH3D8 hexahedral cohesive elements. For convenience, the cohesive layers are modelled with a geometrical thickness of 0.01 mm. Only one element is in the thickness of the cohesive interface. The constitutive thickness used for the computation of the separation of the cohesive elements is equal to 1.0. In this way, the geometrical thickness does not influence the mechanical response of the elements. The geometrical thickness of the cohesive layers is subtracted from the plies to make sure that the total laminate thickness and width are correct.

A mesh sensitivity study was conducted to evaluate the best compromise between computational cost and accuracy of the results. It was found that the results predicted, in terms of the lifetime, were quite insensitive to the mesh size. However, the mesh size is particularly critical for the CZM to predict damage onset and growth. Davila used elements smaller than 0.1 mm to model crack growth over a length of 20 mm in an MMB test [48] and used elements of size 0.03 to avoid noisy results. Leone et al. explored the mesh size requirements to model crack growth using both CDM and CZM approaches. They found that three elements in the crack process zone were necessary for the results to converge [64]. A similar number of elements in the crack process zone was suggested by Turon et al. [65]. The crack process zone length in the case of mode I and mode II crack opening in orthotropic materials was derived in [25]. Adapted to our configuration, the minimum crack process zone was found to be 4.7 mm for mode I. Thus, cohesive elements of size 1 mm or less are sufficient. In order to match the fine mesh around the hole, the element size in the split band was set to 0.1 mm in the finest area. The element length in the delamination interface is 0.3 mm on average, 0.1 mm in the finest area, and 0.5 mm further from the hole. The bulk material is also discretized with an element size of approximately 0.3 mm on average, and down to 0.1 mm in the finest area. A total of 180k elements were created. The model was run on an Intel(R) Xeon(R) W-2155 CPU @ 3.30 GHz, 10 Core(s) in parallel on the 10 CPUs. A computation time of 200–500 s wall-clock time per loading step was obtained, and 250–400 s were used by the Python programme to compute the fatigue damage variables at each integration point. A total of 1.1 GB of memory were required for the simulation. The total simulation time for 25 loading step was 10 h.

3.2.2. Numerical Implementation

The numerical procedure consists of FE loading steps conducted in Abaqus/Standard and fatigue damage steps conducted by a Python programme outside of Abaqus, between each loading step to calculate and extrapolate the fatigue damage over a jump of cycles. During the loading step, the response of the bulk elements is defined by a linear elastic law in the user-defined material subroutine UMAT. The cohesive elements behaviour is also implemented in the UMAT, where the static traction-separation law, including the static damage onset and softening, is defined. For both element types, the material properties used for this loading step are initialised with properties affected by the fatigue damage from the previous steps. This is performed prior to the UMAT by updating the State-Dependent Variables (SDV) in each element using the SDVINI subroutine. At the end of the loading step, the peak stresses are read by the Python programme. The fatigue damage in every element (bulk and cohesive) are then calculated over a given jump of cycles, as per the CDM model presented here for the bulk elements, and as per the CF20 CZM for the cohesive elements. The damage variables are used to update the material properties at the start of the next loading step calculation. The global procedure is illustrated in Figure 9.

3.2.3. Material Properties

The material properties and parameters are presented in

Table 2. Material properties of the GFRP material used for the CDM and CZM models.

Property	Units	Values	References
Ply elastic properties 1,2,3
$E_{11}$	GPa	44.870	[28]
$E_{22}=E_{33}$	GPa	12.130	[28]
$G_{12}=G_{13}$	GPa	3.380	[28]
$G_{23}$	GPa	4.043	Estimated
${\nu }_{12}={\nu }_{13}$	-	0.3	[28]
${\nu }_{23}$	-	0.5	Assumed
Intra-laminar strength 3,4
$X_T$	MPa	1006.30	[28]
$X_C$	MPa	487.00	[28]
$Y_T=Z_T$	MPa	45.95	[28]
$Y_C=Z_C$	MPa	131.9	[28]
$S_{12}=S_{13}=S_{23}$	MPa	49.51	[28]
Intra-laminar fatigue properties 1,2,3,4
${\sigma }_{11,T}^0$	MPa	1006.30	Assumed
${\sigma }_{11,C}^0$	MPa	487.00	Assumed
$a_{11}$	-	0.1	Assumed
${\sigma }_{22,T}^0={\sigma }_{33,T}^0$	MPa	45.95	Assumed
${\sigma }_{22,C}^0={\sigma }_{33,C}^0$	MPa	131.9	Assumed
$a_{22}=a_{33}$	-	0.054	[67]
${\tau }_{12}^0={\tau }_{13}^0={\tau }_{23}^0$	MPa	49.51	Assumed
$a_{12}=a_{13}=a_{23}$	-	0.054	[67]
Inter-laminar properties
${\tau }_I$	MPa	45.95	[28]
${\tau }_{II}={\tau }_{III}$	MPa	49.51	[28]
$G_{Ic}$	N/mm	0.98	[28]
$G_{IIc}=G_{IIIc}$	N/mm	3.71	[28]
$\eta$	-	1.40	Assumed
$K_I$	GPa	12.130	[28]
$K_{II}=K_{III}$	GPa	3.720	[19,48]
CF20 CZM parameters
$ϵ$	-	0.44	Estimated
$\gamma$	-	${10}^7$	[48]
${\eta }_{br}$	-	0.95	[48]
p	-	0.0	[48]
Model-fitting parameters
$p_F$	-	2.0	Calibrated
$p_M$	-	2.2	Calibrated

¹ ${}_{11}$ is in the longitudinal direction. ² ${}_{22}$ and ${}_{33}$ are in the transverse directions. ³ ${}_{12}$, ${}_{13}$ and ${}_{23}$ are in the shear directions. ⁴ $T=tension$, $C=compression$.

. The fibre volume fraction was measured at 54% for this material and process in previous work [28]. Microscopy was conducted to verify that this value was consistent for this batch. Most of the material properties of these laminates were measured in previous work by Perillo and Gagani [28,66,67], while the others were either taken from the literature or assumed. Given the predominance of tensile loads in this study case, and the lack of material data, the fatigue material parameters in compression were taken to be equal to the tensile ones. This had no influence on the model results here, but relevant material properties in compression must be used when simulating load cases with significant compression.

To describe the experimental results, two calibration parameters are used in the model. The first one, $p_F$, is used to determine the S-N curve intercept for the UD longitudinal $0^{°}$ ply. It is applied to offset the intercept from the static strength $X_T$. It is observed experimentally for a $0^{°}$ ply that the S-N curve intercept is not equal to the static strength. For GFRP, it is usually larger [68,69,70]. Due to the lack of values for our material and the difficulty in obtaining reliable values experimentally, the calibration of this material parameter was conducted. The model was calibrated on experiments performed at one load level (30 kN). In the model, the value of $p_F$ governs the catastrophic failure of the specimen. Increasing $p_F$ retards the final failure, reducing $p_F$ advances it. The calibrated value of 2.0 was, thus, determined by iteration until the final failure in the FEA fell within the experimental results at 30 kN.

The second calibration parameter is applied on the matrix strengths and the S-N curve intercepts. This calibration is needed because of the inherent scatter of the matrix strength and the volume effect induced by this scatter. When testing the transverse strength of a plain laminate, the specimen fails at its weakest point. When looking at a bigger laminate, there is a chance that the weakest point is further down the tail of the strength distribution, resulting in a smaller measured strength. This was observed experimentally in [71]. When testing a laminate with a hole, the volume where the damage is induced and develops is smaller than the typical laminate from which the strengths are measured. The material is, thus, locally stronger. To take this effect into account, the authors decided in this study to scale the matrix strength with the calibration parameter $p_M$. It is applied to the matrix strengths and the matrix S-N curves intercepts. A low value of the matrix strengths lead to a full failure of the ${90}^{°}$ ply in the first cycle, while a higher value leads to a gradual failure of the ${90}^{°}$ ply. The value of $p_M$ does not significantly affect the predicted fatigue life. The value of $p_M$ was determined by comparing the residual stiffness evolution in the FEA to the one observed experimentally. The same load level as for the calibration of $p_F$ was used (30 kN) so that the calibration was performed for one single load level. A value of 2.2 was chosen so that the matrix failure in the FEA would occur gradually over the first hundreds of cycles at this load level.

Table 3. A set of matrix strengths and the $0^{°}$ ply S-N curve intercept were calibrated based on the experimental results at one load level.

Property	Unit	Description	Initial Value	Calibrated Value
${\sigma }_{11}^0$	MPa	SN curve intercept 1	1006.30	$p_F·1006.30$
${\sigma }_{22,T}^0={\sigma }_{33,T}^0$	MPa	SN curve intercepts 2	45.95	$p_M·45.95$
${\sigma }_{22,C}^0={\sigma }_{33,C}^0$	MPa	SN curve intercepts 2	131.9	$p_M·131.9$
${\tau }_{12}^0={\tau }_{13}^0={\tau }_{23}^0$	MPa	SN curve intercepts 3	49.51	$p_M·49.51$
$Y_T=Z_T$	MPa	Transverse strength 4	45.95	$p_M·45.95$
$Y_C=Z_C$	MPa	Transverse strength 4	131.9	$p_M·131.9$
$S_{12}=S_{13}=S_{23}$	MPa	Shear strengths 3	49.51	$p_M·49.51$

¹ ${}_{11}$ is in the longitudinal direction. ² ${}_{22}$ and ${}_{33}$ are in the transverse directions. ³ ${}_{12}$, ${}_{13}$ and ${}_{23}$ are in the shear directions. ⁴ $T=tension$, $C=compression$.

summarizes the material parameters, which were scaled using the calibration parameters $p_F$ and $p_M$.

4. Results and Discussion

4.1. Model Behaviour

The results of the calibrated model are plotted along with the experimental results in Figure 10. The failure sequence exhibited by the FEA can, thus, be analysed. Three drops in the stiffness are observed: the first one immediately after the first cycle, the second starting at 100 cycles, and the final catastrophic failure occurring in the range of 200k–400k cycles. To analyse the failure, six stages in the lifetime are defined around the three residual stiffness drops, as displayed in Figure 10. For each of them, the longitudinal stresses and fibre failure are investigated along two paths in the load-bearing $0^{°}$ ply:

• Along the segment [AB], in the width of the specimen,
• Along the segment [AC], along the split band.

The evolution of the axial stress ${\sigma }_{11}$ along these paths in the left part of the $0^{°}$ ply is plotted in Figure 11 at the six stages of interest.

After the first cycle, a small stiffness drop is observed. It is due to the failure of the fibres in the bridge area next to the hole and to the onset of matrix failure in the transverse plies.

A second drop in the residual stiffness occurs between 100 and 300 cycles. It corresponds to the full development and saturation of transverse matrix cracking in the ${90}^{°}$ ply. This was observed by looking at the damage variable $f_{MF}$ between stages 3 and 4. During this stage, the split in the $0^{°}$ ply grows and smears the stress concentration at the hole. The damage variable in the fibres next to the hole increases, but the failure is not initiated. The exposure factor does not yet reach the onset value of 0.8. Thus, this stiffness drop is, solely due to the failure of the ${90}^{°}$ ply. Compared to the gradual stiffness reduction in the experiment, this stage appears quite suddenly in the FEA. This can be explained by the CDM model, which homogenises the matrix properties and does not account for the inherent scatter in the matrix strength. The real failure initiates at the weak spots of the matrix, and the matrix cracks develop much more gradually than the homogenized brittle failure implemented here. This could be better described using discrete CZM bands, as in the work of Yashiro et al. [27], or using a probabilistic description of the strengths, as reported in [24]. However, these approaches are complex and computationally costly. It is assumed that such fidelity in the gradual stiffness reduction would bring little improvement in the prediction of the fibre-governed catastrophic failure, which is the main focus of this study.

Between stages 4 and 5, the residual stiffness decreases with a small slope in a very gradual way. In the FEA, this is due to the onset of fibre damage at the stress concentration next to the hole. The exposure factor reaches the threshold value of 0.8 in some elements and their stiffness is reduced. This happens at the equator-line from the hole and not at the tip of the split. The split creates a stress concentration as well, but smaller than at the hole, and not sufficiently to induce fibre failure there.

After stage 5, catastrophic failure occurs in the FEA by fibre failure propagation from the hole towards the edge of the specimen. When the intact remaining cross-section is too small to withstand the load, the specimen fails. The onset and propagation of the fibre failure is reasonably well-described, but the model does not show a sharp stiffness drop as expected in the structure. This has two causes:

• The number and location of the conducted loading steps and cycle jumps are predefined. The cycle-jump scheme is not able to reproduce a very sharp drop in stiffness because the jumps between the loading steps are too big for the fast degradation and the damage cannot propagate faster than the damage process zone during each cycle jump. This is a compromise between computation time and fidelity in the final failure description.
• For the sake of numerical stability, the stiffness in the failed elements is reduced to 10% of the initial values. The residual stiffness in these elements is obviously not insignificant and a part of the load will still be carried by this region. This artificial effect is particularly pronounced when a large part of the cross-section has developed damage.

The first limitation could be overcome in different ways. First, the number of steps computed close to the final failure could be refined. Such an algorithm must, however, know the time of final failure in advance, which is possible in a descriptive approach, but not in a predictive approach. Computing more cycles all along the lifetime could be a workaround, but would be at the expense of a significant computational cost. Otherwise, an automatic adaptive calculation of the cycle jumps could be implemented, as presented first by Van Paepegem and Degrieck [1] and then in [32,34,37,40,42,43,53]. This is the best way to ensure good agreement between the simulation and the experiments. The size of the maximum admissible damage increment can be adapted and a compromise can be found between the fidelity and the computation time.

The second limitation could be overcome by a modification of the stiffness reduction laws at failure (Figure 2). Which stiffness degradation law should be used to describe the fibre failure is still an open question. The counterpart of this gain in fidelity is an increased difficulty in achieving convergence using an implicit solver. This is why these laws and the value of 10% remaining stiffness were used in this study.

4.2. Numerical Fatigue Life

The experimental fatigue test results are summarized in Table 1. The model was run for the corresponding maximum load values. The normalised residual stiffness was plotted for each simulation (Figure 12) and the failure was determined graphically at the knee point of the final stiffness drop. The knee point was found to occur around a similar value of the normalised residual stiffness of 0.76 for all the specimens. This value was, thus, used as a final failure criterion for the FEA simulations.

The experimental and numerical data points are plotted in Figure 13 and both fitted with Basquin’s S-N curve equations [3] through linear regressions in the log-log space, with N as the dependent variable. Two experimental specimens were excluded from the fit: the run-out specimen at the load level 22 kN, and the specimen at the load level 34 kN. This point was considered to be already in the area where the S-N curve bends towards a flatter slope, which is typical for the short fatigue life of GFRP. Including this point in the linear regression gives a too optimistic (i.e., too flat) S-N curve slope in the high-cycle domain.

The agreement is found to be good, as the numerical S-N curve is well within the experimental scatter band. The value of the experimental slope of −0.11 close to the S-N curve slope of the $0^{°}$ ply used as input of the model was −0.10. The FEA produces a slope of exactly −0.10. Unlike in the experiments, the simulations do not seem to predict any inflexion of the S-N curve at high loads and short lives. This can be explained by the model principle and its inputs. In the model, the mechanism governing the catastrophic failure is the fatigue damage of the $0^{°}$ ply. This damage is ruled by the $0^{°}$ ply S-N curve used as input, which does not include any inflexion itself. The model predicts failure based on the shape of this input S-N curve.

4.3. Split Growth

4.3.1. Numerical Prediction of the Split Length

The length of the split was measured from pictures taken along the lifetime of two specimens at the same load level. The four splits $\{S_1,S_2,S_3,S_4\}$ around the hole were measured, as displayed in Figure 14. The results are shown in Figure 15. Initially, the split in the FEA grows slightly faster than is observed experimentally. It then reaches a plateau, while the splits in the experiments increase rapidly during the end of the lifetime and the specimen’s catastrophic failure.

The initial over-prediction of the model has several origins. First, the intra-laminar matrix crack was modelled using cohesive properties measured on inter-laminar cracks. While quite close, the values of the strength and the energy release rate are usually higher for intra-laminar cracks. This may be due to the significant fibre-bridging which can occur in such cracks. Secondly, it has been shown by Davila in recent work [49] that the CF20 model tends to under-estimate the endurance limit describing a crack growth in mode II. While the onset and first millimeters of splitting are driven by mode I opening, the progression rapidly shifts to a mode-II-dominated failure. Moreover, the mode-mixity is not constant within the crack process zone (CPZ), and the elements further in the CPZ are loaded almost in pure mode II. These considerations should be investigated further to explore whether the model can describe the early split growth with higher fidelity.

Towards the middle and end of the lifetime, the split growth becomes under-predicted by the model. This is most likely due to the softening and failure law implemented in this CDM. Due to convergence limitations, the stiffness of the failed elements is defined as 10% of their initial stiffness. This leads to a significant attenuation of the shear stresses in the split band, even in the failed elements, preventing the split from growing further. The plateau in the split was first believed to be an artifact of the so-called bridge area, which continues to carry some stresses due to the 10% remaining stiffness. To verify this hypothesis, a simulation was conducted without boundary conditions on the bridge area, allowing this surface to be free to move longitudinally along the x-axis. By doing so, the right part of the $0^{°}$ ply is free-hanging and transfers all its load to the left part of the $0^{°}$ ply through shear stresses in the split band. The split growth obtained with this geometry showed no improvement compared to the FEA with the constrained bridge area. Finally, the important split growth observed experimentally happened close to the catastrophic failure of the laminates. At this stage of the lifetime, the simplified model is not able to reproduce the failure mechanisms very precisely because of the fixed cycle-jump increment size and the remaining 10% stiffness in the failed elements. The influence of an adaptive cycle-jump size and more accurate failure behaviour should be explored further.

4.3.2. Influence of the Split on the Predicted Fatigue Life

To assess the influence of the split on the fatigue life and to be able to conclude whether it is necessary to include this failure mechanism in the FEA or not, a model was run without allowing any splitting. Technically, the split band was still present in the FE mesh; however, the interface strengths were set to unrealistically high values to ensure that no damage would onset in the cohesive elements of the split band. The observed fatigue life was slightly shorter than for the reference simulation with splitting. When looking at the stress field and the stresses along the path [AC], higher stresses are observed when splitting is prevented. In this case, the split does not smear the stress concentration at the hole, resulting in a shorter predicted fatigue life. This justifies the need for the inclusion of the split in the analysis.

4.4. Delamination

Pictures of the two specimens were also used to compare qualitatively the development of the inter-laminar delamination between the model and the experiment. Figure 16 presents the pictures taken at different stages of the lifetime and the cohesive damage numerically predicted between the ${90}^{°}$ and $0^{°}$ ply in the FEA. During the experiments, the delamination starts initially at the stress concentration next to the hole. It is believed that the split start simultaneously. The delamination is then observed to develop along the split band. It is believed that the split is developing underneath and that the two failure mechanisms are interacting. The split is essentially a through thickness tunnelling crack in the $0^{°}$ ply, which can drive the delamination at the interface on top of it. But delamination onsets are also observed all over the laminate in bands transverse to the load direction. These are created by the through thickness transverse matrix cracks in the top ${90}^{°}$ ply. This observed failure sequence is well-known and well-described [72]—transverse matrix cracks will induce delamination when reaching a ply interface. Under fatigue loading, these delaminations grow and then coalesce, causing the plies to be almost fully delaminated. This can be clearly seen in the pictures of the test specimens.

In the FEA, the interaction between the delamination and the split is well-captured. This was already demonstrated by Yashiro et al. [27] and is reproduced here. However, the smeared approach of the transverse matrix cracks in the ${90}^{°}$ ply implemented here does not capture this interaction between the cracks and the delamination. Hence, the delamination does not onset anywhere else than above the split and does not grow and coalesce elsewhere. To describe the delamination growth with higher fidelity, discrete or semi-discrete models, such as [26,27], could be used. Analytical numerical models were also recently proposed by Maragoni et al. [73].

5. Conclusions

An FE progressive fatigue model (PFM) is proposed to describe the catastrophic failure of open-hole GFRP laminates. It combines a CDM for intra-laminar failure and a CZM for inter-laminar failure and intra-laminar splitting. The CDM and CZM fatigue models rely on one fibre-dominated and one matrix-dominated S-N curve. They are combined with the Palmgren–Miner damage accumulation rule and Hashin’s residual strength. An empirical softening law was chosen to model the stiffness reduction at failure. In the present implementation, an external cycle-jump scheme with predefined intervals was used to mitigate the convergence issues and to optimise the computational time.

• A set of matrix strength and fibre-dominated S-N curve parameters were calibrated on a first cyclic fatigue experiment at one given load level. After this step, the model was able to describe the fatigue life at other load levels with good agreement: the S-N curve slopes obtained were −0.11 experimentally and −0.10 numerically. It was able to reproduce the S-N curve slope successfully.
• The smeared CDM and the simplistic softening law used for matrix damage make it challenging to reproduce the gradual stiffness degradation with high fidelity. The stiffness drop at an early stage in the model was more sudden than gradual. However, the amount of stiffness loss at catastrophic failure agreed well with the value observed experimentally.
• The description of the split growth during the first stage of the lifetime was satisfying, but the FEA was not able to reproduce the fast split growth at the last stage of the lifetime. This might be due to the simplified softening law and the remaining stiffness in the failed elements. This remaining stiffness was required to achieve convergence but maintained an artificial load bearing capacity in the plies, reducing the shear solicitation on the split band.
• The inter-laminar delamination was described with high fidelity above the split. The interaction between the split and the delamination is well-reproduced by the two cohesive layers inserted in the FEA. However, the smeared CDM modelling of the transverse matrix cracks in the ${90}^{°}$ ply fails to capture the delamination induced by those cracks. This leads to a significant underestimation of the delamination extent and growth along the lifetime.

Despite these limitations, the FE model was successful in reproducing the catastrophic failure of the specimens, using simple inputs and theories available. This was the main focus of this attempt.