**4.1. Stress intensity factor computation**

• Regime II: stable crack growth region—the curve is essentially linear (log-log scale), widely

Regarding the simulation and modeling of crack growth behavior, the computation of the SIFs have been a priority in fracture mechanics, which has given rise to a great diversity of techniques (discussed in Section 4.1). In cases where the SIF range overcomes the threshold value (Δ*K* > Δ *K*th), the velocity and direction of the crack growth should be computed.

Crack propagation rate is usually described by means of a phenomenological law of the type d*a*/d*N* = *f*(Δ*K*). A comprehensive analysis of crack propagation velocity models in fretting fatigue was carried out by Navarro et al. [58] who analyzed nine of the fatigue crack growth models for an aluminum alloy Al7075. This study concluded that the Paris law [59] and the modified SIF model [60] were the most suitable ones for the experimental campaign carried out. Concerning crack orientation criteria, they are generally based on the analysis of the stress and strain fields. The suitability of each criterion mainly depends on the evolution of the stresses and strains along a loading cycle. It should be noted that fretting fatigue usually induces friction between crack faces prone to slip motion during the loading cycle. The problem is therefore nonproportional, and the classical orientation criteria for proportional loading such as the maximum circumferential stress [61] or the minimum of the strain energy density factor *S* [62] among others are not applicable. Giner et al. [63] analyzed the suitability of the nonproportional loading criteria available in the literature for fretting fatigue problem, concluding that the prediction of the crack path observed in the complete contact experiments did not present a good agreement with the models available. Therefore, they developed the criterion of the minimum shear stress range, which is a generalization for nonproportional loading of the so-called criterion of local symmetry well established for proportional loading. The numerical results obtained by this new criterion were in good agreement with the experimental observation.

• Regime III: rapid unstable crack growth leading to catastrophic fracture.

**Figure 4.** Schematic representation of a typical fatigue crack growth rate curve plotted in log-log scale.

known as the Paris regime.

202 Contact and Fracture Mechanics

Many of the initial analytical approaches initially developed for SIFs calculation [64–67] have been now outdated by the versatility offered by numerical methods. In this regard, one of the main methods is the interaction integral [68] through the equivalent domain integral [69] (nowadays implemented in commercial finite element codes such as Abaqus FEA). The interaction integral is an extension of the well-known *J* integral proposed by Rice [70], which is capable of extracting SIFs for each mode separately (*K*<sup>I</sup> , *K*II and *K*III).

The computation of the interaction integral requires first to compute the stresses and strains by means of the FEM. However, the study of the singular problem through the FEM presents several drawbacks. On the one hand, in the classical formulation of the FEM the element edges need to conform to the crack boundaries, which require the use of cumbersome meshing techniques. On the other hand, the shape functions employed are generally of low-order polynomials, leading to the use of very refined mesh to compute reliable stresses and strains around the crack tip. Additionally, fatigue crack requires remeshing techniques to conform to the new crack boundaries.

Once Melenk and Babuška [71] showed that the finite elements could be enriched with additional functions to represent a given function, Möes et al. [72] proposed the eXtended finite element method (X-FEM) as a solution to overcome these issues.

In the X-FEM it is not necessary to have a mesh that conforms to the crack geometry, thus the finite element mesh is independent of the crack shape. To this end, the FEM model is enriched with additional degrees of freedom. On the one hand, the Heaviside function (*H*(**x**) = ±1) is used to introduce discontinuity along the crack faces. On the other hand, the FEM model is additionally enriched with the asymptotic function derived by Irwin, and, therefore, can reproduce the singular behavior of LEFM by the following expression:

$$\{F^\downarrow(\mathbf{x})\} \equiv \sqrt{r} \begin{bmatrix} \sin(\frac{\theta}{2}), \cos(\frac{\theta}{2}), \sin(\frac{\theta}{2})\sin(\theta), \cos(\frac{\theta}{2})\sin(\theta) \end{bmatrix} \tag{8}$$

where *r*, *θ* represents the polar coordinates defined with respect to the local reference system at the crack tip (see **Figure 3**).

Therefore, the classical FEM displacement definition is enriched to obtain the X-FEM approximation to the displacement field, defined as:

$$\mathbf{u}^{\text{X-FEM}}(\mathbf{x}) = \sum\_{i \in \mathbf{u}} \mathbf{u}\_i \, N\_i(\mathbf{x}) + \sum\_{i \in \mathbf{u}} \mathbf{a}\_i \, N\_i(\mathbf{x}) H(\mathbf{x}) + \sum\_{i \in \mathbf{u}} N\_i(\mathbf{x}) \Big(\sum\_{i=1}^4 \mathbf{b}\_i^\dagger F^i(\mathbf{x})\Big) \tag{9}$$

where *I* is the set of all nodes in the mesh, *L* and *K* are the sets of the enriched nodes, *Ni* is the shape function, **u***<sup>i</sup>* is the classical degree of freedom of the FEM and **a***<sup>i</sup>* and **b***<sup>i</sup> l* are respectively the degrees of freedom of the enriched nodes.

Another important aspect of the X-FEM is the geometrical representation of the evolving cracks and the definition of the elements and nodes to be enriched. There are several methods to perform this task and can be divided into two groups: (1) implicit methods, such as the Level Set Method (LSM) [73] or the fast version called the Fast Marching Method (FMM) and (2) explicit methods, such as geometric predicates [74] or other approaches [72]. The suitability of each of the method depends on whether the problem is two or three-dimensional. Although the benefits of the LSM or FMM are not overwhelming in 2D, the use of implicit methods becomes necessary in 3D since the explicit representation can be quite difficult to discretize [75, 76].
