# High Cycle Fatigue Damage Model for Delamination Crack Growth in CF/Epoxy Composite Laminates

Laurent Gornet, Hassan Ijaz

► To cite this version:

Laurent Gornet, Hassan Ijaz. High Cycle Fatigue Damage Model for Delamination Crack Growth in CF/Epoxy Composite Laminates. International Journal of Damage Mechanics, 2011, 20, pp.783 - 807. 10.1177/1056789510374166. hal-01443169

HAL Id: hal-01443169

https://hal.science/hal-01443169

Submitted on 22 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# High Cycle Fatigue Damage Model for Delamination Crack Growth in CF/Epoxy Composite Laminates

LAURENT GORNET AND HASSAN IJAZ

GeM-UMR-CNRS 6183, Ecole Centrale de Nantes – 1 Rue de la Noë, BP 92101, 44321 Nantes Cedex 3, France

**ABSTRACT:** This article presents the development of a fatigue damage model which helps to carry out simulation of the evolution of delamination in the laminated composite structures under cyclic loadings. A classical interface damage evolution law, which is commonly used to predict the static debonding process, is modified further to incorporate fatigue delamination effects due to high cycle loadings. An improved formulation is also presented to incorporate the ‘R’ ratio effects. The proposed fatigue damage model is identified using fracture mechanics tests like double cantilever beam, end-notched-flexure and mixed-mode bending. Then a non-monotonic behavior is used to predict the fatigue damage parameters able to carry out delamination simulations for different mode mixtures. Linear Paris plot behaviors of the above-mentioned fracture mechanics tests are successfully compared with available experimental data on HTA/6376C and AS4/PEEK unidirectional materials.

**KEY WORDS:** fiber-reinforced materials, damage mechanics, delamination, high cycle fatigue, finite element analysis.

## INTRODUCTION

FOR WEIGHT SAVING purposes in aircrafts, trains, ships, and sailing race boats applications, the use of composite materials is no longer limited to secondary structure, but is expending to primary load bearing parts. Due to their laminated nature, composite laminates are prone to delamination failure under static (Davies, 1989; de Morais and Pereira, 2007) and fatigue loadings (Martin and Murri, 1990; Kenane and Benzeggagh, 1997; Hojo and
Gustafson, 1987; Asp et al., 2001). Delamination phenomenon consists of the separation of the adjacent layers of laminated composites. Under fatigue loading conditions, this problem becomes more severe and can cause dramatic reduction of the load carrying capability of the laminated structures.

During the past 20 years, for monotonically applied loading, a lot of work has been carried out at the meso-scale level by authors to model damage mechanism of composite laminates (Allix and Ladevèze, 1992; Chaboche et al., 1997; Allix et al., 1998). Meso scale is strongly connected to the laminated scale, which lies between micro scale (fiber scale) and macro scale (structural scale). A strategy to model laminated composite uses two basic damageable constituents: the layer and the interface. The interlaminar interface, which is a mechanical surface, connects two adjacent layers and depends on the relative orientation of their fibers (Allix et al., 1995). An overview of delamination modeling for static loading can be found in (Tay Te, 2003). Few articles focus on delamination modeling in composite materials under fatigue loadings (Robinson et al., 2005; Tumino and Cappello, 2007; Turon et al., 2007).

In this article, a comprehensive interlaminar interface fatigue damage evolution law is proposed to model the delamination phenomenon under fatigue loadings. The modeling was implemented in the finite element code Cast3M (CEA) through the user subroutine material (UMAT). The double cantilever beam (DCB), end-notched-flexure (ENF) and mixed-mode flexure (MMF) tests were chosen to identify the proposed model for simulations of the crack growth in unidirectional carbon-fiber epoxy-matrix materials.

The proposed model takes its foundations from the classical static damage evolution law proposed by Allix et al. (1992, 1998). The main idea for the fatigue crack growth modeling was first introduced by Peerlings et al. (2000) for metallic parts. The idea of fatigue damage evolution law introduced by Peerlings is also very close to the one proposed by Paas et al. (1993). Robinson et al. (2005) presented the fatigue-driven delamination for the laminated composites using the idea of Peerlings et al. (2000) for the cyclic load, varying between maximum and zero values. In Robinson et al. (2005), fatigue damage evolution is a function of relative displacement of adjacent layers. In this article, the framework of thermodynamics of an irreversible process is used to derive a damageable interface modeling (Allix and Ladevèze, 1992).

The fatigue damage evolution law proposed in this article is a function of thermodynamic forces $Y_{di}$ (also called damage energy release rate) and of critical damage energy release rate $Y_C$ (Allix et al., 1992, 1998). The proposed fatigue damage evolution law permits us to reproduce the linear crack growth rates as obtained by using classical Paris law for fracture mechanics tests (Martin and Murri, 1990; Juntti et al., 1999, Asp et al., 2001). The Paris
law, which depends on two parameters namely B and m, can be expressed as follows:

$$ \frac{da}{dN} = B \left( \frac{\Delta G}{G_C} \right)^m, \qquad (1) $$

where a is the crack length and N the number of cycles. $\Delta G = (G_{max} - G_{min})$ is the cyclic variation of energy release rate with $G_{max}$ and $G_{min}$ representing the maximum and minimum values of energy release rates during the oscillation, respectively. $G_C$ is the fracture toughness of the material, B and m are constants and they are determined experimentally.

In some continuum damage mechanics theories, the failure of interface is taken into account by three damage variables $d_1$, $d_2$, and $d_3$. The delamination crack growth under high cycle fatigue can be considered as a combination of delamination due to the quasi-static loading and due to the cyclic variation of the loading; hence, the total damage evolution for three different modes of failure can be expressed as follows:

$$ \dot{d}_{iT} = \dot{d}_i = \dot{d}_{iS} + \dot{d}_{iF} \quad i = 1, 2, 3 \qquad (2) $$

where the term $\dot{d}_{iS}$ corresponds to delamination growth under static loading and $\dot{d}_{iF}$ is related to the one under fatigue. The details of static and fatigue damage variables will be discussed later.

The article is organized as follows: in 'Review of Static Interface Damage Model' section, classical damage model proposed by Allix et al. (1992, 1998) for the prediction of delamination in laminated composites is recalled. In 'Fatigue Interface Damage Model' section, the proposed fatigue damage model along with simulations and results is presented in detail. In 'Mixed-mode Delamination Criteria' section, a comprehensive criterion for mixed-mode delamination under fatigue is presented and simulation results are successfully compared with the experimental data available. Finally, some concluding remarks are given in the last section.

## REVIEW OF STATIC INTERFACE DAMAGE MODEL

The interface is a surface entity which ensures the transfer of stress and displacement between two adjacent layers as shown in Figure 1. This modeling coupled with the damage mechanics makes it possible to take into account the phenomenon of delamination which can occur during the mechanical loading of structural parts. The relative displacement of one layer to other layer can be written as follows:

$$ U = [U] = U^{+} - U^{-} = U_{1}N_{1} + U_{2}N_{2} + U_{3}N_{3} \qquad (3) $$
<figure><img src="image_2.png" /><figcaption>Figure 1. Interface between plies.</figcaption></figure>

The deterioration of the interface is taken into account by three internal damage variables ($d_1$, $d_2$, and $d_3$). It should be noted that in compression, there will be no damage. The relationship between stress and displacement in orthotropic plane of axis can be expressed as follows:

$$ \begin{pmatrix} \sigma_{13} \\ \sigma_{23} \\ \sigma_{33} \end{pmatrix} = \begin{pmatrix} k_1^0(1-d_1) & 0 & 0 \\ 0 & k_2^0(1-d_2) & 0 \\ 0 & 0 & k_3^0(1-d_3) \end{pmatrix} \begin{pmatrix} U_1 \\ U_2 \\ U_3 \end{pmatrix}, \quad (4) $$

where, $k_1^0$, $k_2^0$, and $k_3^0$ are interface rigidities associated to damage variables in orthogonal directions. The thermodynamic model is built by taking into account the three possible modes of delamination. Three different damage variables can be distinguished according to the three modes of failure. The thermodynamic forces combined with the damage variables associated to the three modes of delamination are (Allix et al., 1998):

$$ Y_{d_3} = \frac{1}{2} \frac{\langle \sigma_{33} \rangle_+^2}{k_3^0 (1 - d_3)^2}, \quad Y_{d_1} = \frac{1}{2} \frac{\sigma_{13}^2}{k_1^0 (1 - d_1)^2}, \quad Y_{d_2} = \frac{1}{2} \frac{\sigma_{32}^2}{k_2^0 (1 - d_2)^2}, \quad (5) $$

where $\langle X \rangle_+$ represents the positive part of X. It is supposed that the three different damage variables corresponding to three modes of failure are very strongly coupled and governed by equivalent single energy release rate function as follows (Allix et al., 1998):

$$ \underline{Y}(t) = \max_{\tau \le t} \left( \left( (Y_{d_3})^\alpha + (\gamma_1 Y_{d_1})^\alpha + (\gamma_2 Y_{d_2})^\alpha \right)^{1/\alpha} \right), \quad (6) $$

where, $\gamma_1$ and $\gamma_2$ are coupling parameters and $\alpha$ a material parameter which governs the damage evolution in mixed mode. The static damage evolution law is then defined by the choice of a material function as follows (Allix et al., 1998):

if $[(d_{3S} < 1) \text{ and } (\underline{Y} < Y_R)]$

then

$$ d_{1S} = d_{2S} = d_{3S} = \omega(\underline{Y}) \quad (7) $$

else

$$ d_{1S} = d_{2S} = d_{3S} = d_c $$
An isotopic damage evolution law is supposed here in order to simulate the delamination process in composite laminates. The damage function $\omega(Y)$ is selected in the form (Allix et al., 1998):

$$ \omega(\underline{Y}) = \left[ \frac{n}{n+1} \frac{\langle \underline{Y} - Y_o \rangle_+}{Y_c - Y_o} \right]^n \quad (8) $$

where $Y_o$ is the threshold damage energy, $Y_c$, critical damage energy, $n$, characteristic function of material, higher values of $n$ corresponds to brittle interface, $Y_R$, damage energy corresponding to $d_c$ failure $Y_R = Y_o + (n+1)/n d_c^{1/n} (Y_c - Y_o) (0 < d_c \le 1)$.

A simple way to identify the propagation parameters is to compare the mechanical dissipation yielded by two approaches of damage mechanics and linear elastic fracture mechanics (LEFM). In the case of pure mode situations when the critical energy release rate reaches its stabilized value at the propagation, it is denoted by $G_c$. Comparison of dissipations between fracture mechanics and damage mechanics approaches leads to (Allix et al., 1998):

$$ G_{IC} = Y_c; \quad G_{IIC} = \frac{Y_c}{\gamma_1}; \quad G_{IIIC} = \frac{Y_c}{\gamma_2}. \quad (9) $$

In order to satisfy the energy balance principle of LEFM, the area under the curve of stress-displacement curve for the whole debonding process (DP) obtained through damage mechanics formulation is set equal to critical energy release rate $G_{iC}$, see Figure 2. Following relations for the modes I, II, and III, critical energy release rates can be written:

$$ G_{IC} = \int_{DP} \sigma_{33} \, dU_3, \quad G_{IIC} = \int_{DP} \sigma_{13} \, dU_1, \quad G_{IIIC} = \int_{DP} \sigma_{23} \, dU_2 \quad (10) $$

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Constitutive model, mode I (right) and mode II (left) (the constitutive behavior of mode III is similar to mode II).</figcaption></figure>
For mixed-mode loading situation, a standard LEFM model is recovered as follows:

$$ \left(\frac{G_I}{G_{IC}}\right)^\alpha + \left(\frac{G_{II}}{G_{IIC}}\right)^\alpha + \left(\frac{G_{III}}{G_{IIIC}}\right)^\alpha = 1 \quad (11) $$

In a general mixed-mode DP, the global fracture energy can be computed as follows:

$$ G_{CT} = G_I + G_{II} + G_{III} \quad (12) $$

## FATIGUE INTERFACE DAMAGE MODEL

### Principle of the Modeling

Some assumptions are made here to simplify the numerical calculation procedure for fatigue delamination modeling. The actual applied cyclic load is oscillating between maximum and minimum values as shown in Figure 3. Hence in the case of high cycle fatigue, the load applied numerically to the structure will be equal to the maximum value of the actual load cycle, see Figure 3. In the case of interface debonding model, Robinson et al. (2005) introduced relative displacement based fatigue damage evolution law.

<figure><img src="image_4.png" /><figcaption>Figure 3. Envelope of applied cyclic load between maximum and minimum values.</figcaption></figure>
However in this article, a fatigue model based on damage energy release rates is proposed as follows:

$$ \dot{d}_F = \frac{\partial d_F}{\partial t} = \begin{cases} g\left(d, \frac{Y}{Y_C}\right) \frac{\dot{Y}(t)}{Y_C} & \text{if } \dot{Y} \ge 0 \text{ and } f \ge 0 \\ 0 & \text{if } \dot{Y} < 0 \text{ or } f < 0 \end{cases} \quad (13) $$

where f is a damage loading function and defines the threshold of fatigue delamination growth. This function can be written in terms of damage energy release rate as $f = \underline{Y} - Y^*$, where $Y^*$ is the threshold damage energy release rate, and damage will grow only and only if $f \ge 0$. This threshold value $Y^*$ is assumed to be zero for all computations done in this article. Here g is a dimensionless function and it depends on damage energy release rate $\underline{Y}$, its critical value $Y_C$, and on total damage.

Since the damage growth defined by Equation (13) is in rate format, it should be integrated over each time increment in the numerical analysis, in order to obtain the damage at the end of the increment. The damage variable at the end of a time increment $\Delta t$ can be written as (Peerlings et al., 2000):

$$ d_F(t + \Delta t) = d_F(t) + \int_t^{t+\Delta t} \dot{d}_F dt = d_F(t) + \sum_{n=N}^{N+\Delta N} P(d, \underline{Y}) \quad (14) $$

Here $t$ and $t + \Delta t$ are the times corresponding to end of cycles $N$ and $N + \Delta N$, respectively. Here $P(d, \underline{Y})$ represents the small variation of damage $d$ within one cycle compared to the total number of cycles and is expressed through the following form:

$$ P(d, \underline{Y}_{\max}, \underline{Y}_{\min}) = \int_{\underline{Y}_{\min}}^{\underline{Y}_{\max}} g\left(d, \frac{\underline{Y}}{Y_C}\right) \frac{d\underline{Y}}{Y_C} = \frac{C \times (1 - R^2)^{1+\beta}}{1+\beta} e^{\lambda d} \left[ \frac{\underline{Y}_{\max}}{Y_C} \right]^{1+\beta} \quad (15) $$

Here $R^2$ can be defined as follows:

$$ R^2 = \frac{\langle \underline{Y}_{\min} \rangle_+}{\underline{Y}_{\max}}, \quad (16) $$

where $\underline{Y}_{\min}$ corresponds to minimum and $\underline{Y}_{\max}$ corresponds to maximum load value during a cycle (envelope of the cyclic load, Figure 3). $\lambda$ is a constant parameter. $\beta$ and $C$ are functions of mode ratios which can also be expressed in more general form as $\beta(\phi)$ and $C(\phi)$. Here $d$ is total damage, and subscript $i$ (indicating any specific mode of failure) is omitted for the
sake of simplicity. For a mode mixture, comprising modes I and II, one can define a local definition of $\phi$ in damage mechanics formulation as: $\phi = Y_{d1}/(Y_{d1} + Y_{d3})$. Where $Y_{d3}$ and $Y_{d1}$ are damage energy release rates for modes I and II, respectively, and are already defined through Equation (5). Similarly, one can also write a global definition of $\phi$ in fracture mechanics formulation as: $\phi = G_{II}/(G_{II} + G_I)$, where $G_I$ and $G_{II}$ are modes I and II energy release rates, respectively.

Hereafter, for the convenience in writing, Equation (15) is rewritten and $\underline{Y}_{max}$ is replaced by $\underline{Y}$:

$$P(d, \underline{Y}) = \frac{C_R}{1 + \beta(\phi)} e^{\lambda d} \left[ \frac{\underline{Y}}{Y_C} \right]^{1+\beta(\phi)} \quad (17)$$

The effect of load ratio 'R,' for loads varying between maximum and minimum values, is taken into account by $C_R$:

$$C_R = C(\phi) \times (1 - R^2)^{1+\beta(\phi)} \quad (18)$$

The sum over the cycle numbers in Equation (14) can be approximated by using numerical integration schemes like trapezoidal rule or Simpson's rule for definite integrals (Hamming, 1987). Here, trapezoidal rule is used by estimating the average of the integrals evaluated at the beginning and end of the increment multiplied by the number of cycles in the increment $\Delta N$.

$$d_F(N + \Delta N) = d_F(N) + \frac{1}{2}[P(d(N + \Delta N), \underline{Y}(N + \Delta N)) + P(d(N), \underline{Y}(N))]\Delta N \quad (19)$$

The delamination crack growth under high cycle fatigue can be considered as a combination of delamination due to the quasi-static loading and due to the cyclic variation of the loading. The variation of damage under cyclic loading is expressed above. Similarly, the variation of static damage evolution as a function of loading cycles can be expressed as follows:

$$\left\{ \begin{array}{l} d_S(N + \Delta N) = d_S(N) + \left[ \frac{n}{n+1} \frac{1}{Y_C - Y_o} \right]^n \left[ \left\langle \underline{Y}(N + \Delta N) - Y_o \right\rangle_+^n - \left\langle \underline{Y}(N) - Y_o \right\rangle_+^n \right] \\ \text{if } \underline{Y}(N + \Delta N) \ge \underline{Y}(N) \end{array} \right. \quad (20)$$

where $\underline{Y}(N + \Delta N)$ and $d_S(N + \Delta N)$ correspond to end of cycles $N + \Delta N$. $\underline{Y}(N)$ and $d_S(N)$ correspond to end of cycles $N$. If $d(N)$ is the total damage at the end of cycles $N$, then the total damage after the increment of number of
cycles $\Delta N$ can be evaluated by combining Equations (19) and (20) for cyclic loading as follows:

$$ \left\{ \begin{array}{l} d(N + \Delta N) = d(N) + \frac{1}{2} \underbrace{[P(d(N + \Delta N), \underline{Y}(N + \Delta N)) + P(d(N), \underline{Y}(N))] \Delta N}_{\text{Fatigue delamination}} \\ \underbrace{\left[ \frac{n}{n+1} \frac{1}{Y_C} \right]^n \left[ \langle \underline{Y}(N + \Delta N) - Y_o \rangle_+^n - \langle \underline{Y}(N) - Y_o \rangle_+^n \right]}_{\text{Static delamination}} \quad \text{if } \underline{Y}(N + \Delta N) \ge \underline{Y}(N) \end{array} \right. \quad (21) $$

Relation (21) is a nonlinear equation in terms of damage variable $d(N + \Delta N)$, because it appears on both sides of the equation. This can be solved iteratively by applying the standard Newton-Raphson method to Equation (21) (Hamming, 1987). One can avoid solving this equation iteratively by replacing $d(N + \Delta N)$ on the right-hand side of the Equation (21), with predictor $d^p$ based on forward Euler step (Peerlings et al., 2000) and defined as,

$$ d^p = d(N) + P(d(N), \underline{Y}(N)) \Delta N \quad (22) $$

In this article, predictor integration scheme is used for all the simulations.

## Identification of the Model

The fatigue damage model presented here is implemented in finite element code in Cast3M (CEA) (Verpeaux et al., 1988). The effectiveness of the fatigue damage model is tested by the finite element simulations of mode I, mode II, and mixed-mode delamination tests.

Two dimensional meshes are used to model the beam arms and interface elements (Beer, 1985) are employed for the modeling of DP. The material used is unidirectional HTA/6376C carbon/epoxy laminate and its properties taken from Asp et al. (2001) are given in Table 1. The experimental results of Asp et al. (2001) are also used for the comparison with the predicted fatigue delamination behavior. The specimen with total length, $L=150$ mm; width, $b=20$ mm; initial crack, $a_0=35$ mm; and thickness, $h=3.1$ mm is used for the simulations. For all the simulations conducted in this article, load is applied in two steps. In the first loading phase, load is applied monotonically to a maximum load point, but this maximum load point should be low enough to avoid the static delamination. Then in the second step, load oscillates between maximum and zero values to simulate the fatigue loading condition.
From the corrected critical energy release rates at propagation (Figure 4) and from the relationship existing between fracture mechanics and damage mechanics, Equation (9), one can find the critical energies $Y_C$ and the coupling coefficient $\gamma_1$. Without any further information on mode III interlaminar fracture, one can have $\gamma_2 = \gamma_1$ (Allix et al., 1998). The identification of $\alpha$, which governs the energy release rate in case of mixed-mode delamination, is done experimentally by performing tests for different mode mixtures and the same is directly used in interface damage modeling. As a general case, normally the value of $\alpha$ is chosen between 1.0 and 2.0 ($1.0 \le \alpha \le 2.0$) (Harper and Hallet, 2008). The identified interface parameters for UD HTA/6376C using Equation (9) are given in Table 2.

<table><caption>**Table 1. Material properties for UD HTA/6376C.**</caption><tr><td>E<sub>11</sub> (GPa)</td><td>120</td></tr><tr><td>E<sub>22</sub> = E<sub>33</sub> (GPa)</td><td>10.5</td></tr><tr><td>G<sub>12</sub> = G<sub>13</sub> (GPa)</td><td>5.25</td></tr><tr><td>G<sub>23</sub> (GPa)</td><td>3.48</td></tr><tr><td>ν<sub>12</sub> = ν<sub>13</sub></td><td>0.30</td></tr><tr><td>ν<sub>23</sub></td><td>0.51</td></tr></table>

Source: Asp et al. (2001).

<figure><img src="image_5.png" /><figcaption>**Figure 4. Energy release rates at propagation for UD HTA/6376C.**</figcaption></figure>
For the model proposed by Allix, the following relationship is derived to calculate the interface rigidities:

$$ \text{with } i \in \{1,2\} $$ $$ k_3^0 = \frac{(2n + 1)^{\frac{(2n+1)}{n}}}{8n(n + 1)Y_C} \sigma_{33}^2, \quad k_i^0 = \frac{\gamma_i(2n + 1)^{\frac{(2n+1)}{n}}}{8n(n + 1)Y_C} \sigma_{3i}^2 \quad (23) $$

Now using Equation (23), the interfacial rigidities in normal and shear directions can be identified. The values for maximum interfacial stresses are assumed as 30 MPa (Robinson et al., 2005). Taking $n = 0.5$ and using Equation (23), the normal and shear rigidities are calculated as $k_0^3 = 9.3 \times 10^3 \text{MPa/mm}$ and $k_1^0 = 2.4 \times 10^3 \text{MPa/mm}$, respectively, for UD HTA/6376C.

Figure 5 presents the variation of normal interfacial stress $\sigma_{33}$ with respect to relative displacement $U_3$, and area under the curve is equal to critical energy release rate $G_{IC}$. Similarly, Figure 6 presents the variation of shear

<table><caption>**Table 2. Interfacial material parameters for UD HTA/6376C.**</caption><thead><tr><th>Interface</th><th>Y<sub>c</sub> (kJ/m<sup>2</sup>)</th><th>γ<sup>1</sup></th></tr></thead><tbody><tr><td>0°/0°</td><td>0.26 ± 0.01</td><td>0.25</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Variation of $\sigma_{33}$ with respect to $U_3$ for UD HTA/6376C.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Variation of $\sigma_{13}$ with respect to $U_1$ for UDHTA/6376C.</figcaption></figure>

interfacial stress $\sigma_{13}$ with respect to relative displacement $U_1$, and area under the curve is equal to critical energy release rate $G_{IIC}$. As a general case, proper values of interfacial stresses and interfacial rigidities are found by making a comparison between numerical and experimental results under static loading conditions.

Interfacial fatigue damage parameters $\lambda, \beta$, and $C$ are directly identified from experimental results for values giving the good fit of numerical results. The influence of three fatigue damage parameters $\lambda, \beta$, and $C$ on Paris plot under mode I fatigue delamination condition is discussed here. The three different values of parameters $\lambda, \beta$, and $C$ are selected in such a way that the simulation results always fall inside the scattered data. Figure 7 shows the influence of $C$ on the crack growth rates as a function of normalized critical energy release rate. From Figure 7, it is clear that varying the value of $C$ strongly affects the cyclic crack growth rates, but the influence on the slope of the curve is not very significant.

Figure 8 presents the effect of $\beta$ on the Paris plot, different values of $\beta$ considerably affect the slope of the curve. However, the variations in $\lambda$ do not seem to affect the slope of the curve significantly, see Figure 9. By taking into account the influence of these fatigue damage parameters on the linear Paris plot behavior, their values can be determined for different modes of failure. The value of $\lambda$ is fixed to 0.5 for all the simulations, and then values of $\beta$ and $C$ are determined to suit the slope of Paris plot for different modes of failure.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Paris plot under pure mode I fatigue delamination for different values of C.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Paris plot under pure mode I fatigue delamination for different values of β.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Paris plot under pure mode I fatigue delamination for different values of $\lambda$.</figcaption></figure>

The delamination toughness values obtained through experiments (Asp et al., 2001) and the associated identified interface and fatigue parameters for UD HTA/6376C are given in Table 3.

## Finite Element Simulations

Fatigue delamination simulations are performed in pure mode I, pure mode II, and mixed mode (for mixed mode, modes I, and II, the components are equal ($\phi = 0.5$)). Two-dimensional meshes comprising of four nodes plane strain elements are used to model the beam arms and interface elements are employed for the modeling of DP. Linear elements of constant size 0.2 mm are used along the length of the arm and four elements are used in thickness direction. No friction effects are considered for the modeling and simulations in this work. For pure mode I, specimen arms are loaded with opposing moments, Figure 10. The opposing moment condition gives the mode I energy release rate that is independent of crack length and therefore fatigue loading at a constant applied moment $M$ results in a constant crack growth rate. Similarly, the loading condition for pure mode II is shown in Figure 11, and that for mixed mode ($\phi = 0.5$) is shown in
<table><caption>**Table 3. Delamination toughness values for UD HTA/6376C and associated Fatigue parameters.**</caption><thead><tr><th>Test method</th><th>G<sub>c</sub> (kJ/m<sup>2</sup>)</th><th>Interface</th><th>Fatigue parameters</th></tr></thead><tbody><tr><td>Mode I<br>φ = 0</td><td>0.26 ± 0.01</td><td>n = 0.5 Y<sub>O</sub> = 0. kJ/m<sup>2</sup><br>Y<sub>C</sub> = 0.26 ± 0.01 kJ/m<sup>2</sup><br>k<sub>0</sub><sup>3</sup> = 9.3 × 10<sup>3</sup> MPa/mm</td><td>λ = 0.5<br>β(φ) = 2.0<br>C(φ) = 6.0 × 10<sup>-4</sup></td></tr><tr><td>Mode II<br>φ = 1.0</td><td>1.002 ± 0.063</td><td>n = 0.5 Y<sub>O</sub> = 0 kJ/m<sup>2</sup><br>Y<sub>C</sub> = 0.26 ± 0.01 kJ/m<sup>2</sup><br>k<sub>0</sub><sup>1</sup> = 2.4 × 10<sup>3</sup> MPa/mm<br>γ<sub>1</sub> = 0.25</td><td>λ = 0.5<br>β(φ) = 2.0<br>C(φ) = 6.0 × 10<sup>-3</sup></td></tr><tr><td>Mixed-mode<br>φ = 0.5</td><td>0.447 ± 0.023</td><td>n = 0.5 Y<sub>O</sub> = 0 kJ/m<sup>2</sup><br>Y<sub>C</sub> = 0.26 ± 0.01 kJ/m<sup>2</sup><br>k<sub>0</sub><sup>3</sup> = 9.3 × 10<sup>3</sup> MPa/mm<br>k<sub>0</sub><sup>1</sup> = 2.4 × 10<sup>3</sup> MPa/mm<br>γ<sub>1</sub> = 0.25 α = 2.0</td><td>λ = 0.5<br>β(φ) = 3.5<br>C(φ) = 8.0 × 10<sup>-2</sup></td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Specimen under pure mode I loading condition.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Specimen under pure mode II loading condition.</figcaption></figure>

Figure 12. For a mode ratio of 50% mode II, the ratio $\rho$ between the two applied moments is calculated as follows (Williams, 1988):

$$ \rho = \frac{1 - \frac{\sqrt{3}}{2}}{1 + \frac{\sqrt{3}}{2}}. \qquad (24) $$

The energy release rate for pure mode I is (Williams, 1988):

$$ G_I = \frac{M^2}{bEI}, \qquad (25) $$
<figure><img src="image_13.png" /><figcaption>Figure 12. Specimen under mixed-mode loading condition.</figcaption></figure>

where the width of specimen is *b*, *E* the longitudinal flexural Young's modulus, and *I* the second moment of area of the specimen's arm.

The energy release rate for pure mode II is (Williams, 1988):

$$G_{II} = \frac{3 M^2}{4 b EI} \qquad (26)$$

In which, $M = cP/2$, see Figure 11. The value of *c* is 30 mm. For mixed-mode case (Williams, 1988):

$$G_I = G_{II} = \frac{3 M^2}{4(1 + \sqrt{3}/2)^2 b EI} \qquad (27)$$

Figure 13 shows the Paris plot behavior for modes I, II, and mixed-mode (50% mode II) delamination results. Simulation results are found in good agreement with experimental results (Asp et al., 2001).

Simulations presented above are based on the assumption that cyclic load is varying between maximum and zero values. However, in many practical situations, load is varying between maximum and minimum values. The effect of this type of load variation can be taken into account by load ratio *R* (Equation (18)).

In order to check the effectiveness of Equations (18), experimental results of Martin and Murri (1990) on delamination growth have been selected for comparison purposes. Martin and Murri (1990) performed fatigue delamination growth experiments for two different load ratios, $R=0.1, 0.5$, on unidirectional AS4/PEEK laminated composite material. The material properties for AS4/PEEK laminate are given in Table 4 (Jen and Lee, 1998).

The proposed fatigue damage evolution law with load ratio *R* effect is tested for pure mode I (DCB) and pure mode II (4ENF) loading conditions. The associated delamination toughness values (Martin and Murri, 1990) and fatigue parameters are given in Table 5. The identification procedure for different values in Table 5 is the same, as explained in 'Principle of the Modeling' section. The geometry of the specimen for DCB test has the following dimensions, $L=140$ mm, $b=25.4$ mm, $a_0=50$ mm, and $h=4.6$ mm, whereas the geometry of the specimen for 4ENF test has the
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Comparison between numerical results and experiments. (Asp et al., 2001).</figcaption></figure>

<table><caption>**Table 4.** Material properties for UD AS4/PEEK laminate.</caption><tr><td>E<sub>11</sub> (GPa)</td><td>140.35</td></tr><tr><td>E<sub>22</sub> = E<sub>33</sub> (GPa)</td><td>9.44</td></tr><tr><td>G<sub>12</sub> = G<sub>13</sub> (GPa)</td><td>5.403</td></tr><tr><td>G<sub>23</sub> (GPa)</td><td>3.48</td></tr><tr><td>ν<sub>12</sub> = ν<sub>13</sub></td><td>0.253</td></tr><tr><td>ν<sub>23</sub></td><td>0.51</td></tr></table>

Source: Jen and Lee (1998).

following dimensions, $L = 102$ mm, $b = 25.4$ mm, $a_0 = 25$ mm, and $h = 4.6$ mm (Martin and Murri, 1990).

The results of simulations of fatigue delamination growth along with experimental results for modes I and II loading conditions are shown in Figures 14 and 15, respectively. It should be noted that Martin and Murri (1990) used British system of units. That is why the same system of units is adopted here for comparison purposes; and for similar reasons the critical energy release rates are also given in both types of system of units in Table 5.
<table><caption>**Table 5.** Delamination toughness values of UD AS4/PEEK and associated fatigue parameters.</caption><thead><tr><th>Test method</th><th>G<sub>c</sub> (kJ/m<sup>2</sup>)</th><th>Interface</th><th>Fatigue parameters</th></tr></thead><tbody><tr><td>Mode I<br>φ = 0.</td><td>1.69–2.47<br>(9.65–14.14 in-lb/in<sup>2</sup>)</td><td>n = 0.5 Y<sub>O</sub> = 0 kJ/m<sup>2</sup><br>Y<sub>C</sub> = 1.69 – 2.47kJ/m<sup>2</sup><br>k<sub>0</sub><sup>3</sup> = 9.3 × 10<sup>3</sup> MPa/mm</td><td>λ = 0.5<br>β(φ) = 2.5<br>C(φ) = 2.0 × 10<sup>-3</sup></td></tr><tr><td>Mode II<br>φ = 1.0</td><td>2.49–3.76<br>(14.2–21.5 in-lb/in<sup>2</sup>)</td><td>n = 0.5 Y<sub>O</sub> = 0 kJ/m<sup>2</sup><br>Y<sub>C</sub> = 1.69 – 2.47kJ/m<sup>2</sup><br>k<sub>0</sub><sup>1</sup> = 2.4 × 10<sup>3</sup> MPa/mm<br>γ<sub>1</sub> = 0.7</td><td>λ = 1.5<br>β(φ) = 2.5<br>C(φ) = 3.0 × 10<sup>-2</sup></td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Comparison between experimental (Martin and Murri, 1990) and numerical results for pure mode I fatigue delamination with different load ratios.</figcaption></figure>

A good agreement is found between numerical and experimental results for the mode I fatigue-driven delamination. While the comparison between numerical and experimental results for mode II loading condition is not as good as for mode I, they still can be considered to be in acceptable range. But one could give more conclusive comments, if more experimental results of fatigue delamination with different *R* would be available for comparison. Similarly, more experimental data on fatigue delamination with different *R* values will also help to further improve the relation for *C*<sub>R</sub>, if needed.
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Comparison between experimental (Martin and Murri, 1990) and numerical results for pure mode II fatigue delamination with different load ratios.</figcaption></figure>

## MIXED-MODE DELAMINATION CRITERIA

### Principle of the Modeling

Experimental results show that the Paris plot behavior can be expressed by Equation (1). Different authors tried to develop a relation between these parameters (*m* and *B*) and mode-mixture *φ*, so that if values of these parameters are known for certain mode mixtures, then the values for others can be predicted. Value of *φ* varies between 0 and 1, where 0 indicates the pure mode I state and 1 represents the pure mode II state. For glass fiber reinforced plastics (GFRP) materials, Kenane and Benzeggagh (1997) showed that variation of these two parameters with respect to *φ* is monotonic. For carbon fiber reinforced plastics (CFRP) materials, Blanco et al. (2004) showed a nonmonotonic type trend and also developed a relation between parameters (*B*, *m*) and mode-mixture *φ*.
Tumino and Cappello (2007) also used the nonmonotonic equation to establish a relation between fatigue damage parameters and mode-mixture $\phi$. In this article, the same approach is used to test the predictive power of proposed fatigue damage evolution law for different mode mixtures.

Expressing Equation (16) in a more general form:

$$ egin{cases} P(d, \underline{Y}) = K(\phi)(1-R)^{\eta(\phi)} e^{\lambda d} \left[ \frac{\underline{Y}}{Y_C} \right]^{\eta(\phi)} & . \\ \\ \text{where } K(\phi) = \frac{C(\phi)}{1+\beta(\phi)}, \quad \eta(\phi) = 1+\beta(\phi) & \end{cases} \quad (28) $$

In order to establish a nonmonotonic relation, values of fatigue damage parameters for three different fracture mechanics test should be known in advance. Hence, for known values of pure modes I and II and for any mode mixture, one can define $K_I$, $K_{II}$, $K_{mix}$, $\eta_I$, $\eta_{II}$, $\eta_{mix}$, respectively. Then, new values of $\eta(\phi)$ and $K(\phi)$ corresponding to any other mode-mixture $\phi$ can be determined by using following equations (Blanco et al., 2004; Tumino and Cappello, 2007):

$$ \eta(\phi) = A_{1\eta}\phi^2 + A_{2\eta}\phi + A_{3\eta} \quad (29) $$

$$ \ln(K(\phi)) = A_{1K}\phi^2 + A_{2K}\phi + A_{3K} \quad (30) $$

where the expressions for $A_{iK}$ and $A_{i\eta}$, $i=1,2,3$ are given below (Tumino and Cappello, 2007).

$$ A_{1\eta} = \frac{\eta_I - \eta_{mix} + (\eta_{II} - \eta_I)\phi}{\phi - \phi^2}, \quad A_{2\eta} = \frac{\eta_{mix} - \eta_I + A_{1\eta}\phi^2}{\phi}, \quad A_{3\eta} = \eta_I \quad (31) $$

$$ egin{cases} A_{1K} = \frac{\ln(K_I) - \ln(K_{mix}) + [\ln(K_{II}) - \ln(K_I)]\phi}{\phi - \phi^2} \\ \\ A_{2K} = \frac{\ln(K_{mix}) - \ln(K_I) + A_{1K}\phi^2}{\phi}, \quad A_{3K} = \ln(K_I) \end{cases} \quad (32) $$

Equations (29)–(32) use the same analogy as used by Tumino and Cappello (2007) in their studies. Here Equations (29) and (30) are used to find the values of $\eta(\phi)$ and $K(\phi)$ for any value of mode mixtures, and then corresponding values of $C(\phi)$ and $\beta(\phi)$ are calculated from the definition of Equation (28). The critical values for the energy release rate are connected to $\phi$ through the following equation (Tumino and Cappello, 2007):

$$ G_C = G_{IC} + (G_{IIC} - G_{IC})\phi^2 \quad (33) $$
In the previous section, the values of fatigue parameters ($C(\phi)$ and $\beta(\phi)$) for pure mode I ($\phi = 0$.), mode II ($\phi = 1$.), and for mode-mixture $\phi = 0.5$ are already found for UD HTA/6376C laminate. The new values of fatigue parameters for any mode-mixture $\phi$ can be found by using Equations (29)–(32). In this study, mode mixture of $\phi = 0.25$ and $\phi = 0.75$ are selected for the simulations in addition to those already described in the section ‘Fatigue Interface Damage Model.’ The loading scheme used is the same as shown in Figure 12. Now the values of $\rho$ and the energy release rate $G$ for the mode-mixture $\phi = 0.25$ can be calculated as follows (Williams, 1988):

$$ \rho(\phi=0.25) = -\frac{1}{5}, \quad G_{I(\phi=0.25)} = \frac{9M^2}{bEI}, \quad G_{II(\phi=0.25)} = \frac{3M^2}{bEI} \quad (34) $$

And similarly for the mode-mixture $\phi = 0.75$, $G$ can be expressed as follows (Williams, 1988):

$$ \rho(\phi=0.75) = \frac{1}{3}, \quad G_{I(\phi=0.75)} = \frac{M^2}{bEI}, \quad G_{II(\phi=0.75)} = \frac{3M^2}{bEI} \quad (35) $$

The trends obtained for $B$ and $m$ as a function of $\phi$ from experiments by Blanco et al. (2004) are found to be in reasonably good agreement with numerical results (Figures 16 and 17). The predicted behavior of fatigue parameters $\beta(\phi)$ and $C(\phi)$, as functions of mode ratio, $\phi$ are given in Figures 18 and 19. The predicted behaviors of Paris plots for mixed mode 0.75 and 0.25 obtained through simulations are given in Figure 20.

<figure><img src="image_17.png" /><figcaption>Figure 16. Paris law coefficient $B$ vs $\phi$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Paris law exponent *m* vs *φ*.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 18.** Prediction of fatigue parameter *β*(*φ*) vs *φ*.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 19.** Prediction of fatigue parameter $C(\phi)$ vs $\phi$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 20.** Predicted behavior of Paris plots for UD HTA/6376C laminate for 75% mode II and 25% mode II.</figcaption></figure>
CONCLUSION

In this article, a comprehensive elastic fatigue damage model based on
damage energy release rate is presented. The effectiveness of the proposed
model has been tested by performing the finite element simulations of dif-
ferent fracture mechanics specimens under cyclic loading conditions. A
slightly modified formulation is also proposed, in order to take into account
the ‘R’ effect. The linear Paris plot behaviors predicted by the proposed
model for pure mode I, pure mode II and for mixed-mode specimens
under fatigue loading condition are found in good agreement with experi-
mental results.

Then, the same fatigue damage model is modified by introducing the
nonmonotonic behavior, for the prediction of fatigue damage parameters
for different mode mixtures. Results generated by this law for Paris plot
parameters are found to be in good accordance with experimental results.

Hence, the proposed fatigue damage model is not only able to reproduce
the linear Paris plot behavior for composite laminates, but also capable of
predicting the fatigue behavior for different mode mixtures.

ACKNOWLEDGMENTS

The authors thank HEC (Higher education commission of Pakistan) and
its collaborating organisation SFERE (Société française d'exportation des
ressources éducatives) for providing financial support.

REFERENCES

Allix, O. and Ladevèze, P. (1992). Interlaminar Interface Modelling for the Prediction of Delamination, *Composite Structures*, **22**: 235-242.

Allix, O., Ladevèze, P. and Corigliano, A. (1995). Damage Analysis of Interlaminar Fracture Specimens, *Composite Structures*, **31**: 61–74.

Allix, O., Ladevèze, P., Gornet, L., Léveque, D. and Perret, L. (1998). A Computational Damage Mechanics Approach for Laminates: Identification and Comparison with Experimental Results. In: Voyiadjis, G.Z., Ju, J.W., Chaboche, J.L. (eds), *Damage Mechanics in Engineering Materials, Studies in Applied Mechanics*, Vol. 46, Elsevier, pp. 481–500.

Asp, L.E., Sjogren, A. and Greenhalgh, E.S. (2001). Delamination Growth and Thresholds in a Carbon/Epoxy Composite Under Fatigue Loading, *Journal of Composite Technology and Research*, **23**: 55–68.

Beer, G. (1985). An Isoparametric Joint/Interface Element For Finite Element Analysis, *International Journal for Numerical Methods in Engineering*, **21**: 585–600.

Blanco, N., Gamstedt, E.K., Asp, L.E. and Costa, J. (2004). Mixed-Mode Delamination Growth In Carbon-Fibre Composite Laminates Under Cyclic Loading, *International Journal of Solids and Structures*, **41**: 4219–4235.
Chaboche, J.L., Girard, R. and Levasseur, P. (1997). On the Interface Debonding Models, *International Journal of Damage Mechanics*, **6**: 220–257.

Davies, P., Cantwell, W., Moulin, C. and Kausch, H.H. (1989). A Study of Delamination Resistance of IM6/PEEK Composites, *Composites Science and Technology*, **36**: 153–166.

de Morais, A.B. and Pereira, A.B. (2007). Application of the Effective Crack Method to Mode I and Mode II Interlaminar Fracture of Carbon/Epoxy Unidirectional Laminates, *Composites Part A*, **38**: 785–794.

Hamming, R.W. (1987). *Numerical Methods for Scientists and Engineers, 2nd edn*, New York, Dover Publications.

Harper, P.W. and Hallet, S.R. (2008). Cohesive Zone Length in Numerical Simulations of Composite Delamination, *Engineering Fracture Mechanics*, **75**: 4774–4792.

Hojo, M. and Gustafson, C.-G. (1987). Delamination Fatigue Crack Growth in Unidirectional Graphite/Epoxy Laminates, *Journal of Reinforced Plastics and Composites*, **6**: 36–52.

Jen, M.R.H. and Lee, C.H. (1998). Strength and Life in Thermoplastic Composites Laminates Under Static and Fatigue Loads. Part I: Experimental, *International Journal of Fatigue*, **20**: 605–615.

Juntti, M., Asp, L.E. and Olsson, R. (1999). Assessment of Evaluation Methods for the Mixed-Mode Bending Tests, *Journal of Composites Technology and Research*, **21**: 37–48.

Kenane, M. and Benzeggagh, M.L. (1997). Mixed-Mode Delamination Fracture Toughness of Unidirectional Glass/Epoxy Composites under Fatigue Loading, *Composites Science and Technology*, **57**: 597–605.

Martin, R.H. and Murri, G.B. (1990). Characterization of Mode I and Mode II Delamination Growth and Thresholds in AS4/PEEK Composites. In: Garbo, S.P. (ed.), *Composite Materials: Testing and Design ASTM STP 1059*, Vol. 9, Philadelphia, American Society for Testing and Materials, pp. 251–270.

Paas, M.H.J.W., Breklemans, W.A.M. and Schreurs, P.J.G. (1993). A Continuum Approach to Brittle and Fatigue Damage: Theory and Numerical Procedures, *International Journal of Solid and Structures*, **30**: 579–599.

Peerlings, R.H.J., Breklemans, W.A.M., de Borst, R. and Geers, M.G.D. (2000). Gradient-Enhanced Damage Modelling of High-Cycle Fatigue, *International Journal of Numerical Methods in Engineering*, **49**: 1547–1569.

Robinson, P., Galvanetto, U., Tumino, D., Bellucci, G. and Violeau, D. (2005). Numerical Simulation of Fatigue-Driven Delamination Using Interface Elements, *International Journal for Numerical Methods in Engineering*, **63**: 1824–1848.

Tay, T.E. (2003). Characterization and Analysis of Delamination Fracture in Composites: An Overview of Developments from 1990 to 2001, *Applied Mechanics Review*, **56**: 1–32.

Tumino, D. and Cappello, F. (2007). Simulation of Fatigue Delamination Growth in Composites with Different Mode Mixtures, *Journal of Composite Materials*, **41**: 2415–2441.

Turon, A., Costa, J., Camanho, P.P. and Davila, C.G. (2007). Simulation of Delamination in Composites under High-Cycle Fatigue, *Composites Part A: Applied Science and Manufacturing*, **38**: 2270–2282.

Verpeaux, P., Charras, T. and Millard, A. (1988). Castem 2000: Une Approche Moderne Du Calcul Des Structures. Fouet, J.M. Ladevèze, P. Ohayon, R. (ed.), pp. 227–261. Available at: http://www-cast3m.cea.fr.

Williams, J.G. (1988). On the Calculation of Energy Release Rates for Cracked Laminates, *International Journal of Fracture*, **36**: 101–119.