This chapter presents a general background and the state of the art of numerical simulation and modeling of fretting phenomenon in terms of wear, fatigue and fracture. First, an introduction of fretting and its implications is exposed. Second, different methodologies for wear modeling and simulation are described and discussed. Afterwards, fatigue and fracture analysis approaches are revised. To that end, multiaxial fatigue parameters are introduced putting an emphasis on the physical basis of the fretting phenomena and the suitability of each model. On the other hand, the propagation phase based on linear elastic fracture mechanics (LEFM) via the finite element method (FEM) and the eXtended finite element method (X-FEM) analysis methods is presented and compared. Finally, different approaches and latest developments for fretting fatigue lifetime prediction are presented and discussed.
- numerical simulation
Fretting phenomena arises when two bodies are in contact subjected to relative movement of small amplitude (0–300 μm), producing damage on the contact surface . Since virtually all machines vibrate, fretting failure can occur in a variety of mechanical components (even the ones that are not intended to move), such as aircraft engine blade housings, ropes, flexible couplings, bearing housings and even orthopedic devices.
It has been reported that up to 50 variables might influence the magnitude and rate of the fretting process ; however, Drobomirski  identified the slip amplitude, the contact pressure and the coefficient of friction as the most influential ones. Different regimes can be determined depending on the slip amplitude, which is related to different failure types :
Stick regime: no slip occurs between the surfaces due to the accommodation of the displacement by elastic deformation. The damage in the surface is very low.
Partial-slip regime: the central zone of the contact interface is motionless or stuck while the outer one is sliding. In this case, there is a minimum wear, and fatigue failure occurs due to crack incubation and growth (fretting fatigue).
Gross-slip regime: the contact interface is in slip regime. The failure mainly occurs due to wear (fretting wear).
Depending on the magnitude of stresses, fretting can cause catastrophic failure of mechanical components. It is noteworthy mentioning that fretting fatigue may reduce the lifetime of a component by half or even more, in comparison to plain fatigue .
Despite the considerable progress made in the understanding of fretting fatigue over the last decades, it is still one of the modern issues for industrial machinery . Accordingly, there is an increasing interest in the use of the finite element method (FEM) to analyze fretting phenomena, since it provides data which currently cannot be obtained through experimental testing or analytical solutions. This chapter presents a general background and the state of the art of numerical simulation and modeling of fretting in terms of wear, fatigue and fracture.
2. Methodologies for wear modeling and simulation
Many wear models have been developed throughout history aimed to describe and predict the phenomena. Those models are extensively reviewed in the literature [7, 8, 9, 10] and can be broadly classified into two main categories, namely, (1) mechanistic models, based on material failure mechanism and (2) phenomenological models, based on contact mechanics.
The most popular and most used among them is the phenomenological Archard and Hirst model :
where represents the worn volume, the wear coefficient (obtained experimentally), the sliding distance and the contact force. It is noteworthy that the wear coefficient depends on the contact force and sliding amplitude, requiring its determination for each test condition.
On the other hand, Fouvry et al.  proposed a new model based on the energy dissipated on the contact surface:
where is the energy wear coefficient and is the dissipated energy of the
Fouvry highlighted the benefits of the energy-based wear model over the classically applied Archard’s approach. This new model requires a unique experimental campaign due to the independence of the energy wear coefficient to the contact force and sliding amplitude. Both models predict the same wear under gross slip condition, since the tangential force follows the equation: where
With the aim to solve wear modeling in more complex configurations and contact conditions, the FEM wear modeling has appeared as a good alternative in the last decade. Contact surface evolves progressively during wear phenomena, which generates an evolution in the stress state. This being so, the strategies for wear simulations consist of updating the geometry at discrete time increments, based on semi-empiric wear models.
McColl et al.  introduced a methodology for numerical wear simulation based on the local implementation of the previously mentioned Archard model (the so-called modified Archard’s model):
where is the incremental wear depth, is the local Archard coefficient, is the local contact pressure and is the relative slip distance increment.
As shown in Figure 1, the simulation consists of an iterative process where the local Archard’s model is solved through contact pressures and sliding distributions obtained by the numerical simulation at discrete time increments. The wear value obtained at each increment is then used to update the worn geometry at the end of the cycle and the process is repeated till a predefined maximum sliding distance is reached.
Since contact problems are nonlinear, the computational demand is considerable, especially in three spatial dimensions where this problem is highly magnified. Among the strategies to optimize wear simulations, the use of the cycle jumping technique should be highlighted. This approach allows to speed up the wear simulation under the assumption that wear remains constant for a small number of cycles. Therefore, a cycle jumping factor multiplying the incremental wear allows using one computational wear cycle to model the material removal of actual cycles .
A further improvement on the computational time was presented by Madge et al.  who programed the spatial adjustment of the contact nodes through the user defined subroutine UMESHMOTION (available on the commercial FE code Abaqus FEA). This subroutine works in an adaptive meshing constrain framework in order to adapt the mesh to the evolving geometries.
Among the several benefits of using UMESHMOUTION subroutine, it should be highlighted that the updating is done incrementally through the fretting cycle, providing more stable results comparing to the updates done at the end of the cycle. Larger cycle jumps can therefore be used, decreasing significantly the computation time. However, the subroutine gives access to the pressure data of only one of the bodies, avoiding the possibility to compute wear on both parts. Cruzado et al. [15, 16, 17] overcame this limitation by transferring the available contact pressure data to the other part by interpolation techniques.
It should be highlighted that recent publications [18, 19] proposed the use of the energy wear approach instead of the Archard’s local equation. Following the previously explained framework proposed by Madge, the energy equation is computed locally as:
As mentioned earlier, Fouvry’s model shows the independency of contact force and sliding amplitude being more versatile than the commonly used Archard’s equation.
3. Fatigue analysis approaches
Material fatigue refers to a progressive degradation of a material caused by loading and unloading cycles. The stress fluctuations suffered over time weakens or breaks the material even at stresses lower than the yielding value. Accordingly, a lot of effort has been directed at developing fatigue-life prediction models.
Fatigue is characterized with a high scatter of the lifetime. Probabilistic approaches are recently arising in the literature to address this problem [20, 21]. However, the majority of the models currently used analyze fatigue in a deterministic way, that is, a structure fails if a given parameter reaches a critical value.
Nowadays, a variety of different approaches for fatigue life prediction exist, such as approaches based on multiaxial fatigue criteria, damage mechanics or micromechanics, which are extensively reviewed in literature [22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33]. The present review focuses on the most widely used classical methodologies, that is, multiaxial fatigue criteria.
Multiaxial fatigue criteria reduce the multiaxial stresses (usually computed by FEM analysis) to an equivalent uniaxial stress state. This way, the results can be compared to an experimental fitting curve obtained from uniaxial fatigue data. A crucial step when selecting a multiaxial criterion is to check whether the simplification from multiaxial stress state to an equivalent uniaxial stress state is acceptable or valid. This task is not simple and requires the detailed study of the evolution of stresses and strains along the loading cycle.
The book published by Socie and Marquis  presents a wide and detailed study about the principal multiaxial parameters, also known as Fatigue Indicator Parameters (FIPs). Those parameters can be broadly classified into three groups: strain-based, stress-based and energy-based FIPs. Strain-based FIPs [34, 35] are generally related with Low-Cycle Fatigue (LCF) where plastic deformation may be predominant. Stress-based FIPs [36, 37] are associated with High-Cycle Fatigue (HCF), where the stresses usually remain in the elastic domain. Finally, energy-based models [38, 39, 40] relate the product of stresses and strains to quantify fatigue life, which generally are applicable to both LCF and HCF regime.
Additionally, fatigue can be categorized into proportional (fixed principal directions along a loading cycle) and nonproportional loading (rotation of the principal directions along a loading cycle). Figure 2 shows the evolution of the stresses at the contact surface along fretting fatigue cycle at three different stages of the loading cycle (maximum, mean and minimum). It can be observed that the normal stress and the shear stress fluctuate along the cycle while the stress remains unaltered. Consequently, principal directions rotate along the cycle generating nonproportional stresses under proportional remote loading. Therefore, the equivalent stress and strain approaches such as Von Mises criterion developed for proportional loading are not applicable in fretting since the problem is highly nonlinear and nonproportional.
For these complex stresses or loading states, other approaches such as the critical plane method are more suitable [41, 42, 43]. The critical plane method has been developed from the experimental observation of nucleation and crack growth under multiaxial loading. The critical plane models include the dominant parameters that govern the type of crack initiation and propagation. An adequate model must be one that estimates correctly both fatigue life and the dominant failure plane. However, several failure modes exist, and there is not a unique parameter that suits all.
A great deal of critical plane based FIPs have been used in the literature to assess fretting fatigue life [44, 45, 46, 47, 48, 49]. Nonetheless, the most popular parameters are the energetic criteria known as Fatemi-Socie (
Under shear loading condition, crack lip surfaces generate frictional forces that reduce stresses at the crack tip, thus increasing fatigue life. However, tensile stresses and strains will separate the crack surfaces, reducing the friction forces. The energetic FIP
where is the maximum range of shear strain on any plane, is the maximum normal stress in that particular plane, is the material yield stress, is a material dependent factor.
Vázquez et al.  recently compared both parameters for the analysis of the initial crack path in cylindrical fretting contact, concluding that the
It should be mentioned that these parameters give a local life prediction and seek to find the hot spot to give the minimum life estimation. However, when high stress gradient events appear, for example, fretting case, an over-estimation of crack nucleation is predicted at the hot spot. Consequently, a nonlocal approach such as the Theory of Critical Distances (TCD)  used extensively in notched fatigue is recommended .
4. Fracture modeling and simulation
Fracture mechanics is the field of mechanics concerned with the study of structures integrity in the presence of cracks. Within this field, there are several approaches, such as the linear elastic fracture mechanics (LEFM), the non linear fracture mechanics (NLFM) or the elasto-plastic fracture mechanics (EPFM) [53, 54, 55]. This chapter focuses on the most widely used one, the LEFM approach.
From a fully elastic point of view, Williams  presented an eigen function expansion method that provides a framework for the description of the stress state near a crack-tip. For each cracked configuration, a sequence of coefficients depending on the geometry and load describes the stress state with respect to the radius of circumference and angle (see Figure 3).
Irwin identified the first physically valid term in this infinite series, the field :
where is the Stress Intensity Factor (SIF) for each of the fracture mode. The use of SIFs assumes that the singular stresses dominate the stress field near the crack front, thus neglecting higher order terms of the Williams series. It can be easily seen that the stress field shows a singularity when
As far as fatigue crack growth behavior is concerned, this is usually described by the relationship between the crack growth length increase per cycle () as a function of the SIF range (). The typical log-log plot of crack growth behavior is shown schematically in Figure 4, which may be divided into three regimes:
Regime I: the near threshold region—below the threshold value of SIF () cracks will not propagate.
Regime II: stable crack growth region—the curve is essentially linear (log-log scale), widely known as the Paris regime.
Regime III: rapid unstable crack growth leading to catastrophic fracture.
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 (), 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 . A comprehensive analysis of crack propagation velocity models in fretting fatigue was carried out by Navarro et al.  who analyzed nine of the fatigue crack growth models for an aluminum alloy Al7075. This study concluded that the Paris law  and the modified SIF model  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  or the minimum of the strain energy density factor
4.1. Stress intensity factor computation
Many of the initial analytical approaches initially developed for SIFs calculation [64, 65, 66, 67] have been now outdated by the versatility offered by numerical methods. In this regard, one of the main methods is the interaction integral  through the equivalent domain integral  (nowadays implemented in commercial finite element codes such as Abaqus FEA). The interaction integral is an extension of the well-known
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  showed that the finite elements could be enriched with additional functions to represent a given function, Möes et al.  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 () 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:
where , 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:
where is the set of all nodes in the mesh, and are the sets of the enriched nodes, is the shape function, is the classical degree of freedom of the FEM and and 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)  or the fast version called the Fast Marching Method (FMM) and (2) explicit methods, such as geometric predicates  or other approaches . 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].
5. Fretting fatigue numerical simulation
Over the years, a significant number of different methodologies for fretting fatigue life prediction have been presented. First models were purely analytical, and numerical models appeared as computers evolved. Nowadays the preferred approaches can be classified into hybrid (combining analytical and numerical methods) or fully numerical methods.
Fretting fatigue life estimation has become an object of interest in the literature of the field. In general, the study of fretting is divided into two stages, the crack initiation () and its subsequent propagation (). Therefore, the sum of the two stages gives a total life prediction:
In those situations in which the initiation phase dominates over propagation, some authors propose to estimate life considering only the initiation stage . Conversely, for cases where crack growth stage dominates the components’ life, the initiation phase is sometimes neglected . Nonetheless, most current approaches combine both phases in order to provide a total life prediction, and this chapter is therefore centered on those combined models.
The number of cycles to initiate a crack is typically obtained with a FIP (the reader is referred to Section 3 for multiaxial criterion details). Crack propagation is subsequently considered using different fatigue crack growth laws, which require the definition of an initial crack length.
It is noteworthy that the development of a combined initiation-propagation model involves two critical steps: (1) the selection of the stress-deformation analysis location and (2) the definition of the initial crack length. On the one hand, due to the high stress gradient that characterizes fretting (see Figure 5(a)), the predicted crack nucleation cycles will increase as the location is further from the contact (see Figure 5(b). On the other hand, longer cracks grow faster, leading to a lower number of propagation cycles (see Figure 5(c)).
Among the different methodologies presented for fretting fatigue lifetime prediction, the simplest approach is to perform the multiaxial analysis at the hot spot (i.e., at the surface), which results in a conservative crack initiation prediction . Once the condition of crack initiation is reached, a predefined crack length is introduced in the numerical model [44, 79]. Alternatively, other authors perform the multiaxial analysis at a certain distance from the surface instead of at the hot spot, obtaining less conservative life estimations [14, 78, 80, 81]. However, the selected distance is usually arbitrary and there is no consensus in the literature.
A further improvement was introduced by Navarro et al.  who presented a nonarbitrary criteria. In this methodology, both the initiation and propagation phases are computed separately for a set of points at different depths, considering the distance from the surface to the point under analysis as the initial crack length. Fretting fatigue life estimation is then calculated for each point () and the minimum of those predictions is considered the estimated life. This methodology combines the FEM for crack initiation analysis and analytical models for crack propagation phase. It has been shown that the use of the X-FEM for crack propagation analysis following the same framework gives a better correlation in fatigue life prediction .
More recently, the use of nonlocal methods such as the TCD (the reader is referred to Section 3) was suggested by some authors in order to define the location of the multiaxial analysis and to set the initial crack length [49, 83].
Despite fretting wear and fretting fatigue being frequently linked, most of the approaches, such as the ones mentioned earlier, neglect wear phenomena. A further improvement is therefore to consider the wear simulation in order to study the effect of wear on fatigue life.
In this regard, one of the most prominent work was presented by Madge et al.  In this work, the 2D numerical simulation of material removal process is performed first under the UMESHMOTION framework (the reader is referred to the Section 2). Then, the multiaxial fatigue analysis coupled with a damage accumulation framework is carried out in order to account for the effect of wear on fatigue lifetime. Finally, the propagation phase is analyzed via submodeling technique, which allows to transfer the stress state of the contact surface from global wear model to crack submodel. It should be noted that this approach does not account for the explicit interaction between the fretting contact and the crack. This study showed that the linkage of wear modeling with fatigue analysis is a key factor to successfully predict the increase in life in the gross sliding regime.
A further improvement of the previous approach was recently published by the present authors  by combining wear, fatigue and fracture phenomena in a single numerical model. Figure 6 shows the flow chart of the coupled numerical analysis approach.
The simulation algorithm is divided into two blocks corresponding to initiation and propagation stages, respectively. The first stage runs under the FEM framework, where the accumulated damage is computed during wear simulation iteratively. This process is repeated until the accumulated damage value reaches the value of 1 (), that is, up to concluding the initiation phase. Afterwards, the propagation is computed through X-FEM code published by Giner et al. . In this stage, the crack propagation is calculated during wear simulation iteratively up to reaching the failure criteria (). The presented model allows the study of the fretting phenomena as a whole, allowing to model explicitly the interaction between the fretting contact and the crack. It is expected that further analysis using this framework will allow a better understanding of the synergies between wear, fatigue and fracture phenomena.
6. Summary and concluding remarks
The complex nature of the physical phenomenon of fretting involves the mixture of different fields of engineering such as fatigue, fracture and tribology. This chapter has introduced the state of the art of the currently available modeling and simulation methods to analyze fretting phenomenon. The benefits and drawbacks of each reviewed technique have been highlighted. Finally, a numerical architecture of coupled wear, fatigue and fracture methodology has been introduced, which allows to analyze the fretting phenomenon as a whole.
Among the main open challenges identified, a predominant role is taken by the need to compare the increasing amount of methodologies to assess fretting lifetime, and it is likely to become the subject of considerable debate in the research community in the near future.