Numerical analysis of three-dimensional thermo-elastic rolling contact under steady-state conditions

In this study, a three-dimensional thermo-elastic model that considers the interaction of mechanical and thermal deformation is developed using a semi-analytic method for steady-state rolling contact. Creepage types in all directions are considered in this model. For verification, the numerical analysis results of shear traction and temperature increase are compared separately with existing numerical results, and the consistency is confirmed. The analysis results include heat flux, temperature increase, contact pressure, and shear traction. Under severe rolling conditions, the thermal effect changes the behavior of the contact interface significantly. Furthermore, the effects of creepage, rolling speed, and conformity under different rolling and creep conditions are investigated.


Introduction
Rolling contacts are widely used in machine components, such as railwheels, rolling bearings, camrollers, and gear systems [1]. In recent years, as awareness of global environmental protection increases worldwide, a reduction in the size and weight of mechanical components has become necessary; this results in severe operating conditions such as high loads and high speeds in rolling contact [2]. In rolling contact interfaces, the basic analysis of relative slip and shear traction is crucial for investigating the fatigue life, traction control, and wear. When rolling contact is severe [1], frictional heat and the resulting temperature increase and thermal deformation become significant issues. Rolling without sliding or rotating is often referred to as "pure rolling". However, even though most mechanical parts operate in a pure rolling state, minimal sliding is included [3]. Johnson [4] argues that the term "pure rolling" is ambiguous because the transmission of tangential forces smaller than the limit friction is not excluded owing to the absence of apparent sliding. Instead of pure rolling, he insisted that the terms "free rolling" when the tangential force is zero, and "tractive rolling" when the tangential force is non-zero are appropriate for describing the rolling motion. Traction rolling is primarily explained via the creep phenomenon. The first analytical solutions of the two-dimensional (2D) rolling contact between identical elastic bodies, including the creep phenomenon, were published by Carter [5]. Nowell and Hills [6,7], Bentall and Johnson [8], and Kalker [9] expanded the existing solutions for 2D steady-state rolling contact between identical elastic cylinders to a solution for two dissimilar elastic cylinders. Johnson [10,11] considered the effects of lateral creep and spin, which could not be explained in Carter's study by extending the 2D rolling contact solution to three-dimensional (3D) cases using an approximate method. Kalker [12][13][14] contributed significantly to the problem of rolling contact by proposing several theories and numerical approaches based on the minimum complementary potential energy principle strategy. CONTACT, a well-known program developed by www Kalker, is widely used in academia and industry as it can solve the problem of transient 3D dissimilar body rolling contact by considering creep and spin in all directions. However, this program requires a significant amount of analysis time. Wang et al. [15] used a semi-analytical method (SAM) to reduce analysis burden. They used the conjugate gradient method (CGM) [16] and the discrete convolution fast Fourier transform (DC-FFT) [17] to solve the defined SAM rolling contact problem; their method incurred a significantly shorter analysis time compared with Kalker's solution. Recently, Xi et al. [1] solved the 2D rolling contact problem using a linear complementarity problem formulation based on the advantage that many numerical algorithms are available and easy to understand and implement. In addition, Xi et al. [18] expanded a 2D model to a 3D model and derived a method to consider creep and spin in all directions; this method was compared and verified with Kalker's solution. However, considering that the analysis results above were obtained on an extremely coarse grid, it is questionable whether solving without using FFT is advantageous.
Owing to improvements in computer performance, the finite element method (FEM) has been applied to the analysis of rolling contact, including effects such as crack propagation [19,20], inelastic behavior [21][22][23][24], rolling contact fatigue [25,26], vibration [27,28], and thermal effects [26,29]. The boundary element method (BEM) was similarly employed to analyze rolling contact problems as a dimension-reduction method [30,31]. However, it is time consuming to resolve rolling contact with creep using the FEM or BEM [15], and most studies did not consider the three possible creepage types-that is, longitudinal, lateral, and spin creepage-simultaneously [18].
The microslip generated in the slip zone of rolling contacts generates frictional heat, which causes a temperature increase and thermal deformation. This thermal effect can result in various problems. Studies regarding the thermal effect of wheel rails indicate that this effect can accelerate thermal fatigue or failure by altering the residual stress and reducing the yield limit of the material [32][33][34][35]. Liao and Lin [36] investigated the change in load distribution and friction torque increase due to the thermal effect under the elasto-hydrodynamic lubrication (EHL) condition of angular contact ball bearings. Similarly, Hao et al. [37] proposed a thermal-fluid-solid coupling model and compared the temperature of rolling bearings with experimental results to validate the coupling model. Recently, Gao et al. [38] proposed a kinematic Hertzian thermal hydrodynamic model for angular contact ball bearings and demonstrated that the sliding behavior caused by insufficient traction on the rolling element significantly increased the bearing temperature. However, in the studies above, the exact thermo-mechanical behaviors within the contact patch were not analyzed. Instead, only the thermomechanical behavior of the system was investigated using a FEM-based solution or a simple EHL equation for shear stress.
In the past two decades, the SAM has been used extensively to solve contact problems associated with frictional heating and has proven to be excellent in terms of efficiency and effectiveness for analysis [39][40][41][42][43][44]. The most representative study pertaining to this mechanical SAM is the analysis of elastoplastic behavior in sliding spherical contact by Chen and Wang [43], where the sliding speed, heat partition, and thermal softening effects were considered. Recently, Zhang and Wang [44] used this thermo-mechanical SAM to investigate the thermo-elastic frictional contact response of a sliding sphere in an elastic-layer-substrate system involving an imperfect interface. Nonetheless, a thermo-mechanical SAM that includes creep has yet to be introduced in rolling contact.
The main objective of this study is to develop a 3D thermo-elastic rolling contact model under a steady state. This model is an improved version of the shear traction calculation method and includes the creep effect introduced in Ref. [15]; therefore, it can account for creepage in all directions. This model was combined with a study pertaining to the thermal effects of contact surfaces in Ref. [43]. The analysis results for shear traction were verified through comparison with the solution of Xi et al. [18], and the analysis results for the temperature increase caused by frictional heating were verified through a comparison with the solution of Tian and Kennedy [45]. Subsequently, the results of heat flux, temperature increase, contact pressure, and shear traction under different rolling and creep conditions were analyzed and compared.

3D rolling contact model
In this study, the variations in the distributions of contact pressure and tangential shear traction were numerically investigated under the effect of thermal deformation caused by friction heat between two spheres during rolling contact at high speeds. The Cartesian coordinate system was used to describe the rolling motion, as shown in Fig. 1. Two identical spheres of radius R were pressed into contact by a normal load W, and a contact region was formed. The z-axis was normal to the contact surface, and the  x y plane was tangential to the contact surface. Linear velocities for the rolling direction   as well as angular velocities  ,  , and  were considered.
The longitudinal, lateral, and spin creepage are denoted by respectively. The assumptions used in the rolling contact analysis are as follows: 1) Two bodies of the same elastic material establish contact with each other. Accordingly, because no dissimilar effect occurs, the combined effect of the normal contact pressure and tangential shear traction is negligible.
2) Surface deformation analyses based on elastic and thermal elastic half-space are valid.
3) The Coulomb friction law is valid. 4) Heat generated between two bodies is distributed and transferred evenly.

5)
Heat propagates only in the slip area, and the stick area is assumed to be adiabatic. 6) The system is in steady state with a constant rolling speed. Figure 2 presents a flowchart of this study. In Section 2.2, each part of the analysis is explained comprehensively.

Contact model
A general dry contact model [4] is derived as follows for completeness: where c A is the real contact area, p is the contact pressure,  is the rigid body approach, and h and The elastic surface deformations at points ( , ) x y by contact pressure   ( , ) p x y can be expressed as follows based on the Boussinesq-Cerruti integrals [4].
where   * ( , ) 1/ G x y E r is the continuous Green's function, * E is the equivalent elastic modulus and   2 2 r x y . Equation (4) can be rewritten in a discrete form for the numerical computations, as follows: where z u p C is the coefficient of influence between the contact pressure and the normal surface deformation. This convolution product can be computed more efficiently using the FFT and inverse FFT.
Herein, the convolution operation is denoted with an asterisk ( * ). The DC-FFT algorithm proposed by Liu et al. [17] was used to efficiently calculate discrete convolutions. The evaluation of the thermal surface deformation t z u is discussed in Section 2.4.

Shear traction
For 3D steady-state rolling, the governing equations are as follows [4]: where  x s and  y s denote the slip velocities between the two surfaces parallel to the x-and y-directions, respectively. x u and y u are the tangential surface deformations, and    /  where  is the friction coefficient. In the stick area, x y s s , and Eqs. (7) and (8) are expressed as follows: Integrating with respect to x yields where ( ) f y and ( ) g y are functions of y and must be resolved through iteration similar to the shear traction. The discrete form of Eq. (12) is as follows: The surface deformations can be calculated using the FFT in terms of the influence coefficients; subsequently, Eq. (13) becomes Based on the assumption that the two objects in this study are of the same material, the dissimilar effect is disregarded, which implies that the shear traction is not a function of the contact pressure, and (14) can be simplified as follows: In this section, the development of the equation for tangential shear traction calculation is the same as the method suggested by Wang et al. [15], except that spin is considered. They employed the CGM to solve Eq. (15) based on the constraints of Eqs. (9) and (10), and the detailed solving process is reported in Ref. [15].

Frictional heat flux and heat partitioning
The concept of heat generation by microsliding and heat transfer for rolling contact is illustrated in Fig. 3. The total heat flux q generated at the contact interface can be calculated by multiplying the tangential traction by the relative sliding speed, as follows: As suggested by Eq. (16), heat is generated in the slip region, whereas the relative sliding speed does not exist in the stick region ; hence, no heat flux occurs. This heat flux is allocated to two bodies in contact and then transferred. The heat partitioning coefficient  is the ratio of the heat transferred to one of the two bodies. Therefore, the equations of heat flux flowing into each body in contact are expressed as shown in Eqs. (17) and (18).
As mentioned above, in this study, it is assumed that frictional heat is equally distributed to the two bodies, resulting in   0.5 .

Temperature increase and thermal deformation
Liu and Wang [42] developed formulations for temperature increase and thermal expansion in halfspace due to irregularly distributed surface heat sources in terms of a frequency response function. They presented this formula for three cases, i.e., transient instantaneous, transient continuous, and steady state, and introduced the following dimensionless values to clarify the following equations:  is the linear thermal expansion coefficient, l is the characteristic length, k is the thermal conductivity, and  is the Poisson ratio. The 2D Fourier transform of the normal thermal deformation and the temperature increase in the half-space for the steady state corresponding to the analysis conditions of this study are as follows: where  

Shear traction
To verify the calculation of shear traction in this study for conditions with different linear and spin creeps, the distribution of shear traction, which was calculated using the current model without considering the thermal effect, was compared with the results of Xi et al. [18] when two spheres of the same size were in rolling contact. The model in this study differs from Xi's model [18]; however, for convenience of comparison, the coordinate system used in Xi's study [18] was used in this study. In Xi's study, the material used was steel, whose E was set to 210 GPa,  to 0. 28 a and 0 b represent the Hertzian contact half widths of rolling and lateral direction, respectively; Xi used a 22 × 22 coarse mesh system within the contact area, whereas in the current model, a 256 × 256 mesh system was used to present more details. The analysis results of the shear traction distribution obtained by the current model for the above five cases are presented in Fig. 4 and consistent with Ref. [18], where the shear traction was normalized by the maximum Hertzian contact pressure 0 p and the x and y coordinates were normalized by 0 a and 0 b , respectively.

Temperature increase
An example of a moving uniform square heat source was solved to verify the calculation of temperature increase in this study. The square heat source was imposed on an area measuring  0 0 2 2 a a . The sliding speed was assumed to be constant in this area. For the steady state, a closed-form solution of the example case is available [46]. Figure 5 shows a comparison of the current solution with the results obtained by Tian and Kennedy [46] when Pe x = 18.8. The comparison, performed along the centerline ( y = 0) in the sliding direction, indicates that the current solution is consistent with that of the existing method.

Thermal deformation effect on rolling contact including creep
To analyze the effect of thermal deformation on rolling contact including creep, it was assumed that creep occurred only in the rolling direction, not in the other direction. In this section, the condition  x = 0.003 is used as a representative example. The analysis results of the total frictional heat, temperature increase, contact pressure, and shear traction are presented for the rolling contact of a half-space. Table 1 lists the Heat partitioning coefficient,  0.5 mechanical and thermal properties of the materials and their settings. The results for the total frictional heat flux, temperature increase, contact pressure, and shear traction are shown in Figs. 6-9, respectively, based on whether thermal deformation was considered. Only the temperature increase was a real-dimensional value, whereas the remaining were expressed as dimensionless values. The contact pressure was normalized by the 0 p , the same parameter used to normalize the shear stress. In Figs. 6-9, the interior of the red solid line represents the initial contact region, whereas the interior of the white solid line represents the initial slip region. The interiors of the dotted line represent the contact and slip regions changed by the thermal deformation. Considering thermal deformation, it was observed that the contact and slip regions reduced. The heat flux inside the contact area was non-uniform. Because no heat flow was present in the stick region, the overall concentration of the heat flux was closer to the trailing edge than the leading edge. In addition, owing to the effect of the contact pressure, the heat flux at the center of the contact was greater than that at the edge of the contact. The uneven distribution of heat flux resulted in an uneven temperature increase and thermal deformation. In areas where the heat flux was high, a greater amount of thermal expansion occurred, resulting in a greater contact pressure in the area. The reduction in the contact region due to thermal deformation was caused by an increase in the local contact pressure. The distribution of shear traction changed with the reduction in the contact region and the change in the contact pressure distribution. Consequently, the slip area was reduced, and the position where the slip commenced shifted toward the trailing edge. Figure 10 shows plots of the total frictional heat flux, temperature increase, contact pressure, and shear traction along the centerline (y = 0) in the rolling direction. When spin was absent, their maximum values appeared along the centerline. The decrease in the slip area due to thermal deformation was the greatest at the centerline. As shown from the heat flux plot in Fig. 10(a), the dimensionless length of the slip in the rolling direction along the centerline changed from 0.78 to 0.5. The distribution of heat flux was concentrated toward the trailing edge via thermal     heat flux in these localized regions, the total amount of heat generated between the two bodies was not proportional to the increase in the maximum heat flux. The total amount of heat was 6.8 W when thermal deformation was not considered, whereas it decreased to 5.6 W when thermal deformation was considered. Even when thermal deformation was considered, no significant difference was observed in the total amount of heat; therefore, the maximum temperature increase did not differ significantly. In contrast to the difference in the maximum value of temperature increase, which was a simple representative value, the temperature increase distribution was concentrated toward the trailing edge owing to the change in the distribution of heat flux due to thermal deformation ( Fig. 10(b)). In the contact pressure plot along the centerline, as shown in Fig. 10(c), a decrease in the contact area and an increase in the local contact pressure were clearly observed. Owing to thermal deformation, the dimensionless x-position of the starting point of the leading edge changed from -1.0 to -0.88; when the dimensionless x was between -0.1 and 1.0, the contact pressure exceeded the initial maximum contact pressure. At this time, the maximum contact pressure was approximately 1.26 times the initial contact pressure, and it appeared near the position where the heat flux in the slip region was maximized, not at the center of the contact. The contact pressure distribution near the maximum pressure exhibited a gentle peak around that point. An interesting observation is that the distribution of contact pressure, reflecting the effect of thermal deformation in a specific dry state, contained pressure peaks, similar to the pressure distribution in the EHL state. As shown in Fig. 10(c), the shear traction distribution shows a change similar to the contact pressure distribution, in which the maximum shear traction increases and the shear traction is concentrated toward the trailing edge. A characteristic change in the shear traction distribution due to thermal deformation is that when thermal deformation is not considered, the curve of shear traction at the boundary between the stick and slip is discontinuous; however, when thermal deformation is considered, the curve becomes relatively smooth and exhibits a continuous shape.

Creepage effect
Section 4.1 shows the general tendency of the thermal deformation effect for a specific creepage ratio, whereas this section presents the effect of thermal deformation in accordance with the creepage ratio. Figure 11 shows the contact pressure and shear traction based on combinations of linear creep at the centerline when spin is absent. In the absence of spin, even if the linear creepage ratios are combined, all values of the total heat flux, temperature increase, contact pressure, and shear traction are maximized in the lateral direction when they are in the centerline and are symmetrical with respect to this line. Herein, the solid line indicates the results based on  x when  y = 0.001, whereas the dotted line indicates the results based on  y when  x = 0.001. As a result of the increase in the thermal deformation effect with the increase in the creepage ratio, the contact region decreased, whereas the maximum contact pressure and maximum shear traction increased. When the creepage ratio was relatively low, the contact pressure changed significantly only at the edge of the contact, particularly near the trailing edge, and the amount of change was not significant. By contrast, when one of the linear creepage ratios exceeded 0.003, a contact pressure greater than the maximum Hertzian contact pressure occurred, and the contact pressure indicated a significant change | https://mc03.manuscriptcentral.com/friction over the total contact area, resulting in a completely different contact pressure distribution than the Hertzian contact pressure distribution. The shear traction curve along the centerline shows a distinct difference between the stick and slip regions. In the slip region, the shear traction curve is parabolic, which is the same as the shape of the contact pressure curve. In the stick region, the shear traction curve is a polynomial of the third or higher order. The value of shear traction was extremely low near the leading edge and increased gradually as it shifted away from the leading edge; meanwhile, it increased rapidly as it approached the maximum value. As the creepage ratio increased, the maximum value of shear traction increased, and the change in the slope of the shear traction curve in the stick region decreased. In both the contact pressure and shear traction curves, the change due to creep in the rolling direction was slightly greater than the effect of creep in the lateral direction. Unlike the case in which only linear creeps were combined, when spin was added, a representative line such as a centerline could not be defined; therefore, a contour map must be used to verify the effect of spin. Figure 12 presents a contour map of the analysis results involving the

Velocity and conformity effect
Because the heat flux is the product of shear traction and sliding speed, even if the slip and friction are high, the heat flux may not be sufficient to induce a thermal deformation effect under slow sliding speeds. Because the sliding speed and rolling speed are proportional, the rolling speed is an important factor that significant affects thermal deformation. In addition to the rolling speed, the shape of the two contact bodies significantly affects thermal deformation. In most mechanical components, the radius of the rolling direction is often uncontrollable because of certain limitations in terms of size. Meanwhile, in the lateral direction, a whole circle is not necessitated, and the shape of the circle is only required at the contacting part; therefore, the radius can be used as a design factor. Hence, we analyzed the results by performing an analysis based on the rolling speed for three conditions: non-conformal, flat, and conformal, with different lateral radii. For each condition, the upper body was a sphere with    , and 165 mm for the non-conformal, flat, and conformal conditions, respectively. When the contact between two bodies becomes extremely conformal, the influence coefficients for the non-conformal condition is no longer applicable. Nevertheless, the analysis results for the conformal condition were derived using the same influence coefficients as those for the non-conformal condition, as the contact in the rolling direction remained non-conformal. If extreme conformal contact does not form in both directions, then this approach will be effective; in fact, it is typically used in contact analyses, such as that for rolling bearings. Figure 13 shows the total frictional heat F Q , maximum temperature increase, stick area and contact area ratios, and ratio of maximum pressure increase for the three  y R conditions based on the rolling speed. In this case, unlike the assumption of a constant creep, as in Sections 4.1 and 4.2, it is assumed that the frictional force in the rolling direction x F is constant at 25 N and the friction force in the lateral direction y F does not exist. The total frictional heat increases with the rolling speed but exhibits a curve with a decreasing slope, and the order of samples with the highest total frictional heat at the same rolling speed, from low to high, is conformal, flat, and nonconformal; this implies that as  y R increases, the total frictional heat increases. The maximum temperature increase was almost linearly proportional to the rolling speed and decreased as  y R increased, as opposed to the total frictional heat sequence. This is because when  y R was large, heat flux was generated in a large area, whereas when  y R was small, heat flux was generated in a narrow area, such that when  y R is small, the concentration of heat flux in the local area will be larger. As shown in Figs. 13(a) and 13(b), a small  y R is advantageous in terms of the total frictional heat, whereas a large  y R is advantageous in terms of the local temperature increase. In Fig. 13(c), A to the real contact area c A after thermal deformation, and the stick area ratio  st is the ratio st A to c A . As the rolling speed increased,  c decreased, and  st increased owing to the greater thermal deformation effect. A larger  y R resulted in a higher decrease rate in  c and a higher increase rate in  st based on the rolling speed; this is because a larger  y R resulted in a more even thermal deformation in the contact area. Finally, we analyzed the maximum contact pressure ratio,   max 0 / p p p based on the rolling speed. The rolling speed and  p indicated an almost proportional relationship, i.e., the larger the  y R , the higher was the contact area reduction ratio; therefore, the slope was smaller. A larger  y R under the same load resulted in a lower maximum contact pressure but a greater thermal deformation effect of the maximum contact pressure. Therefore,  y R must be analyzed sufficiently in terms of design at a high rolling speed.

Conclusions
In this study, a numerical analysis method that considers the thermal effect in steady-state rolling contact including creep was proposed based on an SAM. The method was validated and demonstrated consistency with results of previous solutions. Simulation results for heat flux, temperature increase, contact pressure, and shear traction were obtained and compared under different rolling and creep conditions. In rolling contact including creep, the contact and slip areas reduced owing to the thermal effect. Furthermore, the thermal effect induced an increase in contact pressure and shear traction in localized areas, and the position where slip commenced shifted toward the trailing edge. As creep increased, the thermal effect became more prominent, resulting in a significant change where the distributions of heat flux, temperature increase, contact pressure, and shear traction were locally concentrated. Quantitative investigations revealed that rolling speed significantly affected the thermal effect. When two objects were in a more conformal contact, the total frictional heat increased, whereas the contact increased, resulting in a lower local pressure and shear traction, and an increase in the maximum temperature owing to thermal effects. The proposed method can accurately analyze the behavior of an interface in rolling contact under severe conditions. In further studies, we will improve this model by considering more accurate www.Springer.com/journal/40544 | Friction heat partitioning effects in dry conditions, thermal deformation effects in lubrication conditions through a combination with thermal EHL analysis in the lubrication condition, and transient effects. Finally, we will apply our findings to rolling bearings, railwheels, etc.