#### Abstract

A cyclic cohesive zone model is applied to characterize the fatigue crack growth behavior of a IN718 superalloy which is frequently used in aerospace components. In order to improve the limitation of fracture mechanics-based models, besides the predictions of the moderate fatigue crack growth rates at the Paris’ regime and the high fatigue crack growth rates at the high stress intensity factor levels, the present work is also aimed at simulating the material damage uniformly and examining the influence of the cohesive model parameters on fatigue crack growth systematically. The gradual loss of the stress-bearing ability of the material is considered through the degradation of a novel cohesive envelope. The experimental data of cracked specimens are used to validate the simulation result. Based on the reasonable estimation for the model parameters, the fatigue crack growth from moderate to high levels can be reproduced under the small-scale yielding condition, which is in fair agreement with the experimental results.

#### 1. Introduction

Most of the service failures of aerospace components such as landing gears, turbine blades, fastener joints, and so on are caused by fatigue associated with cyclic loading. The fatigue damage of structural metals normally is associated with the initiation of macrocracks (the nucleation and growth of microcracks), such that a few dominated macrocracks grow until the final unstable fracture happens. By taking into account material degradation instead of one single major crack, continuum damage mechanics provides a sophisticated framework for the fatigue crack initiation analysis of metallic materials [1], while the fracture mechanics-based models are dominated approaches in simulating fatigue crack growth [2].

One of the extensively used models for predicting fatigue crack growth is the Paris’ law [3] which relates the range of the stress intensity factor and the crack growth rate through a power law in double logarithmic coordinates. Generally, the relationship between and for most of metallic materials can be divided into three stages [4], as illustrated in Figure 1. Regime I corresponds to the formation of a fatigue crack at the low , where crack growth cannot be observed if is less than the threshold . Regime II usually maintains the Paris’ law, which exhibits the moderate fatigue crack growth rate (i.e., mm/cycle mm/cycle). The high crack growth rate (i.e., mm/cycle) achieves at Regime III in which asymptotically approaches the fracture toughness . Fatigue life predictions in engineering design have much success with the modifications of the Paris’ law to determine crack growth rates [5, 6].

Since elastic or the small-scale yielding (SSY) condition at crack-tip is assumed in the -based models, many attractive models have also been proposed to simulate fatigue crack growth with considerable plastic deformations [6–10]. One representative model proposed by Dowling and Begley [7] is constructed by replacing with -integral to extend the Paris’ law into the large-scale yielding (LSY) conditions. However, some researchers [11, 12] have pointed out that -integral is invalid upon unloading. Although other researchers [10] found that -integral still maintains its path independence under cyclic loadings, nevertheless, it is a considerable issue how to use -integral correctly in fatigue crack growth simulations.

The cohesive zone model (CZM) is another reliable approach which can describe the gradual degradation of stress-carrying ability at a propagating crack-tip within the framework of the nonlinear fracture mechanics. CZMs were firstly proposed by Dugdale [13] and Barenblatt [14] to consider the presence of process zone at crack-tip, where the local degradation during crack growth is determined by a given constitutive law (cohesive law) to relate the traction with the magnitude of material separation. Except the applications in brittle [15] and composite materials [16], cohesive models have been widely used in simulating the ductile fracture of metallic materials. Tvergaard and Hutchinson [17] introduced a trapeziform cohesive law to investigate crack growth resistance during the fracture process of an elastic-plastic solid. Chen et al. [18] simulated the ductile crack growth of the thick CT specimens made of the pressure vessel steel using a polynomial cohesive law. Uthaisangsuk et al. [19] applied a CZM between martensitic islands and ferritic grains to consider interface debonding in the dual-phase microstructure of DP600 steel. Li et al. [20] studied the ductile fracture of a notched Ti-6Al-4V specimen using a CZM. Several authors [21] summarized classifications and applications of the CZM, which demonstrate the importance of cohesive models in simulating ductile fracture.

For fatigue crack calculations, the material degradation under cyclic loading must be considered; namely, the cyclic cohesive zone model (CCZM) should be applied to simulate fatigue crack growth. So far, there are many existing CCZMs that are able to describe the linear Paris’ region (Regime II) behavior. The first CCZM disinterring the Paris-like behavior of fatigue crack growth for metallic materials was proposed by Nguyen et al. [22]. They simulated fatigue crack growth of a central cracked aluminum specimen and demonstrated that the CCZM is capable of reproducing the Paris’ law. The unloading-reloading hysteresis was applied to account for dissipative mechanisms, which prevents the “shakedown” effect efficiently in the crack growth process. However, they did not mention the predictive ability of the model in Regime III. Maiti and Geubelle [23] introduced a bilinear cohesive law with an irreversible unloading-loading path to assess capacity of the model for predicting fatigue crack growth in polymers. Although they concluded that the developed model can be used in both Regime II and Regime III for an epoxy, they did not show the applications of the proposed model in simulating ductile materials. Roe and Siegmund [12] considered the effects of the interface fatigue crack growth under different loading modes using an irreversible CZM, and while most of the simulation results were performed under the linear Paris’ region, high rates associated with Regime III were not demonstrated. Ural et al. [24] used a damage-based cohesive model to predict the fatigue crack growth behavior in Regime II of an aluminum alloy. However, the application of this model has not been identified in Regime III also. Li and Yuan [25] used a normal stress-dominated cohesive model to simulate fatigue crack growth of cracked and notched specimen and reproduced the linear Paris-like behavior at different loading ratios; nonetheless, the predictive ability in Regime III has not been considered in the CZM. Jimenez and Duddu [26] proposed a strain energy release rate-dependent CZM to model high cycle fatigue delamination of the composite materials under both mode I and mixed mode loadings, and they pointed out that the CZM is suitable for the constant amplitude loading in Regime II and also expected to extended the model to describe all three regimes in fatigue crack simulations. CCZMs can also be applied in analyzing the elastoplastic fatigue cracks by introducing some stress triaxiality-dependent parameters [27, 28]; that is, Jha and Banerjee [28] proposed a triaxiality-dependent cohesive model to simulate the fatigue failure of ductile metals and proved that the CZM is suitable for the geometries with different stress states. Similarly, these models [27, 28] do not mention their applications for fatigue crack growth at the high levels in Regime III.

Although the Paris-like behavior which corresponds to the moderate levels in Regime II can be predicted by using the aforementioned CCZMs [11, 22-27], this is not the advantage of CCZMs, since the Paris regime associated with SSY has already been described successfully according to the -based models. In order to resolve the more challenging task that to simulate the high fatigue crack growth rates in Regime III of the metallic materials, a CCZM was proposed in our previous work [29] in which two damage variables were defined to represent monotonic damage and fatigue damage, respectively. In fact, it is not easy to distinguish these two different damage mechanisms in fatigue crack growth before the final rupture occurs; thus, it is better to describe the material damage evolution uniformly. Meanwhile, it is anticipated that CCZMs have the potential to predict high fatigue crack growth rates in Regime III under the LSY conditions. For the purpose of easing the difficulty level for cohesive modeling of fatigue crack growth at LSY, the predictive ability of CCZMs in simulating high fatigue crack growth rates at SSY should be verified systematically.

Therefore, the objective of this work is to examine a CCZM for predicting high fatigue crack growth rates corresponding to Regime III under the SSY conditions. Meanwhile, the present study attempts to describe the material damage evolution uniformly and identify the effects of the cohesive model parameters on fatigue crack growth under both moderate and high levels. The contents of this paper include the following aspects: in Section 2, the detailed formulations of the CCZM are introduced. In Section 3, the identification of the model parameters is performed firstly, then the numerical results are presented according to the compact-tension (CT) specimens and the single-edge notched (SEN) specimens made of the IN718 superalloy and validated against the experimental data from literatures. Finally, the major conclusions are summarized in Section 4.

#### 2. Cohesive Zone Model

The cohesive model has the ability to consider the ductile rupture of the metallic materials under various loading conditions. For fatigue crack growth simulations, the CCZM has to be able to describe accumulative damage associated with cyclic loading. Assume that the degradation of the original cohesive envelope occurs with damage accumulation, and the constitutive law under fatigue loading is defined as where the normal traction is assumed to be a piecewise function depending upon the normal separation . is the characteristic length of the cohesive law, which is denoted as the cohesive length. denotes the maximum normal separation above which the normal traction vanishes. Following the concept from the continuum damage models, the damage variable stands for the material degradation and represents failure if approaches 1. is the cohesive strength used to identify the peak value of the normal traction and is the initial cohesive strength. The original cohesive law under monotonic loading can be obtained by substituting into (1). The parameter corresponds to the material separation relating to the damage variable and . For both constant and variable amplitude loadings, the damage evolution equation suggested by Roe and Siegmund [12] can be applied to reflect the accumulative damage for the metallic materials, which reads where is the accumulated cohesive length used to scale the increment of the normal separation and . is the cohesive endurance limit. The Macaulay brackets are used to signify nonnegative value of damage evolution, and the Heaviside function states that damage initiates at .

The novel cohesive envelope determined by (1) will degrade and change its shape correspondingly associated with the damage accumulation, and this characteristic is the major difference between the current model and other CCZMs [11, 12, 20, 22–28]. For example, in [25], the combination use of (2) and the cohesive envelope of the exponential form are performed to simulate fatigue crack growth; however, this combination results in the contradictory damage evolution at if a slightly larger value of is provided, and this problem can be solved by using (1). As mentioned above, the similar CCZM has been proposed to investigate the mixed mode fatigue crack growth in our previous simulations [29], the difference between these two models is that the current one does not introduce the monotonic damage, the material degradation under cyclic loadings is described uniformly by (2), and the monotonic damage is determined by the cohesive envelope implicitly. Actually, the combination of (1) and (2) has the capacity to simulate fatigue crack growth at the high levels, the applications need not be distinguished between monotonic damage and fatigue damage as introduced in [29], the key point in such simulations is the form of (1), and more details will be performed in Section 3.

For the numerical implementation of the model under each unloading-reloading cycle, assume that the residual separation exists after a complete unloading cycle as plastic deformations, and the corresponding unloading-reloading path can be denoted by with where is the traction at time related to the separation according to the cohesive stiffness . and are the traction and the separation at time , respectively.

As shown in Figure 2(a), assume that the rising segment () of the loading sequence results in the separation of a material point at the cohesive zone. In order to explain the operating mechanism of the model under cyclic loading in detail, two different cases are considered, namely, (Figure 2(b)) and (Figure 2(c)). Upon unloading, the path () follows (3) without damage evolution. The contact of two cohesive surfaces is constrained by introducing a penalized equation as where is the penalty stiffness. Here, we should mention that this contact algorithm is applied to simulate crack closure explicitly. The crack closure effects are related to the plastic deformation of the bulk material around the crack-tip implicitly under cyclic loading.

**(a)**

**(b)**

**(c)**

The reloading path () specified by (3) deviates from the previous unloading path due to damage evolution. At the subsequent loading segment ( 4), the value of corresponding to reloading may reach the degraded cohesive envelope at . Upon further reloading, the value of moves along the degraded envelope at time ; however, the continued damage accumulation through (2) makes it move along the subsequent degraded cohesive envelope at time . The dot lines in Figures 2(b) and 2(c) are two degraded cohesive envelopes corresponding to and , respectively. In both cases, if the loading level makes the reloading path move along the descending part of the degraded cohesive envelop, there is no damage accumulation that will occur and the material degradation is described by (1) implicitly.

The above CCZM has been implemented based on the user interface UMAT of the commercial code ABAQUS via the cohesive element [30]. Since the numerical simulations are performed according to the ABAQUS/Standard solver, both the updated cohesive traction and the material Jacobian matrix should be provided for constructing the internal nodal force vector and the cohesive stiffness matrix in UMAT. The standard Newton method is applied to solve the global nonlinear equations, and the default values for the solution control parameters are adopted. The global convergence rule is checked at the end of each time increment, such that the tolerance limit on the residual force can be fulfilled.

#### 3. Model Validation

In the present study, the influence of model parameters on fatigue crack growth will be verified systematically; meanwhile, the predictive ability of the model for the high fatigue crack rates in Regime III will be validated experimentally. The material under investigation is a nickel-based superalloy IN718 which is widely used in aeroengine components due to its high yield stress, excellent fatigue resistance, and so on in severe conditions [31]. In modeling a CT specimen with the initial crack, the length is mm and the thickness is mm. The material properties are as follows: Young’s modulus GPa, Poisson’s ratio , the initial yield stress , and the plane strain fracture toughness . The plasticity is applied to describe the plastic behavior of the alloy, and the relationship between true stress and true plastic strain is described by the Holomon’s equation as where the bulk modulus MPa and the strain hardening exponent . We note from the perspective of cohesive modeling for fatigue crack that the isotropic hardening is often assumed to simulate the material response approximately without knowing the exact cyclic deformation behavior of the bulk material [11, 28, 32]. However, the more accurate cyclic plasticity model should be determined experimentally as mentioned in [33]. Due to the symmetry, only the upper half of the specimen is modeled under the plane strain condition, as shown in Figure 3. The whole finite-element model consists of 1725 elements and 1872 nodes. The cohesive zone is predefined along the crack growth path using 240 COH2D4 elements which are uniformly distributed along the ligament of the specimen with the element length of 0.1 mm. The CPE4 elements are used for the rest of the regions.

##### 3.1. Identification of Model Parameters

Based on previous descriptions, the model parameters can be divided into two types. The monotonic model parameters include cohesive strength , cohesive length , and cohesive energy . The cyclic model parameters are the accumulative length and the cohesive endurance limit . Omitting shear stress in the cohesive zone, only the normal traction is active under the mode I crack growth condition. The cohesive energy is determined by integration of the monotonic cohesive envelope as follows: where the form of is represented by (1) with . The monotonic cohesive envelope can be constructed through the initial cohesive strength , the cohesive length , and the cohesive energy . Under the SSY conditions, the cohesive energy is related to the plain strain fracture toughness based on the Irwin’s relation as follows:

According to (8), the cohesive energy is in the present analysis. The cohesive strength is used, which is nearly three times above the initial yield stress of IN718 according to the suggestion in [34], where the corresponding cohesive length . Actually, the relationship between the initial cohesive strength and cohesive length can be attributed to the effect of the initial cohesive stiffness . The choice related to should guarantee that the cohesive zone does not influence the overall compliance before damage initiation and normally a higher value of is used. Since the cohesive energy is the work needed to create a unit area of two mating fracture surfaces, the resulting is obtained according to (1) with .

Identification of the influence of the model parameters is performed using the CT specimen under constant loading amplitudes at the loading ratio , where with and are the maximum and minimum of the uniaxial applied force , respectively. The nominal stress intensity factor is applied to describe the stress field at crack-tip. The range of the stress intensity factor for the CT specimens can be calculated through [35]: where . This formulation is accurate within 0.5% error over the range 0.2 1.

The effect of the accumulative length to the crack growth length curve is indicated in Figure 4(a), where the initial stress intensity factor . At a given crack growth length , the lower value of corresponds to the higher crack growth rate for the same value of . In comparison with the number of loading cycles to the final failure between and , the critical crack growth length of approx. 5.5 mm is achieved under both and at failure. Theoretically, the critical crack length should be reached at also; however, the FE calculation is aborted since the high damage evolution rate associated with is so fast that the mesh resolution of 0.1 mm at crack-tip cannot capture the crack growth rate before the final failure. The effect of on the crack growth curves is shown in Figure 4(b) by comparing two different values of at , where the higher crack growth rate is accompanied by the lower value of at a given crack length. Based on Figure 4, we can conclude that the critical crack length is only determined by the monotonic model parameters together with the level.

**(a)**

**(b)**

The effect of on the versus curve is shown in Figure 5(a), where the data are determined based on the corresponding curves (see Figure 4). In order to examine the effects of cyclic model parameters on fatigue crack growth in the Paris’ region, the numerical results of the lower loading level kN associated with is also displayed. The accumulative length acts like the intercept coefficient which operates in the Paris’ law, that is, . It plays a role in increasing or decreasing crack growth rate in all portions at a certain value of . In contrast, has a tendency to modify the slope of the versus curve as the parameter of the Paris’ law, as shown in Figure 5(b). However, the effects of the cyclic model parameters are also influenced by the monotonic model parameters; namely, the crack growth rate is accelerated inevitably as the level approaches the critical value in Regime III before the final failure. The above results reveal that the high fatigue crack growth rates at high levels can be captured reasonably based on the current model and the uniform damage evolution equation can be applied in such simulations.

**(a)**

**(b)**

##### 3.2. Comparison with Experiments

After the parameter identification, the next the predictive ability of the model will be examined using the CT specimen as displayed in Figure 3(a). The experimental data in validations are obtained by digitizing the results reported by Yamada and Newman [36] and Newman et al. [37], as shown in Figure 6, where the fatigue crack growth tests were performed by using the CT specimens and the SEN specimens at , respectively. The experimental results are fitted according to the Erdogan-Ratwani model [5]. For the CT specimens, the experimental points can be well approximated by the following equation: where is in mm/cycle and is in . For the SEN specimens, all points including the threshold region can be well described by the following equation:

**(a)**

**(b)**

Three constant loading amplitudes, that is, , 18.5 kN, and 20 kN at , are applied to obtain the cyclic cohesive model parameters, where the corresponding initial stress intensity range are , , and , respectively. The fatigue crack growth simulations are performed to fit the experimental data associated with the moderate growth rates. As depicted in Figure 7, the fitting curve constructed by the aforementioned three levels correlates the experiments very well in the linear Paris’ region. The corresponding crack growth versus the applied loading cycle curves are displayed in Figure 8, where the fastest crack growth rate is achieved under the highest level at a given crack length. The determined cyclic model parameters of IN718 are listed in Table 1.

For the fatigue crack growth simulation in Regime III, the constant loading amplitude is applied with at , which corresponds to the initial stress intensity factor range . In order to capture the local response of each material point along the crack growth path, four material points 1, 2, 3, and 4 located at the crack growth path are used to trace the material responses associated with the high propagation rates, as illustrated in Figure 3(c). The distance between two adjoining monitory point is 1.2 mm.

The traction-separation responses during crack growth are indicated in Figure 9 for both points 1 and 4, where the material damage is obtained based on (2) for the two material points. The traction of point 1 approaches 87% of the initial cohesive strength at the first loading cycle while the traction of point 4 is still small; namely, the traction of point 4 is only 10% of at the first loading cycle, as shown in Figure 10. Except for point 1 loses its stress-carrying ability at the 207th cycle, other points have no damage accumulations at the first 548 loading cycles such that their original cohesive stiffnesses still exist, since there are no damage evolutions for these points (see Figure 11). The damage initiation for point 4 starts after 1500 loading cycles, where the maximal value of the traction increases (up to 80% of ) at first and then decreases with the applied loading cycles until the material point is completely damaged (see Figure 10(b)), and this phenomenon is also true for points 2 and 3.

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

During the crack growth process, the distributions of the damage zone, the von-Mises stress, and the equivalent plastic strain under both the moderate and the high levels of are compared in Figure 12. The fatigue crack propagates from length mm to 29.5 mm by measuring the length of the complete damage zone of the two levels of . For two cases, the von-Mises stress ahead of the current crack-tip is greater than the initial yield stress 1196 MPa and the obvious plastic zone is created gradually in the crack wake during crack growth. These results show that the material responses of both the cohesive zone and the surrounding bulk material are simulated reasonably.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

Figure 13 summarizes the numerical and experimental crack growth results under constant amplitude loading. The simulation associated with is the predicting result based on the determined model parameters listed in Table 1. The whole versus curves are assembled through these four different levels. With increment of the crack length, becomes very high and approaches instability limited by the at failure (see Figure 13). As indicated in Figure 1, fatigue crack growth in Regime II is largely influenced by the loading ratio and the frequency in the air; therefore, fatigue crack growth with moderate rates is irrelevant to the geometries of specimens. At high propagation rates as approaches to , there is a good agreement between the prediction and the experiment of the SEN specimens, while a deviation between the simulation and the experiment of the CT specimens exists. Actually, the difference in experiment data also occurs between the CT and SEN specimens in Regime III due to the occurrence of ductile fracture. The ductile fracture process of metallic materials involves nucleation, growth, and coalescence of microvoids, which is greatly influenced by the stress state and may not be described, for example, by (1) simply, and the local constraint effect accompanied by the specimen geometry also affects the crack growth resistance. Nevertheless, based on the cyclic model parameters determined by the experimental data in Regime II, the numerical result captures the accelerated fatigue crack growth phenomenon of the experimental data in Regime III associated with the high levels.

#### 4. Conclusions

In order to overcome the limitations of fracture mechanics-based models, for example, the Paris’ law or the -integral-based model, for predicting elastic-plastic fatigue crack growth rates, the cyclic cohesive model is verified systematically on the fatigue crack growth of the IN718 alloy from moderate to high levels. Since ductile fracture occurs as the level approaches the fracture toughness, the model parameters determined from the experimental data at the moderate levels cannot predict the fatigue crack growth rates at the high levels for all the specimens exactly, and the local constraint effect should be considered during the damage evolution in Regime III.

The simulations reveal that the novel cohesive envelope combined with one damage evolution equation can also be used to capture the accelerated fatigue crack growth behavior at the high levels, and there is no need to distinguish between the monotonic damage and the fatigue damage as introduced in [29]. The crack growth behaviors are largely affected by the accumulative length and the cohesive endurance limit ; meanwhile, the effects of the cyclic model parameters are also influenced by the monotonic model parameters, that is, cohesive energy and the initial cohesive strength as the level approaches the critical value corresponding to the fracture toughness.

After the model parameter identification and the experimental validation, we confirm that the current model has the capacity to simulate the elastic-plastic fatigue crack growth behavior of IN718 at SSY, especially in Regime III. Furthermore, it remains for the next work to extend the application of the CCZM for simulating fatigue crack growth of other alloys at LSY.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This research was made possible through research funding provided by the National Natural Science Foundation of China (no. 11502204), the National Key Research and Development Program of China (no. 2016YFB0701303), the Natural Science Basic Research Plan in Shaanxi Province of China (no. 2017JM1023), and the China Postdoctoral Science Foundation funded project (no. 2017M623234).