Numerical analysis of time-varying wear with elastic deformation in line contact

Wear is an important factor for failures of mechanical components. Current research on wear is mainly focused on experiments while the numerical simulation of wear is hardly used owing to the complexities of the wear process. Explaining the effect of friction on the wear process is important, as it will lead to a deeper understanding of the evolution of wear. This study proposed a numerical method to expound the wear process in the contact between an elastic cylinder and a half-space simulating the ring-block tester. There are two difficulties during the calculation; one is that the contact shapes vary with time, causing the pressure distribution to change simultaneously and the other is the integral equation for calculating the contact pressure under different worn shapes. In the present study, the wear rate was computed using Archard’s law and the wear process was calculated step by step until the specified total sliding distance was achieved. During each step of the calculation, the contact topography was updated. The simulation intuitively reproduced the contact state of change from line to surface contact throughout the wear process. Reasonable agreements on the changes of the wear scar, achieved from experiments and numerical simulations, were obtained.


Introduction
Wear is a continuous process of material loss from friction pairs during contact owing to relative motion. It has been reported that wear failure encompasses roughly 60%−80% of the total failures of mechanical pairs and components [1]. Factors such as temperature, elastoplastic deformations, surface topographies, material properties, and chemistry all contribute to the complex contact conditions. An experimental approach is necessary because of the complexity of the system. Although experimental approaches have provided results consistent with actual situations, it requires considerable time and capital to obtain in situ wear information. Therefore, a numerical prediction method is a good complement to predict the wear evolution.
Friction pairs can be divided into the higher and lower pair according to the contact state. Usually, a pair with surface contact is called as the lower pair because of its lower contact pressure. For a lower pair, the influence of wear on its life expectancy is not significant while for a higher pair that is in the line or the point contact, wear is significant owing to the higher local contact pressure. Numerical methods have been used by many researchers to study wear. Ahmer et al. [2], Hegadekatte et al. [3], and Bortoleto et al. [4] studied wear evolution during point contact based on the pin-on-disc model. Flodin and Andersson [5,6] simulated the wear evolution of spur and helical gears under dry friction conditions. Brauer and Andersson [7] and Hao and Meng [8] studied the wear process in mixed lubrication conditions. There are many formulas to simulate material wear for different working conditions [9]. The Archard's wear model is well known and widely used because of its relatively simple form and consistency with actual situations. Its main idea is that the wear rate is proportional to the normal contact pressure and the emphasis is on determining the pressure distribution between the contact regions. The finite element method (FEM), with the advantage of easy implementation of numerous physical effects, is one of the most popular tools to solve a contact problem [10−12]. The primary drawbacks of these models are that they are inefficient and time-consuming to improve the accuracy of calculation. The Winkler surface model [5,13] and the boundary element method [14] can be applied to speed up the wear simulation. Furthermore, Andersson et al. [15] has used the discrete convolution and FFT method to simulate the wear of a sphere on flat contact pair.
In the present study, the second type of singular integral equation [16,17] was used to compute the contact pressure distribution under various contact topographies. The plastic deformation of asperity was not taken into consideration as this effect was considered indirectly through the use of an experimental value of wear (constant k). Two difficult problems in calculating a wear process were solved. One is the pressure distribution, which influences wear. The contact pressure distribution varies with time because the shapes of the contact surfaces are continuously changed by wearing. The other is the singularity when a numerical method was used to solve the second type of singular integral equation [16]. In this paper, a new calculation model is introduced to obtain the numerical solution. The simulation results were validated through ring-block tester experiments.

Basic equations
A line contact problem can be simplified as the contact between a cylinder and a half-infinite plane. Here, a ring-block friction and wear experimental model is used to illustrate the wear processes (Fig. 1). The ring rotates at a fixed speed under the load L.

Wear calculation
In the relative motion of a friction pair, the wear rate is the function of the normal contact pressure p(x) according to Archard's law. It can be written as follows [18]: where h is the wear depth; s is the relative sliding distance; k is the wear constant, Pa −m ; x is the coordinate and m is the influence index of pressure. m equals to 1 in most cases. The wear constant can be obtained experimentally or numerically. In this study, the influence of temperature is not taken into account, i.e., a linear plot of wear volume as a function of time was considered. This means that the wear rate is constant during sliding.
With discretizing the contact region in the time and space domain, the total number of nodes is N on the contact interface and the total number of time steps is M. If a node is denoted as i and the time step as j, the corresponding contact pressure is ( , ) p i j and the wear depth increment is Δ ( , ) h i j . The total wear depth at the i th node before the j th time step is ( , ) h i j . If the increment of the relative sliding distance between two time steps is s, the relationships between the wear depth ( , ) h i j and the increment of relative sliding distance s are given by The contact pressure distribution should be obtained first in the contact region for each time step to obtain the wear depth increment within the relative sliding distance increment s.

Calculation of contact pressure distribution
Both normal and tangential forces exist over the contact interface when two elastic bodies are in contact. Figure 2 shows the contact between an elastic cylinder  and a half-infinite plane. L is the normal load per unit length, Q is the tangential traction due to friction and n is the rotation speed. E 1 and E 2 and ν 1 and ν 2 are the Young's moduli and the Poisson's ratios of contact bodies, respectively. The initial contact point is set as the origin of the coordinate. a 1 and a 2 (a 1 > 0, a 2 > 0) are the left and right contact width, respectively.
According to the theory of elasticity [16], the relationship between the contact pressure and surface displacement gradient in the interval 1 2 ( ) a x a    can be written as follows: where 1 z u and 2 z u are the surface displacement gradients in the z direction. p and q denote the normal and tangential tractions, respectively.
The normal and tangential tractions on the two elastic bodies are equal in magnitude, but opposite in direction. The two equations in Eq. (3) are summed.
where  is the coefficient of friction and G is shear modulus, a simplified form can be obtained: , and  is a measure of the different elastic constants.
Further, the contact pressure should satisfy the load equilibrium condition, i.e., where L is the normal load per unit length. If the surface displacement gradient is given, the contact pressure distribution in each time step can be solved by using Eqs. (4) and (5).

Calculation of surface displacement gradient
When two elastic bodies are pressed together, the contact of the two bodies can be equivalently represented as the contact between an elastic half-space with the combined elastic modulus, E * , and a rigid cylinder with combined radius R. The equivalent radius and modulus of elasticity are given by In Fig. 3, the contact between a cylinder and an elastic half-space is shown. The surface displacement z ( , ) u i j in the contact region at the i th node on the j th time step is approximated as The surface displacement gradient is obtained by the middle difference method. There is an unknown constant d in Eq. (7), but it does not affect the final  To sum up, the numerical procedure used for wear simulation in this study is as follows. First, the surface displacement gradient is obtained gradually according to the contact topography based on the ring-block test. Subsequently, the contact pressure is computed using the integral equation, which also determined the contact width automatically under the condition of nonnegative normal contact pressure over the contact region. Finally, once the contact pressure on the contact region is determined, Archard's wear law is used to compute the wear depth increment. The aforementioned processes are conducted gradually until the total sliding distance is achieved. In addition, the wear scar and the contact pressure are obtained at the end of each step.

Numerical procedure
In this section, the numerical procedure of wear calculation is described. From the previous analysis, it is vital to solve for the contact pressure correctly during wear simulation. Therefore, the numerical method to obtained the normal contact pressure is explained separately.

Numerical method for solving the contact pressure
The pressure distribution on the contact interface is obtained by solving Eqs. (4) and (5) [17]. As we know, Eq. (4) is the second type of singular integral equation with a Cauchy kernel. Although analytical solutions for singular integral equations are presented in Ref. [19], there are certain problems in using these equations for arbitrary contact shapes, which are hardly described by functions. The finite difference method has been widely used to solve these types of equations. Nevertheless, this method is a local algorithm and the results do not convergence. Refining the mesh will improve the accuracy at the expense of increasing computation time. In this study, the piecewise continuous function method [20] was used to obtain the numerical solution of Eq. (4). The contact area should be normalized first, i.e., supposing     (4) and (5) can be rewritten as follows: According to Ref. [17], the solution of Eq. (8) may be expressed as the product of a bounded function, ( ) g X , and a fundamental function, ( ) w X , i.e., ( ) X   ( ) ( ) w X g X . The fundamental function is given by Therefore, the second term of first formula in Eq. (8) is rewritten as follows: The last integral in Eq. (10) may be solved according to the known result [20]: is a known function. In special cases, 1      , ( ) X  equals to zero. Then the first formula of Eq. (8) converts to The Cauchy singularity is removed because ( ( ) g T  ( ))/( ) g X T X  equals to the derivative of function ( ) g X when T = X. The solution of pressure distribution finally attributes to the solution of ( ) g X from the above analysis. The Gauss-Jacobi quadrature rule is an appropriate choice to solve Eq. (13) for the Jacobi weight function indicated by Eq. (9). At the collocation point i X , Eq. (13) will be discretized as where N is the total number of the integration points, are the weights and abscissas of the standard Gauss-Jacobi quadrature, respectively.
Generally, the selection of the collocation points The values of ( ) g T are obtained by solving the Eqs. (14) and (15). The pressure distribution is obtained by calculating ( ) ( ) ( ) T w T g T   . However, the real contact width is unknown at the beginning, i.e., 1 a and 2 a are unknown parameters when using Eqs. (14) and (15) to compute the contact pressure distribution. In this study, the two unknown parameters are initially assigned two large enough values to ensure they contain a real contact region. Then the redundant contact areas are removed when pressure is less than zero. Finally, the contact region is determined automatically under the condition of non-negative pressure at any node in the contact region through multiple iterations.
It is important that this method does not include the boundary points. In other words, the pressure at these points cannot be obtained directly. It is reported that the pressure distribution between two elastic bodies, whose profiles are continuous through the boundary of the contact area, falls continuously to zero at the boundary according to Johnson's research [16]. Thus, the pressure at the boundary points is set as zero. Figure 4(a) shows a comparison of ( ) p x between the Hertz contact and corresponding numerical results under different friction coefficients, without wear, using the aforementioned numerical method. It is shown that the analytical solution agrees with the numerical solution when the friction coefficient equals to zero. However, the contact region will shift to the left with respect to the initial contact point when 0   Fig. 4(b) are computed according to the method presenting at Ref. [22]. It is shown that the numerical results are in good agreement with the theoretical results. It is shown that the contact width increases monotonically with friction coefficients. The greater the friction coefficient is, the more obvious the deviation of the contact region becomes. The right contact width is only about 69.9% to the left when 1   , which rarely exceeds 1.0 in engineering fields.

Wear calculation
The simulation model presented in this paper is based on the ring-block tester. The wear depth of the ring is negligible compared to the block because the latter shows uniform wear along the circumference. The contact width increases and the contact state is transformed from line to surface contact as well, owing to wear. However, the abscissa of the standard Gauss-Jacobi quadrature is deterministic while the wear widths differ from each other. This will lead to discrete nodes over the contact region that do not coincide with each iteration, as shown in Fig. 5. As a result, the wear data ( , ) h i j cannot be used directly when we calculate the surface displacement at the (j+1) th time step using Eq. (7). To overcome this problem, the following method is introduced. First, the discrete points on the j th time step are mapped onto the (j+1) th time step. Then the wear depth of the corresponding discrete points was obtained by spline interpolation inside the region while on the outside, the wear depth is zero. The surface displacement can be computed using Eq. (7) if the wear depth on the (j+1) th time step is known. Finally, the pressure distribution is obtained from the numerical method described in Section 3.1, and wear depth is obtained from Eq. (2).

Experiment
In this study, the ring and block specimens were used as friction pairs. The materials used in the experiment are AISI 1045 steel and 0Cr18Ni9, which correspond to the block and ring, respectively. The ring surface and one of the block end surfaces composed the sliding line contact. The line-contact tests were performed on a standard friction tester (MRH-3, Yihua Tribology Testing Technology, China). The width and outer diameters of the ring were 13.06 and 49.22 mm, respectively. The length and width of the block were 19.05 and 12.32 mm, respectively. The experimental conditions are presented in Table 1. To omit the influence of temperature on the wear constant and friction coefficient, the sliding velocity should be small.
The variation of the friction coefficient with time is shown in Fig. 6. The friction coefficient was at a constant value of about 0.13 after running-in. This friction coefficient value will be used in the wear simulation process.
The worn scars of the block were measured using a three-dimensional surface profilometer (Talysurf CLI 1000, Taylor Hobson, UK). The wear constant k in Archard's law was experimentally determined. The wear constant k, 9 1 1.245 10 MPa k     , was computed according to the experiment.

Evolution of contact variables
The aforementioned approach was used under the loading conditions given in Table 1. The coefficient of friction and the wear constant were given in Section 4 based on experiments and these values were employed in the wear simulations. Archard's wear law was applied at each sliding step Δs , 500 μm in the simulation procedure, to obtain the wear depth increment Δh . The contact pressure distribution evolves concomitantly with the changes in contact width. Figure 7 shows the numerically predicted evolution of the contact pressure distribution with increasing sliding distance. The peak pressure continuously decreases and finally moves toward uniformity. More details marked by a circle in Fig. 7 indicate that the contact pressure sharply increases at the edge of the contact region. This phenomenon is consistent with the results obtained from the finite element method [8,23]. However, the acting area of the spike pressure is very small. In this study, to simplify the contact problem and improve the efficiency of wear simulation, the pressure spikes at the edges of the contact region are neglected and the pressure is evenly distributed along the contact width when the wear width is 2.5 times the initial contact width.
The variation in the contact width and peak pressure with an increasing sliding distance is demonstrated in Fig. 8. There is a dramatic reduction in peak pressure during the first 20 m, followed by a much more gradual reduction over the subsequent 380 m. After sliding 400 m, the peak pressures reduced to less  than 10% of their initial value. This phenomenon is mainly caused by the change of contact state from line to surface contact. The predicted width also obviously changed as wear continued. There was a rapid increase in the contact width over the first 20 m, followed by a gradual reduction in the rate of the increase. For subsequent sliding, the contact width continued to increase but at a slower rate.

Predicted wear profiles
The numerically predicted worn surface profiles of the block specimens with regard to increasing sliding distances are shown in Fig. 9(a). It is found that as  | https://mc03.manuscriptcentral.com/friction wear progresses, a worn scar develops on the flat surface. The contact moves towards conformity, with a similar radius to the ring. The evolution of the contact pressure in the simulation is shown in Fig. 9(b). It indicates that the maximum contact pressure decreases gradually along with an increasing the sliding distance.
The numerically predicted worn surface profiles of the flat specimens for sliding distances of 75 m, 150 m, 225 m, and 300 m were compared with experimental results and are presented in Fig. 10. It should be noted that the worn profiles in Fig. 10 are the mean values of the three-dimensional scanning contours at each sliding distance. The predicted values of the scar width and the maximum wear depth, together with the corresponding experimental results, are presented in Table 2. For the 75-m case, the scar width is underpredicted by 16%, while for the other three cases, the predicted and experimental results correlate closely. For the 150 and 300-m cases, the maximum wear depths are over-predicted by 39.3% and 18.2%, respectively, while for the other two cases, the experimental results are in agreement with the predicted values.

Effects of friction coefficient on the wear process
As shown in Section 3.1, the contact region will shift with respect to the initial contact point because of the friction coefficient. The degree of deviation depends on the  , a measure for the different elastic constants, and friction coefficient  . The influence of the friction coefficient on the wear process in Section 5.2 is indistinguishable as the material parameters of friction pairs are very similar. The ring is assumed as a rigid body to explain the effects of the friction coefficient on the wear process and the rest of the simulation parameters are the same as in Section 5.1. The evolution of the ratio of right to left contact width with increasing sliding distance is shown in Fig. 11. The contact region is always symmetrical around the initial contact point when the effect of friction force on contact pressure distribution is ignored. However, although the contact area is asymmetric about the initial contact point in the first 100 m, the left and right contact widths are approximately equal eventually for subsequent sliding when 0.5   . . The ring is assumed as a rigid body and the other simulation parameters are the same as in Section 5.1.

Conclusions
In this paper, a numerical method was proposed to analyze mild wear in the initial line contact with elastic deformation. Reasonable agreements for the wear scar between the simulation results and the experiment results were obtained. In addition, the effects of the friction force on the wear procedure were numerically investigated. The main results and conclusions are as follows: 1. The exact contact pressure distribution under different contact topographies with various coefficients of friction was obtained by solving singular integral equations of the second type.
2. The contact region was offset with respect to the initial contact point when taking friction into consideration. The results showed that the deviation in the left and right contact widths depends on  , a measure of the different elastic constants, and friction coefficient  . Nevertheless, the maximum contact pressure and contact width were almost constant for   1, which rarely exceeds 1.0 in engineering fields.
Further research about the relationship between friction and wear is necessary for higher friction coefficients.
3. Line contact was a temporary state during the wear process. The contact state rapidly changed to surface contact from line contact. The maximum contact pressure decreased drastically at the beginning and presented a linear decrease a moment later. The contact region will finally be symmetrical to the initial contact point owing to wear when we take friction into account. However, the contact region always maintains symmetry about the initial contact point when 0   .