Experimental and Numerical Evaluation of Equivalent Stress Intensity Factor Models under Mixed-Mode (I+II) Loading

1. Introduction

The widely accepted Paris rule correlates the stress intensity factor (SIF) with crack length advance per cycle. However, the Paris rule is commonly established under mode I loading. So, when more than one loading mode is present, a crack increment rule under mixed mode and an angle prediction criterion is needed [1], [2]. This is one rule to quantify crack length increment and another to quantify crack direction [3]. In the case of mixed-mode fracture, SIFs appear in orthogonal directions, each contributing to crack growth and crack kinking. Several equivalent SIF (Keq) models are based on different principles [4]. So, it is necessary to identify the Keq models that best account for the combination loading modes [5]. Therefore, when more than one loading is present, the choice of Keq impacts the crack growth rate [6]. Based on previous studies [7], [8] suggesting a significant decrease in the toughness of metals when mixed mode failures occur, this paper tests three Keq models using two different numerical approaches in combination with three fatigue crack growth rules. Numerical results are validated with experimental data from the literature at ambient temperature.

Recently, the importance of Keq in fatigue crack growth has been acknowledged in Aluminum samples [9] and both Aluminum and steel specimens [10]. Additionally, numerical studies [11] investigated the effects of mode-mixity of surface cracks, observing a difference in crack growth rate under mixed mode in steel. Moreover, [12] reproduced FCG in a double cantilever sample of low-carbon steel. Further research by [13] tested the Tanaka model against numerical simulations. Plate samples were tested under modes I and III with different Keq criteria [14]. On the other hand, [15] compared various Keq models using the Paris rule and found an acceptable agreement for the Tanaka and Irwin models. Furthermore, [16] explored various crack growth models considering load amplitude. Additionally, [17] focused on the Tanaka Keq model combined with the Paris rule, finding a good correlation but lacking comparison to other models. Consequently, the literature review revealed no consensus on the Keq model that best describes FCG.

The structural integrity of an element is characterized by how it experiences failure. Which is categorized into three fracture modes: mode $I$, where the stresses act perpendicular to the crack, generating the separation of the opposing faces; mode II, characterized by shear stresses causing in-plane crack edge displacement; and mode III, involving out-of-plane stress, usually associated with elements subjected to torsion [4]. In Linear elastic fracture mechanics (LEFM), the stress field ahead of the crack tip is described by the SIF in each direction. Additionally, LEFM provides tools to assess the lifespan of mechanical elements with flaws [13]. Moreover, a crack is a flaw in a solid, which can be produced by a stress gradient, corrosion, plastic deformation, fatigue [4], or thickness reduction [18], which requires adequate stress and time [13]. A body subjected to alternating loads, below its elastic limit, presents a localized plastic phenomenon known as fatigue, which causes the crack to grow in length and change its direction if more than one stress direction is involved [1], [2]. Mixed mode loading may also occur due to multiaxial loads, a stress concentrator, or a crack. Furthermore, the fracture toughness (Kc) measures the material's resistance to crack propagation [19] as a parameter for assessing failure in cracked bodies. However, in practical case scenarios, a body is subjected to multiaxial loads, so it is necessary to consider a mixed-mode fracture [20]. Here is where the models to quantify the combination of SIF under mixed modes come into play [5].

Moreover, FCG is usually conducted using analytical methods for simple loading cases, numerical methods [11], [12], [13], a combination of both, or even with non-parametric models [21]. A review of some of them can be found in [22]. The choice of the Keq may impact crack growth prediction. In this study, we compare the finite element method implemented in ANSYS's (Separating Morphing and Adaptive Remeshing Technology (SMART)) with the dual boundary element method (BEM), which is implemented and described elsewhere [23] for FCG under mixed-mode (I and II).

In this paper, a comparison is made between three Keq models combined with three FCG models, each covering different zones of the Paris plot. There has been a lack of studies concurrently exploring multiple Keq models with different crack growth models to assess the combination that best matches experimental results. This study seeks to fill this void by providing a comprehensive analysis that evaluates the performance of different Keq and crack growth model combinations under varying experimental conditions.

2. Background

This section provides the formulations used on the paper.

Equivalent SIF models

Despite the wide variety of models available, this paper evaluates three, each with a distinct formulation. Some models are experiment-based, while others are based on numerical data fitting.

The Tanaka criterion [24], [25] is derived from the Weertman dislocation model, which determines that crack propagation occurs when strains in the plastic zone reach the critical value. Equation (1) represents the Tanaka model, which has demonstrated a close correlation between experimental and numerical results, particularly in for steel specimens [26].

$${\mathrm{K}}_{\mathrm{Tan}\mathrm{aka}}={\left({\mathrm{K}}_{\mathrm{I}}^4{+8\mathrm{K}}_{\mathrm{II}}^2+\frac{\left({8\mathrm{K}}_{\mathrm{III}}^4\right)}{1-\mathrm{v}}\right)}^{0.25}$$

(1)

where Ki is the SIF, and v is the Poisson Module. Moreover, Asaro's analysis of stresses at the crack tip derived from the consideration of the release of potential energy [27] is represented by Equation (2). Furthermore, due to the influence of the Poisson ratio, this model depends on whether plane stress or plain conditions if mode III is present.

$${\mathrm{K}}_{\mathrm{Asaro}}={\left({\mathrm{K}}_{\mathrm{I}}^2{+\mathrm{K}}_{\mathrm{II}}^2\right)}^{0.5}$$

(2)

Pook [28] proposed the criterion shown in Equation (3) for a beam subjected to three-point bending. Although the tests on the specimen initially anticipated a dominant mode I, the crack tip rotation occurred until eventually aligned with the specimen's side.

$${\mathrm{K}}_{\mathrm{Pook}}=\frac{0{.83\mathrm{K}}_{\mathrm{I}}+\sqrt{0{.4489\mathrm{K}}_{\mathrm{I}}^2{+3\mathrm{K}}_{\mathrm{II}}^2}}{1.5}$$

(3)

So, one can see that an appropriate equivalent SIF model defines the accuracy of the lifespan estimation [29].

Crack growth models

Understanding crack growth is essential for assessing the lifespan of structural components. This section examines three models.

The Paris rule, shown in Equation (4), establishes a relationship between the crack growth rate and the equivalent SIF [30]. The constants C and m depend on the material, geometry, and test conditions. However, this model is only valid for stable crack growth and does not account for the load inversion rate [31].

$$\frac{\mathrm{da}}{\mathrm{dN}}=\mathrm{C}\cdot {\left({\mathsf{\Delta }\mathrm{K}}_{\mathrm{eq}}\right)}^{\mathrm{m}}$$

(4)

The Klesnil-Lukas model, represented by Equation (5), introduces the concept of the fracture threshold, Kth, which is particularly useful in cases where the applied load either does not cause crack growth or barely exceeds the threshold value. The constants C and m remain consistent with those in the Paris rule [32].

$$\frac{\mathrm{da}}{\mathrm{dN}}=\mathrm{C}\cdot {(\mathsf{\Delta }\mathrm{K}}_{\mathrm{eq}}^{\mathrm{m}}{-\mathsf{\Delta }\mathrm{K}}_{\mathrm{th}}^{\mathrm{m}})$$

(5)

A more detailed model, the modified Forman-Newman, can describe the three regions of crack growth by incorporating the threshold and fracture toughness, Kc, as depicted in Equation (6).

$$\frac{\mathrm{da}}{\mathrm{dN}}=\frac{\mathrm{C}{\left({\mathsf{\Delta }\mathrm{K}}_{\mathrm{eq}}\right)}^{\mathrm{m}}{\left({\mathsf{\Delta }\mathrm{K}-\mathsf{\Delta }\mathrm{K}}_{\mathrm{th}}\left(1-0.82\mathrm{R}\right)\right)}^{\mathrm{p}}\left(1+\frac{0{.82\tan }^{-1}\left(\mathsf{\eta }\mathrm{R}\right)}{\mathsf{\eta }}\right)}{{\left(1-\frac{{\mathsf{\Delta }\mathrm{K}}_{\mathrm{eq}}}{\left(1-\mathrm{R}\right){\mathrm{K}}_{\mathrm{c}}}\right)}^{\mathrm{q}}}$$

(6)

where the additional constants p, q, and η are adjusted to experimental results, and R represents the load reversal ratio.

Figure 1 illustrates the comparison between the Paris, Klesnil, and modified Forman-Newman models for R = 0.1. The Klesnil model exhibits conservative behavior in Region I. As it enters Region II, its trend converges with the Paris rule. Finally, the modified Forman-Newman model predicts less conservatively, with higher growth rates, and can predict unstable growth in Region III. In all models, the lifespan of each crack increment Δa or increment in cycles ΔN can be estimated by separating variables in Equation (4). In the case scenario of mixed-mode loading, the crack increment is a function of Keq or ΔKeq for alternating loads [2]. Integrating the Paris rule yields Equation (7), incorporating established parameters and introducing the variable T to accommodate additional terms alongside the constants C and m as described in the Klesnil-Lukas and Forman-Newman models.

$$\underset{0}{\overset{\Delta \mathrm{a}}{{\displaystyle \int }}}\frac{\left(\mathrm{da}\right)}{\mathrm{C}\cdot \mathrm{T}\cdot {\left({\mathrm{K}}_{\mathrm{eq}}\right)}^{\mathrm{m}}}=\underset{0}{\overset{\Delta \mathrm{N}}{{\displaystyle \int }}}\mathrm{dN}=\Delta \mathrm{N}$$

(7)

Moreover, fully defining the crack path under mixed-mode requires establishing the crack kinking angle $\mathsf{\theta }$, as shown in Equation (8). The model is derived from the maximum tangential stress (MTS) criterion [13]. Recent evaluations of several models based on experimental data concluded that, in some cases, the computational cost could be high despite reducing the error in prediction [2].

$$\mathsf{\theta }=2\arctan \left(\left(\frac{{\mathrm{K}}^1}{{4\mathrm{K}}^2}\right)+\mathrm{sgn}\left(\frac{1}{4}\right)\sqrt{{\left(\frac{{\mathrm{K}}^1}{{\mathrm{K}}^2}\right)}^2+8}\right)$$

(8)

3. Materials and Methods

A modified C(T) specimen data, as shown in Figure 2, was retrieved from the literature [33], where SIFs were experimentally estimated for different crack lengths. The sample was made of low-carbon steel, and had a material composition detailed in Table 1. The specimen was drilled to modify the stress fields ahead of the crack tip, inducing mixed modes I and II. Sample dimensions are shown in Figure 2a.

The modified C(T) specimen initially developed a horizontal crack, which later curved due to the influence of the drilled hole. The SIF ranges in modes I and II were found by fitting displacement fields to Williams' displacement model in an infinite plate using the digital image correlation (DIC) technique, as extensively detailed in [33]. So, calculations were conducted for six crack lengths, designated as 0, a, b, c, d, and e, respectively, as illustrated in Figure 2b, with numerical values in Table 2. Extensive details about SIF calculations are available in [33].

Once the experimental results were available for evaluating the numerical models, obtaining the crack growth constants was necessary, as depicted in Table 3. For the Paris rule, tests with R = 0.1 are required. As for the Forman and Modified Forman-Newman models, data for low-carbon plain steels can be found in the literature [34].

So, to assess the equivalent SIF models, both the Finite Element Method (FEM) and the Dual Boundary Element Method (DBEM) [35] were employed. ANSYS introduced the SMART methodology, which allows for dynamic mesh adjustments during the simulation progression [9]. This technology has been recently assessed and successfully compared to experimental results [36]. The geometry, as illustrated in Figure 2a, was modeled in ANSYS, using quadratic tetrahedral meshing, with an influence sphere of 20 mm around the crack tip and element size of 0.5 mm, this resulted in a mesh compromising 127305 nodes and 91640 elements shown in Figure 3a. Boundary conditions consisted of free cylindrical support in the radial direction and an average bearing load of 5420 N. The analysis configuration involved 15 steps, aligning with the experimental crack length.

On the other hand, the DBEM simplifies 2D problems to one dimension, facilitating meshing and reducing computational costs [35]. The GID software was used for the geometric modeling, as shown in Figure 3b, which displays the mesh nodes. Boundary conditions were set for a plane stress problem, resulting in 458 nodes and 229 elements. The lower left corner was fixed in both directions, and the upper left corner was constrained in the x-direction. Elements of 0.5 mm in length were used at the crack tip; and the mesh was refined by implementing 30 elements around the drilled hole, which the crack tip approached as it grew. The implementation of the DBEM method was carried out using an in-house code described somewhere else [23]. It was noted that the computation time was shorter than with FEM. Each simulation lasted approximately 18 seconds, unlike FEM simulations, which took about 15.37 minutes using the same computer.

4. Results

Parametric models

Table 4 displays the SIF values obtained using the equivalent SIF models previously described. Employing the presented crack growth models is feasible based on the results obtained. Since the crack increment Δa is known, Equation (7) can be applied to determine the estimated number of cycles N. This procedure was executed for each growth model. As a result, graphs illustrate the calculated number of cycles for each interval about the experimentally recorded cycles and the percentage error.

The results are presented in five intervals, each representing 0-a, a-b, and so forth.

Paris rule

In this case, a conservative prediction is obtained (Figure 4). The error decreases after the third interval once the curvature of the crack becomes noticeable. The SIF model with the lowest mean percentage error, 40.9, corresponds to the one presented by Tanaka.

Modified Forman Newman

When evaluating the Forman-Newman model (Figure 5), the prediction for the first interval diverges from the conservative trend. However, subsequent intervals exhibit a deviation from the experimental results. At the point of the highest growth rate, the prediction tends to be conservative. Notably, both projections are significantly close in the fourth interval. Beyond this point, as the crack curve assumes higher values, indicating the onset of unstable crack growth. However, the number of experimental cycles remains higher. This phenomenon could be attributed to the yielding at the crack tip due to the material ductility.

Klesnil

Finally, the growth model that presented the lowest average error percentage was Klesnil. This behavior can be attributed to the influence of the fatigue threshold, ${\mathrm{K}}_{\mathrm{th}}$, indicating that its predictions are more accurate for cracks propagating in Regions I and II. When comparing the material's toughness with the recorded ${\mathrm{K}}_{\mathrm{I}}$ values, it becomes evident that the material is in the transition zone towards Region II, where the growth rate tends to stabilize. Moreover, by using Paris's constants, Klesnil incorporates a specific R. In Figure 6, a conservative trend in cycle predictions is observable, consistent with the previous models, and it satisfactorily aligns with the recorded experimental value in the fourth interval.

The results obtained from the growth models reveal that Asaro's SIF model demonstrated the lowest average percentage error for each of the proposed crack growth models, as shown in Figure 7.

The Klesnil model yields intermediate results compared to the other models during the initial intervals. As the growth rate increases, all three models tend to linearize. Upon closer examination of the Klesnil model, its proximity to K_th might be the cause. It is noteworthy that both this model and the Modified Forman-Newman model only consider mode I in their adjustment parameters, as there is no specific methodology to calculate them in mixed mode. However, studies suggest that the initial crack propagation energy is higher when the mixed mode is present, especially when the presence of K_II is noticeable and K_I decreases [37].

Figure 8 compares the most suitable combination of the crack growth model and equivalent SIF (Asaro-Klesnil) with the experimentally obtained data. The results tend to be slightly conservative, with a difference of 10% in the total number of cycles. This trend can be attributed to various phenomena that smooth the crack tip, either due to the material's chemical or morphological composition or additional overload conditions [38].

Numerical FEM

Following the acquisition of the results from the two numerical methods, the comparison is depicted in Figure 9 alongside the experimental results. Figure 9a shows a similarity in the increase of the SIF concerning the crack size between the FEM simulation and the experimental results. This factor is related to the ratio between the crack size and the specimen width. However, In Figure 9b, this trend does not persist, possibly due to the relatively low values attained by ${\mathsf{\Delta }\mathrm{K}}_{\mathrm{II}}$. In FEM, results are computed concerning the immediately preceding interval, requiring the establishment of a new coordinate system at each evaluated point. During the final measured intervals, the crack exhibits its highest curvature about the initial measurement point. However, when measured from the immediately preceding point, the curvature is minimal. These results are reflected in the values recorded by FEM, although the simulation shows a trajectory similar to experiments [2].

Figure 9c displays the results of the equivalent SIF concerning the crack size. FEM results are higher compared to the experimental ones, as the values of ${\mathsf{\Delta }\mathrm{K}}_{\mathrm{I}}$ were also high. However, accelerated growth is observed after reaching 35% of the total crack size. The results associated with the SIF and crack size are similar in the numerical methods and experimental results. Figure 9d illustrates that the simulations yield a conservative result regarding the relationship between the crack size and the number of cycles. This is because the SIF values are close, and this difference is reflected in the growth rate, especially in the growth constants, as shown in Figure 9e. The experimental growth rate was slower; this is evident in the first reported value for the equivalent SIF, which is practically equal but experimentally corresponded to 109,000 cycles, while in the FEM simulation, it was reached at 18,100 cycles.

Numerical BEM

The results of the BEM numerical method are split into two groups: the results with the Paris rule (BEM-Paris) and Klesnil (BEM-Klesnil). In Figure 9a, the reported values of ΔK_I are higher than the experimentally reported ones; however, they exhibit a similar trend to the FEM values. Figure 9a shows that the reported values for ΔK_II are near zero. However, the experimental values measured the horizontal deformation of the specimen, while neither of the two numerical methods accurately estimated this value; in the case of the programs, none calculates it this way. Additionally, the values of ΔK_I and ΔK_II are equal in the case of BEM-Paris and BEM-Klesnil. In Figure 9c, the trend exhibits a behavior similar to Figure 9a because ΔK Asaro predominates over the influence of ΔK_I. In Figure 9d, in the case of BEM-Paris, the results are equal to those reported by FEM, which also employs the same growth model, and when determining the crack propagation, the Keq is the same.

In contrast, in the BEM-Klesnil combination, the crack propagation rate is slower at the beginning; this may reflect some plasticity ahead of the crack tip [39]. In this case, the reported values are close to the experimental ones. It is important to emphasize that for this figure, the implemented growth constants were those of the mode I reported in Table 3. Finally, in Figure 9e, where the growth rate is presented concerning the equivalent SIF, the results of BEM-Paris are the same as those of FEM for the reasons discussed earlier. A slow initial growth in BEM-Klesnil can be observed.

Additionally, some of the points overlap with the experimental ones. However, the general experimental growth continues to present a lower rate. In this case, new constants are proposed for mode I/II, for which a potential adjustment was made, obtaining the constants C and m reported in Table 5.

Simulations were rerun but implementing the constants reported in Table 5. In this case, the values of ΔK_I and ΔK_II are minimal, so we omitted these graphs for practicality. Figure 10a illustrates that all the results, except BEM-Klesnil, overlap. However, as the propagation progresses, the results converge as Region II develops. Figure 10b demonstrates that both numerical methods present acceptable predictions, even slightly conservative compared to the experimental results. If specific constants are known for mixed mode, it is possible to use the Paris rule. In this case, the Klesnil model predicts a conservative number of almost double cycles. This is because, in this case, the constants already consider slow growth. This model is appropriate if the specimen presents slow propagation and only the Paris constants are known for mode I.

5. Discussion

The experimental data acquired through DIC was processed and compared with numerical methods FEM and BEM see Figure 9a,b. The results unveiled a consistent trend of higher ${\mathrm{K}}_{\mathrm{I}}$ values in FEM, aligning with other studies [40]. However, it has been shown FEM cannot reproduce conditions such as fatigue-induced plasticity, crack roughness, debris, friction, crack flank interlocking, and the microstructural effect on the crack path [41]. Their impact on LEFM parameters can be seen when selecting the experimental data to work with. For example, choosing the full field displacement may capture the plasticity in the very near region potentially leading to higher K_I, or may yield negative K_I values if the fields are in compression. Alternatively, choosing points behind the crack tip may pick the effect of crack roughness and crack flank interlocking [42], which may shadow K_II. Of course, these scenarios are out of the scope, but they are worth mentioning.

Other researchers [43] have encountered similar results, addressing them by implementing a theta-method in FEM simulations to calculate ${\mathrm{K}}_{\mathrm{I}}$ values and computing the crack opening displacement factor (COD) using DIC data or the J integral to avoid the plasticity fields around the crack tip [42]. However, in our experience, the J integral implementation for mixed-mode fields has been a challenge. Similarly, in the study by [44], they resolved it through a hybrid experimental-computational approach. The FEM robustness and RPIM (Reduced-Order Proper Orthogonal Decomposition Integration Method) computations were experimentally validated, through Thermal Stress Analysis (TSA), with the achieved SIF solution. Such an approach appears suitable for cases of mixed-mode loading and could serve as a potential solution for this issue.

Regarding ${\mathrm{K}}_{\mathrm{II}}$, the simulations did not show a specific trend, potentially due to their low magnitudes. However, a consistent pattern emerged, showing higher values for both ${\mathrm{K}}_{\mathrm{I}}$ and ${\mathrm{K}}_{\mathrm{II}}$ in FEM [45], as reported in Figure 9d. For future studies, it is essential to recognize that achieving convergence between experimental and computational results may require the implementation of alternative algorithms.

The sample's thickness and the assumption of plain strain conditions influence the accuracy of computational predictions. Despite the sample's relatively small thickness (only 8.7 mm), the assumptions regarding plain strain conditions influenced the computational results. In the BEM simulations, the superficial calculation led to lower values than experimental findings [46], see Figure 8c. This discrepancy suggests that the idealized conditions assumed in the computational model may not fully represent the real-world behavior of the material.

Examining the relationship between crack size and the number of cycles, the discrepancy reaches 100% at its peak. This can be attributed to crack tip-induced plasticity, as the material is a low-carbon steel. According to the COD formulation, the plasticity mechanisms governing crack closure [47], are proportional to the SIF, but inversely proportional to r^0.5, significantly influencing the experimental growth slope [48] as shown in Figure 9e.

A significant finding was the decrease in growth constants for ${\mathrm{K}}_{\mathrm{Asaro}}$ relative to those presented in Table 3. This implies a potential reduction in mixed-mode Paris constants. This trend is consistent with findings from [49], where changes in constants were noted based on the severity of the applied mixed mode. However, a definitive trend was not identified, attributing the variation to shifts in force distribution throughout the body as the mixed mode becomes more pronounced. Nonetheless, the underlying cause of this phenomenon remains unclear, whether it arises from present failure modes or plasticity effects. Nevertheless, it is evident that alterations in loading conditions, such as changes in mode mixity, can lead to corresponding adjustments in growth constants, as highlighted before [50].

6. Conclusion

Three fatigue crack growth rules (Paris, Klesnil-Lukas, and Modified Forman-Newman) were tested with experimental data from the literature combined with three equivalent SIF models (Asaro, Pook, and Tanaka). It was found that the equivalent Stress Intensity Factor (SIF) model proposed by Asaro exhibits the lowest percentage error in predicting cycle differentials among the evaluated growth models, with a mean error within the range of 18.4 to 41.2 %. Asaro, Tanaka, and Pook's models exhibit close results, originating from maximum stress theoretical principles, resulting in conservative predictions Experimental ${\mathsf{\Delta }\mathrm{K}}_{\mathrm{I}}$ and ${\mathsf{\Delta }\mathrm{K}}_{\mathrm{II}}$ values compared to ASTM E647-23 are notably lower in the initial three intervals (where mostly mode I is present), leading to slower crack propagation. This phenomenon may be attributed to localized plasticity.

If specific constants are known for mixed mode, it is possible to use the Paris rule. When not available, the Klesnil fatigue crack growth model showed the smallest range of error in predicting cycle differentials, ranging from 18.4% to 22.1%. The total cycle count calculation using the Asaro equivalent SIF model showed slightly conservative behavior, with a satisfactory adjustment of 10%. The Klesnil growth model offered closer predictions due to its ability to model regions I and II of crack propagation. Observing the material's propagation threshold, SIF values indicated that the data reside at the entry limit of Region II.

DIC-acquired and processed experimental data were compared against commercial FEM and BEM numerical methods using SIF as a benchmark. The comparison revealed a higher ${\mathsf{\Delta }\mathrm{K}}_{\mathrm{I}}$ value in FEM, agreeing with other studies [37]. Differences were attributed to idealized conditions in the simulations, notably impacting ${\mathsf{\Delta }\mathrm{K}}_{\mathrm{II}}$ values without a specific trend. Notably, BEM's over-the-surface calculations led to lower values, as shown in Figure 9c.