Fretting Fatigue Behavior under Tension–Bending Mixed-Mode Loading

Abstract

The mixed-mode loading fretting fatigue caused by the complex geometry of components and combinations of boundary conditions is a common failure mechanism in engineering components, which can dramatically reduce fatigue life. In this paper, a cylinder-on-flat numerical model was established to investigate tension–bending mixed-mode fretting fatigue. The finite element method in conjunction with two criteria, plane parameters McDiarmid (MD) and Smith–Watson–Topper (SWT), were used to evaluate the effects of mode angle, oblique loading, and stiffness ratio on the contact width, the maximum equivalent stress of the specimen, the surface stress, the fretting damage initiation location, and the extent of the damage initiation. The results indicate that the extent of fretting damage increases with the mode angle, and the characterization parameters are sensitive to smaller mode angles. The contact width, peak surface stress, maximum damage parameters, and damage initiation location can be effectively adjusted by the stiffness ratio. The findings may provide insights into fretting fatigue behavior under complex loading conditions, potentially contributing to enhanced structural safety and reliability for tension–bending mixed-mode loading.

1. Introduction

Fretting can be defined as a small amplitude slip caused by cyclic loading. It generally occurs between two nominally fastened bodies, causing surface damage initiation as well as material fatigue fracture [1]. Fretting fatigue is a major source of failure in many engineering components such as bolted joints [2,3,4], riveted joints [5,6], wire ropes [7,8,9], and dovetail joints [1,9]. It is generally accepted that a slip band is produced at the edge of the contact surface in the fretting fatigue process. The crack can be initiated and grow due to the high gradient of cyclic stress there, which can severely reduce the fatigue life of the structure [10,11,12].

Researchers have adopted three single-mode models according to the loading conditions of the structure to simulate fretting fatigue phenomena, i.e., tension fretting fatigue [13,14,15], torsional fretting fatigue [16], and bending fretting fatigue [17,18,19]. These three single modes of fretting were set up in a series of experiments to study the effects of multiple variables (such as loading [20,21], material properties [22,23], and environmental conditions [24,25]) on fretting damage. The results of the experiment provide guidance for enhancing the fretting fatigue resistance of materials. With an in-depth understanding of the damage evolution of fretting fatigue, numerical analysis has also been developed. The finite element method with multiaxial fatigue theory and fracture mechanics is combined to analyze and predict the initiation and propagation of cracks and fatigue life in fretting fatigue [26,27,28,29], and the results show a good correlation with experimental data.

However, due to the complex geometry and combinations of the boundary conditions of the components, the mixed-mode fretting fatigue that is subjected to mixed loading conditions is more common in engineering applications. Few studies have been found regarding mixed-mode loading conditions in fretting damage. Zhu et al. [27] developed a novel experimental apparatus to simulate combined fretting conditions. The apparatus integrated tangential fretting and radial fretting loading modes, allowing the investigation of dynamic behavior, contact conditions, and damage processes under combined loads by varying the maximum load and inclination angles. The results revealed that material damage was governed by the competition between the two fretting loading modes, closely related to experimental parameters such as normal force, inclination angle, and material properties. Contact surface observations showed asymmetric fretting scars, with some tests exhibiting eccentricity between the center of the adhered region and the overall scar pattern. Wang et al. [26] explored the influence of the torsion angle on the multiaxial fretting fatigue behavior of steel wires using a custom-built tension–torsion combined multiaxial fretting fatigue test rig. They found that increasing the torsion angle led to an increase in wire deflection angle, wear scar size, friction coefficient, and maximum crack depth, indicating that larger torsional loads would exacerbate fretting fatigue damage in wires. Wang et al. [30] conducted tests on steel wires using a custom-built tension–torsion fretting corrosion fatigue test apparatus to investigate the effects of tensile stress ratio and amplitude on fretting fatigue behavior. The results demonstrated that increasing the tensile stress ratio and amplitude led to an overall increase in dissipated energy, friction coefficient, maximum wear depth, crack propagation rate, lifetime impact, fractal dimension, and surface complexity of the wires, suggesting that higher tensile stresses accelerate combined fretting fatigue damage. Peng et al. [31] studied the bi-directional composite fretting wear behavior of Zr-4 alloy tubes under different phase differences, simulating the normal and radial composite loading conditions experienced by fuel assemblies in pressurized water reactors. Utilizing a self-developed multimode fretting wear test rig, they conducted bidirectional combined fretting wear tests at room temperature and 300 °C. The results revealed that the temperature variation and phase difference changes significantly influenced fretting wear characteristics, with a 90° phase difference exhibiting the most severe wear, and the 300 °C tests showing intensified wear and damage compared to room temperature.

It is noteworthy that the aforementioned studies mainly focused on fretting damage under normal–radial or tension–torsion combined loading conditions. However, in engineering applications, components often need to withstand the effects of combined tensile-bending mixed-mode loads. For example, in aerospace engines, blade roots experience combined tension and bending stresses, leading to fretting fatigue at the contact interface with the disk, which can cause blade fracture and engine failure [32,33]. Similarly, railway axles, subjected to tensile and bending stresses from the vehicle weight and track irregularities, can develop fatigue cracks at the axle-bearing contact area, potentially resulting in axle breakage and derailments [11]. Wind turbine main shafts face similar issues, with varying wind loads causing fretting fatigue at the shaft-bearing interface, impacting turbine operation and requiring maintenance [34]. Additionally, in the bridge cables of suspension and cable-stayed bridges, combined tensile and bending stresses from deck loads and wind can induce fretting fatigue at the cable–anchor interface, posing risks to bridge integrity [35]. Despite the importance of these conditions, there are few reports investigating fretting fatigue behavior under tension–bending mixed-mode loading conditions. Under this complex loading, relative slip in the contact regions can trigger fretting fatigue crack initiation, potentially leading to component failure. This not only accelerates fatigue crack nucleation and propagation but also influences contact stress distributions through the evolution of wear scar morphology, thereby altering crack propagation paths. Consequently, investigating tension–bending mixed-mode fretting fatigue is crucial for accurately predicting and assessing the service life of critical components, enhancing their reliability and safety in practical engineering applications.

In this article, the finite element method in conjunction with two critical plane parameters is employed to investigate fretting fatigue under tension–bending mixed-mode loading conditions. The organization of this paper is as follows: Section 2 establishes a cylindrical pad-on-flat contact model to study tension–bending mixed-mode fretting fatigue, where dimensional analysis is also performed to reduce the number of independent variables. Section 3 presents the numerical model and its details. The results and discussion regarding the effects of oblique load, stiffness ratio, and mode angle are provided in Section 4. Finally, the main conclusions of the present work are drawn in Section 5.

2. Problem Formulation

2.1. Model Description

The numerical model of cylindrical pad-on-flat contact (see Figure 1) is established to evaluate the fretting fatigue behavior of the material under tension–bending mixed-mode loading. Regarding the load setting, two constant normal loads, P, are applied at the top of each cylindrical pad, and the cyclic oblique load, F, is applied to the right end of the flat specimen to simulate fretting cycles. The stress ratio of the cyclic load F is 0.1, and it is applied in a direction that has an angle (defined as the mode angle θ) with the x-axis. Here, θ = 0° represents tension-mode fretting fatigue, θ = 90° indicates bending-mode fretting fatigue, and 0° < θ < 90° signifies a tension–bending mixed-mode fretting fatigue condition. The left end of the specimen is fixed, and the horizontal displacement of the pad is zero. Both the specimen and the pad have the same thickness (the z-directional size), t, and the coefficient of friction between them is μ. The magnitude of F and θ should ensure that the tangential stress, q(x), distributed along the contact surface is less than the frictional force provided by the normal stress, p(x), i.e., q(x) < μp(x), and thus the contact remains in a partial slip regime. Under the action of a small relative oscillatory motion, the fretting fatigue phenomenon occurs at the edge of the contact surface.

2.2. Dimensional Analysis

Dimensional analysis [28] is a powerful tool to analyze the relationship between variables in science and engineering. This method is conducive to simplifying the obtainment of the physical phenomena either from calculations or experiments by reducing the number of independent variables in functions. For the problem investigated in this paper, the most important variable stress, σ_ij, which characterizes fretting damage, can be derived by other variables using general expression f as follows:

(1)${\sigma }_{ij}=f(P,F,\theta ;E_1,E_2,{\nu }_1,{\nu }_2,\mu ;R,t,x,y),$

where E₁, ν₁, E₂, and ν₂ are the elastic modulus and Poisson’s ratio of the specimen and the pad, respectively. R is the radius of the cylindrical pad. Among the twelve governing parameters, two of them, namely E₁ and t, have been selected as the independent variables, and other variables are dependent variables whose dimensions can be given by:

(2)$\begin{array}{l}\left[{\sigma }_{ij}\right]=\left[E_2\right]=\left[E_1\right]\\ \left[P\right]=\left[F\right]=\left[E_1\right]{\left[t\right]}^2\\ \left[R\right]=[x]=\left[y\right]=\left[t\right]\\ \left[{\nu }_1\right]=\left[{\nu }_2\right]=\left[\mu \right]=\left[\theta \right]=1.\end{array}$

By applying the Π-theorem to dimensional analysis, a dimensionless stress function expression is obtained:

(3)$\frac{{\sigma }_{ij}}{E_1}=\prod (\frac{P}{E_1{t}^2},\frac{F}{E_1{t}^2},\theta ;\frac{E_2}{E_1},{\nu }_1,{\nu }_2,\mu ;\frac{R}{t},\frac{x}{t},\frac{y}{t}).$

For the problem studied in this paper, the focus is on the tension–bending mixed-mode loading, thus the parameters P, ν₁, ν₂, μ, and R are set to constant. Consequently, the dimensionless stress function can be simplified as follows:

(4)$\frac{{\sigma }_{ij}}{E_1}=\prod (\frac{F}{E_1{t}^2},\theta ;\frac{E_2}{E_1};\frac{x}{t},\frac{y}{t}).$

Based on the above dimensional analysis, the number of governing parameters is reduced, and three important dimensionless variables, $\frac{F}{E_1{t}^2}$, θ, and $\frac{E_2}{E_1}$, are obtained. Defining $\frac{F}{E_1{t}^2}$ as λ and $\frac{E_2}{E_1}$ as γ, λ is referred to as the dimensional oblique load, γ is the stiffness ratio, where $\gamma <1$ represents a stiff specimen pressed by a compliant pad, and $\gamma >1$ signifies that a compliant specimen is pressed by a stiff pad. $\gamma =1$ means the specimen matches the pad.

3. Numerical Model and Method

3.1. Finite Element Model

In order to investigate the fretting behavior of the cylindrical pad-on-flat specimen under tension–bending mixed-mode loading, a 2D plane strain finite element model (shown in Figure 2) is developed by employing ABAQUS software (version 6.14) along with MATLAB standard code. ABAQUS is used for pre-processing and solving the FEM model, and MATLAB code is used to perform the critical plane method to predict the fretting damage initiation.

As is indicated in Figure 2, the left side of the specimen is set to a fixed constraint, and the pad is restricted to move just in a vertical direction. Normal load P is applied to the top of both cylindrical pads using a multi-point constraint (MPC). Oblique force F applied to the right end of the specimen is resolved into the tension load, Fcosθ, and the bending load, Fsinθ, to simulate the tension–bending mixed-mode loading.

Four nodes, plane strain, and a reduced integration quadrilateral element (CPE4R) are used for the analysis. The master–slave algorithm is adopted to calculate the contact between the cylindrical pad and the specimen. Al2024-T351 with elastic modulus E₁ = 74.1 GPa is used for the specimen, and the constant dimensionless parameters of the model are shown in Table 1. The material and geometric data of the model are derived from Ref. [13]. The cylindrical pad has the same Poisson’s ratio as the specimen. The maximum Von Mises equivalent stress in all the cases is lower than the yield stress, so the elastic constitutive relation is implemented to simulate the mechanical properties of the material.

The load sequence is shown in Figure 3, where the loads are applied in three steps, and the time increments in each step are small enough. Normal pressure, P, is applied at the first step to establish contact and remains constant in the subsequent steps. The bending load, Fsinθ, and the tension load, Fcosθ, each with a stress ratio of 0.1, are simultaneously loading in the second step and unloading in the third step. According to the numerical observation, the equivalent stress of the specimen under tension–bending mixed-mode loading increased to the maximum at step 2. The node with the maximum stress is found at the leading edge of the contact interface (the contact edge, which is towards the cyclic load, as shown in the upper dashed rectangle of Figure 2), which is consistent with the fretting damage initiation location observed in many experiments [9,13,18,21]. Therefore, the nodes on the upper contact interface of the specimen are defined as the path for stress analysis and calculating damage parameters.

Table 2 shows the sixteen groups of numerical experimental data designed for simulating the tension–bending mixed-mode fretting fatigue. Based on the dimensional analysis, the experimental data are tabulated into three categories, λ, θ, and γ. Four levels of oblique loads, seven stiffness ratios, and seven mode angles are designed to investigate the effect of governing parameters on mixed-mode fretting fatigue.

Due to the high stress gradient in the potential zone of fretting damage, an appropriate mesh size should be selected to obtain accurate stress in this zone while considering the calculation efficiency. The convergence study is carried out to obtain a suitable mesh size under the loading case of Exp. 4. Figure 4 shows the maximum Von Mises equivalent stress, σ_eff,max, at each step for different mesh sizes. It reveals that the stress is not sensitive to mesh size under the single pressure load (the blue line in Figure 4) but is sensitive in the loading and unloading steps. As the mesh density increases, the value of σ_eff,max first increases, and then tends to converge. It is noteworthy that when the mesh size changes from 15 μm to 10 μm, a significant variation in the σ_eff,max of step 3 is observed. This is because the refinement reaches a critical threshold, where the mesh begins to more accurately capture the finer details of the stress gradients. This phenomenon is common in mesh convergence studies, where certain mesh sizes reveal a more accurate representation of the physical phenomena being modeled. When the mesh size is 8 μm, the maximum equivalent stress relative error between two adjacent mesh density models drops to 1%. Figure 5 shows the relationship between the shear stress (in step 2) and the mesh sizes. It indicates that the accurate stress is obtained with the refined mesh. When the mesh size is 8 μm, it is sufficient to capture the precise stress distribution in the contact area. Considering the results from Figure 4 and Figure 5 together, a mesh size of 8 μm is deemed to have achieved mesh convergence and is suitable for the mixed-mode fretting problem simulation.

Although the convergence of the numerical model is confirmed by Figure 4 and Figure 5, it is necessary to compare the FEM solutions with the analytical solutions to validate the model. Since no analytical solutions for tension–bending mixed-mode fretting contact have been reported in the literature, the FEM results (the mesh size of the model is 8 μm) can only be compared with the analytical solution for the contact problem derived by Hertz [29], as shown in Figure 6. The diagram shows the normal contact stress, p(x), of the FEM and the analytical solution have almost coincident profiles, and both reach the peak value at the contact center. The percentage error of peak contact stress is less than 2%, which is because the half-space assumption is not satisfied by the boundary conditions of the FEM model.

3.2. Fretting Damage Initiation Prediction Method

High-gradient multiaxial and non-proportional stress produced by fretting fatigue will cause the fatigue life of the component to be lower than the fatigue limit of the material, which presents a great challenge to predict the fretting damage initiation. Several multiaxial fatigue criteria such as the critical plane approach, continuum damage mechanics approach, and stress invariance approach have been developed [36,37,38,39] and successfully applied to the numerical prediction of fretting fatigue initiation. Among them, the critical plane approach, which is based on historical stress–strain parameters and agrees well with many experimental observations, has been extensively used to predict the crack initiation location, angle, and life of fretting fatigue. The key to the critical plane approach is to determine the maximum fatigue damage plane of the material based on the selected fatigue parameters. The maximum value of the damage parameter in conjunction with some material constants is used to calculate the fatigue life and takes the peak value point as the crack initiation position. Different categories of damage parameters are suitable for different fatigue mechanisms (high-cycle fatigue or low-cycle fatigue) and different failure modes, and the classification of the categories is based on different combinations of mechanical parameters.

Tension–bending mixed-mode fretting fatigue is numerically assessed by two critical plane parameters—McDiarmid (MD) [40] and Smith–Watson–Topper (SWT) [41]—to predict the fretting damage initiation location and to estimate the extent of the fretting damage. McDiarmid used two shear fatigue limits: one for the case where cracks grow parallel to the surface (Case A), and another for the case where cracks grow towards the surface (Case B). In our model, the torsional fatigue limit ${\tau }_{f-1}$ was used, following the assumption of the original MD model in which cracks initiate parallel to the surface (Case A), which is consistent with the mechanism of fretting fatigue crack initiation. The stress-based criterion parameter MD is suitable for materials whose failure is dominated by tensile-mode cracking, while the energy-based criterion parameter SWT is suitable for materials whose failure is dominated by shear-mode cracking. The value of these two fatigue parameters can be calculated following Equations (5) and (6):

(5)$MD=\frac{\Delta {\tau }_{\max }}{2}+\left(\frac{{\tau }_{f-1}}{2{\sigma }_u}\right){\sigma }_{n,\max }$

(6)$SWT={\left({\sigma }_{\max }\frac{\Delta \epsilon }{2}\right)}_{\max }$

where $\frac{\Delta {\tau }_{\max }}{2}$ is the maximum shear stress amplitude, and ${\sigma }_{n,\max }$ is the maximum normal stress perpendicular to the direction of maximum shear stress amplitude. Material constants ${\tau }_{f-1}$ and ${\sigma }_u$ are the torsional fatigue limit and ultimate tensile strength, respectively, and their values are derived from the literature [13] for Al2024-T351. The ${\sigma }_{\max }$ is the maximum normal stress on a plane of the material point, and $\Delta \epsilon$ is the maximum strain range on the same plane. The predicted fretting damage initiation life can be calculated by the damage parameter combined with Basquin’s equation, and the extent of the component under fretting fatigue is inversely related to the value of the fretting damage parameter [1,12,13].

4. Results and Discussion

In this section, the results are explained first for the contact state of the fretting contact zone. Then, the maximum Von Mises equivalent stress and contact surface stress of the specimen are analyzed. Finally, the critical plane method combining two multiaxial damage parameters, MD and SWT, is used to analyze the effect of the governing parameters on the fretting damage initiation. The results of Exp. 1 to Exp. 7 are used to analyze the effect of the mode angle θ. The results of Exp. 8 to Exp. 10 in combination with Exp. 4 are used to study the effect of oblique loads λ, and the results of Exp. 11 to Exp. 16 in combination with Exp. 4 are used to analyze the effect of stiffness ratio γ. The analysis of these two parameters is based on the mode angle of 45°, and this is the most typical load setting for tension–bending mixed-mode fretting fatigue. Furthermore, to present the computational results more distinctly, certain data are displayed in normalized form, with normalization based on the computational results from Exp. 4.

4.1. Contact Analyses

The experimental results of many researchers have indicated that the fretting fatigue crack initiation at the contact edge and subsequently rapid propagation is the main reason for reducing the fatigue life of structures [13,18,42]. Table 3 shows that the stick–slip ratio of the contact zone decreased with the increase in the dimensionless oblique load, while the contact width remained almost constant. The reasons can be explained as follows. Since the oblique loading, λ, increased, its normal component, Fsinθ, and tangential component, Fcosθ, increased as well. The tangential traction of the contact surface, which causes the slip zone to expand, increases with the increase in Fcosθ, and the magnitude of Fsinθ is much smaller than the normal force, P; the resultant force in the normal direction is slightly reduced, so the contact width is almost unchanged. From Table 4, it is evident that with the increase in stiffness ratio, the contact width decreases. This trend aligns with the findings of fretting fatigue experiments conducted by Lee et al. [43] using different mating materials and is consistent with the Hertz contact theory [29].

Table 5 shows the effect of the mode angle on the contact zone, which reveals that the mode angle has no obvious effect on the contact width, while it has the most significant effect in adjusting the stick–slip ratio of the contact zone. The ratio of the stick–slip width decreases with the increase in the mode angle, and the decreased gradient is large when the mode angle is between 0° and 45°. The contact condition of θ = 75° is close to the result of θ = 90°. It means that the slip regime is mainly affected by the mode angle under the tension-dominated (0° < θ < 45°) tension–bending mixed-mode fretting fatigue.

4.2. Stress Analysis

The maximum Von Mises equivalent stress of the specimen σ_eff,max, the distribution of normal stress p(x), the horizontal component of stress σ_xx, and the tangential stress q(x) on the contact surface are examined in this section to understand the fretting fatigue behavior of the material under tension–bending mixed-mode loading. The σ_eff,max produced in any loading case is less than the yield stress, σ_y, for the material (σ_y =310 Mpa [13]), as shown in Figure 7, and it demonstrates that the linear elastic constructive relation is reasonable. Figure 7 also shows that the σ_eff,max increases with the increase in the dimensionless stiffness ratio, γ, and oblique load, λ, but the slope of the curve related to the mode angle θ decreases rapidly. The curve of σ_eff,max vs. θ grows rapidly when θ < 45° but slowly increases when θ > 45°, as shown in Figure 7c, which indicates that the risk of material yielding is the highest for the bending-dominated (45° < θ < 90°) mixed-mode loading and it can be conservative and simple to estimate the bending-dominated fretting fatigue by the bending fretting fatigue strength.

Figure 8 shows the dimensionless normal stress, $p(x)/E_1$, distribution on the contact surface. As shown in Figure 8a,c, the peak value and the curve profile of the normal stress are almost the same for different oblique loads and mode angles, but the location of the peak normal stress (i.e., contact center) moves towards the trailing edge of the contact as the two variables increase. These phenomena indicate that the peaks of the contact stress and the contact width are dominated by the normal load, P, but the position of the contact center and the slip regime of the contact zone will be affected by the mode angle and oblique load. Unlike λ and θ, the peak value of contact pressure increases with the increase in stiffness ratio γ, but the position of the contact center remains constant (shown in Figure 8b), which indicates that the variable, γ, can function in adjusting the preload between two contact bodies (such as bolts and other fasteners).

Figure 9 demonstrates the effects of θ, λ, and γ on the dimensionless surface stress component ${\sigma }_{xx}/E_1$. As shown in Figure 9, ${\sigma }_{xx}/E_1$ is compressive at most of the contact surface but transforms into tensile stress at the leading edge, and the fretting fatigue crack propagation is dominantly by the high gradient of tensile stress [39]. As expected, when a larger oblique load is applied, the peak of tensile stress is higher, and its position moves toward the contact center, as shown in Figure 9a. Different stiffness ratios on the surface stress component, ${\sigma }_{xx}/E_1$, are shown in Figure 9b, in which the peak tensile stress increases as the stiffness ratio increases, and the changing γ from 0.4 to 1 has a more powerful adjustment capability on ${\sigma }_{xx}/E_1$ than that from 1 to 2.5. This indicates that when the stiff specimen is pressed by the compliant pad, the crack propagation driving force is effectively reduced by adjusting the stiffness ratio, γ. With the increase in the mode angle, the peak of tensile stress increases, and its position moves toward the contact center, as shown in Figure 9c. This trend is obvious in the tension-dominant mode, but it is not pronounced in the bending-dominant mode, and the stress curves almost coincide with θ = 75° and θ = 90°, which indicates that the mode angle should be controlled in the range of 0° ≤ θ ≤ 45° to reduce the driving force of fretting fatigue crack propagation.

Figure 10 shows the distribution of dimensionless shear stress on the contact surface. As shown in Figure 10a, as the oblique load increases, the shear stress gradient in the vicinity of the contact edge is reduced, where the peaks increase and positions move closer to the contact center, which is related to the stick–slip ratio. If λ continues to be applied, the fretting regime will change from partial slip to gross slip, and the mode of fretting damage will also change from fretting fatigue to fretting wear. The two peaks of the shear stress increase with the increase in the stiffness ratio as well, but the distance between them decreases, as shown in Figure 10b, which can be explained by the narrowing of the contact width. When the stiffness of the pad increases, the variable amplitude of the peak at the leading edge is smaller than that at the trailing edge, and the position of the maximum shear stress is shifted from the leading edge to the trailing edge. It is seen in Figure 10c that the $q(x)/E_1$ vs. θ curve shows that peaks appear on both sides of the contact center, and the peak of the leading edge is slightly larger than that of the trailing edge. The shear stress increases as the mode angle increases, and its two peaks move towards the trailing edge, but the stress profiles of θ = 75° and θ = 90° almost coincide.

4.3. Damage Initiation Prediction

4.3.1. Damage Initiation Location

The prediction of the damage initiation location of the MD and SWT parameters for all cases is shown in Figure 11, Figure 12 and Figure 13 and Table 6, Table 7 and Table 8. The graphs of the two parameters are not completely consistent, but the general profiles are similar. For all the numerical experimental results, the two critical plane parameters are smoothly distributed in the contact center region, while for the contact edge, the damage parameters increase suddenly, which indicates that the damage is highly concentrated in this region. Although the damage parameter peaks appear at both edges of the contact, the damage parameter value at the leading edge is the highest, which indicates that the damage is likely to be initiated at the leading edge of the contact. This phenomenon is consistent with many experimental results [13,15,18,20,38,42,44], and the high gradient damage parameter that is in the vicinity of the leading edge is due to the high stress concentration that is caused by the cylindrical geometry of the fretting pad.

As shown in Figure 11 and Table 6, with the increase in oblique load, the damage initiation location moves closer to the contact center. This is consistent with the findings of Han et al.’s dovetail joint fretting fatigue experiments, where the crack initiation location shifted away from the loading end as the fatigue load increased [45]. Figure 12 and Table 7 indicate that a higher stiffness ratio causes the fretting damage initiation location to move toward the trailing edge of the contact. This shift is due to the reduction in contact width with the increasing stiffness ratio, resulting in the crack initiation location being closer to the contact center.

For all the cases, compared to the MD parameter, the predicted position of the SWT parameter is slightly closer to the leading edge. From Figure 13 and Table 8, it can be seen that the damage initiation positions predicted by the two fretting damage parameters, MD and SWT, are both on the leading edge of the contact surface, and the positions are very close. With the increase in the mode angle, θ, the damage initiation position is closer to the contact center, but this tendency is less obvious in the bending-dominated mode, and the damage initiation positions predicted by the two parameters are almost the same when θ = 75° and θ = 90°.

4.3.2. The Extent of Damage Initiation

During the fretting damage initiation process, we care not only about the damage initiation location but also the extent of the damage initiation. The fretting damage initiation extent of the specimen can be represented by the maximum damage parameter, and its effectiveness has been verified in many experiments [9,37,46,47].

In order to clearly analyze the relationship between the control variable and the extent of fretting damage initiation, the normalized maximum damage parameters were separately extracted and plotted, as shown in Figure 14. Despite the two damage parameters being based on different mechanical parameters, their prediction trends given by different dimensionless variables for the extent of fretting damage initiation are generally consistent. For the relationship between the inclined load and the maximum damage parameter, the result of λ = 1.594 × 10⁻⁴ (Exp. 9) was used for normalization, as shown in Figure 14a. The positive correlation between the oblique load λ and the damage initiation extent is shown in the figure, which reveals that it is necessary to properly control the range of cyclic load to improve the fretting fatigue damage initiation life of the component. Figure 14b illustrates that the maximum damage parameters increase with increasing stiffness ratio, and the relationship between them is almost centrally symmetric. This shows that the stiffness ratio is a very effective parameter for adjusting the tension–bending mixed-mode fretting fatigue initiation life of the component, and a reasonable design of the stiffness ratio between two contact bodies will improve the fretting fatigue resistance of the structure. Figure 14c shows that with the increase of the mode angle, θ, the damage initiation extent of the specimen increases, but the damage extent for pure bending mode (θ = 90°) is lower than that of θ = 75°. This phenomenon demonstrates that the fretting damage extent of the specimen under tension–bending mixed-mode loading is more serious than that under single-tension or single-bending mode loading. It can also be seen from Figure 14c that the extent of fretting damage is more sensitive at smaller mode angles (tension dominant mode), which is consistent with the bridge-type pad fretting experimental results [44].

5. Conclusions

This paper presents and numerically analyzes the tension–bending mixed-mode fretting fatigue behavior of a cylindrical pad-on-flat model. The effects of oblique loading, stiffness ratio, and mode angle on the contact width, maximum equivalent stress in the specimen, surface stress, fretting damage initiation location, and extent of damage initiation are analyzed. The finite element method, in combination with the critical plane parameter approach, has been employed to assess fretting damage initiation. Our findings can inform the design and manufacturing processes of fretting-resistant components, enabling engineers to better understand fatigue behavior under mixed-mode fretting loading conditions. By identifying critical stress locations and dominant modes of fatigue loading, these insights facilitate targeted design improvements to enhance the structural integrity and durability of aerospace and mechanical components. The key results can be summarized as follows:

• As the dimensionless oblique load increases, the stick–slip ratio decreases, and the contact center and the location of the damage initiation move toward the contact trailing edge. Both the maximum equivalent stress and the extent of fretting damage are positively related to the oblique load.

• The contact width, the peak surface stress, and the extent of fretting damage increase with the increase in the stiffness ratio. The position of the maximum shear stress will be transformed from the leading edge to the trailing edge when the stiffness ratio is greater than unity.

• When 0° ≤ θ ≤ 45°, the stick–slip ratio, surface stress, maximum equivalent stress, and the extent of damage are very sensitive to the mode angle. However, when 45° ≤ θ ≤ 90°, the parameter gradient is not obvious.

• The damage parameters are concentrated at the contact edge, but the leading edge is more severe. When bending-dominant (45° ≤ θ ≤ 90°), the maximum value of the damage parameters is similar, and the damage extent for the pure bending mode (θ = 90°) is lower than that of θ = 75°.

Author Contributions

X.Z.: Conceptualization; data curation; formal analysis; investigation; methodology; resources; software; validation; visualization; and writing—original draft. X.C.: Conceptualization; funding acquisition; project administration; supervision; and writing—review and editing. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All data generated and analyzed in this research are included in this published article.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Schematic of cylindrical pad-on-flat tension–bending mixed-mode fretting fatigue model.

Figure 2. Sketch of the finite element model.

Figure 3. Loading sequence of the finite element model.

Figure 4. The maximum equivalent stress, σeff,max, at each step for different mesh sizes.

Figure 5. Dimensionless shear stress, [Forumla omitted. See PDF.], at the contact surface for different mesh sizes.

Figure 6. Comparison between FEM result and Hertz analytical solution of the contact pressure.

Figure 7. Finite element model variations in maximum equivalent stress σeff,max with different control variables: (a) dimensionless oblique load λ, (b) stiffness ratio γ, and (c) mode angle θ.

Figure 8. Variations in dimensionless normal stress [Forumla omitted. See PDF.] with different control variables: (a) dimensionless oblique load λ, (b) stiffness ratio γ, and (c) mode angle θ.

Figure 9. Variations in dimensionless surface stress component [Forumla omitted. See PDF.] with different control variables: (a) dimensionless oblique load λ, (b) stiffness ratio γ, and (c) mode angle θ.

Figure 10. Variations in dimensionless shear stress q(x)/E1 with different control variables: (a) dimensionless oblique load λ, (b) stiffness ratio γ, and (c) mode angle θ.

Figure 11. Variations in damage parameters with dimensionless oblique load λ: (a) MD and (b) SWT.

Figure 12. Variations in damage parameters with stiffness ratio γ: (a) MD and (b) SWT.

Figure 13. Variations in damage parameters with mode angle θ: (a) MD and (b) SWT.

Figure 14. Variations in normalized maximum damage parameters with different control variables: (a) dimensionless oblique load λ, (b) stiffness ratio γ, and (c) mode angle θ.

Table 1

Constant dimensionless parameters of the contact model.

P/E1t2	ν 1	ν 2	μ	R/t
$5.405\times {10}^{-4}$	0.33	0.33	0.65	14.016

Table 2

Numerical experimental data in this paper.

Exp. No.	λ (10−4)	θ	γ	Exp. No.	λ (10−4)	θ	γ
1	1.063	0°	1	9	1.594	45°	1
2	1.063	15°	1	10	1.860	45°	1
3	1.063	30°	1	11	1.063	45°	0.4
4	1.063	45°	1	12	1.063	45°	0.5
5	1.063	60°	1	13	1.063	45°	0.67
6	1.063	75°	1	14	1.063	45°	1.5
7	1.063	90°	1	15	1.063	45°	2
8	1.328	45°	1	16	1.063	45°	2.5

Table 3

Contact sizes for different oblique loads (θ = 45°, γ = 1).

λ (10−4)	1.063	1.328	1.594	1.860
Contact width (x/t)	0.259	0.259	0.260	0.260
Stick zone width (x/t)	0.162	0.143	0.132	0.132
Ratio of stick–slip width	1.680	1.231	1.038	1.039

Table 4

Contact sizes for different stiffness ratios (θ = 45°, $\lambda =1.063\times {10}^{-4}$).

	γ = 0.4	γ = 0.5	γ = 0.67	γ = 1	γ = 1.5	γ = 2	γ = 2.5
Contact width (x/t)	0.341	0.316	0.289	0.259	0.236	0.225	0.217
Stick zone width (x/t)	0.232	0.211	0.188	0.162	0.142	0.131	0.123
Ratio of stick–slip width	2.126	2.022	1.867	1.680	1.515	1.394	1.312

Table 5

Contact sizes for different mode angles ($\lambda =1.063\times {10}^{-4}$, γ = 1).

	θ = 0°	θ = 15°	θ = 30°	θ = 45°	θ = 60°	θ = 75°	θ = 90°
Contact width (x/t)	0.259	0.259	0.259	0.259	0.259	0.259	0.259
Stick zone width (x/t)	0.243	0.214	0.186	0.162	0.143	0.132	0.132
Ratio of stick–slip width	15.510	4.785	2.570	1.680	1.231	1.038	1.039

Table 6

Damage initiation locations for different oblique loads (θ = 45°, γ = 1).

λ (10−4)	1.063	1.328	1.594	1.860
MDmax location (x/t)	0.109	0.104	0.099	0.093
SWTmax location (x/t)	0.113	0.109	0.105	0.100

Table 7

Damage initiation locations for different stiffness ratios (θ = 45°, $\lambda =1.063\times {10}^{-4}$).

	γ = 0.4	γ = 0.5	γ = 0.67	γ = 1	γ = 1.5	γ = 2	γ = 2.5
MDmax location (x/t)	0.142	0.131	0.121	0.109	0.101	0.097	0.095
SWTmax location (x/t)	0.145	0.135	0.124	0.113	0.105	0.101	0.098

Table 8

Damage initiation locations for different mode angles ($\lambda =1.063\times {10}^{-4}$, γ = 1).

	θ = 0°	θ = 15°	θ = 30°	θ = 45°	θ = 60°	θ = 75°	θ = 90°
MDmax location (x/t)	0.127	0.121	0.114	0.109	0.105	0.102	0.102
SWTmax location (x/t)	0.113	0.123	0.118	0.113	0.110	0.107	0.107

References

1. Hu, C.; Wei, D.S.; Wang, Y.R.; Shi, L. Experimental and Numerical Study of Fretting Fatigue in Dovetail Assembly Using a Total Life Prediction Model. Eng. Fract. Mech.; 2019; 205, pp. 301-318. [DOI: https://dx.doi.org/10.1016/j.engfracmech.2018.08.001]

2. Benhamena, A.; Amrouche, A.; Talha, A.; Benseddiq, N. Effect of Contact Forces on Fretting Fatigue Behavior of Bolted Plates: Numerical and Experimental Analysis. Tribol. Int.; 2012; 48, pp. 237-245. [DOI: https://dx.doi.org/10.1016/j.triboint.2011.12.008]

3. Juoksukangas, J.; Lehtovaara, A.; Mäntylä, A. Experimental and Numerical Investigation of Fretting Fatigue Behavior in Bolted Joints. Tribol. Int.; 2016; 103, pp. 440-448. [DOI: https://dx.doi.org/10.1016/j.triboint.2016.07.021]

4. Jiménez-Peña, C.; Talemi, R.H.; Rossi, B.; Debruyne, D. Investigations on the Fretting Fatigue Failure Mechanism of Bolted Joints in High Strength Steel Subjected to Different Levels of Pre-Tension. Tribol. Int.; 2017; 108, pp. 128-140. [DOI: https://dx.doi.org/10.1016/j.triboint.2016.11.014]

5. Han, L.; Chrysanthou, A.; O’Sullivan, J.M. Fretting Behaviour of Self-Piercing Riveted Aluminium Alloy Joints under Different Interfacial Conditions. Mater. Des.; 2006; 27, pp. 200-208. [DOI: https://dx.doi.org/10.1016/j.matdes.2004.10.014]

6. Zhao, L.; He, X.C.; Xing, B.Y.; Lu, Y.; Gu, F.S.; Ball, A. Influence of Sheet Thickness on Fatigue Behavior and Fretting of Self-Piercing Riveted Joints in Aluminum Alloy 5052. Mater. Des.; 2015; 87, pp. 1010-1017. [DOI: https://dx.doi.org/10.1016/j.matdes.2015.08.121]

7. Cruzado, A.; Urchegui, M.A.; Gómez, X. Finite Element Modeling and Experimental Validation of Fretting Wear Scars in Thin Steel Wires. Wear; 2012; 289, pp. 26-38. [DOI: https://dx.doi.org/10.1016/j.wear.2012.04.018]

8. Cruzado, A.; Urchegui, M.A.; Gómez, X. Finite Element Modeling of Fretting Wear Scars in the Thin Steel Wires: Application in Crossed Cylinder Arrangements. Wear; 2014; 318, pp. 98-105. [DOI: https://dx.doi.org/10.1016/j.wear.2014.06.019]

9. Shi, L.; Wei, D.S.; Wang, Y.R.; Tian, A.M.; Li, D. An Investigation of Fretting Fatigue in a Circular Arc Dovetail Assembly. Int. J. Fatigue; 2016; 82, pp. 226-237. [DOI: https://dx.doi.org/10.1016/j.ijfatigue.2015.07.025]

10. Juuma, T. Torsional Fretting Fatigue Strength of a Shrink-Fitted Shaft with a Grooved Hub. Tribol. Int.; 2000; 33, pp. 537-543. [DOI: https://dx.doi.org/10.1016/S0301-679x(00)00102-X]

11. Gürer, G.; Gür, C.H. Failure Analysis of Fretting Fatigue Initiation and Growth on Railway Axle Press-Fits. Eng. Fail. Anal.; 2018; 84, pp. 151-166. [DOI: https://dx.doi.org/10.1016/j.engfailanal.2017.06.054]

12. Li, X.; Yang, J.W.; Li, M.H.; Zuo, Z.X. An Investigation on Fretting Fatigue Mechanism under Complex Cyclic Loading Conditions. Int. J. Fatigue; 2016; 88, pp. 227-235. [DOI: https://dx.doi.org/10.1016/j.ijfatigue.2016.03.030]

13. Szolwinski, M.P.; Farris, T.N. Observation, Analysis and Prediction of Fretting Fatigue in 2024-T351 Aluminum Alloy. Wear; 1998; 221, pp. 24-36. [DOI: https://dx.doi.org/10.1016/S0043-1648(98)00264-6]

14. Walyekar, A.A.; Leonard, B.D.; Sadeghi, F.; Jalalahmadi, B.; Bolander, N. An Experimental Study and Fatigue Damage Model for Fretting Fatigue. Tribol. Int.; 2014; 79, pp. 183-196. [DOI: https://dx.doi.org/10.1016/j.triboint.2014.06.006]

15. Jin, O.; Mall, S. Effects of Slip on Fretting Behavior: Experiments and Analyses. Wear; 2004; 256, pp. 671-684. [DOI: https://dx.doi.org/10.1016/S0043-1648(03)00510-6]

16. Peng, J.F.; Jin, X.; Xu, Z.B.; Zhang, J.; Cai, Z.B.; Luo, Z.P.; Zhu, M.H. Study on the Damage Evolution of Torsional Fretting Fatigue in a 7075 Aluminum Alloy. Wear; 2018; 402, pp. 160-168. [DOI: https://dx.doi.org/10.1016/j.wear.2018.02.008]

17. Peng, J.F.; Liu, J.H.; Cai, Z.B.; Shen, M.X.; Song, C.; Zhu, M.H. Study on Bending Fretting Fatigue Damages of 7075 Aluminum Alloy. Tribol. Int.; 2013; 59, pp. 38-46. [DOI: https://dx.doi.org/10.1016/j.triboint.2012.06.016]

18. Ebara, R.; Fujimura, M. Fretting Fatigue Behaviour of Ti-6Al-4V Alloy under Plane Bending Stress and Contact Stress. Tribol. Int.; 2006; 39, pp. 1181-1186. [DOI: https://dx.doi.org/10.1016/j.triboint.2006.02.007]

19. Peng, J.F.; Zhu, M.H.; Cai, Z.B.; Liu, J.H.; Zuo, K.C.; Song, C.; Wang, W.J. On the Damage Mechanisms of Bending Fretting Fatigue. Tribol. Int.; 2014; 76, pp. 133-141. [DOI: https://dx.doi.org/10.1016/j.triboint.2013.12.018]

20. Abbasi, F.; Majzoobi, G.H. Effect of out-of-Phase Loading on Fretting Fatigue Response of Al7075-T6 under Cyclic Normal Loading Using a New Testing Apparatus. Eng. Fract. Mech.; 2018; 188, pp. 93-111. [DOI: https://dx.doi.org/10.1016/j.engfracmech.2017.08.010]

21. Massingham, M.; Irving, P.E. The Effect of Variable Amplitude Loading on Stress Distribution within a Cylindrical Contact Subjected to Fretting Fatigue. Tribol. Int.; 2006; 39, pp. 1084-1091. [DOI: https://dx.doi.org/10.1016/j.triboint.2006.02.025]

22. Han, Q.N.; Rui, S.S.; Qiu, W.H.; Ma, X.F.; Su, Y.; Cui, H.T.; Zhang, H.J.; Shi, H.J. Crystal Orientation Effect on Fretting Fatigue Induced Geometrically Necessary Dislocation Distribution in Ni-Based Single-Crystal Superalloys. Acta Mater.; 2019; 179, pp. 129-141. [DOI: https://dx.doi.org/10.1016/j.actamat.2019.08.035]

23. Su, Y.; Han, Q.N.; Zhang, C.C.; Shi, H.J.; Niu, L.S.; Deng, G.J.; Rui, S.S. Effects of Secondary Orientation and Temperature on the Fretting Fatigue Behaviors of Ni-Based Single Crystal Superalloys. Tribol. Int.; 2019; 130, pp. 9-18. [DOI: https://dx.doi.org/10.1016/j.triboint.2018.09.006]

24. Liao, J.P.; Wu, X.Q.; Tan, J.B.; Tang, L.C.; Qian, H.; Xie, Y.C. Fretting Corrosion Fatigue of Alloy 690 in High-Temperature Pure Water. Corros. Sci.; 2018; 133, pp. 423-431. [DOI: https://dx.doi.org/10.1016/j.corsci.2018.01.037]

25. Majzoobi, G.H.; Hojjati, R.; Soori, M. Fretting Fatigue Behavior of Al7075-T6 at Sub-Zero Temperature. Tribol. Int.; 2011; 44, pp. 1443-1451. [DOI: https://dx.doi.org/10.1016/j.triboint.2011.03.021]

26. Wang, X.R.; Wang, D.G.; Zhang, D.K.; Ge, S.R.; Araújo, J.A. Effect of Torsion Angle on Tension-Torsion Multiaxial Fretting Fatigue Behaviors of Steel Wires. Int. J. Fatigue; 2018; 106, pp. 159-164. [DOI: https://dx.doi.org/10.1016/j.ijfatigue.2017.09.021]

27. Zhu, M.H.; Zhou, Z.R.; Kapsa, P.; Vincent, L. An Experimental Investigation on Composite Fretting Mode. Tribol. Int.; 2001; 34, pp. 733-738. [DOI: https://dx.doi.org/10.1016/S0301-679x(01)00061-5]

28. Barenblatt, G.I. Scaling, Self-Similarity, and Intermediate Asymptotics: Dimensional Analysis and Intermediate Asymptotics; Cambridge Texts in Applied Mathematics Cambridge University Press: Cambridge, UK, 1996.

29. Johnson, K.L. One Hundred Years of Hertz Contact. Proc. Inst. Mech. Eng.; 1982; 196, pp. 363-378. [DOI: https://dx.doi.org/10.1243/PIME_PROC_1982_196_039_02]

30. Wang, D.; Zhu, Z.; Song, D. Effects of Tensile Stress Ratio and Amplitude on Tension-Torsion Fretting-Corrosion-Fatigue Behaviors of Non-Perpendicularly Crossed Steel Wires. Eng. Fail. Anal.; 2020; 117, 104839. [DOI: https://dx.doi.org/10.1016/j.engfailanal.2020.104839]

31. Peng, J.; Ding, S.; He, Z.; Li, B.; He, J.; Liu, J.; Zhu, M. Study on Bi-Directional Composite Fretting Wear Characteristics of Zr-4 Alloy Tube with Different Phase Differences under Temperature Conditions. Tribol. Int.; 2024; 191, 109139. [DOI: https://dx.doi.org/10.1016/j.triboint.2023.109139]

32. Golden, P.J.; Nicholas, T. The Effect of Angle on Dovetail Fretting Experiments in Ti-6Al-4V. Fatigue Fract. Eng. Mater. Struct.; 2005; 28, pp. 1169-1175. [DOI: https://dx.doi.org/10.1111/j.1460-2695.2005.00956.x]

33. Qu, Z.; Liu, K.C.; Wang, B.Z.; Chen, Z.Y. Fretting Fatigue Experiment and Finite Element Analysis for Dovetail Specimen at High Temperature. Appl. Sci.; 2021; 11, 9913. [DOI: https://dx.doi.org/10.3390/app11219913]

34. Rauert, T.; Herrmann, J.; Dalhoff, P.; Sander, M. Fretting Fatigue Induced Surface Cracks under Shrink Fitted Main Bearings in wind Turbine Rotor Shafts. Procedia Struct. Integr.; 2016; 2, pp. 3601-3609. [DOI: https://dx.doi.org/10.1016/j.prostr.2016.06.449]

35. Guo, T.; Liu, Z.X.; Correia, J.; de Jesus, A.M.P. Experimental Study on Fretting-Fatigue of Bridge Cable Wires. Int. J. Fatigue; 2020; 131, 105321. [DOI: https://dx.doi.org/10.1016/j.ijfatigue.2019.105321]

36. Chen, X.J.; Zhu, X.D. Fretting Fatigue Analysis of FGM Coating under out-of-Phase Loadings. Tribol. Int.; 2019; 134, pp. 165-177. [DOI: https://dx.doi.org/10.1016/j.triboint.2019.01.049]

37. Bhatti, N.A.; Wahab, M.A. A Numerical Investigation on Critical Plane Orientation and Initiation Lifetimes in Fretting Fatigue Under out of Phase Loading Conditions. Tribol. Int.; 2017; 115, pp. 307-318. [DOI: https://dx.doi.org/10.1016/j.triboint.2017.05.036]

38. Bhatti, N.A.; Pereira, K.; Wahab, M.A. A Continuum Damage Mechanics Approach for Fretting Fatigue under out of Phase Loading. Tribol. Int.; 2018; 117, pp. 39-51. [DOI: https://dx.doi.org/10.1016/j.triboint.2017.08.009]

39. Martínez, J.C.; Useche, L.V.V.; Wahab, M.A. Numerical Prediction of Fretting Fatigue Crack Trajectory in a Railway Axle Using XFEM. Int. J. Fatigue; 2017; 100, pp. 32-49. [DOI: https://dx.doi.org/10.1016/j.ijfatigue.2017.03.009]

40. McDiarmid, D.L. A General Criterion for High Cycle Multiaxial Fatigue Failure. Fatigue Fract. Eng. Mater. Struct.; 1991; 14, pp. 429-453. [DOI: https://dx.doi.org/10.1111/j.1460-2695.1991.tb00673.x]

41. Smith, K.N.; Watson, P.; Topper, T.H. Stress- Strain Function for the Fatigue of Metals. J. Mater.; 1970; 5, pp. 767-778.

42. Navarro, C.; Vázquez, J.; Domínguez, J. Nucleation and Early Crack Path in Fretting Fatigue. Int. J. Fatigue; 2017; 100, pp. 602-610. [DOI: https://dx.doi.org/10.1016/j.ijfatigue.2016.12.028]

43. Lee, H.; Jin, O.; Mall, S. Fretting Fatigue Behavior of Ti-6Al-4V with Dissimilar Mating Materials. Int. J. Fatigue; 2004; 26, pp. 393-402. [DOI: https://dx.doi.org/10.1016/j.ijfatigue.2003.07.004]

44. Farrahi, G.H.; Minaii, K.; Bahai, H. Fretting Fatigue Behavior of 316L Stainless Steel under Combined Loading Conditions. Int. J. Fatigue; 2019; 128, 105206. [DOI: https://dx.doi.org/10.1016/j.ijfatigue.2019.105206]

45. Han, Q.N.; Lei, X.S.; Su, Y.; Rui, S.S.; Ma, X.F.; Cui, H.T.; Shi, H.J. In-Situ Observation and Finite Element Analysis of Fretting Fatigue Crack Propagation Behavior in 1045 Steel. Chin. J. Aeronaut.; 2021; 34, pp. 131-139. [DOI: https://dx.doi.org/10.1016/j.cja.2020.10.024]

46. Hojjati-Talemi, R.; Wahab, M.A.; Giner, E.; Sabsabi, M. Numerical Estimation of Fretting Fatigue Lifetime Using Damage and Fracture Mechanics. Tribol. Lett.; 2013; 52, pp. 11-25. [DOI: https://dx.doi.org/10.1007/s11249-013-0189-8]

47. Bhatti, N.A.; Pereira, K.; Wahab, M.A. Effect of Stress Gradient and Quadrant Averaging on Fretting Fatigue Crack Initiation Angle and Life. Tribol. Int.; 2019; 131, pp. 212-221. [DOI: https://dx.doi.org/10.1016/j.triboint.2018.10.036]