A high-cyclic elastic fatigue damage model for carbon fibre epoxy matrix laminates with different mode mixtures

This paper presents the development of a fatigue damage model able to carry out simulation of evolution of delamination in the laminated composite structures under cyclic loadings. A classical interface damage evolution law, which is commonly used to predict static debonding process, is modified further to incorporate fatigue delamination effects due to high cycle loadings. The proposed fatigue damage model is identified using fracture mechanics tests DCB, ENF and MMB. Then a non-monotonic behaviour is used to predict the fatigue damage parameters able to carry out delamination simulations for different mode-mixtures. Linear Paris plot behaviour of the above mentioned fracture mechanics tests are successfully compared with available experimental data on HTA/6376C and AS4/PEEK unidirectional materials.


Introduction
For weight saving purposes in aircrafts, trains and ships, the use of sandwich structures made of composite laminate skins are no longer limited to secondary structure, but is expending to primary load bearing parts. Due to their laminated nature, composite laminate skins are prone to delamination failure under static and fatigue loadings [1,2]. This phenomenon consists of the separation of the adjacent layers of laminated composites. Delamination growth under fatigue loadings causes total collapse of the load-bearing properties of composite structures.
During the last twenty 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 [3][4][5][6]. Meso-scale is strongly connected to the laminate scale which lies between micro scale (fibre 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 fibres [5]. Few papers focus on delamination modelling under fatigue loadings in composite laminates [7][8][9].
In this paper a comprehensive interlaminar interface fatigue damage evolution law is proposed to model the delamination phenomena under high cycle fatigue loadings. When the fatigue failure occurs above 10 4 cycles, it is usually called High-cycle fatigue [9]. A state of the art of fatigue behaviour response of composite laminates is presented in [10].
The modelling was implemented in the finite element code Cast3 M (CEA, http://www-cast3m.cea.fr). The DCB, ENF and MMB tests were chosen to identify the proposed model for simulation of the crack growth in unidirectional carbon-fibre epoxymatrix materials. The proposed model takes its foundations from the classical static damage evolution law proposed by Allix and co-workers [5,6]. Main ideas for the fatigue crack growth modelling were first introduced by Peerlings et al. [11] and Pass et al. [12] for metallic parts. Robinson et al. [7] presented fatigue driven delamination for the laminated composites using the idea of [11] for the cyclic load, varying between maximum and zero values. In Ref. [7], the fatigue damage evolution is a function of relative displacement at the interface. However in the present work a damage energy release rate based formulation is proposed for the fatigue delamination.
The proposed fatigue damage evolution law requires the identification of certain fatigue parameters by comparing with experimental results. A non-monotonic behaviour [8,13] is used to predict the fatigue damage parameters for different mode-mixtures. This prediction method requires the fatigue damage parameters for three different states of loading to be known as a priori, i.e. for pure mode I, mode II and for mixed-mode. Once the parameters for three different loading conditions are known then the predictions can be made for fatigue delamination for different mode-mixtures /.
The proposed fatigue damage evolution law permits to reproduce the linear crack growth rates as obtained by using classical Paris law [14] for fracture mechanics tests. The Paris law can be described by Eq. (1) which relates the delamination increment per cycle with the cyclic variation of the energy release rate with two parameters B and m that can be obtained from experiment.
where a is the crack length, N is the number of cycles, DG = (G max À G min ) is the cyclic variation of the energy release rate with G max and G min represent the maximum and minimum values of energy release rates during the oscillation. G c is the fracture toughness of the material, B and m are constants and are determined experimentally.
In Damage Mechanics theory, 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 the combination of delamination due to the quasi-static loading and due to the cyclic variation of the loading [7], hence total damage evolution for three different modes of failure can be expressed as follows: Where the term d iS corresponds to delamination growth under static loading and d iF is related to fatigue loading. The paper is organised as follows, in Section 2, the classical damage model proposed by Allix and co-workers [5,6] for the prediction of delamination in laminated composites is recalled. In Section 3, the proposed fatigue damage model is presented and simulations results are given in Section 4. In Section 5, a comprehensive criterion for mixed-mode delamination under fatigue is presented and simulation results are successfully compared with available experimental data. Finally some concluding remarks are given in Section 6.

Review of static interface damage model
The interface is a surface entity which ensures the transfer of stress and displacements between two adjacent layers. This modelling coupled with 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 where N 1 , N 2 and N 3 represents the orthotropic directions of the interface. The deterioration of the interface is taken into account by three internal damage variables (d 1 , d 2 and d 3 ). The relationship between stress and displacement in orthotropic plane of axis can be expressed as: where, k where hXi + 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 are governed by equivalent single energy release rate function as follows [5]: where c 1 and c 2 are coupling parameters and a is 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: The damage function is selected in the form: where Y O is threshold damage energy, Y C is critical damage energy, n is Characteristic function of material (higher values of n corresponds to brittle interface). If Y O is initialized to any value x then one will have linear response of stress vs interface displacement and damage will be zero during this phase. Damage ''d'' will start to grow once the value x of Y O is reached, see Fig. 1. Y R is energy corresponding to rupture, i.e. d = 1.0 and can be expressed as: 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 denoted by G C , comparison of dissipations between Fracture Mechanics and Damage Mechanics approaches lead to [5,6]: 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 , and the following relations for the mode I, mode II and mode III critical energy release rates can be written: The area under the curve (as shown in Fig. 1 for mode I only) will always equal to critical energy G c for any mode of delamination. For mixed-mode loading situation, a standard LEFM model is recovered as: In a general mixed-mode debonding process, the global fracture energy can be computed as follows: 3. Fatigue interface damage model Before going into details of fatigue delamination interface modelling, some assumptions are made here, in order to simplify the numerical calculation procedure. The actual applied cyclic load is oscillating between zero and maximum value as shown in Fig. 2 (the case where load is varying between maximum and minimum values will be discussed later). Hence in 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 Fig. 2.
In case of interface debonding model, Robinson et al. [7] introduced relative displacement based fatigue damage evolution law. However, according to the thermodynamic formulation, a new fatigue interface model is proposed. The fatigue damage evolution law is based on damage energy release rate. The damage fatigue part is proposed as follows: 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 where Y Ã is threshold damage energy release rate and damage will grow only and only if f P 0. For all computations done in this article this threshold value Y Ã is assumed zero. Here g is a dimensionless function and depends on damage energy release rate Y, on its critical value Y C and on the total damage itself d (here the subscript i, representing different modes of failure, is omitted for the sake of simplicity). This function is selected of the following form: where k is a constant parameter. b and C are function of mode ratio and can also be expressed in more general form as b(/) and C(/). For a mode-mixture, comprised of mode I and mode II, one can write a local definition of / in Damage Mechanics formulation as, Where Y d3 and Y d1 are damage energy release rates for mode I and mode II and are already defined through Eq. (5). Similarly one can also write a global definition of / in fracture mechanics formulation as, / = G II /(G II + G I ). Where G I and G II are mode I and mode II energy release rates. Since the damage growth defined by Eq. (13) is in rate format, it should be integrated over each time increment on the numerical analysis in order to obtain the damage at the end of increment. The damage variable at the end of a time increment Dt can be written as [11]: Here Pðd; YÞ is the evolution of damage within one cycle due to fatigue and can be written as: Here the value of Y corresponds to the peak during a cycle, i.e. envelope of the cyclic load. The effect of load ratio ''R'', for loads varying between maximum and minimum values, can be taken into account by rewriting the above relation in more general form: Here ''R'' can be defined as the ratio of square root of minimum and maximum value of square root of Y during a cycle.Then the sum over the cycle numbers in Eq. (15) can be approximated by using numerical integration schemes like trapezoidal rule or Simpson's rule for definite integrals [15]. 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 DN. If t and t + Dt are the times corresponding to end of cycles N and N + DN respectively then: The delamination crack growth under high cycle fatigue can be considered as the 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: In this article predictor integration scheme is used for all the simulations.

Simulations and results of fatigue interface damage model
The fatigue damage model presented here is implemented in finite element code in Cast3 M (CEA). 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 comprised of 4 nodes plane strain elements are used to model the beam arms and interface elements [16] are employed for the modelling of debonding process, Fig. 3. It should be highlight that Finite element predictions using beam elements for arms lead to equivalent results. The material used is unidirectional HTA/6376C carbon/epoxy laminate and its properties taken from Ref. [14] are given in Table 1. The experimental results obtained in [14] are also used for the comparison with predicted fatigue delamination behaviour. 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 and is in accordance with specimen used in experiments. 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 a numerical point of view, this method is close to the experimental loading procedure [21] and hence avoids numerical instability. The delamination toughness values obtained through experiments [14] and associated fatigue parameters are given in Table 2.
Fatigue delamination simulations are performed in pure mode I, pure mode II and for mixed-mode (for mixed-mode, mode I and mode II components are equal (/ = 0.5)). For pure mode I, specimen arms are loaded with opposing moments, Fig. 4a. 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 conditions for pure mode II and is shown in Fig. 4b and for mixed-mode (/ = 0.5) is shown in Fig. 4c. For a mode ratio of 50% mode II, the ratio q between the two applied moments is as follows [17]: The energy release rate for pure mode I [17] Where the width of the specimen is b, E is the longitudinal flexural Young's modulus and I is the second moment of area of the specimen's arm. The energy release rate for pure mode II [17] G II ¼ 3 4 In which M = cP/2, see Fig. 4b. The value of c is 30 mm. For mixed mode case [17] Fig . 5 shows the Paris plot behaviour for mode I, mode II and mixedmode (50% mode II). Simulation results are found in good agreement with experimental results [14]. Fig. 5 also shows scattering of experimental data. The parameters (k, b and C) are selected for simulations to give the best trend of scattered data. 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, Eq. (17).
In order to verify the efficacy of Eq. (17) experimental results of Martin and Murri [18] on delamination growth have been selected  Table 1 Material properties for UD HTA/6376C [14].  Table 2 Delamination toughness values for UD HTA/6376C and associated Fatigue parameters.
Test method G C (kJ/m 2 ) [14] Interface for comparison purposes. Authors in Ref. [18] 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 has been taken from Ref. [19] are given in Table 3. The proposed fatigue damage evolution law with load ratio, R, effect is tested for pure mode I (DCB) as also done in [18]. The associated material toughness values as quoted in [18] and associated fatigue parameters are given in Table 4. The geometry of the specimen for DCB test is with following dimensions, L = 140 mm, b = 25.4 mm, a 0 = 50 mm and h = 4.6 mm, [18]. The results of simulations of fatigue delamination growth along with experimental results for mode I are shown in Fig. 6 and are found in good agreement with experi-mental results. It should be noted that the authors in Ref. [18] used British system of units, that's 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 4. Mode I quasi-static tests are classically performed following the ISO 15024 but no standardised test procedures exist for fatigue [21].

Mixed-mode delamination criteria
Experimental results show that the Paris plot behaviour can be expressed by Eq. (1). In Eq. (1), parameters B and m depend on the material and on the fracture mode. Different authors tried to develop a relation between these parameters and mode-mixture /, so that if values of these parameters are known for certain mode-mixtures then values for other 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 GFRP materials, Kenane and Benzeggagh [20] showed that variation of    [14]. Table 3 Material properties for UD AS4/PEEK laminate [19].  Table 4 Delamination toughness values of UD AS4/PEEK and associated fatigue parameters.
these two parameters with respect to / is monotonic. For CFRP materials, Blanco et al. [13] showed a non-monotonic type trend and also developed a relation between parameters (B, m) and mode-mixture /. Tumino and Cappello [8] also used the non-monotonic equation to establish a relation between fatigue damage parameters and mode-mixture /. In this article the same approach is used to test the proposed fatigue damage evolution law for different mode-mixtures. Expressing Eq. (17) in general form In order to establish a non-monotonic relation, fatigue damage parameters values for three different fracture mechanics test should be known in advance. Hence for known values of pure mode I and mode II and for any mode-mixture, one can define K I , K II , K mix , g I , g II , g mix respectively. Then new values of g(/) and K(/) corresponding to any other mode-mixture / can be determined by using following equations [8,13]: where the expressions for A iK and A ig , i = 1, 2, 3 are given below [8].
Here Eqs. (26) and (27)  The critical values for the energy release rate are connected to / through following equation [8]:

Simulations and results of mixed-mode delamination criteria
In the previous section, the values of fatigue parameters (C(/) and b(/)) for pure mode I (/ = 0.), mode II (/ = 1.) and for mode mixture / = 0.5 are already found. The new values of fatigue parameters for any mode-mixture / can be found using Eqs. (26)-(29). In the present work, mode-mixture of / = 0.25 and / = 0.75 are selected for the simulations in addition to those already described in the previous section. The loading scheme used is the same as shown in Fig. 4c. Now the values of q and the energy release rate G for the mode mixture / = 0.25 can be calculated as follows [17]: q ð/¼0:25Þ ¼ À ; G Ið/¼0:25Þ ¼ 9M 2 bEI ; G IIð/¼0: And similarly for the mode mixture / = 0.75 can be expressed as follows [17]: q ð/¼0:75Þ ¼ The trends obtained for B and m as a function of / from experiments by Blanco et al. [13] are found in reasonably good agreement with numerical results, Fig. 7a and b. The predicted behaviours of Paris plots for mixed mode 0.75 and 0.25 obtained through simulations are given in Fig. 8.

Conclusion
In this article a comprehensive elastic fatigue damage model including ''R'' effect, based on damage energy release rate is pre- sented. The effectiveness of the proposed model has been tested by performing the finite element simulations of different fracture mechanics specimens under cyclic loading conditions. Pure and mixed modes finite element simulations are performed on unidirectional HTA/6376C and AS4/PEEK carbon fibres. The linear Paris plot behaviours 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 experimental results. The proposed fatigue damage model is not only able to reproduce the linear Paris plot behaviour for composite laminates but is also capable of predicting the fatigue behaviour for different mode mixtures. Experimental scatter of fatigue data is governed both by manufacturing imperfections and the material variability. This variability should be introduced in the parameters of the damage modelling.   [21].