**1. Introduction**

Since the end of the Second World War, the failure of materials under stresses even lower than the yield stresses has gained significant attention [1]. Ensuring the stability of critical structures while establishing their safe working condition is vital. In most industries, the accurate estimation of both crack path and fatigue life are crucial in terms of reliability. In various applications, such as aerospace manufacturing and the aviation industry, experimental studies are necessary for fatigue analysis, but, because of high costs, precise computational methods are required for crack propagation analysis to predict the direction of crack growth and fatigue lifetime in both static and dynamic loading conditions [2]. The failure is related to (a) the presence of flaws such as interfaces and cracks, and (b) the nature of fluctuating loads. Cracks tend to initiate and propagate when subjected to fluctuating loads until a point when the structure does not bear the load that contributes to complete failure. These cracks are considered fatigue cracks and the expected life is one of the major parameters to determine the safety of the structure. This is computed by adding the number of loading cycles required to nucleate the fatigue cracks that lead to failure. The calculation of the growth rate of the cracks is usually based on the relation between the range of the stress intensity factors (SIFs) and the cracks' geometry.

The extended finite element method (XFEM) is an alternate way to predict the SIFs using computational methods. In general, the initiation and propagation of cracks must be associated to the SIFs in a complicated state [3–6]. The extended finite element method proposed by Belytschko and Black in 1999 [7] has been widely used in recent studies. It is based on the standard finite element framework and uses a special displacement feature to allow discontinuities to occur, overcoming the need to re-mesh continuously throughout the crack tip expansion process. To evaluate the SIFs, XFEM was used to perform crack

**Citation:** Fageehi, Y.A. Fatigue Crack Growth Analysis with Extended Finite Element for 3D Linear Elastic Material. *Metals* **2021**, *11*, 397. https://doi.org/10.3390/met11030397

Academic Editor: Ricardo Branco

Received: 2 February 2021 Accepted: 19 February 2021 Published: 1 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

growth analysis without updating the mesh [8]. Extensive work was undertaken to develop efficient models to evaluate the fatigue crack growth (FCG) and fatigue life in order to overcome fatigue failures. There are several proposed experimental models, but they become prohibitive both in terms of cost and time. An effective way to reduce the laboratory work, time, and costs is to incorporate a simulation methodology that involves numerical analysis and use of the ANSYS APDL.19.2 extended finite element method. Many fatigue crack problems identified in the literature to date use different computational approaches in simulating simple and complex two- and three-dimensional geometries [2,9–15].

This work used the ANSYS APDL 19.2 XFEM to precisely predict the mixed-mode stress intensity factors along with the associated fatigue life for a modified four-point bending beam and a cracked plate with three holes. In particular, three methods have been widely used to illustrate the fatigue assessment of materials: the fracture mechanics method developed by Paris and Erdogan [16], the strain-life method independently proposed by Coffin [17], and the stress-life (SN) method proposed by Wöhler [18]. The first approach, by which the crack tip can be described separately by the SIFs, was employed in this study for predicting fatigue life. The second approach is suitable in the lower cycle fatigue range, whereas the third, SN approach estimates the time spent to initiate and grow a crack until the component breaks into parts, which requires stress results from a linear static analysis. The main motivation for this work was to make a significant contribution to the use of ANSYS as an alternative tool for simulating fatigue crack propagation problems during mixed-mode loading and to monitor the trajectory of crack growth in cases of the presence of holes in the geometry.

### **2. Mixed-Mode Fatigue Life Evaluation Procedure Using ANSYS**

ANSYS (version 19.2, Ansys, Inc., Canonsburg, PA, USA), can model three kinds of cracks: arbitrary, semi-elliptical, and pre-meshed. The pre-meshed crack method requires a crack front employed by the Smart Crack growth analysis tool, whereby the stress intensity factor is the criterion of failure. The node sets that were rendered were distributed to the front, top, and bottom of the crack. The latest feature presented in ANSYS is the Smart Crack growth mesh-based tetrahedron, which adds the pre-meshed crack requirement after completion, enabling the selection of the type of crack growth. The sphere of influence process can be used in refining the mesh around the crack tip about the geometric edge that passes through the thickness. The geometric regions to be described are the crack tip, the crack top, and the bottom surfaces of the crack; each of these regions is associated with a node set to be used for analysis. The ANSYS software considered mixed-mode loading where the maximum circumferential stress criterion was implemented. The following are the formulas for the direction angle of crack propagation in ANSYS [19,20]:

$$\theta = \cos^{-1}\left(\frac{3K\_{II}^2 + K\_I\sqrt{K\_I^2 + 8K\_{II}^2}}{K\_I^2 + 9K\_{II}^2}\right) \tag{1}$$

where *KI* = max *KI* during cyclic loading and *KII* = max *KII* during cyclic loading.

In ANSYS Mechanical APDL 19.2, by using XFEM, crack growth simulation was restricted to region II of the typical fatigue crack growth graph, which can be represented as:

$$\frac{da}{dN} = \mathbb{C}(\Delta K\_{c\eta})^m \tag{2}$$

From Equation (2), fatigue lifecycles may be predicted for crack increments as:

$$\int\_0^{\Delta a} \frac{da}{\mathcal{C}(\Delta \mathcal{K}\_{eq})^m} = \int\_0^{\Delta N} dN = \Delta N \tag{3}$$

The equivalent range of the stress intensity factor formula is determined as follows [19]:

$$
\Delta K\_{eq} = \frac{1}{2} \cos \frac{\theta}{2} [\left(\Delta K\_I (1 + \cos \theta)\right) - 3\Delta K\_{II} \sin \theta] \tag{4}
$$

where Δ*KI* = the stress intensity factor range in mode I loading and Δ*KI I* = the stress intensity factor range in mode II loading.

Based on numerical analysis, there are many methods formulated for evaluating the stress intensity factors. The interaction integral technique is usually the most accurate method that has the ability of estimating *KI* and *KII* separately. ANSYS proposes two methods for computing SIFs: the displacement extrapolation method (DEM) and the interaction integral method (IIM). The second method was adopted because it is numerically easier to implement and has better precision with fewer mesh requirements. This approach uses the domain integral approach [21] where an auxiliary field is used to separate *KI* from *KII*, as this ability is missing in the domain integral itself. The energy release rate in terms of the mixed-mode stress intensity factors *KI*, *KII*, and *KIII* was proposed as [21,22]:

$$J(s) = \frac{K\_I^2 + K\_{II}^2}{E^\*} + \frac{1+\nu}{E} K\_{III}^2 \quad E^\* = \begin{bmatrix} \frac{E}{(1-\nu^2)} \text{ Plane strain} \\ \frac{E}{E} \text{ Plane stress} \end{bmatrix} \tag{5}$$

The superimposed state Equation (5) becomes:

> *I*

$$\begin{split} I^S(s) &= \frac{1}{E^\*} \left[ (K\_I + K\_I^{aux})^2 + (K\_{II} + K\_{II}^{aux})^2 \right] + \frac{1+\nu}{E} (K\_{III} + K\_{III}^{aux})^2 \\ &= J(s) + J^{aux}(s) + I(s) \\\\ I(s) &= \frac{1}{\Gamma^\*} (2K\_I K\_I^{aux} + 2K\_{II} K\_{II}^{aux}) + \frac{1+\nu}{\Gamma} (2K\_{III} K\_{III}^{aux}) \end{split} \tag{6}$$

*E*

Here, superscript  $(\mathbf{S})$  denotes the superimposed state;  $\mathbf{J}(\mathbf{s})$  is the domain integral for the actual state;  $\mathbf{J}^{\text{max}}(\mathbf{s})$  is the domain integral for the auxiliary state; and  $\mathbf{I}(\mathbf{s})$  is an integral with interacting actual and auxiliary terms.

 *I I*

By setting *Kaux I*= 1 and *Kaux II*= *Kaux III*= 0, Equation (6) yields:

$$K\_I = \frac{E^\*}{2} I(s) \tag{7}$$

 *III*

By setting *Kaux I I* = 1 and *Kaux I* = *Kaux III* = 0, and selecting *Kaux III* = 1, *Kaux I* = *Kaux I I* = 0 gives the relationships between *KII* and *KIII*:

$$K\_{II} = \frac{E^\*}{2} I(s) \tag{8}$$

$$\mathbb{K}\_{III} = \mu \operatorname{I}(\mathbf{s}) \tag{9}$$

where *E* and *μ* are the modulus of elasticity and the modulus of rigidity, respectively.
