A Numerical Model for Solving Three-Dimensional Rolling Contact Problems with Elastic Coating Layers

In this work, a numerical model is proposed for three-dimensional rolling contact problems with one or two elastic layers, and the tangential contact solution is emphasized. Previous works on this topic have mostly been two-dimensional, in which only longitudinal creepage has been considered. With the three-dimensional model presented in this work, all possible creepages, such as the longitudinal, lateral and spin creepages are taken into account. In order to improve the calculation efficiency, the conjugate gradient method and the FFT technique are employed. The influence coefficients for displacement and stress are obtained from the corresponding frequency response functions. The numerical results are validated against existing results and good agreement can be found. The effects of the different layers’ thicknesses and elastic moduli under different creepage combinations on the traction distribution and stick/slip results are investigated. It can be seen that by adjusting the layer parameters the traction and stick/slip results can be modified significantly, and it may, therefore, be very useful information for improving the rolling contact fatigue and mitigating wear problems in various mechanical systems.

The corresponding influence coefficient matrices for pressure and displacement C u i q The corresponding influence coefficient matrices for traction and displacement P The normal force (N) d x , d y The mesh sizes with respect to the x-and y-axis S x , S y The relative slip velocities in the x-and y-direction x The creepage ratio in rolling direction y The creepage ratio in lateral direction The spin creepage ratio (rad/ mm) q x , q y The traction results in x and y direction (mm) The friction coefficient

Introduction
Rolling contacts are widely found in mechanical systems, e.g. wheel and rail contact interfaces, cam and roller systems, gears, and rolling element bearings. The calculation of traction distribution and the stick/slip zones for rolling contact problems is significant for understanding traction control, lifespan prediction, wear and rolling contact fatigue. A lot of research has been conducted about this topic, and some sophisticated numerical models have been proposed by several researchers, such as CONTACT and FASTSIM algorithms proposed by Kalker [1] and the revised algorithms proposed by Vollebregt [2], the finite element method by Li [3,4], the boundary element method by Abascal [5,6], semi-analytical method by Wang [7] and the linear complementarity problem formulations by Xi [8,9]. In order to improve and to adjust the contact interface behaviors without changing the bulk material, thin layers, usually in the form of surface coatings are often used. For example, the ceramic coatings, alloy/metal coatings and composite coatings [10], are often utilized. Moreover, sometimes the layers presented within the contact interfaces are unavoidable. For example, different third layers, such as leaves, dust and dirt, are often involved in wheel and rail contact problems [11].
To solve the rolling contact problems with elastic layers, many numerical studies have been proposed to obtain the traction distribution and the stick/slip zones. Among them, Meierhofer et al. proposed a model for calculating the traction characteristic for the wheel and rail contact with a third body layer [12]. The normal contact solution was simplified using the Hertz solution. A Winkler similar method was used for the tangential contact solution, and the layer deformation was assumed to be described by the independent brushes' deformation. Guler et al. [13][14][15] focused on the solution of rolling contact problems with graded coatings. In their work, the relationship between the surface displacement gradients and the traction was derived by several integral equations, and the Goodman approximation was adopted to decouple the deriving singular integral equations. The solution of the integral equations was conducted using the Gauss-Chebyshev integration method. A similar method was also used in a recent work of Alinia et al. to analyze the stress for rolling contact problem with functionally graded material coating [16].
Most of the reports mentioned above are focused on twodimensional rolling contact problems involving coating layers only. However, it is surprising that so little work has been done for three-dimensional ones with elastic coating layers. Here, it should be noted that, the three-dimensional rolling contact problems under viscoelastic layer conditions are considered by some researchers [17][18][19]. For viscoelastic materials, the imperfect elasticity plays a significant role for rolling traction force, that is different from the origin of elastic rolling traction forces. In our topic, for rolling contact problems of elastic materials, as important factors, creepages indicate the velocity differences within the rolling contact. For three-dimensional rolling contact problems, there are three kinds of creepages as per the different directions. Especially, the longitudinal creepage denotes the velocity difference in the rolling direction, the lateral creepage indicates the velocity difference in the lateral direction and the spin creepage denotes the angular velocity difference along the normal direction of the contact surfaces. Lateral and spin creepages are normally involved besides the longitudinal one, and they are also important for mechanical performance. For example, when wheels are negotiating with curved rails, squeal noise is often generated as the result of the lateral force [20,21]. Moreover, asymmetrical wear profiles can also be found for the curved track because of the spin creepage [22,23]. To understand and to explain these behaviors, it is required to take the effects of lateral and spin creepages into account. While for two-dimensional rolling contact conditions, only longitudinal creepage, or in other words slide-to-roll ratio (SRR), should be considered.
In this work, a three-dimensional rolling contact numerical model for point contact with one or two elastic layers are considered. The influence coefficients of the stress and the displacement for the contact with layers are obtained from the frequency response functions using the fast Fourier transform method. For calculating the normal contact and the tangential contact of the rolling contact problems, the conjugate gradient method is utilized during the iteration processes. Moreover, to speed up the calculation, the FFT technique is also employed. All possible creepages, i.e., the longitudinal, lateral and spin creepages are taken into account in the current rolling contact model. The results of traction distribution and the stick/slip zones are presented with various coating elastic modulus and thickness under different creepage combinations.

Problem Formulation
The schematic diagram of the current three-dimensional rolling contact problem is shown in Fig. 1, and , and denote the angular velocities about x, y and z axes, respectively. In this work, each layer is homogeneous and linear elastic, and the interfaces are bounded perfectly. The elastic modulus for layer 1 and layer 2 are E 1 and E 2 , respectively and v 1 and v 2 and h 1 and h 2 are the corresponding Poisson's ratios and the layer thicknesses.

The Normal Contact Solution
For the normal contact solution, the complementarity conditions for the normal gap g z and the corresponding normal pressure p are employed, and they can be given by The normal gap g z can be described by In which h z and z are the initial separation and the rigid displacement between the contact surfaces, while u z is the displacement and can be calculated by where C u z p is the corresponding influence coefficient matrices (ICMs) and the derivation will be provided later. The normal force P can be calculated by Here, d x and d y are the mesh sizes with respect to the xand y-axis, respectively.

The Tangential Contact Solution
To solve the tangential contact problem, firstly the governing equations for three-dimensional rolling contact problems are given by [1] In which S x and S y are the relative slip velocities in the x-and y-direction, and x , y and denote the longitudinal, lateral and spin creepages, respectively, while u x and u y are the corresponding deformations in the x-and y-direction.
An important characteristic, for rolling contact problems is that there are both stick and slip zones within the contact region. As per the report of Wang [7], mathematically, the basis to determine the stick and slip zones in this work can be expressed by For the slip zone. where q x and q y are the traction results in x and y direction, is the friction coefficient and p is the corresponding normal pressure.
In order to solve the current rolling contact model efficiently, the method proposed by Wang et al. [7] is adopted in this work, in which the conjugate gradient method is employed during the iteration processes. However, it should be noted that in Wang's work the lateral and spin creepages are not considered, and in this work, the effects of these creepages on the traction distribution and the stick/slip zones are investigated. Therefore, based on Wang's work, some necessary processes are adjusted to consider the lateral and spin creepages. To guarantee the completeness of the content in this work, the main methodology is also provided.
If there is no slip in the whole contact region, both S x and S y are zero. Therefore, the governing equation for 3D rolling contact problems can be expressed by By integrating both sides with respect to x, we have In which f (y) and g(y) are functions of y , after discretization they can be given by To obtain the results of u x ij and u y ij , the ICMs to relate the deformation with the traction results are required. In contrast from the conditions without layers, it is hard to obtain the usable ICMs directly from the analytical solution, while for the layered condition, the ICMs can be derived from the corresponding frequency response functions [24]. In the work of Liu [25], the solution details for derivation of the ICMs with coating layers has been provided. Therefore, u x and u y can be calculated by .
In which IFFT denotes inverse Fourier transform, C u x q x indicates the influence coefficients to describe the relationship between traction q x and deformation u x , and similar patterns for C u x q y , C u x p , C u y q x , C u y q y and C u y p , while the symbol ^ is the discrete Fourier transform and the symbol * denotes the convolution calculation. For more details on the derivation of the ICMs, the reader is referred to the works of Liu, Cai and Sullivan [25][26][27].
In term of the Papkovich-Neuber potentials φ and ( 1 , 2 , 3 ) , stress and displacement u can be given by In which k denotes the kth layer, and by double Fourier transform for x and y the stress and deformation can be expressed by In which both symbol ≈ and FT xy denote double Fourier transform, and the solutions can be expressed by are unknown coefficients needed to be solved.
There are several boundary conditions needed to acquire the solution, for the surface they can be given by For the layer interfaces the boundary conditions can be given by It can be seen that for the single layer condition there are 9 unknowns, while for the double layer conditions 15 unknowns are required to be solved. By using these equa- can be solved, and the results can be found in [26,27] for one and two layers.
When u x and u y are obtained, the residual vectors r x and r y can be expressed by To obtain the sum of squares of r x and r y in the stick zones, we have If r n r n−1 < 1 , the conjugate directions in both stick and slip zones can be calculated by In which t n x and t n y are the conjugate directions for the nth iteration, and if r n r n−1 ≥ 1 the old conjugate directions are employed.
The step length for the nth iteration α n can be calculated by The traction distribution q x and q y can be given by r n In which is the friction coefficient. If the convergence condition is arrived, the results of q x and q y can be output. Otherwise, the iteration from Eq. 9 should be conducted again until the convergence condition is reached. A flowchart is provided in Fig. 2.

Validation
To validate the normal pressure results from the current numerical model, dimensionless results of normal pressure against contact radius along the x-axis are provided in Fig. 3. In the following results, the sphere radius is set as 50 mm, the elastic modulus is 210 GPa and the Poisson ratio is 0.3. The normal load is set as 100 N. The potential contact area is selected as 2 × 2mm 2 , and 64 × 64 mesh density is selected.
As shown in Fig. 3a, the results with one layer are given. The layer thickness parameters are varied from 0.25a to a,  Fig. 3b. The present results for a single layer with a thickness of a and the results of two layers, both with a thickness of 0.5a, are in good agreement with the results in Yu's work [28]. Moreover, it can be seen that for layered conditions, the contact radius is increasing with the decreasing elastic modulus of layers. For compliant layers, i.e., the layer elastic modulus is smaller than the substrate material, the thicker the layer thickness, the larger the contact radius. Due to the larger contact area for compliant layer conditions, the maximum contact pressure results are much smaller by comparing with the hard one or no layer conditions. While for stiff layers, i.e. larger layer elastic modulus, the contact radius is decreasing with increasing layer thickness.
In order to validate the current results for traction distribution, the Linear Complementarity Problem (LCP) method for solving rolling contact problems is adopted [9]. Due to the high computation cost limitation of the LCP method, a coarse mesh, i.e., 22 × 22 , is adopted. By the application of the FFT technique a much finer mesh may be used, and for the numerical solutions presented here, a 64 × 64 mesh was considered to yield representative results. Similar to the study of Nyqvist [29], for the depth direction the resolution of z can be coarser and set as 0.1a if the substrate stress results are required.
In this work, the relative tolerance for the normal contact is set as 1e-5, and for the tangential it is set as 1e-7. The relaxation factor for the normal contact is initiated from 64 and is decreasing with iterations but not smaller than 0.5. All the numerical simulations were performed on a computer with an i7-7820HQ CPU (2.9 GHz) and 16 G RAM, and the running time for each case is normally less than 6 s. Different creepages combinations, including longitudinal, lateral and spin creepages, are employed. It should be noted that in this work the creepage parameters are seen as input parameters, and for each case the values are fixed in all iteration. As shown in Fig. 4a and b, pure longitudinal creepage condition, i.e., x = 0.002 , is compared. In Fig. 4c and d, both longitudinal and lateral creepages are applied, specifically, x = 0.002 and y = 0.001 . While for Fig. 4e and f, spin creepage is involved, i.e., x = 0.002 and = 0.01 . Lastly, as given in Fig. 4g and h, all possible creepages, i.e. longitudinal, lateral and spin creepages are employed in the numerical simulations, i.e., x = 0.002 , y = 0.001 and = 0.01 . The slip zones are indicated by the filled circles on the corresponding grids. From the results presented in Fig. 4, it can be seen that the present results for all conditions agree well with the ones obtained with the LCP method.

Traction Distribution Results with Layers Under Pure Longitudinal Creepage
The traction distribution results under pure longitudinal creepage are provided in this section. Firstly, the traction distribution results along the x-axis with a single layer are given in Fig. 5, in which the layer thicknesses are varied from 0.25a to a and the elastic moduli are selected as 0.5 E s and 2 E s . The longitudinal creepages from 0.02 to 0.1 are applied for each case. From Fig. 5a to 5c, it can be seen that for the compliant layer conditions the ratio of slip zones to the whole contact zones is increasing with increasing layer thickness, while the ratio of slip zones is decreasing with increasing layer thickness for the stiff layer conditions as shown from Fig. 5d to 5f. By comparing the compliant layer results with the stiff layer ones, it can be seen that the slip zones are more predominated for hard ones with the same longitudinal creepage conditions. In Fig. 6, the traction results with two layers under pure longitudinal creepage conditions are provided. Similar to the results for a single layer, larger slip zone ratios can also be found for the results of two hard layers, and smaller ones for the results of two compliant layers. Moreover, for the compliant layer conditions, slip zones also show more aggressive behavior for thick layers, while for the hard ones, more slip zones can be found for thinner layers.
From the results of Fig. 5 and 6, it can be seen that the longitudinal creepage is a significant factor to determine the traction distribution and the stick/slip zones. Moreover, by employing suitable layer parameters the desired traction distribution and stick/slip results can be obtained. These derived behaviors agree well with the conclusions proposed by Guler et al. [13,15], in their works the traction distribution results for rolling contact problems are provided, but it should be noted that the graded elastic coatings are employed by them, which can be seen as a combination of multi-layered elastic coatings [30].

Traction distribution results with layers under longitudinal, lateral and/or spin creepages
The main purpose of this work is to provide the traction distribution and the stick/slip zones results with layers under different creepage conditions. In contrast to the results under pure longitudinal creepages, the results with lateral and spin creepages the traction distribution results are provided for the whole contact area. The mesh density parameters are the same as used in the previous section. The traction distribution results for a single layer are shown in Fig. 7 and Fig. 8, while in Fig. 9 and Fig. 10 the results of two layers are given. In order to show the input parameters more clearly for the following results from Fig. 7 to Fig. 10, the elastic modulus and the thickness of each layer for all the numerical simulations are summarized in Table 1, while the creepage parameters are provided in Table 2. Figure 7a to 7f show the traction distribution results with a single layer, and the layer thickness is 0.5a and the corresponding elastic modulus is 0.5E s . The results for the creepage parameters are x = 0.002 , y = 0.001 and x = 0.002 , y = 0.002 , can be found in Fig. 6a and 6b, respectively. The traction results are almost symmetric with respect to the x-axis, and the slip zones are located at the trailing edges. For larger lateral creepage, the slip zones become larger. The results with the combination of longitudinal and spin creepages are given in Fig. 7c and 7d, in which x = 0.002 , = 0.005 and x = 0.002 , = 0.01 are employed separately. Moreover, Fig. 7e and 7f show the results when the longitudinal, lateral and spin creepages, i.e. x = 0.002 , y = 0.001 , = 0.005 and x = 0.002 , y = 0.001 , = 0.01 are applied. Due to the involving spin creepage, the symmetrical behavior of the traction distribution is ruined. By comparing with the results of Fig. 7c with 7d, more stick zones can be found for smaller spin creepage condition, and the slip zones are located at the top part of the trailing end. Similar behavior can also be found in Fig. 7e and 7f. The spin creepage can be seen as the combined effects of longitudinal creepage and lateral creepage, and they are not fixed values but dependent on specific location. In the current conditions, at the area of x < 0 and y > 0 , the spin creepage can enhance the traction and therefore to make more grids to reach the corresponding traction limits. While in the area of x < 0 and y < 0 , the effects of the longitudinal and the lateral creepages on the traction are neutralized partly by the spin creepage, and some original slip zones are taken over by stick domain.
As shown in Fig. 7g to 7l, a stiff layer, i.e. E 1 = 2E s , with 0.5a thickness is adopted. By comparing the results with the ones as provided in Fig. 7a to 7f, similar patterns for each case can be found, and it also can be seen that the contact areas are significantly smaller than the ones with a compliant layer, while for the same creepage conditions the slip zones account for larger ratio within the contact area. Figure 8a to 8f illustrate the results with a single compliant layer of a and E 1 = 0.5E s , by comparing the results with the ones from Fig. 7a to 7f, the effects of compliant layer thickness on the traction distribution and stick/slip results can be obtained. It can be noted that for the compliant layer larger traction results can be obtained with thinner layer under same creepage conditions, the slip portions are increasing for thicker compliant layer conditions. While the opposite behaviors can be found from the results as given by Fig. 8g to 8l and Fig. 7g to 7l for stiff layer conditions. Traction distribution results with two compliant layers of 0.25a and E 1 = 0.5E 2 = 0.25E s are illustrated in Fig. 9a to Fig. 9f, while the results with two stiff layers of 0.25a and E 1 = 2E 2 = 4E s are presented in Fig. 9g to Fig. 9l. Here, 0.25a is the thickness for each layer and therefore the total thickness is 0.5a for each calculation. For compliant layer conditions, similar to the results of a single layer, for the conditions with both longitudinal and lateral creepages, the slip zones are located at the trailing edge and show   Fig. 10a to f, traction distribution results with two compliant layers of 0.5a and E 1 = 0.5E 2 = 0.25E s are shown, and the ones with two stiff layers of 0.5a and E 1 = 2E 2 = 4E s are given in Fig. 10g to l. It is interesting to note that for compliant layers, the ratios of slip zones are larger when the thicker layers are adopted by comparing the results of Fig. 10a to f with the results of Fig. 9a to f. While for the stiff layers, smaller ratios of slip zones can be found for thicker layers and this behavior can also be seen in the corresponding results of Figs. 8 and 9.
As can be seen from the above results, different layer thickness and elastic modulus parameters are employed in the numerical simulation under different creepage combinations. Opposite behaviors can be found for the conditions with compliant coating layers and hard coating layers. On one hand, for compliant coating layers, larger contact area can be achieved, and at the same time the mean contact pressure results are smaller. Under the same creepage combinations, for thicker compliant layers, the slip zone ratios are increasing. On the other hand, for hard coating layers, the contact area results are smaller and larger contact pressure results can be found. The slip zone ratios are decreasing with decreasing layer thickness. Therefore, it can be concluded that the traction distribution and stick/slip zone results can be modified using different coating layers.

Conclusions
In this work, the rolling contact problems with one or two elastic layers are considered. The influence coefficient relationship between the displacement and the stress is obtained from the corresponding frequency response functions, and we extend this method to solve the rolling contact problems with elastic layers. An efficient numerical model using the conjugate gradient iteration method and the FFT technique has been proposed. All possible creepages, such as the longitudinal creepage, lateral creepage and spin creepage, are considered in the current model to consider about the rolling contact problems with elastic coating layers.
Several numerical simulations have been carried out and as per the results, the following remarks can be derived: • Larger contact area and smaller normal pressure results can be obtained with compliant coating layers, while smaller contact area and larger normal pressure results can be achieved using hard coating layers.
• The traction distributions can be significantly adjusted by different coating layers with the same creepage combinations, and it is found that the slip zones are more predominated for hard coating layers. • For compliant coating layers, the slip zone ratios are increasing with increasing layer thickness. • For hard coating layers, the slip zone ratios are increasing with decreasing layer thickness.
As for future work, we will improve the numerical model to include surface roughness, and to incorporate three or more coating layers.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.