Improvement of a Cohesive Zone Model for Fatigue Delamination Rate Simulation

Abstract

The cohesive zone model (CZM) has found wide acceptance as a tool for the simulation of delamination in composites and debonding in bonded joints and various implementations of the cohesive zone model dedicated to fatigue problems have been proposed in the past decade. In previous works, the authors have developed a model based on cohesive zone to simulate the propagation of fatigue defects where damage acts on cohesive stiffness, with an initial (undamaged) stiffness representative of that of the entire thickness of an adhesive layer. In the case of a stiffness that is order of magnitude higher than the previous one (for instance, in the simulation of the ply-to-ply interface in composites), the model prediction becomes inaccurate. In this work, a new formulation of the model that overcomes this limitation is developed. Finite element simulations have been conducted on a mode I, constant bending (constant G)-loaded double cantilever beam (DCB) joint to assess the response of the new model with respect to the original one for varying initial stiffness K₀ and cohesive strength σ₀. The results showed that the modified model is robust with respect to changes of two orders of magnitude in initial stiffness and of a factor of two in σ₀.

1. Introduction

Composite and hybrid metal/composite structures are nowadays present not only in the aerospace industry, but thanks to continuous performance improvement and cost reduction, also many more industrial fields are approaching the use of multimaterial structural elements. This requires, in turn, extensive use of adhesive bonding and a more sophisticated capability to simulate and predict the strength of bonded connections where, for this purpose, analytical methods are being progressively integrated or replaced by finite element analysis (FEA). In engineering applications, it is well established that fatigue is the root cause of many structural failures. In the case of bonded joints, fatigue life is related to the initiation and propagation of defects starting at the free edges of joining regions or other features, such as through-thickness holes. In the case of composite or metal/composite joints, fatigue can start also from defects at the same locations cited above, with the difference that the crack may either run into the adhesive or become a delamination crack. In particular in the case of damage tolerant or fail safe design, it is necessary to know how cracks, or defects in general, propagate during the service life of a component.

Originally proposed by Barenblatt [1], the cohesive zone model is extensively used for the prediction of fracture propagation under quasi-static conditions along the interfaces. Considering a bi-linear cohesive law (see for example Figure 1) in a quasi-static simulation, the interface behaves linearly up to δ₀. Above this value, the stiffness of the interface is progressively reduced in order to yield the descending branch of the law until the critical opening δ_c is attained. In this phase, the unloading follows the dashed line. The cohesive zone model was further developed in a discrete damage zone model and integrated with the extended finite element method (XFEM) by Wang et al. [2] in order to simulate the crack propagation in an arbitrary direction.

In the literature, a certain number of works can be found dealing with the simulation of the fatigue crack propagation using the cohesive zone model, many of them being recently reviewed in [3,4]. When dealing with fatigue, crack growth may occur also sub-critically, that is at a value of stress lower than that shown by the solid line in Figure 1. Among cyclic cohesive models, some were applied to the modelling and prediction of fatigue in bulk, ductile materials [5,6,7,8,9,10,11,12,13,14,15,16] or quasi-brittle polymers [17] while, at the same time, several others were developed with reference to applications involving failure at the interfaces, such as debonding in adhesive joints or delamination and matrix cracking (followed by delamination) in polymer composites [18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37].

When modelling fatigue, an important point is the integration of the damage (hence crack growth) rate. Most of the models for bulk, ductile or quasi-brittle materials [6,7,8,9,10,11,12,13,14,15,16,17] as well as some for interfaces [18] rely on a full incremental solution of each individual load cycle, eventually adopting some kind of extrapolation scheme in order to reduce the computation time (see for example [8]). Others adopt instead a load envelope strategy (simulation is performed applying the maximum load of the cycle without any unloading-reloading) [19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37] to save computation time. In these types of model, damage evolution equations are formulated in terms of damage rate per load cycle and they generally include external load parameters (like load ratio) in the damage evolution equation. They are, therefore, defined as non-constitutive [38] to differentiate them from models ([5,6,7,8,9,10,11,12,13,14,15,16,17] and others reported in [38]) where a constitutive traction-displacement behavior was formulated with evolution equations for the internal variables accounting for the possibility of cyclic damage accumulation. The model developed by the authors in [27,31,36] belongs to the non-constitutive kind, and therefore the analysis of the literature will be restricted in the following to the other models of this kind.

The cohesive zone damage evolution under fatigue loading is accounted for in different ways. One first way is to establish a priori a phenomenological law [19,20,22,23,24] where fatigue damage, D_f, sums up with the quasi-static damage, D_s, defined as,

$${\mathrm{D}}_{\mathrm{s}}=\frac{{\mathsf{\delta }}_{\mathrm{c}}}{\mathsf{\delta }}\frac{\mathsf{\delta }-{\mathsf{\delta }}_0}{{\mathsf{\delta }}_{\mathrm{c}}-{\mathsf{\delta }}_0}$$

(1)

to yield the overall effect of fatigue loading to the cohesive zone. In this case, the fatigue damage evolution law contains parameters that have to be adjusted to calibrate the numerical model with experimental results, usually by trial and error, with possible limitations on the simulation of different conditions.

A different approach to fatigue damage was formulated in [21], where a link between damage and fracture mechanics was established assuming that the crack growth rate dA/dN is equal to the sum of the damaged area growth rates of i-th elements in the cohesive zone, i.e.,

$$\frac{\mathrm{dA}}{\mathrm{dN}}={\displaystyle \sum }_{{\mathrm{A}}_{\mathrm{CZ}}}\frac{{\mathrm{dA}}_{\mathrm{d}}{}^{\mathrm{i}}}{\mathrm{dN}}$$

(2)

coming to the following relationship for cyclic damage:

$$\frac{\mathrm{dD}}{\mathrm{dN}}=\frac{1}{{\mathrm{A}}_{\mathrm{CZ}}}\frac{{\left[\left(1-\mathrm{D}\right){\mathsf{\delta }}_{\mathrm{c}}+{\mathrm{D}\mathsf{\delta }}_0\right]}^2}{{\mathsf{\delta }}_0{\mathsf{\delta }}_{\mathrm{c}}}\frac{\mathrm{dA}}{\mathrm{dN}}$$

(3)

while the quasi-static component of damage is still given by Equation (1) and the overall damage rate is of course the sum of the quasi-static and fatigue rates. In this way, the fatigue damage rate can be determined directly from the experimental fatigue delamination rate dA/dN. The last, in fact, can be expressed as a function of the cyclic strain energy release rate, ∆G, by means of a Paris-like equation

$$\frac{\mathrm{dA}}{\mathrm{dN}}=B\Delta G^d$$

(4)

where ΔG = (1 − R²)G_max and G_max is maximum value of G of the cycle. However, in this model the damage variable D representing the loss of stiffness, i.e., $\mathrm{K}=\left(1-\mathrm{D}\right){\mathrm{K}}_0$, does not coincide with the density of microcracks on a representative interface element, $\mathrm{D}={\mathrm{A}}_{\mathrm{d}}/{\mathrm{A}}_{\mathrm{e}}$, as it may be deduced from the application of the effective stress concept and strain equivalence principle of continuum damage mechanics (see [39]). Rather, the density of microcracks on a representative interface element is related to the damage variable representing the loss of cohesive strength $\mathsf{\sigma }=\left(1-\overline{\mathrm{D}}\right){\mathsf{\sigma }}_0$ (linear softening). This implies also that $\overline{\mathrm{D}}$ is coincident with the ratio of the energy dissipated during the damage process, Ξ, to the critical energy release rate, G_c.

$$\overline{\mathrm{D}}=\frac{{\mathrm{A}}_{\mathrm{d}}}{{\mathrm{A}}_{\mathrm{e}}}=\frac{\Xi }{{\mathrm{G}}_{\mathrm{c}}}=1-\left(1-\mathrm{D}\right)\frac{\mathsf{\delta }}{{\mathsf{\delta }}_0}$$

(5)

The same kind of approach was adopted in [25,26,28,29,30,32], although coming to a slightly different formulation of the fatigue damage rate:

$$\frac{\partial {\mathrm{D}}_{\mathrm{f}}}{\partial \mathrm{N}}=\frac{1-{\mathrm{D}}_{\mathrm{s}}-{\mathrm{D}}_{\mathrm{f},\mathrm{u}}}{{\mathrm{A}}_{\mathrm{fat}}}\frac{\mathrm{dA}}{\mathrm{dN}}$$

(6)

in 25, 26, 28 and 29, where D_f,u accounts for unwanted fatigue damage, i.e., fatigue damage that occur before the value of δ in Figure 1, and,

$$\frac{\partial {\mathrm{D}}_{\mathrm{f}}}{\partial \mathrm{N}}=\frac{1-{\mathrm{D}}_{\mathrm{s}}}{{\mathrm{A}}_{\mathrm{eff}}}\frac{\mathrm{dA}}{\mathrm{dN}}$$

(7)

in [32], while [30] exploited Equation (3). In these works, an effort was also put on the development of alternative ways to identify the values of A_CZ (A_eff) and ∆G with respect to [21], where A_CZ was estimated using mode I Rice’s closed-form equation and simulations were performed for geometries where ∆G was independent of crack length. In particular, [32] defined A_eff as the portion of the area on the crack tip element subjected to fatigue damage, while [30] adopted a modified form of Rice’s solution that accounts for the difference in A_cz between mode I and mode II loading. In both works, the simulation was done on geometries where G was dependent on the crack length.

In [25,26,28,29], A_CZ was evaluated numerically by performing a quasi-static analysis before going on to the fatigue damage step. The quasi-static cohesive zone length, A_CZ,s was extracted as the sum of the areas of elements where 0 < D < 1 and, since fatigue crack growth (FCG) occurs at a value G < Gc, A_CZ was finally estimated as:

$${\mathrm{A}}_{\mathrm{CZ}}=\left(\frac{\mathrm{G}}{{\mathrm{G}}_{\mathrm{c}}}\right){\mathrm{A}}_{\mathrm{CZ},\mathrm{s}}$$

(8)

Concerning the calculation of the fatigue damage rate, A_fat in Equation (6) is the portion of A_CZ where damage develops subcritically, that is at a stress lower than the one dictated by the quasi-static traction-displacement behavior. The value of A_fat was taken equal to half of A_CZ based on the results of simulation examples. The value of A_eff in Equation (7) corresponds instead to the undamaged length of the crack tip element when a new loading is applied to the model. For example, if D_s = 0.4 at the crack tip cohesive element after a quasi-static loading, the element will enter the fatigue loading step with a 60% (A_eff = 0.6A_e) of residual area. In [30], a closed-form formulation of A_CZ was adopted alike [21] but depending in this case on the mode mixity.

Regarding ∆G, the authors of [25,26,28,29,30] used the instantaneous (i.e., at the current increment) integrated traction-displacement response of cohesive elements, where the subcritical portion (see Figure 1) was assumed to be vertical. It is worth underlining that these models calculate ∆G pointwise within A_CZ, and therefore the fatigue damage rate varies pointwise accordingly. In other words, the experimental fatigue crack growth rate dA/dN evaluated at the scale of the specimen is used to model local damage processes. The authors of [32] choose to concentrate fatigue damage into the crack tip element; therefore, the value of G was evaluated in the same way as in [25,26,28,29,30] but only at the crack tip element. Since in this way G increased with delamination growth contradicting experimental results, the value considered was calculated as a weighted average of the values of G at the crack tip and at the element just ahead of the crack tip. In [33] the authors recognized that the “instantaneous” G may vary with fatigue degradation and mesh refinement, while the value of G at the cohesive element failure is a better estimate. The strain energy release rate at failure was, therefore, extracted from elements in the wake of the numerical crack front, assuming that the variation of G between consecutive cohesive elements is small.

All of the models mentioned above consider the development of damage starting from the softening branch of the cohesive law. Since fatigue damage is likely to begin below the quasi-static stress threshold σ₀, or it is possible that it develops in absence of a pre-existing crack (i.e., cohesive element would be little stressed in this case), a condition for crack initiation was added in [26]. Here a fatigue damage initiation parameter d_fi was accumulated, depending on the stress level of cohesive elements with respect to fatigue stress-life behaviour. In order to identify correctly the initiation region, a failure index was also introduced [28] and assessed [29].

The definition of damage of Equation (5) and the link between damage mechanics and fracture mechanics represented by Equation (2) were adopted also in [34,35]. They however defined an equivalent δ that incorporates the cumulative effects of the static (D_s) and fatigue (D_f) damage on the material stiffness, as the variable to be continuously recorded during the analysis. In this way the damage, independently of its origin, follows Equation (1) where D_s is replaced by D. The resulting fatigue damage evolution is:

$$\frac{\partial {\mathrm{D}}_{\mathrm{f}}}{\partial \mathrm{N}}=\frac{{\mathrm{r}}_{\mathrm{w}}}{{\mathrm{A}}_{\mathrm{CZ}}}\frac{{\mathsf{\delta }}_0{\mathsf{\delta }}_{\mathrm{c}}}{{\mathsf{\delta }}^2}\frac{\mathrm{dA}}{\mathrm{dN}}$$

(9)

where r_w represents the relative weight of each Integration Point (IP) in the Newton–Cotes integration scheme chosen by the authors. The three IPs of each cohesive element used here, have therefore relative weights of (1, 4, 1), respectively. With a few algebraic manipulations, it can be shown that Equation (9) is equivalent to Equation (5) except the presence of the weight r_w. Again, G is evaluated at each point within A_CZ implying that the macroscopic fatigue crack growth behavior still holds at the mesoscale.

Also the model developed in [27,31,36] made use of the same link between damage and fracture mechanics established in [21], i.e., Equation (2). As in other works cited before, a focus was done on the improvement of the precision of the model, namely [27]:

-

evaluation of G as a whole model value using the contour integral;

-

numerical evaluation of A_CZ increment-by-increment during the simulation;

-

introduction of a threshold for fatigue crack growth, δ_th, that can be lower than the quasi-static stress threshold δ₀.

Mixed-mode I/II criteria were then implemented in [31] and three-dimensional crack fronts were managed in [36]. The effects of the density of microcracks on a representative interface element, A_d/A_e was related to the loss of stiffness, i.e., D = 1 − K/K₀, and the fatigue damage evolution law was developed on this basis resulting:

$$\frac{{\mathrm{dD}}_{\mathrm{f}}}{\mathrm{dN}}=\frac{1}{{\mathrm{A}}_{\mathrm{CZ}}}\frac{\mathrm{dA}}{\mathrm{dN}}=\frac{1}{{\mathrm{A}}_{\mathrm{CZ}}}B\Delta G^d$$

(10)

while quasi-static damage follows Equation (1). Stiffness-based damage is one of the damage measures found in the literature; other measures for definition of damage include [38]: (i) current traction; (ii) current separation; (iii) energy dissipated during the monotonic loading process up to the current state. All of these types were found in the models analysed previously (stiffness was used also in [19,20]). In [38] it has been illustrated, for the case of the bilinear cohesive law, that these damage measures can be uniquely converted one another. Among the different measures, stiffness-based damage is quite sensitive to the value of initial stiffness, i.e., the higher K₀, the steeper the damage evolution. In [40] a study was performed in order to define bound values of the initial stiffness for a cohesive zone model. The sensitivity of the damage evolution with respect to the initial stiffness is quite clear looking at the graphical representation of Equation (1) in Figure 2 for increasing values of K₀. This implies that, using this measure, a high value of damage may not necessarily mean a high loss of strength. In [38] it was also pointed out that an energetic damage variable may be an appropriate quantity to interpret the damage states but the formulation of cohesive laws with this measure is not recommended. Instead, the use of the separation-type damage variable is more effective and efficient for this purpose.

The development of the model [27] was done focusing on bonded joints. In finite element modelling of bonded joints, the adhesive layer, although not infinitely thin, is generally not introduced explicitly in the finite element model because of modelling and computation time convenience. The adhesive layer behavior and damage is therefore embedded into the cohesive zone, that has the purpose not only to model decohesion but also the stiffness of the bondline which may not be as high as usually assumed in CZM. For instance, a typical figure would be a 0.2 mm thick layer of an epoxy adhesive having an elastic Young modulus equal to 2000 MPa, that yields a stiffness K₀ = 2000/0.2 = 10⁴ MPa/mm. This value is 2–3 order of magnitude lower than that used typically to model nominally zero-thickness interfaces with CZM. Looking at Figure 2, damage evolution with such a value of K₀ may not be drastically steep as with order of magnitudes higher values. A stiffness-based damage variable is also used in the commercial finite element program ABAQUS, where the model [27,31,36] was implemented by means of the embedded user subroutines USDFLD and URDFIL. These motivations can justify from the numerical side the use of a stiffness-based damage measure.

Damage in the adhesive layer can develop in the form of microcavities and crazes as a result of the presence of filler particles and of the nature of the polymer matrix. Therefore, a fracture occurs in several cases in a quasi-brittle fashion. For this reason, it is believed that the use of a stiffness-based damage can be justified in this case also from a physical point of view. This might not be the same for very thin (ideally infinitely thin) interfaces such as metal-ceramic joints or, depending on the retained resin content after curing, ply-to-ply interface in composite laminates. In those cases, very high values of stiffness can be expected, that can be hardly evaluated in practice; therefore, working with a definition of damage based on the loss of interface (cohesive) strength might be more appropriate from the physical standpoint.

Recently, Reference [41] compared the performance of six models [19,21,25,31,33,37]. The conclusions drawn by this benchmark study are that all the models, with a few exceptions [37,33] are influenced by parameters not directly linked to the damage rate model, which significantly reduce their capabilities to perform effective simulations on a broad range of parameters. In particular concerning [31], the dependence of simulation results on the initial stiffness K₀ and on the cohesive strength σ₀ were indicated as the main responsible of the inaccuracy. It is worth underlining that in [31] and related papers, the initial stiffness was always set to K₀ = 10⁴ MPa/mm, as the target of the application was on adhesive joints. When the authors of [41] used this value in the simulations with the model of [31], they found a good agreement between theoretical and simulated FCG rate.

The objectives of this work is therefore to develop and validate a new formulation of the model [27] that would not be sensitive to changes in K₀ and σ₀ and, therefore, demonstrate that a stiffness-based damage measure is suitable also in the case of stiff interfaces. The steps are:

-

to identify the root causes of the poor performance of the model at high values of K₀;

-

to review critically the implementation strategy;

-

to propose a revised version of the model;

-

to assess the response of the revised model for varying K₀ and σ₀.

2. Model Implementation into Abaqus™

The model [27] was implemented into the software Abaqus™ 6.13 (Dassault Systèmes, Vélizy-Villacoublay, France) by programming the user subroutines USDFLD and URDFIL. The USDFLD subroutine allows to define a field-variable dependent material behaviour, such as the dependence of stiffness on damage, while URDFIL is used to post-process results during the analysis, in this case to obtain G from stress and opening of cohesive elements and to update the number of elapsed cycles. A field variable dependence through USDFLD can be easier to implement than a user-defined material behaviour by the UMAT subroutine, although at the expense of the introduction of an explicit solution dependence since USDFLD provides access to material point quantities only at the beginning of the increment. Therefore, the accuracy of the results may depend on the increment size, that is in this case on the value of ∆D_max set by the user.

The analysis is divided in four steps (see Figure 3):

• ramp-up until ${\mathrm{G}}_{\mathrm{th}}=\frac{\Delta {\mathrm{G}}_{\mathrm{th}}}{1-{\mathrm{R}}^2}$ and evaluate the fatigue damage threshold δ_th as the value of δ at the crack tip cohesive element;
• ramp down to zero force and remove damage possibly developed in the previous step;
• ramp up to the maximum load of the cycle;
• simulation of fatigue phase in subsequent increments (j-th), with the procedure described in the following, while keeping the load constant along the increments; ΔG^j is evaluated using the contour integral over a path surrounding the cohesive zone. In this step, both static ad fatigue damages are considered.

The workflow of the algorithm used for the fatigue crack growth simulation during step 4 includes the following two phases, that are repeated sequentially every increment j of the analysis step:

• Solution phase (USDFLD subroutine)
  For every integration point i of cohesive elements:
  • Get σ^j,i, δ^j,i
  • Initialize D^j,i = D^j−1,i
  • Update cohesive law: ${\mathsf{\delta }}_0{}^{\mathrm{j},\mathrm{i}}=\frac{{\mathsf{\delta }}_0{\mathsf{\delta }}_{\mathrm{c}}}{{\mathsf{\delta }}_0+\left({\mathsf{\delta }}_{\mathrm{c}}-{\mathsf{\delta }}_0\right)\left(1-{\mathrm{D}}^{\mathrm{j},\mathrm{i}}\right)}$ (see Figure 4)
  Update field variable FV = D^j,i and store it
  Loop over integration points of cohesive elements
• Post-processing phase (URDFIL subroutine)
  Evaluate A_cz^j
  Evaluate G^j using contour integral; ∆G^j = G^j(1 − R²)
  Initialize N^j = N^j−1
  For every integration point i:
  • Impose tentative damage increment $\Delta {\mathrm{D}}^{\mathrm{j},\mathrm{i}}=\min \left\{\Delta {\mathrm{D}}_{\max };1-{\mathrm{D}}^{\mathrm{j},\mathrm{i}}\right\}$
  • (∆D_max is a user-defined value)
  Calculate $\Delta {\mathrm{N}}^{\mathrm{j},\mathrm{i}}=\frac{\Delta {\mathrm{D}}^{\mathrm{j},\mathrm{i}}}{\mathrm{B}{\left(\Delta {\mathrm{G}}^{\mathrm{j}}\right)}^{\mathrm{d}}}{\mathrm{A}}_{\mathrm{CZ}}{}^{\mathrm{j}}$ Loop over integration points of cohesive elements
  For the entire model:
  • Find $\Delta {\mathrm{N}}^{\mathrm{j}}{}_{\min }={\min }_{\mathrm{i}}\left\{\Delta {\mathrm{N}}^{\mathrm{j},\mathrm{i}}\right\}$
  • Update damage N^j = N^j + $\Delta {\mathrm{N}}^{\mathrm{j}}{}_{\min }$ and store it for increment j + 1
  • For every integration point i:
  • Calculate $\Delta {\mathrm{D}}^{\mathrm{j},\mathrm{i}}=\Delta {\mathrm{N}}^{\mathrm{j}}{}_{\min }\frac{\mathrm{B}{\left(\Delta {\mathrm{G}}^{\mathrm{j}}\right)}^{\mathrm{d}}}{{\mathrm{A}}_{\mathrm{CZ}}{}^{\mathrm{j}}}$
  • Update damage D^j,i = D^j−1,i + $\Delta {\mathrm{D}}^{\mathrm{j},\mathrm{i}}$ and store it for increment j + 1
  Loop over integration points of cohesive elements
  Damage is shared between subroutines in a Fortran COMMON block.

3. Model Performance Checkout

3.1. Modelling

The fatigue crack growth rate dependence on model parameters is tested in the present work on a mode I double cantilever beam (DCB) joint geometry loaded with a constant bending moment, Figure 5, that yields a nominally constant G = M²/EI for increasing crack length (E is the adherends elastic modulus and I is the second moment of area of the beam section). The analysis is two-dimensional and only the upper cantilever was modelled due to symmetry. The cantilever is modelled with four-node, reduced integration plane stress elements, with an average mesh size of 1 mm (length direction) and 0.25 mm (thickness direction), respectively. The debonding interface instead, is represented through four-node cohesive elements. In order to mesh the process zone with enough accuracy (at least 4–5 elements enclosed in the cohesive process zone) even in the case of a very high cohesive stiffness, the cohesive element size has to be set to a low value, as well as it does for the maximum damage increment ∆D_max. A preliminary study was done where four combinations of cohesive element size (0.2 or 0.02 mm) and ∆D_max (0.2 or 0.02) were tested, yielding that a cohesive element size of 0.02 mm and a ∆D_max = 0.02 is a reasonable trade-off to keep the accuracy of the solution as high as possible compatibly with the simulation time. A rigid kinematic coupling is set between the cantilever and the cohesive zone.

The stiffness of cohesive elements is dependent on the current value of damage according to the procedure described in Section 2. Cantilevers are assigned isotropic elastic constants E = 70 GPa and ν = 0.3. The CZM and FCG parameters that hold for all simulations are reported in

Table 1. Cohesive zone model (CZM) and fatigue crack growth (FCG) parameters that are not varied in the simulations, from 21.

Parameter	Value
Gc [N/mm]	0.26
∆Gth [N/mm]	0.06
σ0 [MPa]	30
δc[mm]	0.0173
B [mm/cycle × (N/mm)−d]	4.443
d	5.4

. Fatigue is simulated with a load ratio R = 0 (i.e., ∆G = maximum G of the cycle).

The CZM parameters that are varied in the simulations are reported instead in

Table 2. CZM parameters that are varied in the simulations.

Parameter	Value
K0 [MPa/mm]	104	105	106
	30	30	30
σ0 [MPa]		50
		70

. The values were chosen in the following way:

-

K₀ = 1 × 10⁴ MPa/mm is the initial stiffness adopted in all the works done previously by the authors, representing a 0.2 mm thick adhesive layer of a polymer with Young’s modulus equal to 2000 MPa. Other values were set one- and two-order of magnitude greater than that, in order to represent much more thin layer such ply-to-ply interface in composite laminates;

-

load levels corresponding to G/G_c = 0.25, 0.5 and 0.75, respectively, were simulated for each of the values of K₀;

-

the influence of σ₀ has been assessed only in the case of an intermediate value of initial stiffness to limit the number of simulations (see Table 2).

The introduction of δ_th instead of δ₀ as a fatigue threshold has the effect of moving backward or forward the transition point between quasi-static and subcritical damage evolution, alike K₀ and σ₀. Since the extent of this effect depends, at the same time, on the values assumed by K₀ and σ₀, for the sake of simplicity a detailed analysis has not been carried out regarding δ_th.

3.2. Results

The FCG rate predicted by the model for the parameter combinations in Table 2 is summarized in Figure 6. The simulated crack growth varied from about 1 mm to 5 mm depending on the FCG rate. A linear regression of crack length vs. number of cycles has been done in order to extrapolate the FCG rate and to evaluate the correlation of the simulation with the theoretical trend represented by Equation (4). At a value of K₀ = 10⁴ MPa/mm the model yields results in line with the theoretical value, showing just a small deviation at the highest value of G/G_c. For increasing values of K₀, as shown in 41, an increasing deviation from the theoretical trend is found that may become very large. This confirms also why in previous works of the authors, where simulations were done with a low value of stiffness, representative of an adhesive layer (10⁴ MPa/mm), a good correspondence between theoretical and simulated FCG rate was always found.

This behavior is consistent with the dependence of the transition point between quasi-static and subcritical fatigue damage with the value of K₀ represented in Figure 7. The higher the initial stiffness and strain energy release rate, the higher is the amount of damage accumulated following the quasi-static cohesive law, i.e., non-subcritically.

Not surprisingly, the opposite trend is found in Figure 8 by increasing σ₀, where the difference between simulation and theory decreases. Thus, the higher σ₀ the lower the non-subcritical damage accumulation.

4. Modification of the Model and Validation

4.1. Modification

From the analysis conducted in Section 3, the model yields accurate predictions of the FCG rate under conditions of low values of K₀ and G/G_c and high values of σ₀, that is when damage accumulates mainly under subcritical conditions. Therefore, a term accounting for fatigue damage developing under “critical” (i.e., quasi-static) stress-opening conditions is foreseen.

The model [27] makes use of the Abaqus™ USDFLD subroutine for an easier implementation of a damage-dependent cohesive stiffness, with respect to the use of UMAT (User-defined material). However, as anticipated in Section 2, the cohesive element stiffness is not updated by USDFLD during the increment and the stress-opening distribution in the cohesive zone may therefore overshoot the quasi-static cohesive law. An example of this behavior is shown in Figure 9. The overshooting of the quasi-static limit corresponds to an unintended damage increment, ∆D*^j,i, that is not accounted for in the procedure described in Section 2. Since the overshooting of the cohesive law tends to be more pronounced at high values of K₀ and/or G/G_c a high discrepancy between theoretical and simulated FCG rate can be expected under these circumstances, as seen in Figure 6. A high stiffness leads also to a more pronounced overshooting of the cohesive law especially at the end of the elastic phase, always because of the explicit solution dependence introduced by the USDFLD subroutine.

The idea of a correction for unintended quasi-static damage was developed in [37], due in that case to the mesh-dependency of opening profile and traction field in the cohesive zone, and implemented by a novel predictor-corrector approach. The predictor step performed an approximate prediction of the total damage increment of ΔN and the corrector step handled the unintended evolution of quasi-static damage in terms of an adjustment of the number of loading cycles.

For the sake of simplicity, a first-guess, direct estimate of the damage increment ∆D*^j,i is done referring to the scheme illustrated in Figure 10:

∆D*^j,i = D*^j,i − D^j,i

(11)

where,

$${\mathrm{D}}^{*\mathrm{j},\mathrm{i}}=\frac{{\mathsf{\delta }}_{\mathrm{f}}\left({\mathsf{\delta }}^{\mathrm{j},\mathrm{i}}-{\mathsf{\delta }}_0\right)}{{\mathsf{\delta }}^{\mathrm{j},\mathrm{i}}\left({\mathsf{\delta }}_{\mathrm{c}}-{\mathsf{\delta }}_0\right)}\mathrm{if}{\mathsf{\delta }}^{\mathrm{j},\mathrm{i}}>{\mathsf{\delta }}_0{}^{\mathrm{j},\mathrm{i}}$$

(12)

$${\mathrm{D}}^{*\mathrm{j},\mathrm{i}}={\mathrm{D}}^{\mathrm{j},\mathrm{i}}\mathrm{if}{\mathsf{\delta }}^{\mathrm{j},\mathrm{i}}\le {\mathsf{\delta }}_0{}^{\mathrm{j},\mathrm{i}}$$

(13)

In other words, ∆D*^j,i represents an additional fatigue damage increment that is accumulated non-subcritically. The condition ${\mathsf{\delta }}^{\mathrm{j},\mathrm{i}}\le {\mathsf{\delta }}_0{}^{\mathrm{j},\mathrm{i}}$ corresponds to subcritical fatigue crack growth where the stress lies below the cohesive law. The field variable where damage is stored for the following increment is then updated to the value of D*^j,i.

Similarly to the corrector step in [37], the increment of damage ∆D*^j,i is then related to an increment of number of cycles ∆N*^j. In this case, it has been done by defining an equivalent damaged area by the knowledge of ${\mathrm{A}}_{\mathrm{e}}{}_{\mathrm{i}}$, the effective area of each of the i-th elements lying on A_CZ:

$${\mathrm{A}}_{\mathrm{d}}{}^{*\mathrm{j}}={\displaystyle \sum }_{\mathrm{i}\in {\mathrm{A}}_{\mathrm{CZ}}}{\mathrm{A}}_{\mathrm{e}}{}_{\mathrm{i}}\Delta {\mathrm{D}}^{*\mathrm{j},\mathrm{i}}$$

(14)

and dividing this latter by the fatigue crack growth rate, yielding,

$$\Delta {\mathrm{N}}^{*}{}^{\mathrm{j}}={\left(\frac{{\mathrm{A}}_{\mathrm{d}}{}^{*}}{\mathrm{dA}/\mathrm{dN}}\right)}^{\mathrm{j}}=\frac{{\sum }_{\mathrm{i}\in {\mathrm{A}}_{\mathrm{CZ}}}{\mathrm{A}}_{\mathrm{e}}{}_{\mathrm{i}}\Delta {\mathrm{D}}^{*\mathrm{j},\mathrm{i}}}{\mathrm{B}{\left(\Delta {\mathrm{G}}^{\mathrm{j}}\right)}^{\mathrm{d}}}$$

(15)

that has to be added to ∆N^j_min (see Section 2) to have ∆N^j = ∆N^j_min + ∆N*^j. The workflow is therefore:

• Solution phase (USDFLD subroutine)
  For every integration point i of cohesive elements:
  • Get σ^j,i, δ^j,i
  • Initialize D^j,i = D^j-1,i
  • Update cohesive law: ${\mathsf{\delta }}_0{}^{\mathrm{j},\mathrm{i}}=\frac{{\mathsf{\delta }}_0{\mathsf{\delta }}_{\mathrm{c}}}{{\mathsf{\delta }}_0+\left({\mathsf{\delta }}_{\mathrm{c}}-{\mathsf{\delta }}_0\right)\left(1-{\mathrm{D}}^{\mathrm{j},\mathrm{i}}\right)}$
  • Check for IPs exceeding quasi-static cohesive stress limit:
  • if ${\mathsf{\delta }}^{\mathrm{j},\mathrm{i}}>{\mathsf{\delta }}_0{}^{\mathrm{j},\mathrm{i}}$ update D^j,i to D*^j,i = $\frac{{\mathsf{\delta }}_{\mathrm{f}}\left({\mathsf{\delta }}^{\mathrm{j},\mathrm{i}}-{\mathsf{\delta }}_0\right)}{{\mathsf{\delta }}^{\mathrm{j},\mathrm{i}}\left({\mathsf{\delta }}_{\mathrm{c}}-{\mathsf{\delta }}_0\right)}$
  • Update field variable FV = D*^j,i and store it
  Loop over integration points of cohesive elements
• Post-processing phase (URDFIL subroutine)
  Evaluate A_cz^j
  Evaluate G^j using contour integral; ∆G^j = G^j(1 − R²)
  Initialize N^j = N^j-1
  For every integration point i:
  • Impose tentative damage increment $\Delta {\mathrm{D}}^{\mathrm{j},\mathrm{i}}=\min \left\{\Delta {\mathrm{D}}_{\max };1-{\mathrm{D}}^{*\mathrm{j},\mathrm{i}}\right\}$
  • Calculate $\Delta {\mathrm{N}}^{\mathrm{j},\mathrm{i}}=\frac{\Delta {\mathrm{D}}^{\mathrm{j},\mathrm{i}}}{\mathrm{B}{\left(\Delta {\mathrm{G}}^{\mathrm{j}}\right)}^{\mathrm{d}}}{\mathrm{A}}_{\mathrm{CZ}}{}^{\mathrm{j}}$
  Loop over integration points of cohesive elements
  For the entire model:
  • Find $\Delta {\mathrm{N}}^{\mathrm{j}}{}_{\min }={\min }_{\mathrm{i}}\left\{\Delta {\mathrm{N}}^{\mathrm{j},\mathrm{i}}\right\}$
  For every integration point i:
  • Calculate $\Delta {\mathrm{D}}^{\mathrm{j},\mathrm{i}}=\Delta {\mathrm{N}}^{\mathrm{j}}{}_{\min }\frac{\mathrm{B}{\left(\Delta {\mathrm{G}}^{\mathrm{j}}\right)}^{\mathrm{d}}}{{\mathrm{A}}_{\mathrm{CZ}}{}^{\mathrm{j}}}$
  • Calculate ∆D*^j,i = D*^j,i − D^j,i
  • Calculate $\Delta {\mathrm{N}}^{*}{}^{\mathrm{j}}=\frac{{{\displaystyle \sum }}_{\mathrm{i}\in {\mathrm{A}}_{\mathrm{CZ}}}{\mathrm{A}}_{\mathrm{e}}{}_{\mathrm{i}}\Delta {\mathrm{D}}^{*\mathrm{j},\mathrm{i}}}{\mathrm{B}{\left(\Delta {\mathrm{G}}^{\mathrm{j}}\right)}^{\mathrm{d}}}$
  • Update damage D^j,i = D*^j,i + ∆D^j,i and store it for increment j + 1
  Loop over integration points of cohesive elements
  For the entire model:
  • Update no. of cycles ∆N^j = ∆N^j_min + ∆N*^j and store it for increment j + 1

The results of this modification will be assessed in the next section.

4.2. Results

The fatigue crack growth rate obtained by adding ∆N*^j to ∆N^j_min is shown in Figure 11 as a function of G/G_c. It is evident that now the FCG rate is coherent with the theoretical value (Equation (4)) independent of K₀. Also the influence of G/G_c visible in Figure 6 is now absent. This is a confirmation that the dependence of the FCG rate on the value of initial stiffness in Figure 6 is related to the explicit solution dependence introduced by the use of the USDFLD subroutine, rather than on the choice of stiffness as a measure of damage.

Concerning the effect of σ₀, the results are summarized in Figure 12 for K₀ = 10⁵ MPa/mm. Also in this case the FCG rates respect the theoretical trend with a limited underestimation at the lowest value of G/G_c.

Globally, it can be said that the modification proposed in Section 4.1 yields good results over a broad range of cohesive initial stiffness (K₀) and cohesive strength (σ₀) for a fixed value of the quasi-static critical value of strain energy release rate.

5. Conclusions

The model developed in [27], has been reviewed critically concerning damage evolution and implementation strategy. A modification to improve the accuracy of the model in yielding the theoretical FCG rate has been proposed and the response has been assessed by performing simulations varying the following model parameters:

• Initial stiffness K₀
• Cohesive strength σ₀

The results showed that the modified model is robust with respect to changes of two orders of magnitude in initial stiffness. The influence of cohesive strength σ₀ was tested in the case of an intermediate value of stiffness K₀ = 10⁵ MPa/mm and a limited difference between theoretical and numerically simulated FCG rate was found also in this case, for changes of a factor of two in σ₀.