Fretting wear mechanical analysis of double rough surfaces based on energy method

A fretting wear model of a rough surface that conforms to the actual situation is established to accurately reveal the wear mechanism of the connection structure. In the ABAQUS software, the UMESHMOTION subroutine and the energy dissipation model are used to simulate the fretting wear of double rough surfaces. The new model, a single rough surface model, and a smooth model are compared to analyze their differences. In addition, the influence of surface roughness, material, and friction coefficient on the fretting wear of rough surfaces is systematically explored through finite element simulation. The results show that the model's reliability has been verified through Hertz's theory and experiments. The stress and wear of the contact surface are more realistically reflected by the double roughness model. Besides, with the increase of surface roughness and material rigidity and the decrease of friction coefficient, the wear of the double rough surface model becomes more severe. The research work provides a theoretical basis for the design and performance prediction of the connection structure.


Introduction
Mechanical equipment is widely used in practical projects like transportation, power transmission equipment, bridge engineering, navigation, and aerospace. They are usually assembled by several parts in different ways, such as bearings, bolts, rivet connections, spline fits, etc.. 1 Fretting is a small amplitude (less than 100 μm, generally 2∼20 μm) reciprocating motion between two contact surfaces caused by vibration, which usually occurs in these seemingly close-fitting mechanical connections. The fretting will cause wear of the contact surface, accelerate the sprouting, and expand cracks. It will also cause the member to bite and loosen so that the fatigue life of the member is seriously reduced. 2,3 Fretting wear has become the main cause of failure for connection parts. Therefore, understanding the fretting wear mechanism of the connecting structure is extremely important.
At present, experimental testing and finite element simulation are two main methods to study the phenomenon of fretting wear. 4 Compared with the experimental method, the finite element simulation can produce a much more detailed set of results. Besides, it is very convenient and time-saving, suiting to explore the potential mechanism of fretting wear. As the first numerical model, the Archard model gives the relationship between the amount of wear and the relative slip distance. 5 Then, this model is used in finite element simulation and verified by experimental data in Mccoll et al.. 6 Subsequently, the finite element model was extended to more general problems in fretting wear. For example, an optimized finite element model which gives more accurate results of fretting wear was established by Madge et al., 7 where the UMESHMOTION subroutine and ALE adaptive mesh technology in ABAQUS are used. Kasarekar et al. 8 pointed out that surface roughness is a nonnegligible factor when calculating wear and proposed a model considering the effect of surface roughness on the evolution of wear under partial slip. Leonard et al. 9 studied the effect of surface roughness in the gross slip regime and found that roughness only acts at the beginning of the test Shen et al. 10 combined the Archard model and a discrete quasi-static model to analyze the dynamic wear process of bearings, where complex nonlinear wear processes were captured. Tang et al. 11 established a two-dimensional cylindrical/planar contact finite element model to study the fretting wear behavior of zirconium alloys. They found that the wear characteristics under partial sliding and full sliding conditions are inconsistent.
The Archard model has a limitation; the influence of the friction coefficient on the wear rate under the alternating sliding state 12 is not considered. In order to resolve this issue, a wear calculation model based on energy dissipation was proposed by Fouvry and Sauger et al., 13,14 where the relationship between wear volume and accumulated dissipated energy is determined by the energy wear coefficient. In addition, experimental studies have shown that the energy dissipation method is superior to the Archard method. 15,16 Paulin et al. 17 developed a finite element model based on an energy equation verified by experimental results. After that, the energy method is applied to practical problems. Zhang et al. 18 developed a finite element model considering the local deformation of the thread surface to simulate the self-loosening phenomenon of bolts under lateral loads. Jiang et al. 19 studied the electrical contact durability and electromagnetic induction failure under fretting wear. The surface roughness has a great influence on the wear, fatigue, temperature rise, [20][21][22][23][24] and friction characteristics 25 of the connecting structure. However, the effect of roughness is not considered in most finite element models. 4,16,17,26 For the characterization of rough surfaces, Pérez-Ràfols and Almqvist 27 proposed a method that allows to specify the power spectrum and the height probability distribution to generate a realistic surface topography. Bonari et al. 28 developed a computational method for simulating rough surface contact problems in finite elements. The method can be used to study wear problems considering surface roughness. To study temperature rise distribution in fretting contact, Qin et al. 29 employed the Weierstrass-Mandelbrot function (W-M function) to approximate the rough contact surface. In their finite element model, the two rough surfaces are characterized by a rigid, smooth plane and an equivalent rough surface. Pereira et al. 30 developed a multi-scale method to study the effect of roughness on fretting wear, but the rough contact surface is replaced by a smooth surface. The above study did not sufficiently consider the roughness of both surfaces, so the actual fretting wear mechanism could not be accurately revealed.
In this paper, a finite element model of fretting wear considering two rough contact surfaces is developed. The fractal curve generated by the W-M function is used to simulate the rough contact surface. The dynamic simulation of fretting wear is realized through the energy dissipation model and the UMESHMOTION subroutine. The difference between the rough-, the smooth-, and single rough surface models is analyzed through comparison. The significant factors that may affect the fretting wear are systematically investigated, including the roughness, material properties, and friction coefficient.

Characterization of surface topography
The fractal theory has been found useful to characterize rough surfaces because the contour curves of rough surfaces have self-similarity. 31 The curve generated by the W-M function has the characteristics of continuous and non-differentiation, and its expression is given as 32,33 where Z(x) is the height of the curve, x is the corresponding coordinate, G is the characteristic scale coefficient, D is the fractal dimension (1 < D < 2), γ is a fixed value larger than 1 (usually 1.5), n is a spatial frequency coefficient, γ n is the spatial frequency, which is equal to the inverse of the surface wavelength (γ n = 1/λ n ), and n 1 is an ordinal number corresponding to the lowest cut-off frequency, which is related to the sampling length L, γ n 1 = 1/L. For the case L = 0.6 mm and γ = 1.5 and the characteristic scale factor G = 4.35 × 10 −12 m, the fractal curve is shown in Figure 1. Then, the surface roughness Ra corresponding to the fractal dimension is obtained by substituting the curve value into equation (2). 34 The calculated results of surface roughness Ra are 0.8 μm, 0.5 μm, and 0.2 μm when the fractal dimension D is 1.52, 1.57, and

Ra
Finite element software modeling The sample is composed of a semi-cylindrical body and a flat, 35,36 as shown in Figure 2(a). Table 1 shows the material parameters set by the finite element model. Among them, the material of the semi-cylindrical body is alloy steel, ductile iron (QT800-2), or gray cast iron (HT200), and the material of the flat is alloy steel. 37 The radius of the semi-cylindrical body is R = 6 mm, and the length and width of the flat are 12 mm and 6 mm, respectively. At the center area of the upper and lower specimens, rough contour curves are introduced separately to form two rough surfaces in contact, as shown in Figure 2(b). The length of the rough profile curve is 0.4 mm. The rough contour curve is fitted by using a Python script to import the coordinate values of the curve in Figure 1 into the ABAQUS software. Figure 2(c) shows the loading history of the model, which can be divided into two stages. The uniform load P with a magnitude of 20 MPa is applied in step 1 and remains unchanged in the subsequent analysis steps. A periodic tangential displacement with an amplitude of 3 μm is applied in step 2 and circulates five times in a complete analysis step. Static generic analysis steps are used in the model.
Before performing mechanical analysis on the model, meshing is first required. A four-node plane strain element (CPE4) is used to discretize the model. 38 It is crucial to note that in order to obtain the contact data accurately, the mesh of the rough contact area needs to be refined. The mesh near the rough contact area is refined with a size of about 1.5 μm × 1.5 μm. In addition, the setting of contact attributes and interactions in finite element simulation is also essential. The contact form of the smooth surfaces is a positive point contact. However, the rough surface consists of many asperities, its contact form is multiple point contact. The contact form of a single asperity is either side contact or positive contact. In the contact properties, the tangential behavior is set to isotropic Coulomb friction. Lagrange multipliers are used to define tangential constraints. The hard contact method is used to define the normal behavior of the contact pair, and the master-slave interaction algorithm is adopted by contact discretization in the rough finite element model. Among them, the lower surface of the semicylindrical body is defined as master, and the top surface of the flat is defined as slave.
In the force analysis process in Figure 2, the bottom of the flat is fully constrained. The normal load P is applied on the top surface of the semi-cylindrical body to bring the two rough surfaces into contact. A reference point is set at the center of the top surface of the semi-cylindrical body. The reference point and the top surface are bound together through the multi-point constraint (MPC). A periodic horizontal displacement is applied at the reference point to make the two parts slide relative.

Calculation of fretting wear
The prerequisite for accurate calculation of fretting wear is to obtain a suitable numerical model. The Archard model is a classic wear calculation model; however, it is    the surfaces, which can be defined as where V is the wear volume, α is the wear coefficient, E d is the energy consumed by the wear, Q is the friction shear force, and S is the displacement amplitude of the sample. The wear depth of each node is calculated from the corresponding shear stress. In the process of a complete cyclic load, the wear depth at x is where h(x) is the wear depth, q(x) i is the shear force at the position x of the i incremental step, Δs(x) i is the relative slip distance at the location x in the i incremental step, and T is the number of incremental steps.
In the simulation of fretting wear, the amount of wear generated in one cycle is very small. Therefore, one incremental step is used instead of the cyclic load of ΔN times. This step not only ensures accuracy but also reduces simulation time. The acceleration of the wear simulation is realized. After acceleration, the wear depth at the location x is expressed as The number of accelerations ΔN used in the model is 1000, 18 and the wear coefficient α is 3.33 × 10 −8 MPa −1 . 4 In addition, variation of the surface profile leads to the change of the nodal stress and strain, which will affect the accuracy of the final simulation results. For this reason, the grid near the rough surface was set as an adaptive grid (ALE) to make it automatically adjust The UMESHMOTION subroutine is called at the end of each incremental step to obtain data and calculate the amount of wear. At the same time, the grid nodes are updated. Therefore, the accuracy of the wear is guaranteed by keeping the mesh in a proper shape at all times.

Validation of the model
The developed model was validated by comparing the finite element simulation results with the cylinder/plane wear experiment. 6 Here, for facilitating verification, surface roughness is not considered. The material used in the literature is high-strength alloy steel, and the material parameters are the same as the model in this paper. The mesh size of the contact area is 10 μm, the friction coefficient is 0.8, the amplitude of the displacement load is S = 25 μm, the normal load is F = 185 N, and the other parameters are consistent with the previous ones.
Before verifying the wear model, the model's accuracy before wear should be guaranteed first The contact pressure at the end of step 1 was obtained. The correctness of the model before wear was verified by comparing it with the analytical solution of the Hertzian stress distribution. Under a normal load F, the Hertzian contact halfwidth a is given as 19 and the stress distribution is expressed as where E* = [(1-μ 1 )/E 1 + (1-μ 2 )/E 2 ] −1 , R = (1/R 1 + 1/ R 2 ) −1 , P 0 = [FE*/(πR)] 1/2 , E 1 and E 2 represent the elastic modulus of the two materials, μ 1 and μ 2 represent the Poisson's ratio of the two materials, and R 1 and R 2 represent the radius of curvature of the contact surface of the two samples. For the cylinder/plane contact type, the radius of curvature R 2 of the plane is taken to be infinite in the above formula. Figure 3(a) shows the comparison between the finite element results and the Hertz theory. It can be seen that the finite element results agree with the theoretical results, which in turn verifies the accuracy of the model before wear occurs. Finite element results on both sides of the contact area slightly deviate from the transverse coordinates of the Hertz model, which is caused by the mesh size.
The wear depth after 18,000 cycles is compared with the experimental data in the literature, as shown in Figure 3(b). The wear width and maximum wear depth of the finite element model are consistent with the experimental results. The same UMESHMOTION subroutine and finite element model are used in the model with a rough surface and the previous verification model. Therefore, the rough model of fretting wear is reliable. From the finite element results, it can also be found that the contour of the contact surface after wear is smooth, which is different from the actual situation. The critical effect of rough surfaces on fretting wear has been emphasized again.

Effect of rough surface on wear
Smooth surface, single rough surface, and rough surface models were established to explore the influence of rough surface on fretting wear. In the single rough surface model, only the lower contact surface is rough. For the rough surface model, both contact surfaces are rough. The roughness of the rough surface is Ra = 0.2 μm. Other parameters are the same for all three models.
The distribution of contact pressure and stress of the three models are shown in Figure 4(a)-(c), respectively. It can be found that the distribution of contact pressure is discrete in both the single rough model and the rough model. The larger contact pressure is in the middle area. The distribution of contact pressure is continuous in the smooth model, and its maximum is at the edge of the contact area. This result is because the contact surface produces stress concentration at the locations where the surface profile changes, leading to a sudden change in the contact pressure. The contact pressures of the three models are shown in Figure 4(d). It can be seen that the value of the contact pressure decreases with the number of rough surfaces. The rough model has the largest contact pressure, and the smooth model has the smallest contact pressure, which means the rough model has the smallest contact area, and the smooth model has the largest contact area. Figure 5(a) shows the changes in the contact area of the three models. It can be seen that the contact area is inversely proportional to the number of rough surfaces. The contact area for the smooth model is much larger than other models, and it increases at the largest speed compared to other models when the normal load is applied. In addition, the contact area constantly increases as the wear progresses.
As shown in Figure 5(b), the value of shear stress oscillates between positive and negative values, which means that the direction of the shear stress on the asperities is not consistent. Asperities exist on rough surfaces, and their size and shape are different. The protruding asperities came into contact after the normal load was applied. At this time, some contacts are face-to-face contacts, and more are side contacts. The elastic action leads to a distribution of the shear stress oscillating between positive and negative values after the load is applied. For the smooth model, there are no asperities between the surfaces. However, the direction of the shear stress is also inconsistent due to the sudden changes in the contours at both ends of the contact area and the effect of elasticity. In addition, the shear stress of the smooth model is symmetrical, and its value in the middle region is zero, which indicates a partial slip state. From a macro point of view, the absolute value of the shear stress increases as the number of rough surfaces increases. Figure 5(c) shows the relative slip distance of the three models. It can be seen that the direction of relative slip is opposite at both ends of the contact area of the three models. The smooth model is in a mixed slip state because the middle of the contact area is sticky, and the two ends have a relative slip. It can be seen from the inset that the rough model has fewer asperities with relative slip than the single rough model. This result also shows that the rough model has a smaller contact area.
The contact is discontinuous when rough surfaces are present. The wear depth of the three models is shown in Figure 5(d). It can be found that the rough model has the greatest depth of wear compared with other models. The location with the greatest wear depth is the same as the single rough model. The reason is because the same fractal curve is used for all the rough surfaces, and the most prominent position of the contour is the same.
The wear depth in the middle area of the smooth model is zero because it is in a mixed slip state. In this state, the shear stress and relative slip distance of the smooth model are both zero in the middle area. Compared with the shear stress, the relative slip distance of the three models is of a smaller order of magnitude. According to the wear depth equation (6), it can be concluded that the rough model has a larger wear depth, which is consistent with the trend in Figure 5(d).

Effect of surface roughness on wear
For the rough model, the roughness of the lower contact surface is fixed at Ra = 0.2 μm. The roughness of the upper contact surface is Ra = 0.2 μm, Ra = 0.5 μm, and Ra = 0.8 μm, and other conditions are the same. The fretting wear under different surface roughness is analyzed, as shown in Figures 6-7. Figure 6(a)-(c) show the distribution of contact pressure and stress under different roughness, respectively, and the comparison of contact pressure is shown in Figure 6(d). It can be seen that the contact pressure increases as the roughness increases. Within the same contact length, the distribution of contact pressure is denser when the roughness is smaller. The reason is because the asperities have different distributions in the horizontal and vertical directions under different roughness. When the roughness increases, the height difference of the micro protrusions in the vertical direction increases. Therefore, the number of micro protrusions in contact decreases, resulting in reducing the contact area, as shown in Figure 7(a). When the contact area decreases, the contact pressure increases under the same normal load.
The contact pressure becomes larger as the roughness increases. At this time, the friction coefficient does not change, then the friction between the two surfaces will become larger. Therefore, the shear stress increases under the same displacement, as shown in Figure 7(b). Figure 7(c) shows the relative slip distance under different surface roughness. It can be seen that the relative slip distance also increases with the increase of roughness. The relative slip distance on both sides of the contact zone is greater compared to the middle region. It indicates that the wear is most severe at both ends of the contact area, exhibiting the same trend as shown in Figure 5(d). The comparison of wear depth under different roughness is shown in Figure 7(d). It can be seen that the wear depth increases with roughness. Compared with the middle area, the wear depth on both sides of the contact area is greater, consistent with the law of relative slip distance. The number of incremental steps T does not change because the number of wear cycles does not change. The shear stress and relative slip distance increase with the increase of roughness. Based on the wear equation (6) of the energy model, it can be inferred that the wear depth increases with increasing roughness, which is consistent with the simulation results. Therefore, the roughness between the contact surfaces of the connection structure should be reduced as much as possible to reduce the adverse effects of fretting wear.

Effect of materials on wear
In addition, the effect of materials on the fretting wear is systematically studied, as shown in Figures 8-9. The material of the lower specimen of the rough model is defined as alloy steel. The upper test pieces are defined as alloy steel, QT800-2, and HT200, respectively. Other conditions are the same. Figure 8(a)-(c) show the contact pressure and stress under different materials, respectively, and the comparison of contact pressure is shown in Figure 8(d). It can be found that the contact pressure shows the same trend with the change of Young's modulus of the upper specimen material, that is, the contact pressure increases with the increase of Young's modulus. The distribution of contact pressure does not change while its magnitude varies. The relationship between the contact area over time is shown in Figure 9(a). It can be seen that the contact area under the three materials fluctuates up and down, although the macroscopically continues to increase. Among them, HT200 has the largest contact area, and alloy steel has the smallest contact area. The above phenomenon occurs because the roughness and the contact position do not change. The magnitude of Young's modulus represents the stiffness of the material, and it affects the deformability of the material. The deformability of the material decreases as the rigidity of the material increases. Therefore, the contact area decreases as the rigidity of the material increases, and the contact pressure presents the opposite trend. Figure 9(b) shows the relationship between the shear stress and the material. The distribution of shear stress does not change under the same rough profile; however, the position of the maximum shear stress is changed. The maximum shear stress of alloy steel is near position x = 0.05, QT800-2 is near x = 0.22, and HT200 is near x = 0.45. In addition, the maximum shear stress increases with Young's modulus. The change of the shear stress is due to the different resistance to deformation of the material and the change of the contact area. Materials that are more susceptible to deformation can promote the occurrence of relative slip. Therefore, the relative slip distance is the largest when the material is HT200, as shown in Figure 9(c). The shear stress determines the size of the wear depth when the material is the same. Therefore, the wear depth is the largest when the material is alloy steel, as shown in Figure 9(d). In practical applications, materials with a smaller Young's modulus should be selected to slow down the effect of fretting wear.

Effect of the friction coefficient on wear
The coefficient of friction is another critical parameter that affects fretting wear. The friction coefficients of the rough models were set to μ = 0.4, μ = 0.6, and μ = 0.8, respectively, while the other parameters concerning roughness and material were the same. The wear under different friction coefficients is analyzed, as shown in Figure 10 and Figure 11. Figure 10(a)-(c) show the contact pressure and stress under different friction coefficients, respectively, and the comparison of contact pressure is shown in Figure 10(d). It can be seen that the distribution of contact pressure does not change for different friction coefficients. The size of the contact pressure changes slightly with the change of the friction coefficient. Besides, many asperities of varying sizes make up the rough surface. Fretting occurs on partial, discontinuous, and contacting asperities. This effect makes the actual contact area very small. The contact area changes slightly with the change of the friction coefficient, as shown in Figure 11(a). Small changes in the total contact area are allocated to each asperity. As a result, the change in the contact area of each asperity will be smaller. Therefore, the contact pressure changes little. Figure 11(b) shows the shear stress under different friction coefficients. It can be found that the shear stress decreases with the decrease of the friction coefficient. The location of the maximum shear stress is changed. The increase in the friction coefficient means that the friction between the surfaces becomes larger. As a result, some contact parts change from a slipping state to an adhesive state, and other contact parts are affected. Therefore, the change of the friction coefficient greatly influences the shear stress.
The increase in friction prevents relative slip between two contact surfaces. Therefore, the relative slip distance decreases as the friction coefficient increases, as shown in Figure 11(c). According to the wear depth equation (6) of the energy model, it can be known that the wear depth is reduced. Figure 11(d) shows the wear depth under different friction coefficients, which is consistent with the theoretical results. In addition, the position of the maximum wear depth also changes with the change of the friction coefficient because the position of the maximum shear stress has changed.

Conclusion
A finite element model of rough contact was developed to study the fretting wear of mechanical connections better. The stress results of the rough surface model are compared with other models. The influence of roughness, material parameters, and the friction coefficient on the fretting wear characteristics under rough contact is investigated. The main conclusions are as follows.
1. The rough model can simulate real fretting wear and calculate wear more accurately. The wear depth and contact pressure increase with the number of rough surfaces. However, the contact area decreases with the increase in the number of rough surfaces in the model.  Hertz contact pressure (MPa) q(x) i shear force at the position x of the i incremental step Q friction shear force R equivalent radius of curvature (mm) R 1 radius of curvature of the upper specimen (mm) R 2 radius of curvature of the lower specimen (mm) Ra surface roughness (μm) Δs(x) i relative slip distance at the location x in the i incremental step S displacement amplitude of the sample (μm) T number of incremental steps V wear volume Z(x) height of fractal curve (mm) Α wear coefficient Γ fixed value larger than 1 (usually 1.5) γ n spatial frequency λ n surface wavelength μ coefficient of friction μ 1 Poisson's ratio of the upper specimen μ 2 Poisson's ratio of the lower specimen