On the Convergence of Stresses in Fretting Fatigue

Fretting is a phenomenon that occurs at the contacts of surfaces that are subjected to oscillatory relative movement of small amplitudes. Depending on service conditions, fretting may significantly reduce the service life of a component due to fretting fatigue. In this regard, the analysis of stresses at contact is of great importance for predicting the lifetime of components. However, due to the complexity of the fretting phenomenon, analytical solutions are available for very selective situations and finite element (FE) analysis has become an attractive tool to evaluate stresses and to study fretting problems. Recent laboratory studies in fretting fatigue suggested the presence of stress singularities in the stick-slip zone. In this paper, we constructed finite element models, with different element sizes, in order to verify the existence of stress singularity under fretting conditions. Based on our results, we did not find any singularity for the considered loading conditions and coefficients of friction. Since no singularity was found, the present paper also provides some comments regarding the convergence rate. Our analyses showed that the convergence rate in stress components depends on coefficient of friction, implying that this rate also depends on the loading condition. It was also observed that errors can be relatively high for cases with a high coefficient of friction, suggesting the importance of mesh refinement in these situations. Although the accuracy of the FE analysis is very important for satisfactory predictions, most of the studies in the literature rarely provide information regarding the level of error in simulations. Thus, some recommendations of mesh sizes for those who wish to perform FE analysis of fretting problems are provided for different levels of accuracy.


Introduction
Fretting happens when two contacting surfaces are normally loaded and subjected to small amplitude oscillatory relative movement. This amplitude generally varies from 5 to 100 µm [1], but it can be as low as, or even below, 1 µm [2]. Due to its cyclic characteristic and the high stresses gradient in the vicinity of contact, fretting may lead to unexpected failure due to fretting fatigue, being responsible for the premature failure of many common mechanical assemblies, such as bolted joints, shrink-fitted shafts, and dovetail joints. As a consequence, it has been an important research topic that has been vastly studied in the literature [3][4][5][6].
In order to evaluate the effects of different variables (surface finishing, coefficient of friction, normal load, relative slip amplitude, among others) on the characteristics of fretting, different laboratory tests are generally used. One of the most common is a cylinder-on-pad configuration, as illustrated in Figure 1. In this set-up, two cylindrical pads are maintained in contact with a flat specimen through the application of a constant clamping or normal force, F. The specimen is fixed at one end and the other end is subjected to an oscillatory bulk stress σ axial . On application of the bulk stress, the compliance springs transmit an oscillatory tangential force, Q, at the pads. Generally, the tangential load |Q| is smaller than the product of the normal load, F, by the coefficient of friction µ and the contact is divided into two regions: A stick zone and a slip region. In the early 1970s, Nishioka and Hirakawa [7] had already used this configuration to study the effects of slip amplitude in the fatigue strength of specimens. Even in recent research, this test set-up is still very common. For instance, Pierres et al. [8] proposed a combined numerical and experimental approach to simulate fretting fatigue crack growth of 2D and 3D configurations. A similar methodology was used by Luke et al. [9], however, they were interested in simulating crack initiation using different damage parameters and they used laboratory tests to validated their predictions.
Materials 2016, 9,639 2 of 15 the premature failure of many common mechanical assemblies, such as bolted joints, shrink-fitted shafts, and dovetail joints. As a consequence, it has been an important research topic that has been vastly studied in the literature [3][4][5][6].
In order to evaluate the effects of different variables (surface finishing, coefficient of friction, normal load, relative slip amplitude, among others) on the characteristics of fretting, different laboratory tests are generally used. One of the most common is a cylinder-on-pad configuration, as illustrated in Figure 1. In this set-up, two cylindrical pads are maintained in contact with a flat specimen through the application of a constant clamping or normal force, F. The specimen is fixed at one end and the other end is subjected to an oscillatory bulk stress σaxial. On application of the bulk stress, the compliance springs transmit an oscillatory tangential force, Q, at the pads. Generally, the tangential load |Q| is smaller than the product of the normal load, F, by the coefficient of friction μ and the contact is divided into two regions: A stick zone and a slip region. In the early 1970s, Nishioka and Hirakawa [7] had already used this configuration to study the effects of slip amplitude in the fatigue strength of specimens. Even in recent research, this test set-up is still very common. For instance, Pierres et al. [8] proposed a combined numerical and experimental approach to simulate fretting fatigue crack growth of 2D and 3D configurations. A similar methodology was used by Luke et al. [9], however, they were interested in simulating crack initiation using different damage parameters and they used laboratory tests to validated their predictions. Under fretting conditions, the fatigue limit of a material may be shortened by up to 50% [10]. It is known that, in this case, the crack growth phase is significantly different from plain fatigue propagation phase, due to the influence of contact stresses distributions on the crack and vice versa [11]. This contactcrack interaction is particularly important for cracks' length smaller than the magnitude of the contact zone dimension [12] and must be taken into account. For a cylinder-on-plane configuration, the stress and strain field in the specimen can be analytically estimated by a combination of the normal pressure distribution p(x) (due to the normal force, F) and surface traction q(x) (due to the tangential and bulk loads, Q and σaxial, respectively). However, these solutions are valid under a series of conditions, such as infinite and idealized bodies, elastic material properties, and loading conditions, among others. In addition, the stress field near the contact region is variable, multiaxial and non-proportional [13], which provides extra complexity to the phenomena.
Fretting fatigue is a complex phenomenon due the stick-slip zone at the contact interface. This complex phenomenon is not well understood and a recent research report [14] has questioned the applicability of the analytical solution (Cattaneo-Mindlin problem) to the stick-slip problems. In the analytical solution, the superimposing of shear stress due to normal load and due to fatigue load is a linear approximation and ignores the effect of interaction between both loads. Furthermore, recent laboratory measurements [15,16] indicated that the transition from 'static' to 'dynamic' friction (stickslip) can be described by classical fracture mechanics singular solutions of shear cracks, rather than by Coulomb law. This motivates us to investigate whether or not stress singularity takes place at the stickslip zone in fretting conditions. Under fretting conditions, the fatigue limit of a material may be shortened by up to 50% [10]. It is known that, in this case, the crack growth phase is significantly different from plain fatigue propagation phase, due to the influence of contact stresses distributions on the crack and vice versa [11]. This contact-crack interaction is particularly important for cracks' length smaller than the magnitude of the contact zone dimension [12] and must be taken into account. For a cylinder-on-plane configuration, the stress and strain field in the specimen can be analytically estimated by a combination of the normal pressure distribution p(x) (due to the normal force, F) and surface traction q(x) (due to the tangential and bulk loads, Q and σ axial , respectively). However, these solutions are valid under a series of conditions, such as infinite and idealized bodies, elastic material properties, and loading conditions, among others. In addition, the stress field near the contact region is variable, multiaxial and non-proportional [13], which provides extra complexity to the phenomena.
Fretting fatigue is a complex phenomenon due the stick-slip zone at the contact interface. This complex phenomenon is not well understood and a recent research report [14] has questioned the applicability of the analytical solution (Cattaneo-Mindlin problem) to the stick-slip problems. In the analytical solution, the superimposing of shear stress due to normal load and due to fatigue load is a linear approximation and ignores the effect of interaction between both loads. Furthermore, recent laboratory measurements [15,16] indicated that the transition from 'static' to 'dynamic' friction (stick-slip) can be described by classical fracture mechanics singular solutions of shear cracks, rather than by Coulomb law. This motivates us to investigate whether or not stress singularity takes place at the stick-slip zone in fretting conditions. Numerical methodologies have become an interesting option to evaluate stresses at contact and its impact on fretting fatigue lifetime. In this regard, the finite element analysis (FEA) has been widely used over the past few decades. For instance, McVeigh and Farris [17] used finite element analysis to study the influence of the bulk loading σ axial on the contact stresses distributions, and compared the results with analytical approximations, validating the latter. Tur et al. [18] treated the problem considering the effects of plasticity on the contact stress distribution for a Titanium material and analyzed the impact of plastic deformations on the size of the stick zone and peak stresses. They concluded that the plastic zone started at the trailing edge (the edge of the largest slip zone) and that the effects of contact stresses decayed rapidly as the distance from the contact increased.
The focus of this paper is to recognize the existence of stress singularity at the stick-slip zone in fretting fatigue conditions using FEA. In order to do that, a finite element model of a fretting test configuration (cylindrical pad and flat specimen) was created and stresses at the contact interface were monitored and compared with analytical solutions for different mesh sizes and fretting contact conditions. The paper is organized in the following way. Firstly, the analytical solutions of the contact stresses used as references in this study are described in Section 2. Then the finite element models are constructed and details of them are provided in Section 3. Finally, the results are presented and discussed in Section 4 and conclusions are drawn in Section 5.

Analytical Solutions
In this section, we first present the Hertzian solutions for the pressure distribution at the contact interface of a cylinder and a flat surface under normal load. Then, we consider the effect of combined normal and tangential loads, and, finally, we shall present solutions for the effect of bulk stresses on fretting fatigue conditions.

Hertzian Solutions for the Pressure Distribution
As discussed by Johnson [19], the contact pressure distribution, p(x), due to the normal clamping force, F, between the elastic pad and elastic specimen, can be calculated analytically if the following contact conditions hold:
Small strains at contact region; 3.
Bodies can be approximated as a semi-infinite elastic half-space near the contact zone; 4.
In this case, the contact pressure, p(x), is elliptical at a distance, x, from the center of the contact zone (see Figure 1) and is given by [19]: where p max is the maximum contact pressure at the center of the contact; R is the combined curvature; and E˚is the combined modulus of elasticity. Both R and E˚can be defined as: where E i , for i = 1,2 are the Young's Modulus and ν i , for i = 1,2 are the Poisson's ratio for the first and second bodies, respectively. The flat specimen can be considered as a cylinder with an infinitely large radius R 1 = 8 and the combined curvature, R, becomes equal to the radius of the surface of the pad R 2 .
Considering that contact should occur only inside the loaded area, and, also, the fact that all contact regions must be in compression, the semi-contact width, a, and the applied load, F, are related by: where t is the thickness of cylinder pad. The elastic deformation of the surfaces results in a rectangular contact region of area equal to 2aˆt.

Solutions for Combined Normal and Tangential Loads
When studying fretting, it is necessary to consider, not only the normal loading condition, but also the effect of the tangential frictional force, Q. The Coulomb friction law can be used to model the contact shear traction, q(x), at an arbitrary position, x, as a function of the normal contact pressure, p(x), and the coefficient of friction, µ. If Q is smaller than the product of µ and the normal load, F, the contact region will be divided into two different zones: Stick and slip, in which the width of the stick zone is denoted by c. In this case, the contact shear traction can be seen as combination of a pressure distribution and two superposed shear tractions, q'(x) due to p(x) and q"(x) due to Q, as shown in Figure 2. Considering that contact should occur only inside the loaded area, and, also, the fact that all contact regions must be in compression, the semi-contact width, a, and the applied load, F, are related by: where t is the thickness of cylinder pad. The elastic deformation of the surfaces results in a rectangular contact region of area equal to 2a × t.

Solutions for Combined Normal and Tangential Loads
When studying fretting, it is necessary to consider, not only the normal loading condition, but also the effect of the tangential frictional force, Q. The Coulomb friction law can be used to model the contact shear traction, q(x), at an arbitrary position, x, as a function of the normal contact pressure, p(x), and the coefficient of friction, μ. If Q is smaller than the product of μ and the normal load, F, the contact region will be divided into two different zones: Stick and slip, in which the width of the stick zone is denoted by c. In this case, the contact shear traction can be seen as combination of a pressure distribution and two superposed shear tractions, q'(x) due to p(x) and q''(x) due to Q, as shown in Figure 2. The complete expression for the shear traction q(x) can be written as [12,19]: where 1 .

Effect of Bulk Load σaxial on Contact Shear Traction
According to Hills and Nowell [12], the contact shear traction presented above can be adjusted for the presence of bulk stresses σaxial. This causes an eccentricity to the solution presented in Section 2.2, and for the case of negative tangential load, it can be written as [12]: The complete expression for the shear traction q(x) can be written as [12,19]: where c a " b 1´Q µF .

Effect of Bulk Load σ axial on Contact Shear Traction
According to Hills and Nowell [12], the contact shear traction presented above can be adjusted for the presence of bulk stresses σ axial . This causes an eccentricity to the solution presented in Section 2.2, and for the case of negative tangential load, it can be written as [12]: where c a " b 1´Q µF and e " aσ axial 4µp max .  Figure 3 shows a typical normalized shear traction distribution for fretting fatigue specimen using Equation (6). Note that, based on this distribution, it is possible to determine the size of the stick and slip zones and also the peak values of shear stresses. For this paper, we monitored two peak values of shear tractions q(x 1 ) and q(x 2 ), at the leading edge and at trailing edge sides (the edge of the largest slip zone [20]), respectively. where 1 and . Figure 3 shows a typical normalized shear traction distribution for fretting fatigue specimen using Equation (6). Note that, based on this distribution, it is possible to determine the size of the stick and slip zones and also the peak values of shear stresses. For this paper, we monitored two peak values of shear tractions q(x1) and q(x2), at the leading edge and at trailing edge sides (the edge of the largest slip zone [20]), respectively.

Effect of Bulk Load σaxial on Subsurface Stresses
In his literature review, Mutoh [4] mentioned studies showing that fretting fatigue cracks, which propagate to material final ruptures, originate in the edge of the contact area (x = a), while small arrested cracks are initiated near the maximum shear traction q(x2). Other research [12,21,22] has also pointed out that the principal crack initiates near the trailing edge (x = a). The reason for that may be related to the contribution of the principal stress σxx in the stress state at the contact interface. As discussed by Szolwinski and Farris [23], studies showed that the sharp peak in tangential stresses σxx,max, at trailing edge of the contact region (see Figure 4), might play a significant role on fretting fatigue crack initiation. There are analytical solutions for subsurface elastic stresses, σxx, as function of x for a given normal and tangential loads (F and Q) and coefficient of friction, μ, in the slip zone [12,19,24]. For instance, Szolwinski and Farris [24] provided an analytical solution for the stress distribution, σxx, treating the problem as a superposition of individual stress components, caused by the normal pressure distribution and surface tractions, q'(x) and q''(x).

Effect of Bulk Load σ axial on Subsurface Stresses
In his literature review, Mutoh [4] mentioned studies showing that fretting fatigue cracks, which propagate to material final ruptures, originate in the edge of the contact area (x = a), while small arrested cracks are initiated near the maximum shear traction q(x 2 ). Other research [12,21,22] has also pointed out that the principal crack initiates near the trailing edge (x = a). The reason for that may be related to the contribution of the principal stress σ xx in the stress state at the contact interface. As discussed by Szolwinski and Farris [23], studies showed that the sharp peak in tangential stresses σ xx,max , at trailing edge of the contact region (see Figure 4), might play a significant role on fretting fatigue crack initiation. where 1 and . Figure 3 shows a typical normalized shear traction distribution for fretting fatigue specimen using Equation (6). Note that, based on this distribution, it is possible to determine the size of the stick and slip zones and also the peak values of shear stresses. For this paper, we monitored two peak values of shear tractions q(x1) and q(x2), at the leading edge and at trailing edge sides (the edge of the largest slip zone [20]), respectively.

Effect of Bulk Load σaxial on Subsurface Stresses
In his literature review, Mutoh [4] mentioned studies showing that fretting fatigue cracks, which propagate to material final ruptures, originate in the edge of the contact area (x = a), while small arrested cracks are initiated near the maximum shear traction q(x2). Other research [12,21,22] has also pointed out that the principal crack initiates near the trailing edge (x = a). The reason for that may be related to the contribution of the principal stress σxx in the stress state at the contact interface. As discussed by Szolwinski and Farris [23], studies showed that the sharp peak in tangential stresses σxx,max, at trailing edge of the contact region (see Figure 4), might play a significant role on fretting fatigue crack initiation. There are analytical solutions for subsurface elastic stresses, σxx, as function of x for a given normal and tangential loads (F and Q) and coefficient of friction, μ, in the slip zone [12,19,24]. For instance, Szolwinski and Farris [24] provided an analytical solution for the stress distribution, σxx, treating the problem as a superposition of individual stress components, caused by the normal pressure distribution and surface tractions, q'(x) and q''(x). There are analytical solutions for subsurface elastic stresses, σ xx , as function of x for a given normal and tangential loads (F and Q) and coefficient of friction, µ, in the slip zone [12,19,24]. For instance, Szolwinski and Farris [24] provided an analytical solution for the stress distribution, σ xx , treating the problem as a superposition of individual stress components, caused by the normal pressure distribution and surface tractions, q'(x) and q"(x).
Although the addition of the bulk stress σ axial brings some extra complexity to the problem, there are still some simplified equations to estimate stresses at contact. McVeigh and Farris [17] adjusted the analytical solution from Szolwinski and Farris [24] by adding bulk stress in the distribution of σ xx . Szolwinski and Farris [23], based on the work done by McVeigh and Farris [17], provided a simplified equation to estimate the maximum peak stress σ xx,max as:

Finite Element Model: Cylinder Pad on Flat Specimen
A parametric 2D finite element model was created in ABAQUS ® and an analysis of the fretting cycle was performed, aiming to study the model response to different mesh sizes. Three values of coefficients of friction were considered (0.3, 0.85 and 2.0). These variable values of coefficient of friction (COF) allowed us to study different configurations of stick-slip regions and, therefore, to simulate different fretting scenarios.
The model details, such as geometry, material properties, mesh details, boundary conditions and loading history, are presented here. Two FE models were developed and their dimensions and boundary conditions are shown in Figure 5. The models were composed of only two parts: A pad and a specimen, which represents half of the experimental set-up, due to its symmetry. In order to check the influence of different geometries, the radius of the pad was also variable in those models, and two values were chosen: 50 mm and 10 mm. Both parts were made of aluminum 2420-T3, having material properties which are summarized in Table 1. We did not consider any plasticity effect in this study, only an elastic material response. Stress analysis was carried out by applying a normal load (F = 543 N) and oscillatory axial and reaction stresses to the specimen, reflecting a fretting cycle. Although the addition of the bulk stress σaxial brings some extra complexity to the problem, there are still some simplified equations to estimate stresses at contact. McVeigh and Farris [17] adjusted the analytical solution from Szolwinski and Farris [24] by adding bulk stress in the distribution of σxx. Szolwinski and Farris [23], based on the work done by McVeigh and Farris [17], provided a simplified equation to estimate the maximum peak stress σxx,max as:

Finite Element Model: Cylinder Pad on Flat Specimen
A parametric 2D finite element model was created in ABAQUS ® and an analysis of the fretting cycle was performed, aiming to study the model response to different mesh sizes. Three values of coefficients of friction were considered (0.3, 0.85 and 2.0). These variable values of coefficient of friction (COF) allowed us to study different configurations of stick-slip regions and, therefore, to simulate different fretting scenarios.
The model details, such as geometry, material properties, mesh details, boundary conditions and loading history, are presented here. Two FE models were developed and their dimensions and boundary conditions are shown in Figure 5. The models were composed of only two parts: A pad and a specimen, which represents half of the experimental set-up, due to its symmetry. In order to check the influence of different geometries, the radius of the pad was also variable in those models, and two values were chosen: 50 mm and 10 mm. Both parts were made of aluminum 2420-T3, having material properties which are summarized in Table 1. We did not consider any plasticity effect in this study, only an elastic material response. Stress analysis was carried out by applying a normal load (F = 543 N) and oscillatory axial and reaction stresses to the specimen, reflecting a fretting cycle.    The master-slave algorithm in ABAQUS ® was used to describe the contact behavior and the Lagrange multiplier formulation was used to define the tangential behavior of the contact pair. The surface-to-surface and finite sliding options were used to define the contact interaction.
A 2D quadrilateral, 4-node (bilinear), plane strain, reduced integration element (CPE4R) was used in both models. Different mesh sizes were considered at the contact interface and increased as the distance from the contact region increased. In order to create a fine mesh at the contact region, the models were partitioned and the edges were seeded. The values of the mesh element size along the contact region varied according to the following list: 20, 10, 5, 2.5, 1.25, 0.625 and 0.3125 µm. Details of the seeding used to generate the mesh and also of the model partition dimensions are shown in Figure 6. The partition dimensions were dependent on the radius of the pad, being calculated based on the semi-contact width, a, from Equation (4). An illustration of one of the meshes used in this study is also presented in Figure 6. The master-slave algorithm in ABAQUS ® was used to describe the contact behavior and the Lagrange multiplier formulation was used to define the tangential behavior of the contact pair. The surface-to-surface and finite sliding options were used to define the contact interaction.
A 2D quadrilateral, 4-node (bilinear), plane strain, reduced integration element (CPE4R) was used in both models. Different mesh sizes were considered at the contact interface and increased as the distance from the contact region increased. In order to create a fine mesh at the contact region, the models were partitioned and the edges were seeded. The values of the mesh element size along the contact region varied according to the following list: 20, 10, 5, 2.5, 1.25, 0.625 and 0.3125 μm. Details of the seeding used to generate the mesh and also of the model partition dimensions are shown in Figure 6. The partition dimensions were dependent on the radius of the pad, being calculated based on the semi-contact width, a, from Equation (4). An illustration of one of the meshes used in this study is also presented in Figure 6. Due to the symmetry of the problem, the bottom of the specimen (representing the axial centerline of the specimen) was restricted from vertical movement in the y direction (Uy = 0). The sides of the pads were restricted from horizontal movement in the x direction (Ux = 0) and the Multiple point constraints (MPC) tie constraint was also used at the top surface of the pad to guarantee that it would not rotate due to the applied concentrated load, F.
The effect of the compliance spring and tangential load Q were modeled as a cyclic reaction stress (σreaction). This reaction stress is obtained as: where b is the specimen width (b = 10 mm); t is the specimen thickness (t = 4 mm). The values of Q and σaxial are obtained from experimental data (see Table 2). For this study, they are taken from the experimental set-up FF1 in Reference [25]. Due to the symmetry of the problem, the bottom of the specimen (representing the axial centerline of the specimen) was restricted from vertical movement in the y direction (U y = 0). The sides of the pads were restricted from horizontal movement in the x direction (U x = 0) and the Multiple point constraints (MPC) tie constraint was also used at the top surface of the pad to guarantee that it would not rotate due to the applied concentrated load, F.
The effect of the compliance spring and tangential load Q were modeled as a cyclic reaction stress (σ reaction ). This reaction stress is obtained as: where b is the specimen width (b = 10 mm); t is the specimen thickness (t = 4 mm). The values of Q and σ axial are obtained from experimental data (see Table 2). For this study, they are taken from the experimental set-up FF1 in Reference [25]. Table 2. Values of maximum and minimum σ reaction and σ axial , based on data from experimental test FF1 from Reference [26].

Steps σ axial [MPa] σ reaction [MPa] Q [N]
Step In order to simulate fretting fatigue conditions, the loads are applied in three steps (see Figure 7), with adaptive time steps in ABAQUS ® . In the first loading step, the top pad was pressed against the specimen surface by a normal load F = 543 N and this compressed condition was held constant until the end of the cycle. Then, both axial and reaction maximum stresses were applied to the sides of the specimen (values are presented in Table 2). Finally, in the third loading step, both axial and reaction minimum stresses were applied. In order to simulate fretting fatigue conditions, the loads are applied in three steps (see Figure 7), with adaptive time steps in ABAQUS ® . In the first loading step, the top pad was pressed against the specimen surface by a normal load F = 543 N and this compressed condition was held constant until the end of the cycle. Then, both axial and reaction maximum stresses were applied to the sides of the specimen (values are presented in Table 2). Finally, in the third loading step, both axial and reaction minimum stresses were applied.

Results and Discussion
In order to recognize a singularity's presence, the methodology presented by Sinclair [27] was adopted. Accordingly, the element size in the models was successively systematically halved for a sequence of seven analyses and the magnitude of maximum stress values was examined. The following stress components were monitored at the maximum axial loading condition (end of loading step 2): The contact shear traction peak at trailing edge side q(x1) and at leading edge side q(x2) (see Figure 3) and the peak tangential stress in the x direction, σxx,max (see Figure 4). The influence of the mesh size on the values of the ratios between stick and slip zones sizes (c/a) is also considered here. The slip zone size, c, is obtained by measuring the position in the contact that have non-zero values of slip and the contact width, a, is obtained by the position in the x direction of the edges of the contact region, both calculated from ABAQUS ® .
The results of various stress components and for the ratios between stick and slip zones sizes (c/a) are presented in Table 3, for different values of coefficient of friction and different radius of cylindrical pad. FEA results were also compared with analytical solutions (Equations (6) and (7), presented in Sections 2.3 and 2.4, respectively). The values of shear traction at trailing and leading edges seem to converge to the analytical solution for all values of coefficient of friction. Regarding the values of peak tangential stress σxx,max, they seem to converge, but to a different value than the estimate from Equation (7).

Results and Discussion
In order to recognize a singularity's presence, the methodology presented by Sinclair [27] was adopted. Accordingly, the element size in the models was successively systematically halved for a sequence of seven analyses and the magnitude of maximum stress values was examined. The following stress components were monitored at the maximum axial loading condition (end of loading step 2): The contact shear traction peak at trailing edge side q(x 1 ) and at leading edge side q(x 2 ) (see Figure 3) and the peak tangential stress in the x direction, σ xx,max (see Figure 4). The influence of the mesh size on the values of the ratios between stick and slip zones sizes (c/a) is also considered here. The slip zone size, c, is obtained by measuring the position in the contact that have non-zero values of slip and the contact width, a, is obtained by the position in the x direction of the edges of the contact region, both calculated from ABAQUS ® .
The results of various stress components and for the ratios between stick and slip zones sizes (c/a) are presented in Table 3, for different values of coefficient of friction and different radius of cylindrical pad. FEA results were also compared with analytical solutions (Equations (6) and (7), presented in Sections 2.3 and 2.4, respectively). The values of shear traction at trailing and leading edges seem to converge to the analytical solution for all values of coefficient of friction. Regarding the values of peak tangential stress σ xx,max , they seem to converge, but to a different value than the estimate from Equation (7). This is reasonable, since this equation provides only an approximate value of σ xx,max .
Note that the non-dimensional parameter (c/a) also converges on the analytical solution for all values of coefficient of friction and pad radius. To examine convergence in Table 3, the relative error between FE and analytical solutions was considered. If those errors do not decrease with successively refined analysis, divergence occurs and the presence of a singularity is detected.

Stress Singularity Check: Influence of Mesh Size on Stress Components
In order to analyze the influence of mesh size on the contact shear traction, the analytical solution was chosen as a reference value. The relative error between FE and analytical solutions (e rel,an ) was calculated as: where φ i max is the maximum variable output (contact shear stress, maximum tangential stress or ratio between stick slip zone sizes) in the ith model and φ a max in the analytical solution (see Table 3). Higher coefficient of friction implies stronger gradients in the stress distribution, and singularities are expected to happen for higher values of coefficient of friction. The relative error between FE and analytical solutions for the contact shear traction stress component for different coefficients of friction and pad radius are presented in Figure 8. The results show that the error is decreasing as the mesh size reduces, independent of the value of coefficient of friction and pad radius. Thus, the analysis is converging, even if only slowly, and no singularity was found for any of the tested loading conditions, pad radius, and coefficients of friction.
for higher coefficients of friction, the rate of convergence reduces and it is necessary to use relative fine meshes to guarantee reasonable results. For instance, for coefficient of friction equal to 2.0, a mesh size of 1.25 μm is enough to guarantee that the relative error on the shear traction peak is smaller than 10% for all analyzed cases. However, for the same coefficient of friction and a mesh size of 5 μm, the error can increase to almost 40%, for the contact shear traction peak at leading edge. The dependence of the rate of convergence on the coefficient of friction can be further investigated by analyzing the contact shear traction at contact interface. As can be seen in Figure 9, for the case of high coefficient of friction, the contact shear traction distribution has very sharp peaks at both leading and trailing edges, justifying the necessity of a very fine mesh to accurately capture those steep gradients. It is clear that for a higher coefficient of friction, a very fine mesh size is required in order to achieve convergence. This is due to the fact that the value of the friction coefficient affects the contact stress distribution and the larger the coefficient of friction, the steeper the stress gradient. The value of coefficient of friction also affects the stick-slip zone size, which is an important parameter to determine a suitable mesh size as explained in Section 4.2. Moreover, it can also be seen that the rate of convergence is dependent on the coefficient of friction for both cases of pad geometry. As different values of coefficient of frictions represent different loading conditions (various sizes of stick zone in comparison with the contact dimension), one might conclude that the rate of convergence of the solution depends on the loading condition. For the smallest coefficient of friction, a relative coarse mesh (around 20 µm) at the contact is sufficient for obtaining reasonable accurate shear stresses at contact, with relative error smaller than 10% for all analyzed cases. However, for higher coefficients of friction, the rate of convergence reduces and it is necessary to use relative fine meshes to guarantee reasonable results. For instance, for coefficient of friction equal to 2.0, a mesh size of 1.25 µm is enough to guarantee that the relative error on the shear traction peak is smaller than 10% for all analyzed cases. However, for the same coefficient of friction and a mesh size of 5 µm, the error can increase to almost 40%, for the contact shear traction peak at leading edge.
The dependence of the rate of convergence on the coefficient of friction can be further investigated by analyzing the contact shear traction at contact interface. As can be seen in Figure 9, for the case of high coefficient of friction, the contact shear traction distribution has very sharp peaks at both leading and trailing edges, justifying the necessity of a very fine mesh to accurately capture those steep gradients. It is clear that for a higher coefficient of friction, a very fine mesh size is required in order to achieve convergence. This is due to the fact that the value of the friction coefficient affects the contact stress distribution and the larger the coefficient of friction, the steeper the stress gradient. The value of coefficient of friction also affects the stick-slip zone size, which is an important parameter to determine a suitable mesh size as explained in Section 4.2. As discussed before, the peak stress σxx,max, seems to converge to a different value than the estimate from Equation (7). Therefore, in order to study the convergence of the results of the, instead of considering the analytical solution as a reference, the maximum stresses between two subsequent mesh refinements were used to calculate the relative error erel as: (10) where is the maximum variable output (contact shear stress, maximum tangential stress or ratio between stick slip zone sizes) in the ith model and in the (i + 1)th model. The results of the relative error between two consecutive mesh sizes for the maximum tangential stress are presented in Figure 10. The results also show that the error reduces with successively refined analysis, independent of the value of coefficient of friction or pad radius. Thus, once again, one may conclude that the analysis is converging and no singularity's presence was found.
Additionally, some comments regarding the convergence rate can be made. For a mesh size of 0.625 μm, the relative error is around 5% for all analyzed scenarios, and results can be considered to be good [28]. As for the shear traction component, the convergence rate of the maximum tangential stress depends upon the coefficient of friction. Again, for the smallest coefficient of friction, the convergence rate is the highest. This dependency may be further investigated by checking the tangential stress distribution at contact, as shown in Figure 11. As discussed before, the peak stress σ xx,max , seems to converge to a different value than the estimate from Equation (7). Therefore, in order to study the convergence of the results of the, instead of considering the analytical solution as a reference, the maximum stresses between two subsequent mesh refinements were used to calculate the relative error e rel as: where φ i max is the maximum variable output (contact shear stress, maximum tangential stress or ratio between stick slip zone sizes) in the ith model and φ i`1 max in the (i + 1)th model. The results of the relative error between two consecutive mesh sizes for the maximum tangential stress are presented in Figure 10. The results also show that the error reduces with successively refined analysis, independent of the value of coefficient of friction or pad radius. Thus, once again, one may conclude that the analysis is converging and no singularity's presence was found.
Additionally, some comments regarding the convergence rate can be made. For a mesh size of 0.625 µm, the relative error is around 5% for all analyzed scenarios, and results can be considered to be good [28]. As for the shear traction component, the convergence rate of the maximum tangential stress depends upon the coefficient of friction. Again, for the smallest coefficient of friction, the convergence rate is the highest. This dependency may be further investigated by checking the tangential stress distribution at contact, as shown in Figure 11. As mentioned in Reference [29], increasing the pad radius causes a reduction on the peak contact pressure and increases the contact width. Therefore, for the same loading conditions, the contact pressure distribution has a steeper gradient for pads with smaller radius. As discussed by the authors of [12,19], and presented in Section 2, the analytical distribution of shear stress at contact can be seen as a superposition of contact pressure distribution and two shear tractions ( Figure 2). Thus, it is expected that the gradient of the distribution of the tangential stresses at the trailing edge is higher for the model with smaller pad radius. The smaller contact width for smaller pad radius also implies higher peak values of tangential stresses in a smaller area ( Figure 11). Therefore, for the same loading conditions, a finer mesh is necessary to properly to capture this changes in the model with a smaller pad radius.
It can also be observed in Figure 11 that the peak values and, therefore, the gradients of the distribution of the tangential stresses are smaller for low coefficients of friction. Consequently, the convergence rate, for the model with a 10 mm pad radius, is slower than for the model with a 50 mm pad radius, as can be seen in Figure 10. This impact of geometry on convergence rate is expected, as the smallest radius implies smaller contact region (Equation (4)) for the same loading condition. It also implies higher peaks of tangential stresses in a smaller area. Thus, for the same level of accuracy, a finer mesh is required in the model with a smaller pad radius.

Fretting Fatigue Convergence Map
As pointed out by Ainsworth and Oden [30], although aware of the existence of numerical errors, the analyst is seldom interested in quantifying them. In fretting fatigue, the quality of a simulation is generally assessed by visual comparison between finite element results and analytical solutions [18,31,32], and information regarding the error is rarely provided [22]. The accuracy of the fretting contact stress  As mentioned in Reference [29], increasing the pad radius causes a reduction on the peak contact pressure and increases the contact width. Therefore, for the same loading conditions, the contact pressure distribution has a steeper gradient for pads with smaller radius. As discussed by the authors of [12,19], and presented in Section 2, the analytical distribution of shear stress at contact can be seen as a superposition of contact pressure distribution and two shear tractions ( Figure 2). Thus, it is expected that the gradient of the distribution of the tangential stresses at the trailing edge is higher for the model with smaller pad radius. The smaller contact width for smaller pad radius also implies higher peak values of tangential stresses in a smaller area ( Figure 11). Therefore, for the same loading conditions, a finer mesh is necessary to properly to capture this changes in the model with a smaller pad radius.
It can also be observed in Figure 11 that the peak values and, therefore, the gradients of the distribution of the tangential stresses are smaller for low coefficients of friction. Consequently, the convergence rate, for the model with a 10 mm pad radius, is slower than for the model with a 50 mm pad radius, as can be seen in Figure 10. This impact of geometry on convergence rate is expected, as the smallest radius implies smaller contact region (Equation (4)) for the same loading condition. It also implies higher peaks of tangential stresses in a smaller area. Thus, for the same level of accuracy, a finer mesh is required in the model with a smaller pad radius.

Fretting Fatigue Convergence Map
As pointed out by Ainsworth and Oden [30], although aware of the existence of numerical errors, the analyst is seldom interested in quantifying them. In fretting fatigue, the quality of a simulation is generally assessed by visual comparison between finite element results and analytical solutions [18,31,32], and information regarding the error is rarely provided [22]. The accuracy of the fretting contact stress As mentioned in Reference [29], increasing the pad radius causes a reduction on the peak contact pressure and increases the contact width. Therefore, for the same loading conditions, the contact pressure distribution has a steeper gradient for pads with smaller radius. As discussed by the authors of [12,19], and presented in Section 2, the analytical distribution of shear stress at contact can be seen as a superposition of contact pressure distribution and two shear tractions ( Figure 2). Thus, it is expected that the gradient of the distribution of the tangential stresses at the trailing edge is higher for the model with smaller pad radius. The smaller contact width for smaller pad radius also implies higher peak values of tangential stresses in a smaller area ( Figure 11). Therefore, for the same loading conditions, a finer mesh is necessary to properly to capture this changes in the model with a smaller pad radius.
It can also be observed in Figure 11 that the peak values and, therefore, the gradients of the distribution of the tangential stresses are smaller for low coefficients of friction. Consequently, the convergence rate, for the model with a 10 mm pad radius, is slower than for the model with a 50 mm pad radius, as can be seen in Figure 10. This impact of geometry on convergence rate is expected, as the smallest radius implies smaller contact region (Equation (4)) for the same loading condition. It also implies higher peaks of tangential stresses in a smaller area. Thus, for the same level of accuracy, a finer mesh is required in the model with a smaller pad radius.

Fretting Fatigue Convergence Map
As pointed out by Ainsworth and Oden [30], although aware of the existence of numerical errors, the analyst is seldom interested in quantifying them. In fretting fatigue, the quality of a simulation is generally assessed by visual comparison between finite element results and analytical solutions [18,31,32], and information regarding the error is rarely provided [22]. The accuracy of the fretting contact stress calculations is of significant importance, as these stresses impact directly on the crack propagation phase. Therefore, estimating errors for those stresses is of great interest to ensure an accurate analysis.
Aiming to help researchers to easily determine the required element size for their finite element analysis for a given stick-slip ratio and desired accuracy, a 'fretting fatigue convergence map' was produced and is presented in Figure 12. This map was constructed by plotting the stick-slip ratio (c/a) against the element size in the contact zone for different numerical accuracies (1%, 2% and 5%) and may be used as a reference for choosing the element size in FEA of fretting fatigue (cylinder on plane configuration). calculations is of significant importance, as these stresses impact directly on the crack propagation phase. Therefore, estimating errors for those stresses is of great interest to ensure an accurate analysis. Aiming to help researchers to easily determine the required element size for their finite element analysis for a given stick-slip ratio and desired accuracy, a 'fretting fatigue convergence map' was produced and is presented in Figure 12. This map was constructed by plotting the stick-slip ratio (c/a) against the element size in the contact zone for different numerical accuracies (1%, 2% and 5%) and may be used as a reference for choosing the element size in FEA of fretting fatigue (cylinder on plane configuration).

Conclusions
In this paper, we investigate the singularity's presence in fretting fatigue stresses distributions at a contact interface. Different scenarios were studied: Three different coefficients of friction (replicating different loading conditions) and two different pad geometries (radius of 50 mm and 10 mm). For the considered loading conditions and coefficient of frictions, we could not find any indications that a singularity happens as the mesh becomes denser, our results all converged as the element size reduced. Additionally, the convergence rate of the finite element models was discussed as a function of the coefficient of friction. We have confirmed that the convergence rate is smaller for higher coefficients of friction. This means that, for a fixed element size, the level of error in the analysis depends on the loading condition. Therefore, it is recommended that the analyst perform a mesh convergence study for each of his loading condition of interest, as it may impact the accuracy of results. Considering all scenarios that we have studied, a choice of element size of 0.625 μm at contact provided the smallest relative error for all variables, being around, or even smaller than, 5% and producing satisfactory results. A 'fretting fatigue convergence map' was also constructed, providing information on the required element size for a specific stick-slip ratio and different levels of accuracy.

Conclusions
In this paper, we investigate the singularity's presence in fretting fatigue stresses distributions at a contact interface. Different scenarios were studied: Three different coefficients of friction (replicating different loading conditions) and two different pad geometries (radius of 50 mm and 10 mm). For the considered loading conditions and coefficient of frictions, we could not find any indications that a singularity happens as the mesh becomes denser, our results all converged as the element size reduced. Additionally, the convergence rate of the finite element models was discussed as a function of the coefficient of friction. We have confirmed that the convergence rate is smaller for higher coefficients of friction. This means that, for a fixed element size, the level of error in the analysis depends on the loading condition. Therefore, it is recommended that the analyst perform a mesh convergence study for each of his loading condition of interest, as it may impact the accuracy of results. Considering all scenarios that we have studied, a choice of element size of 0.625 µm at contact provided the smallest relative error for all variables, being around, or even smaller than, 5% and producing satisfactory results. A 'fretting fatigue convergence map' was also constructed, providing information on the required element size for a specific stick-slip ratio and different levels of accuracy.