A simulation method for fatigue-driven delamination in layered structures involving non-negligible fracture process zones and arbitrarily shaped crack fronts

TL;DR

This paper introduces a new 3D cohesive zone-based computational method for simulating fatigue-driven delamination in layered structures, accurately predicting crack growth without fitting parameters and validated against experimental data.

Contribution

The paper presents a novel 3D simulation method that accounts for non-negligible fracture process zones and arbitrarily shaped crack fronts, validated with experimental benchmarks.

Findings

Accurately predicts fatigue delamination in 3D structures.

Does not require fitting parameters, using experimental inputs.

Validated against experimental benchmark cases.

Abstract

Most of the existing methods for fatigue-driven delamination are limited to two-dimensional (2D) applications or their predictive capabilities have not been validated in three-dimensional (3D) problems. This work presents a new cohesive zone-based computational method for simulating fatigue-driven delamination in the analysis of 3D structures without crack migration. The method accurately predicts fatigue propagation of non-nelgigible fracture process zones with arbitrarily shaped delamination fronts. The model does not require any kind of fitting parameter since all the input parameters are obtained experimentally from coupon tests. The evaluation of the energy release rate is done using two new techniques recently developed by the authors (the growth driving direction and the mode-decomposed J-integral) leading to an accurate prediction of delamination propagation under mixed-mode and non-self-similar growing conditions. The new method has been implemented as a UEL for Abaqus and validated against an experimental benchmark case with varying crack growth rate and shape and extension of the fracture process zone.

Figures

Tables

Table 1. Table 1: Dependencies and partial derivatives of the variables in the system using the CZM presented in [ 16 , 17 ] .

Dependencies	Partial derivatives
${\lambda }_o\left({\mu }_o\right)$	$\frac{\partial {\lambda }_o}{\partial {\mu }_o}=\frac{1}{K}$	(14)
${\lambda }_c({\mu }_o,{𝒢}_c)$	$\frac{\partial {\lambda }_c}{\partial {\mu }_o}=-\frac{2{𝒢}_c}{{\mu }_o^2}$	(15)
	$\frac{\partial {\lambda }_c}{\partial {𝒢}_c}=\frac{2}{{\mu }_o}$	(16)
${\mu }_o\left(B\right)$	$\frac{\partial {\mu }_o}{\partial B}=\frac{\eta \left({\tau }_{so}^2-{\tau }_{Io}^2\right)B^{\eta -1}}{2{\mu }_o}$	(17)
${𝒢}_c\left(B\right)$	$\frac{\partial {𝒢}_c}{\partial B}=\eta \left({𝒢}_{sc}-{𝒢}_{Ic}\right)B^{\eta -1}$	(18)
${\omega }_d({\lambda }_o,{\lambda }_c,\lambda )$	$\frac{\partial {\omega }_d}{\partial {\lambda }_o}=\frac{1}{2}K{\lambda }_c\frac{{\lambda }_o^2-2{\lambda }_c{\lambda }_o+{\lambda }_c\lambda }{{\left({\lambda }_o-{\lambda }_c\right)}^2}$	(19)
	$\frac{\partial {\omega }_d}{\partial {\lambda }_c}=\frac{1}{2}K{\lambda }_o^2\frac{{\lambda }_o-\lambda }{{\left({\lambda }_o-{\lambda }_c\right)}^2}$	(20)
	$\frac{\partial {\omega }_d}{\partial \lambda }=\frac{1}{2}K{\lambda }_o{\lambda }_c\frac{1}{\left({\lambda }_c-{\lambda }_o\right)}$	(21)
${𝒟}^e({\omega }_d,{𝒢}_c)$	$\frac{\partial \left(\frac{{\omega }_d}{{𝒢}_c}\right)}{\partial {\omega }_d}=\frac{1}{{𝒢}_c}$	(22)
	$\frac{\partial \left(\frac{{\omega }_d}{{𝒢}_c}\right)}{\partial {𝒢}_c}=\frac{-{\omega }_d}{{𝒢}_c^2}$	(23)

Table 2. Table 2: Laminate elastic properties of the validation material.

Laminate properties
$E_{11}$: Longitudinal Young’s modulus	154	GPa
$E_{22}=E_{33}$: Transversal Young’s modulus	8.5	GPa
$G_{12}=G_{13}$: Shear modulus in the longitudinal planes	4.2	GPa
$G_{23}$: Shear modulus in the transversal plane	3.04	GPa
${\mu }_{12}={\mu }_{13}$: Poisson’s coefficient in the longitudinal planes	0.35	-
${\mu }_{23}$: Poisson’s coefficient in the transversal plane	0.4	-

Table 3. Table 3: Static interlaminar fracture properties of the validation material.

Property	Value	Units
${𝒢}_{Ic}$	0.305	N/mm
${𝒢}_{IIc}$	2.77	N/mm
$\eta$	2.05	-
${\tau }_{Io}$	32.6	MPa
${\tau }_{IIo}$	98	MPa

Table 4. Table 4: Paris’ law-based fitting parameters (Equation ( 30 )) used in 2D simulations.

Property			Units
Mode I $\Phi$=0
	$R=0.1$	$R=0.5$
$p$	8.39	13.8	-
$A$	6.45E-2	8.26E2	mm/cycle
Mode II $\Phi$=1
	$R=0.1$	$R=0.3$
$p$	3.62	4.17	-
$A$	7.03E-1	7.41E-1	mm/cycle
Mixed-mode $R$=0.1
	$\Phi =0.5$	$\Phi =0.75$
$p$	5.95	4.14	-
$A$	1.26	3.96E-1	mm/cycle

Table 5. Table 5: Analysis and model settings of 2D simulations.

Solver	Newton-Raphson
	Linear static
Increment	Adaptive increment size
Element type in beams	C3D8 (Abaqus 6.12)
Elements in beam thickness	4
Integration of cohesive element	2 x 2 Newton-Cotes quadrature
$L_e$: Cohesive element length	0.2 mm
$\Delta a_t$: Crack increment target	0.2 mm
$K$: Penalty stiffness	1E5 N/mm3
ACI [34]:	On

Table 6. Table 6: Analysis and model settings of 3D simulations.

Solver	Newton-Raphson
	Linear static
Increment	Adaptive increment size
Element type in beams	SC8R (Abaqus 6.12)
Elements in beam thickness	1
Elements in reinforcement thickness	1
Integration of cohesive element	2 x 2 Newton-Cotes quadrature
$L_{X_1}$: Cohesive element length	0.36810 mm
$L_{X_2}$: Cohesive element width	0.21605 mm
$\Delta a_t$: Crack increment target	0.36 mm
$K$: Penalty stiffness	1E5 N/mm3
ACI [34]:	On

Table 7. Table 7: Fatigue properties used in 3D simulations.

Property	Value	Units
$p_I$	8.39	-
$A_I$	6.45E-2	mm/cycle
$p_{II}$	3.62	-
$A_{II}$	7.03E-1	mm/cycle
$p_m$	-6.20	-
$A_m$	1.28E5	mm/cycle
${𝒢}_{Ith}$	8.54E-2	N/mm
${𝒢}_{IIth}$	8.29E-2	N/mm
${\eta }_2$	-	-
$R$	0.1	-

Equations

$$\begin{array}{rl}{\tau }_j& =\left(1-{\mathcal{D}}^K\right)K{\delta }_j\text{for}j=1,2\\ {\tau }_3& =\left(1-{\mathcal{D}}^K\right)K{\delta }_3-{\mathcal{D}}^{K}K⟨-{\delta }_3⟩\end{array}$$

$$\lambda =\sqrt{{\left({\delta }_I\right)}^2+{\left({\delta }_s\right)}^2}$$

$${\delta }_I=⟨{\delta }_3⟩,{\delta }_s=\sqrt{{\left({\delta }_1\right)}^2+{\left({\delta }_2\right)}^2}$$

$${\mathcal{G}}_c={\mathcal{G}}_{Ic}+\left({\mathcal{G}}_{sc}-{\mathcal{G}}_{Ic}\right)B^{\eta }$$

$$B=\frac{{\delta }_s^2}{{\delta }_I^2+{\delta }_s^2}$$

$${\mu }_o=\sqrt{{\left({\tau }_{Io}\right)}^2+\left[{\left({\tau }_{so}{)}^2-({\tau }_{Io}\right)}^2\right]B^{\eta }}$$

$${\tau }_{so}={\tau }_{Io}\sqrt{\frac{{\mathcal{G}}_{sc}}{{\mathcal{G}}_{Ic}}}$$

$${\lambda }_o=\frac{{\mu }_o}{K},{\lambda }_c=\frac{2{\mathcal{G}}_c}{{\mu }_o}$$

$${\mathcal{D}}^e=1-\frac{\left(1-{\mathcal{D}}^K\right){\lambda }_{\mathcal{D}}}{{\lambda }_o}$$

$${\lambda }_{\mathcal{D}}=\frac{{\lambda }_o{\lambda }_c}{{\lambda }_c-{\mathcal{D}}^K\left({\lambda }_c-{\lambda }_o\right)}$$

$$\frac{\text{d}{\mathcal{D}}^e}{\text{d}N}=\left(\frac{\partial {\mathcal{D}}^e}{\partial B}\frac{\partial B}{\partial a}+\frac{\partial {\mathcal{D}}^e}{\partial \lambda }\frac{\partial \lambda }{\partial a}\right)\frac{\text{d}a}{\text{d}N}$$

$$\begin{array}{rl}\frac{\partial {\mathcal{D}}^e}{\partial B}& =\frac{\partial \left(\frac{{\omega }_d}{{\mathcal{G}}_c}\right)}{\partial {\omega }_d}\left(\frac{\partial {\omega }_d}{\partial {\lambda }_o}\frac{\partial {\lambda }_o}{\partial {\mu }_o}+\frac{\partial {\omega }_d}{\partial {\lambda }_c}\frac{\partial {\lambda }_c}{\partial {\mu }_o}\right)\frac{\partial {\mu }_o}{\partial B}+\left(\frac{\partial \left(\frac{{\omega }_d}{{\mathcal{G}}_c}\right)}{\partial {\omega }_d}\frac{\partial {\omega }_d}{\partial {\lambda }_c}\frac{\partial {\lambda }_c}{\partial {\mathcal{G}}_c}+\frac{\partial \left(\frac{{\omega }_d}{{\mathcal{G}}_c}\right)}{\partial {\mathcal{G}}_c}\right)\frac{\partial {\mathcal{G}}_c}{\partial B}\\ \frac{\partial {\mathcal{D}}^e}{\partial \lambda }& =\frac{\partial \left(\frac{{\omega }_d}{{\mathcal{G}}_c}\right)}{\partial {\omega }_d}\frac{\partial {\omega }_d}{\partial \lambda }\end{array}$$

$$\begin{array}{rl}\frac{\partial {\mathcal{D}}^e}{\partial B}& =\frac{\eta \left({\mathcal{G}}_{sc}-{\mathcal{G}}_{Ic}\right)B^{\left(\eta -1\right)}\lambda }{K{\lambda }_c{\lambda }_o\left({\lambda }_o-{\lambda }_c\right)}\\ \frac{\partial {\mathcal{D}}^e}{\partial \lambda }& =\frac{1}{{\lambda }_c-{\lambda }_o}\end{array}$$

$$\begin{array}{rl}\frac{\partial B}{\partial a}& \approx \frac{\partial B}{\partial x_1}\\ \frac{\partial \lambda }{\partial a}& \approx \frac{\partial \lambda }{\partial x_1}\end{array}$$

$$\begin{array}{r}\frac{\partial B}{\partial x_1}=\frac{\partial B}{\partial {\delta }_j}\frac{\partial {\delta }_j}{\partial x_1}\\ \frac{\partial \lambda }{\partial x_1}=\frac{\partial \lambda }{\partial {\delta }_j}\frac{\partial {\delta }_j}{\partial x_1}\end{array}$$

$$\begin{array}{rl}\frac{\partial B}{\partial {\delta }_1}& =\frac{2{\delta }_1⟨{\delta }_3{⟩}^2}{{\lambda }^4};\frac{\partial B}{\partial {\delta }_2}=\frac{2{\delta }_2⟨{\delta }_3{⟩}^2}{{\lambda }^4};\frac{\partial B}{\partial {\delta }_3}=-\frac{2{\delta }_s^2⟨{\delta }_3⟩}{{\lambda }^4}\\ \frac{\partial \lambda }{\partial {\delta }_1}& =\frac{{\delta }_1}{\lambda };\frac{\partial \lambda }{\partial {\delta }_2}=\frac{{\delta }_2}{\lambda };\frac{\partial \lambda }{\partial {\delta }_3}=\frac{⟨{\delta }_3⟩}{\lambda }\end{array}$$

$$\frac{\partial {\delta }_j}{\partial x_1}=R_{ji}\left[-\frac{\partial N_{ik}}{\partial {\eta }_{\alpha }},\frac{\partial N_{ik}}{\partial {\eta }_{\alpha }}\right]{\left[R_{1p}\frac{1}{2}\frac{\partial N_{pm}}{\partial {\eta }_{\alpha }}\left(C_m^{-}+C_m^{+}+Q_m^{-}+Q_m^{+}\right)\right]}^{-1}\left[Q_k^{-},Q_k^{+}\right]$$

$$R=\sqrt{\frac{{\mathcal{G}}_{min}}{{\mathcal{G}}_{max}}}$$

$$\Phi =\frac{{\mathcal{G}}_s}{{\mathcal{G}}_I+{\mathcal{G}}_s}$$

$$\frac{\text{d}a}{\text{d}N}=\{\begin{array}{rl}& A{\left(\frac{{\mathcal{G}}_{max}\left(1-R\right)}{{\mathcal{G}}_c}\right)}^p\text{for}{\mathcal{G}}_{th}<{\mathcal{G}}_{max}<{\mathcal{G}}_c\\ & 0\text{for}{\mathcal{G}}_{max}\le {\mathcal{G}}_{th}\end{array}$$

$$\begin{array}{rl}p& ={\Phi }^2\left(p_s-p_I-p_m\right)+\Phi p_m+p_I\\ \text{log}\left(A\right)& ={\Phi }^2\text{log}\left(\frac{A_s}{A_m{A}_I}\right)+\Phi \text{log}\left(A_m\right)+\text{log}\left(A_I\right)\end{array}$$

$${\mathcal{G}}_{th}={\mathcal{G}}_{Ith}+\left({\mathcal{G}}_{sth}-{\mathcal{G}}_{Ith}\right){\Phi }^{{\eta }_2}$$

$${\mathcal{G}}_{max}=J_I+J_{II}+J_{III}$$

$$\Phi =\frac{J_{II}+J_{III}}{J_I+J_{II}+J_{III}}$$

$$\Delta N=\frac{\Delta a_t}{{\left(\frac{\text{d}a}{\text{d}N}\right)}_{max}}$$

$${\tau }_{Io}={\tau }_{so}\sqrt{\frac{{\mathcal{G}}_{Ic}}{{\mathcal{G}}_{sc}}}$$

$$\frac{\text{d}a}{\text{d}N}\approx \frac{\Delta n_{\mathcal{D}}{L}_e}{\Delta N}$$

Full text

††journal: Composites Part A

1 Introduction

The two alternatives to deal with interlaminar fatigue damage in aircraft design are non-growth criterion and the damage tolerance approach [1]. Some of the certification processes of aircraft composite structures rely on extensive, and very costly, fatigue testing aimed at ensuring that no growth occurs during the service life. Testing is performed, from small components to large substructures, to evaluate the number of cycles required to make a crack grow perceptibly as a function of the load. On the other hand, recent methods for assessing the certification of aircraft composite structures [2, 3, 1], allow for a “slow growth” criterion that ensures that any initial delamination inherently present in the structure will grow up to a critical size that can cause the structure collapse during service. If advanced computational tools capable of predicting fatigue damage growth behavior with accuracy are available, an inspection schedule that enables detection of damage before it becomes critical can be reliably programmed. Consequently, a damage tolerance approach can be adopted, enabling lighter and more efficient aircraft designs.

One of the most recent reviews of simulation methods for delamination propagation [4] splits them into two categories: linear elastic fracture mechanics (LEFM) and cohesive zone model (CZM) based methods. Another very recent review [5], adds the categories of stress/strain and extended finite element method (XFEM) based models. However, the stress-based expressions are phenomenological models which can be rewritten in terms of LEFM parameters, and most of the XFEM approaches resort to LEFM to define the load. Therefore, attention will be focused on LEFM and CZM based simulation methods.

In general terms, the LEFM based methods [6, 7] directly apply any Paris’ law-based expression [8] for the crack growth rate, $\textrm{d}a/\textrm{d}N$, evaluated using virtual crack closing technique (VCCT) [9, 10], or any other technique for determining the energy release rate or stress intensity factor. The strength of such methods is that, since they use a phenomenological expression for the crack growth rate, any experimental evidence of the effect of the load ratio, $R$, and/or the mode mixity, $\Phi$, may easily be included in the simulation (provided that the method used in the simulation for characterizing the load is capable of measuring $R$ and $\Phi$). The weaknesses are the assumption of negligible fracture process zone (FPZ), which is not always appropriate, and the need for a pre-existing flaw. An additional limitation related to the VCCT is the requirement of orthogonality of the mesh with the delamination front [11]. In contrast, the CZM approach accounts for a non-negligible FPZ, allows to model crack initiation from a pristine interface surface and avoids the need for re-meshing.

Current cohesive zone models are formulated within the framework of damage mechanics to ensure irreversible crack propagation [12, 13, 14, 15, 16, 17, 18, 19]. In interface damage mechanics, a damage parameter, $\mathcal{D}$, acts to prevent the restoration of the previous cohesive state between the interfacial surfaces. Most of the existing fatigue formulations are extensions of the cohesive laws for quasi-static loading accounting for an additional criterion for damage development due to fatigue loading.

The existing CZMs for fatigue delamination can be divided into two main groups [4]: the loading-unloading hysteresis models [20, 21, 22, 23, 24, 25, 26], which simulate the whole load cycle to compute the evolution of the damage variable, and the envelope load models [27, 28, 29, 30, 31, 32, 33], which only model the maximum cyclic load. One of the main advantages of the loading-unloading hysteresis models is that they are capable to model variable loading spectra. In return, as they simulate cycle by cycle, these models become computationally unfeasible for high-cycle fatigue simulation, being only used in low-cycle fatigue applications. Conversely, in envelope load methods, the number of cycles is discretized and the damage variable, $\mathcal{D}$, is updated for each increment in cycles, $\Delta{N}$. The damage at a given number of cycles, $\mathcal{D}(N_{n}+\Delta N)$ is determined by integration of the damage rate, $\textrm{d}\mathcal{D}/\textrm{d}N$. Thus, the envelope load methods are more efficient for simulating high-cycle fatigue.

In [28], the authors incorporated a link between the damage rate, $\textrm{d}\mathcal{D}/\textrm{d}N$, and the crack growth rate, $\textrm{d}a/\textrm{d}N$, which was later followed in [29, 30, 31, 34, 35, 33]. By using this approach, any influence of the load ratio, $R$, and/or the mode mixity, $\Phi$, is included in the formulation by means of a phenomenological model for the crack growth rate, $\textrm{d}a/\textrm{d}N=f(\mathcal{G}_{max},R,\Phi...)$. A challenge related to these CZMs is the accurate calculation of the envelope energy release rate, $\mathcal{G}_{max}$, which is needed to evaluate the Paris’ law-based expression for $\textrm{d}a/\textrm{d}N$.

The envelope energy release rate, $\mathcal{G}_{max}$, can be computed using the $J$-integral evaluated at the cohesive zone interface [36]; i.e. integrating the quantities from the cohesive zone model all across the length of the cohesive zone. This approach, used in the methods presented in [29, 37, 34], leads to an accurate evaluation of the energy release rate, $\mathcal{G}$. However, since no formulations for the computation of the $J$-integral in 3D modeling of delamination growth, using a cohesive zone model approach, were available at the time they were formulated, these methods were only applicable to 2D simulations. To overcome this limitation, the calculation of the envelope energy release rate, $\mathcal{G}_{max}$, was done by integration of the traction-separation historical response of the most opened point in the cohesive zone, in the methods proposed in [31, 35, 38]. This implicitly assumes that the historical energy dissipation of the point at the crack tip is comparable to the macroscopic energy release rate, thus, assuming self-similar crack growth. Another approach, used in [28] and [30], is replacing the energy release rate, $\mathcal{G}$, by the total specific work, $\omega_{tot}$, in the Paris’ law-based expression for $\textrm{d}a/\textrm{d}N$. The total specific work, $\omega_{tot}$, is the area under the quasi-static law using instantaneous local information at each material point (c.f. Figure 2). This energy measure, which is a field quantity of the cohesive zone, circumvents the need for the more computationally expensive $J$-integral evaluation of the envelope energy release rate, $\mathcal{G}_{max}$ [39].

In addition, some of the previous methods make use of the length of the cohesive zone, $l_{cz}$, in the link between the crack growth rate, $\textrm{d}a/\textrm{d}N$, and the damage rate, $\textrm{d}\mathcal{D}/\textrm{d}N$. The length of the cohesive zone, $l_{cz}$, is a parameter that can be either approximated analytically [28, 29, 37] or related to that extracted from quasi-static simulations [30]. Note that the use of fitting parameters that are not possible to be determined experimentally or need to be adjusted depending on the problem make the predictive character and accuracy of the methods questionable. Others [31, 33], make use of the length associated to an integration point in the crack propagation direction, $l_{e}$, instead. In this case, the propagation direction must be determined in advance. On the other hand, the methods in [34, 38] rely on the computation of the spatial derivatives of certain quantities in the CZM formulation along the propagation direction. However, the lack of efficient formulations for the identification of the crack propagation direction prevented these models of being applicable to 3D analysis of delamination.

In [32], a benchmark study on six of the aforementioned computational methods [34, 28, 27, 40, 30, 31] is performed. The accuracy of the predicted crack growth rate is studied and compared. Also, the sensitivity to different quasi-static material parameters and method-related fitting parameters are analyzed. It is shown that the method in [34] is superior in terms of robustness and accuracy. Moreover, it does not include any fitting parameter. It is worth mentioning that, although some of the most recent methods [35, 38, 33] were not benchmarked in [32], either they are based on any of the previous models or they are not formulated for their applicability to 3D structures.

In this work, a model capable of predicting delamination growth under both quasi-static and fatigue loading in 3D analysis is presented. It is based on the link between the crack growth rate, $\textrm{d}a/\textrm{d}N$, and the damage rate, $\textrm{d}\mathcal{D}/\textrm{d}N$ presented in [34]. An efficient methodology to determine the growth driving direction in 3D analysis [41] is used. In addition, the model uses a 3D CZ $J$-integral formulation [42] for calculating the mode-decomposed energy release rates. This 3D CZ $J$-integral formulation takes into account the current loading state over the entire cohesive zone. The model verification is done by simulating the same tests that provided the phenomenological expression for the crack growth rate, $\textrm{d}a/\textrm{d}N$, as input under different mode mixities, $\Phi$, and load ratios, $R$. Finally, the method is validated by comparison of the results obtained from an experimental benchmark test on a partially reinforced double cantilever beam (DCB) specimen with varying crack growth rate and front shape [43, 44].

2 Cohesive zone model for quasi-static and fatigue-driven delamination in 3D structures

The point of departure for the formulation developed in this work is the simulation method for high-cycle fatigue-driven delamination presented in [34], which was formulated for its applicability to 2D structures. The method in [34] is, in turn, an extension of the CZM for quasi-static loading proposed in [16, 17].

The quasi-static cohesive model is formulated in the framework of damage mechanics. A thermodynamical representation of the irreversible process of interlaminar fracture is done by defining an energy-based damage variable, $\mathcal{D}^{e}$. It corresponds to the fraction of dissipated specific energy, $\omega_{d}$, to the energy necessary to create a unit of new surface under quasi-static loading conditions, which is the fracture toughness, $\mathcal{G}_{c}$ (c.f. Figure 2). The fatigue model is based on an envelope load approach, in which the energy-based damage rate, $\textrm{d}\mathcal{D}^{e}/\textrm{d}N$, is directly related to the crack growth rate, $\textrm{d}a/\textrm{d}N$, avoiding making use of any fitting parameter. This is achieved by maintaining the quasi-static relationship between tractions, $\boldsymbol{\tau}$, and displacement jumps, $\boldsymbol{\delta}$, during fatigue propagation [34].

2.1 Quasi-static damage model

The CZM in [16, 17] describes the cohesive behavior of a band of micro-cracked material ahead of the crack tip capable of transferring stress until complete separation. The distribution of cohesive tractions, $\boldsymbol{\tau}(\boldsymbol{\delta})$, is a function of the separations between crack faces, $\boldsymbol{\delta}$, referred to as displacement jumps. Defining a local Cartesian coordinate system on the delamination mid-surface such that $\boldsymbol{e_{1}}$ and $\boldsymbol{e_{2}}$ are the local tangential directions and $\boldsymbol{e_{3}}$ is the local normal direction (c.f. Figure 1), the constitutive relation is defined as:

[TABLE]

where $\mathcal{D}^{K}\in[0,1]$ is a scalar damage variable reducing the penalty stiffness, $K$, introduced in the numerical implementation of the CZM in a finite element framework.

The evolution of the stiffness degrading damage variable, $\mathcal{D}^{K}$, is governed by an equivalent one-dimensional cohesive law (c.f. Figure 2) to ensure irreversibility under changing mixed-mode loading conditions. The equivalent one-dimensional displacement jump is defined as:

[TABLE]

where $\delta_{I}$ is the displacement jump associated to mode I and $\delta_{s}$ is the shear sliding resultant of the displacement jumps associated to any combination of mode II and mode III.

[TABLE]

The process of fracture initiates when the stress reaches the interlaminar strength, $\mu_{o}$. Then, the cohesive law describes a stress-softening behavior. With increasing displacement jump, $\lambda$, cohesive traction, $\mu$, decreases. The amount of energy dissipated per unit of newly created crack area in the opening process (from 0 to the critical displacement jump, $\lambda_{c}$) is the fracture toughness, $\mathcal{G}_{c}$. The fracture toughness, $\mathcal{G}_{c}$, for a given mode-mixity, is determined by [45]:

[TABLE]

where subscripts $Ic$ and $sc$ denote the pure mode I and shear quasi-static critical values, respectively, and $\eta$ is a mode interaction parameter determined by fitting Equation (4) to a batch of experimentally measured fracture toughnesses under different mode-mixities. Note that Equation (4) was originally formulated in [45] for the mode I-II interpolation of the fracture toughness and, thus, $\mathcal{G}_{IIc}$ was used instead of $\mathcal{G}_{sc}$. In the applied CZM [16, 17], modes II and III are not disregarded and the shear properties are treated as mode II ($\mathcal{G}_{sc}=\mathcal{G}_{IIc}=\mathcal{G}_{IIIc}$). $B$ is the local displacement-based mode-mixity and is defined in terms of the displacement jump as:

[TABLE]

The mixed-mode intepolation of the interlaminar strength, $\mu_{o}$, is done by:

[TABLE]

where $\tau_{Io}$ and $\tau_{so}$ are the tensile and shear interfacial strengths, where $\tau_{so}=\tau_{IIo}=\tau_{IIIo}$. A mode-independent penalty stiffness, $K$, is used which sets the interfacial shear strength dependent on the other material properties [17]:

[TABLE]

The onset, $\lambda_{o}$, and propagation, $\lambda_{c}$, of delamination in terms of the equivalent one-dimensional displacement jump are expressed as:

[TABLE]

Finally, the stiffness degrading, $\mathcal{D}^{K}$, and the energy-based, $\mathcal{D}^{e}$, damage variables are related by:

[TABLE]

where $\lambda_{\mathcal{D}}$ is the equivalent one-dimensional displacement jump associated to the current damage state, defined as:

[TABLE]

2.2 Fatigue damage rate model

Since $\mathcal{D}^{e}$ depends on the local displacement-based mode-mixity, $B$, and the equivalent one-dimensional displacement jump, $\lambda$, the rate of energy-based damage, $\frac{\textrm{d}\mathcal{D}^{e}}{\textrm{d}N}$, is deduced by applying the chain rule [34] as:

[TABLE]

Each of the partial derivatives in Equation (11) are addressed in the following. The derivative of the energy-based damage with respect to the mode mixity, $\frac{\partial\mathcal{D}^{e}}{\partial B}$, and the derivative of the energy-based damage with respect to the equivalent one-dimensional displacement jump, $\frac{\partial\mathcal{D}^{e}}{\partial\lambda}$, are derived once again applying the chain rule:

[TABLE]

The partial derivatives in Equation (12) are solved for the particular application of the CZM from [16, 17] and listed in Table 1. The expressions for $\frac{\partial\mathcal{D}^{e}}{\partial B}$ and $\frac{\partial\mathcal{D}^{e}}{\partial\lambda}$, obtained after substituting the partial derivatives in Equation (12) by the expressions listed in Table 1 read:

[TABLE]

In a self-similar crack growth, $\frac{\partial B}{\partial a}$ and $\frac{\partial\lambda}{\partial a}$ can be interpreted as the slopes of the local mode-mixity and the equivalent one-dimensional displacement jump along the direction of crack propagation (for further details on this assumption, the reader is referred to the original paper [34] describing the 2D method). Thus, if $x_{1}$ is defined as the crack growth direction coordinate, the rates $\frac{\partial B}{\partial a}$ and $\frac{\partial\lambda}{\partial a}$ may be approximated by:

[TABLE]

Then, by application of the chain rule, $\frac{\partial B}{\partial x_{1}}$ and $\frac{\partial\lambda}{\partial x_{1}}$ read:

[TABLE]

where the terms $\frac{\partial B}{\partial\delta_{j}}$ and $\frac{\partial\lambda}{\partial\delta_{j}}$ are solved as:

[TABLE]

Note that, in order to compute $\frac{\partial\delta_{j}}{\partial x_{1}}$ of Equation (25), the local Cartesian coordinates $(e_{1},e_{2},e_{3})$ must be rotated around $e_{3}$ (c.f. Figure 1) so that a new local Cartesian coordinate system $(x_{1},x_{2},x_{3})$ is defined, where $x_{1}$ is the direction of crack propagation and tangent to the mid-surface, $x_{3}$ is the normal direction to the mid-surface and $x_{2}$ is perpendicular to $x_{1}$ and $x_{3}$. Therefore, the propagation direction must be identified. In the framework of LEFM, the propagation direction is usually understood as the direction normal to the crack front. Conversely, if a CZM approach is used, the crack front is not defined by a line, but there is a band of damaged material ahead of the crack tip, and thus the normal direction to a front line does not apply. In this case, the concept of growth driving direction (GDD) [41] can be used. The GDD can be defined as the gradient vector field of the total specific work, $\omega_{tot}$, over the fracture toughness, $\mathcal{G}_{c}$, with respect to the coordinates tangent to the mid-surface (c.f. Figure 3). It can be efficiently evaluated at any point within the cohesive zone using information available at finite element level. Moreover, the formulation presented in [41] also addresses the calculation of the slopes $\frac{\partial\delta_{j}}{\partial x_{1}}$, although the main steps are introduced in the following for completeness. By assuming that the curvature of the interface within the element domain is small, the derivatives $\frac{\partial\delta_{j}}{\partial x_{1}}$ read:

[TABLE]

where $R_{ji}$ is the transformation tensor that relates the global coordinate $i$ to the local crack coordinate $j$. $\frac{\partial N_{ik}}{\partial\eta_{\alpha}}$ is the derivative of the shape function matrix with respect to the $\alpha$-isoparametric coordinate of the finite element, in which subscript $k$ runs from 1 to the number of degrees of freedom of each of the surfaces of the cohesive element (lower and upper). $C_{m}^{-}$ and $C_{m}^{+}$ are the global coordinates of the nodes at the lower and upper surfaces, respectively. $Q_{m}^{-}$ and $Q_{m}^{+}$ are the nodal displacements of the lower and upper surfaces and subscript $m$ runs from 1 to the number of degrees of freedom of each of the surfaces of the cohesive element.

The last term in Equation (11), $\textrm{d}a/\textrm{d}N$, is the crack growth rate, which can be expressed by any variant of the Paris’ law [8] based on the energy release rate, $\mathcal{G}$. Thus, any experimentally observed effect on the crack growth rate, such as the effect of the mode mixity, $\Phi$, or the load ratio, $R$, can be included in the fatigue damage model. In this work, the crack growth rate, $\textrm{d}a$/$\textrm{d}N=f(\mathcal{G}_{max},R,\Phi)$, is assumed to depend on the maximum cyclic energy release rate, $\mathcal{G}_{max}$, the load ratio, $R$, defined as:

[TABLE]

and the energy-based mode mixity, $\Phi$:

[TABLE]

where subscript $s$ denotes the shear component of the energy release rate.

The following Paris’ law-based expression is proposed for the crack growth rate, $\textrm{d}a$/$\textrm{d}N$, though any other expression that incorporates, for instance, the energy release threshold [46, 47, 48, 49] could be used. For a given load ratio, $R$, the Paris’ law-like expression reads:

[TABLE]

where the exponent, $p$, and coefficient, $A$, are mode-dependent parameters determined by [50]:

[TABLE]

$p_{I}$ and $A_{I}$ are the parameters for pure mode I, $p_{s}$ and $A_{s}$ are the parameters for shear mode, and $p_{m}$ and $A_{m}$ are mode interpolation parameters. Note that, in the original formulation presented in [50], the pure mode II parameters, $p_{II}$ and $A_{II}$, where used instead of $p_{s}$ and $A_{s}$. Likewise in the quasi-static model [16, 17], modes II and III are not disaggregated in the current implementation ($p_{s}=p_{II}=p_{III}$ and $A_{s}=A_{II}=A_{III}$). $\mathcal{G}_{th}$ in Equation (30) is the energy release rate threshold below which no propagation occurs. Its dependence with the mode mixity is assumed to follow a Benzeggagh-Kenane based [45, 28] expression:

[TABLE]

where subscripts $Ith$ and $sth$ denote the pure mode I and shear threshold values, respectively, and $\eta 2$ is an experimentally determined mode interaction parameter.

The maximum cyclic energy release rate, $\mathcal{G}_{max}$, and the energy-based mode mixity, $\Phi$, are computed using the 3D CZ $J$-integral formulation presented in [42]. According to [42], the integration paths cross the CZ following the GDD trajectory (c.f. Figure 3). The computation of the $J$-integral is done using the current traction-displacement jump field over the entire integration path. Thus, any variation in the cohesive quantities due to an abrupt change in the loading scenario or geometry is captured at the current time. Each point within the cohesive zone belongs to a single integration path; and each integration path is understood as a crack propagating in the GDD at a velocity of $\textrm{d}a/\textrm{d}N$. In other words, all the integration points in the cohesive zone along a certain GDD are attributed the same $\textrm{d}a/\textrm{d}N$, which is a scalar measure described by a Paris’ law-like expression (Equation 30). The model enforces that the speed at which the cohesive zone advances with the number of cycles is equal to $\textrm{d}a/\textrm{d}N$.

Using the mode-decomposed 3D CZ $J$-integral formulation, the maximum cyclic energy release rate, $\mathcal{G}_{max}$, is computed as the sum of the $J$-terms:

[TABLE]

and the energy-based mode mixity, $\Phi$, defined in Equation (29) is evaluated as:

[TABLE]

Note the difference between the displacement jump-based mode mixity, $B$, expressed in Equation (5) as a local quantity, and the energy release rate-based mode mixity, $\Phi$, which is extracted from the mode-decomposition of the $J$-integral.

Finally, for the cycle jump strategy, a target for the maximum increment in crack length per solution sub-step, $\Delta a_{t}$, is set. Based on this, the increment in cycles is determined as:

[TABLE]

where $\left(\frac{\text{d}a}{\text{d}N}\right)_{max}$ is the instantaneous maximum crack growth rate in the model.

The formulation described in this section has been implemented as a user-defined cohesive element in the commercial finite element code Abaqus. The implemented procedure for the calculations of the simulation method for fatigue-driven delamination is summarized in the flowchart shown in Figure 4.

3 Comparison with experimental results

In order to evaluate the performance of the presented method, two different studies are performed. First, the capability of accurately reproduce the phenomenological expression used as input in the model for the crack growth rate is assessed in Section 3.1: Verification. This is done through 2D simulations of the same specimens that provided the Paris’ law-based expression under different mode mixities, $\Phi$, and load ratios, $R$. Secondly, the capability of the method to accurately predict delamination growth in problems that cannot be simplified to 2D models is assessed in Section 3.2: Validation. The benchmark case presented in [43, 44], which exhibits change in crack front shape and growth rate, is used for this purpose.

All the specimens used in these studies were made of unidirectional (UD) Carbon Fibre Reinforced Polymer (CFRP) prepreg plies of 0.184 mm nominal thickness stacked at 0∘. AERNNOVA Engineering, the company supplying the material, provided the elastic properties for the unidirectional laminate material (Table 2). The static interlaminar fracture properties are listed in Table 3. The interlaminar fracture toughness were measured following the standards recommendations [51, 52, 53] under 4 different mode mixity conditions: 0% ($\mathcal{G}_{Ic}$), 50%, 75% and 100% ($\mathcal{G}_{IIc}$). The nominal values are presented in [43]. The Benzeggagh-Kenane [45] mode interpolation parameter, $\eta$, in Equation (4) has been fitted to the results and included in Table 3. Moreover, the interfacial shear strength, $\tau_{IIo}$, was measured following the ASTMD2344 standard. The obtained nominal value was 98 MPa. Then, the interfacial tensile strength, $\tau_{Io}$, is expressed as a function of $\tau_{so}$, $\mathcal{G}_{Ic}$ and $\mathcal{G}_{sc}$ (c.f. Equation (7)) [17]:

[TABLE]

3.1 Verification

Crack growth rate curves under mode I loading conditions were obtained using the multi-specimen testing methodology described in [54] with $R=0.1$ and $R=0.5$. Under mode II and mixed-mode loading conditions, the methodology described in [55] and [56] was used instead. This consisted of varying the applied cyclic displacement in order to increase the energy release rate range covered by a single test. Two different load ratios ($R=0.1$ and $R=0.3$) were tested under mode II conditions and two different mode mixities ($\Phi=0.5$ and $\Phi=0.75$) were tested with $R=0.1$. The Paris’ law-based fitting parameters (according to Equation (30)) from all the mentioned crack growth rate curves are listed in Table 4.

To validate the capability of the model to reproduce the phenomenological expressions used as inputs, these tests have been simulated with 2D FE models assuming plane strain conditions. The specimens were made of 16 layers of unidirectional CFRP and were 25 mm wide. The details of the FE models and analyses are given in Table 5. Note that due to fine mesh discretization, a standard 2x2 Newton-Cotes quadrature rule has been used to numerically integrate element quantities as opposed to more advanced integration schemes [57, 18], which are needed for coarser meshes. The results from the simulations are shown in figures 5 (mode I), 6 (mode II) and 7 (mixed-mode) together with the experimental data and the fitted Paris’ law curves. The crack growth rate in the simulations is determined at every converged increment as:

[TABLE]

where $\Delta n_{\mathcal{D}}$ is the increment in number of damaged elements, $L_{e}$ is the element length and $\Delta N$ is the increment in number of cycles. The results from the simulations accurately reproduce the crack growth rate obtained from experimental tests for all the mode mixities, $\Phi$, and load ratios, $R$, analysed.

Note that the threshold region is not explored in the experimental tests. Therefore, the energy release rate threshold is set to $\mathcal{G}_{th}=0$ in the simulations without influencing the results in figures 5 (mode I), 6 (mode II) and 7 (mixed-mode).

3.2 Validation

The capabilities of the simulation method for predicting propagation of fatigue-driven delamination are validated against the experimental benchmark case presented in [43, 44]. The test specimen is a double cantilever beam (DCB) built by stacking 16 unidirectional CFRP plies at 0∘ with an initial straight delamination front made by a Teflon insert at the midplane. Two reinforcements, made of 8 plies of the same material and orientation, are bonded on top and bottom sides of the specimen in order to promote curved delamination front. The benchmark case generates a rich phenomenology of crack advance in terms of varying crack growth rate and crack front shape.

The test consisted on four loading steps (c.f. Figure 8):

1. Quasi-static step from the initial unloaded position until 5 mm of the prescribed displacement at a loading rate of 1 mm/min.

2. Fatigue step with a maximum cyclic displacement of 5 mm and $R=0.1$ during 420,000 cycles.

3. Quasi-static step until 10 mm of the prescribed displacement at a loading rate of 1mm/min.

4. Fatigue step with a maximum cyclic displacement of 10 mm and $R=0.1$ during 10,000 cycles.

The experimental tests were stopped at the programmed stops indicated with black arrows in Figure 8 to perform X-ray inspection of the delamination front location. A batch of three specimens (numbered 1, 2 and 3) were tested following the described loading sequence. The reader is referred to [43] for further details on the employed test methodology and data reduction techniques.

Exploiting $X_{2}$-symmetry, only one half of the specimen was modeled to reduce the required computational resources. The dimensions and boundary conditions of the partially reinforced DCB specimen model are shown in Figure 9. The simulation and model settings of the simulations are listed in Table 6. Finally, the fatigue properties used as input parameters for the user-defined cohesive elements were obtained from [43] and listed in Table 7. The $p_{m}$ and $A_{m}$ were obtained by fitting Equation (31) to the $p$ and $A$ parameters for pure mode I, mode II, 50% and 75% mode-mixities with $R=0.1$ listed in Table 5. The energy release rate threshold, $\mathcal{G}_{th}$, was not mode-interpolated using Equation (32) because there was no mode-mixity data available that could be used to determine the mode interaction parameter, $\eta_{2}$. Instead, a linear interpolation was used. This is not affecting the results since the loading conditions are very lose to pure mode I loading. The pure mode values, $\mathcal{G}_{Ith}$ and $\mathcal{G}_{IIth}$, were obtained from [43].

The force-displacement relationship of the 3 specimens [43, 44] is plotted in Figure 10 together with the simulation results. Note that during fatigue steps 2 and 4, the maximum cyclic displacement is constant at 5mm and 10 mm, respectively. The simulation method accurately predicts the mechanical response observed in the experiments.

In figures 11 and 12, the position of the delamination fronts of specimens 1, 2 and 3 is plotted for every stop done during the test. Due to the presence of micro-cracks, the fracture zone might contain some contrast liquid. Thus, it is hard to ensure whether the extracted delamination front location from the X-ray pictures corresponds to the beginning or the end of the fracture process zone. For this reason, both the 1 and 0-damage isolines delimiting the cohesive zone are plotted in figures 11 and 12 for the comparison with results from experiments.

It can be observed that the delamination front shape and growth rate with the number of cycles is reproduced with high accuracy by the proposed simulation tool. In figures 11.a-11.h, the delamination front shape changes from initial straight to curved front shape as the delamination approaches the stiffened region. Under fatigue loading, the front is arrested at the reinforcement edges due to the supplied stiffness. It is observed in Figure 11.h that the delamination front of specimens 1, 2 or 3 did not surpass the reinforcement edges within 410,000 cycles. Delamination arrestment is also captured by the simulation. Then, during the quasi-static increment of displacement in Step 3, the delamination front surpasses the reinforcement edge and keeps advancing under the reinforcement during fatigue loading in Step 4. During fatigue Step 4, the crack front shape changes once again.

The historical evolution of the growth driving direction (GDD) and the mode I component of the $J$-integral are plotted in figures 13 and 14. The mode II and III components of the $J$-integral are not plotted since they are negligible under these loading conditions. The results are evaluated within the cohesive zone and projected on the mid-surface. It can be visually inspected that the values for the GDD correspond to the angle of the normal to the delamination front with the global $X_{1}$-axis. More accurate demonstrations of the capabilities of the GDD formulation are provided in [41]. In Figure 13.a, it can be observed that the computed $J$-integral during quasi-static propagation (Step 1) is equal to the mode I fracture toughness, $\mathcal{G}_{Ic}$. Then, during Step 2 (fatigue loading with constant maximum cyclic displacement), the crack front is advancing with a $J$-value lower than the fracture toughness, $\mathcal{G}_{Ic}$, and $J$ is reducing as the crack front advances (c.f. figures 13.b and 13.c). The same is repeated for steps 3 and 4: in Figure 14.d, the $J$-integral during quasi-static propagation again equals the mode I fracture toughness, $\mathcal{G}_{Ic}$, whereas, in figures 14.e and 14.f, $J$ decreases with crack propagation. The accuracy of the 3D CZ $J$-integral in the computation of the mode-decomposed energy release rate is further demonstrated in [42].

4 Summary and conclusions

In this work, a new simulation method for high-cycle fatigue-driven delamination has been proposed. The method is applicable to 3D analysis of delamination in layered composite structures. A cohesive zone model (CZM) approach is used. The model is capable to account for a non-negligible fracture process zone and arbitrarily-shaped delamination fronts avoiding the need of using re-meshing techniques. Moreover, the model is based on a Paris’ law-based expression, as most of the current phenomenological models for representing crack growth rate observations.

Recent developments by the authors, like the formulation for the efficient identification of the growth driving direction (GDD) and the computation of the mode-decomposed $J$-integral in 3D analysis (3D CZ $J$-integral), have been integrated to this numerical tool. First, the GDD formulation is used to compute the spatial derivatives of some quantities of the cohesive zone in the direction of crack advance. This enhances the accuracy and robustness of the method as it avoids the introduction of any fitting parameter. Secondly, the 3D CZ $J$-integral is used to evaluate the mode-decomposed energy release rate that feeds the expression for the crack growth rate. By applying this formulation, a macroscopic measure of the energy release rate is obtained, which takes into account the current traction-displacement jump field of the entire cohesive zone. To the authors knowledge, this has not been previously addressed in 3D analysis using a CZM approach.

The method has been implemented in an FE framework as a user-defined cohesive element in Abaqus. The capability to reproduce the Paris’ law curve under different loading conditions of mode mixity and load ratio has been demonstrated. Moreover, in order to assess its predictive capabilities in specimens that cannot be simplified to 2D representations, it has been evaluated against a benchmark case with varying crack growth rate and front shape. In all cases, the documented discrepancy between experimental and simulation data is very small.

The method proposed is a key step in providing a simulation tool that contributes to reduce the amount of experimental testing needed in the dimensioning of real layered components and structures undergoing fatigue delamination. Further work on this method is ongoing and embraces enhancing the model capabilities to include fiber bridging effects and R-curve behavior on d$a$/d$N$ [59, 60]. The formulation presented is compatible for its use with other cohesive zone law shapes that are better suited to model delamination with large scale bridging and that can be obtained from experiments [19, 61]. In addition, due to the fact that the method is based on the concept of the GDD, the formulation may, in future works, be enhanced to include in-plane anisotropy with respect to fracture properties, i.e. $\mathcal{G}_{c}$ and d$a$/d$N$, which may depend on the off-axis angle between the lamina orientation and the crack growth direction [62].

5 Acknowledgements

This work has been partially funded by the European Union by the financial supports of ERANet AirTN 01/2013 under the project entitled “Methodology to design composite structures resistant to intra- and interlaminar damage (static & fatigue)- MERINDA” and by the Spanish Government (Ministerio de Economia y Competitividad) under contract TRA2015-71491-R, cofinanced by the European Social Fund. In addition, this project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 763990.

6 Data availability

The processed data required to reproduce these findings are available to download from

http://dx.doi.org/10.17632/v7bgzx9gmw.1 and http://dx.doi.org/10.17632/87r49xbrp3.3, open-source online data repositories hosted at Mendeley Data [44, 58].