Issue 
Mechanics & Industry
Volume 19, Number 4, 2018



Article Number  405  
Number of page(s)  14  
DOI  https://doi.org/10.1051/meca/2018033  
Published online  01 November 2018 
Regular Article
Semianalytical transient solutions of fully thermomechanical coupled problem for sliding contact
School of Mechanical Engineering, Yanshan University,
066004
Qinhuangdao, PR China
^{*} email: ysupfsun@sina.com
Received:
20
December
2017
Accepted:
31
July
2018
The semianalytical transient solutions of fully thermomechanical coupling for sliding contact have been deduced by using the planepoint heat source equivalent method, the thermoelastic displacement potential and the central difference method. To improve the calculated efficiency, the discrete convolution and fast Fourier transform are applied. The effects of coupled term, sliding velocity, frictional coefficient and initial temperature on temperature rise and the influence of thermal behavior on stress are studied. Results indicate that the coupled term can significantly affect the temperature rise for high velocity and it can be ignored for low velocity; the temperature rise is sensitive to the sliding velocity and the frictional coefficient; the Mises stress on the contact patch is decreased when considering thermal effect.
Key words: Sliding contact / fully thermomechanical coupling / semianalytical transient solutions / temperature rise / stress
© AFM, EDP Sciences 2018
1 Introduction
The dissipative mechanical energy due to the frictional work during the sliding contact process is transformed into frictional heating, which leads to temperature rise in the contact bodies. This thermal phenomenon exists widely in engineering applications and can significantly affect the mechanical properties of contact bodies. The knowledge of the thermoelastic stress field in sliding solid bodies is essential for accurate fatigue life and strength analysis because the total stress is comprised of the thermal stress and elastic stress. The thermal and mechanical analyses should be fully coupled and conducted simultaneously in sliding contact.
An overview of thermoelastic mechanics has been summarized in Johnson's book [1] and Yan's book [2]. The instantaneous point heat source solution and the heat source method proposed by Jaeger [3] are commonly used to solve various complicated heat transfer problems in tribology and manufacturing processes [4,5]. According to the heat source method, Hou [5] developed the comprehensive solutions for stationary/moving plane heat sources of various shaped. Levin [6] suggested that Jaeger's solution for the quasisteadystate problem of moving heat source is inconsistent with the character of source intensity, he gave the improved fundamental solutions of various heat source (such as point, line and circular) by exponential transformation of moving coordinates, and the comparison between Jaeger's solution and fundamental solution was conducted to explain the issue he pointed out. Bosman [7] developed a semianalytical transient temperature model by discretizing arbitrary shaped heat source into uniform square heat source, the method adopted here actually applied the idea of calculus. This solution is appropriate for microcontact and macrocontact. However, it can be only used to solve the surface temperature of bodies and cannot be used to analyze the temperature below the surface. Even so, his work is quite a useful reference in terms of temperature analysis. The separation of variables is also an effective feasible method for exploring the analytical solution of heat conduction equation [8–10]. Kuo et al. [9] investigated the 3D heat conduction problem with a rectangularshaped moving plane heat source by applying the method of separation of variables, the temperature rise expressions derived are similar to the form of solutions obtained through the heat source method. Nevertheless, the solutions of moving heat source presented by Araga et al. [9] and Belghazi et al. [10] are not convenient to use due to the complicated construction. Hamraoui et al. [11] presented an exact analytical solution of temperature in the form of series, but each term of expression is computed impossibly in actual application, as a result, his solution is still approximate. In addition, the solution is valid for the case of uniform heat source, steady state and cylindrical coordinate, which leads to a limited application scope.
Kennedy et al. [12], a brief review of the earlier work in terms of contact temperature, pointed out that although analytical solutions of contact temperature are readily available, most of them do not consider the potentially important factors such as surface roughness, convective cooling or finite body geometry. To select the simple, accurate, and easytouse expressions that can be used to predict contact temperatures, Kennedy et al. [13] considered a variety of models for temperature rise, and the predictions of those models were compared with finite element thermal models of sliding contacts, the evaluations show that the models suggested by Alilat et al. [14] and Laraqi et al. [15] gave a more accurate representation of nominal steadystate temperature rise. Bansal et al. [16] studied the accuracy of approaches presented by Blok [17], Jaeger [3] and Tian and Kennedy [18] for the estimation of maximum and average temperatures at the interface for the case of sliding contacts, it was found that the accuracy of formulas developed by Tian and Kennedy is the highest, and the typical errors of the method presented by Blok less than 3%, whereas Jaeger's hypothesis leads to the accuracy that has typical errors of around 6%. Based on this research, Tian and Kennedy's formulas are the first choice, followed by Blok's method, and Jaeger's method at last for calculating the maximum and average temperature rise. Wu et al. [19] used the formulas developed by Tian and Kennedy when analyzing the characteristics of contact temperature under threebody dry friction conditions.
In addition to the analytical and semianalytical models mentioned above, the numerical method is prevalent method for studying the contact temperatures recently, of which the finite element method is the most widely used, such as the thermal analysis of automotive brake, gear drives, etc. [20–23]. However, the foregoing literatures involving thermal do not consider the effects of thermomechanical coupling. Research manifested that the thermomechanical effect plays a significant role in the variation of contact temperature and pressure in some cases [24].
In view of the mathematical complexity in thermomechanical problems, the finite element method has been widely used to investigate the contacts which involve frictional heating [25–28]. The finite element thermomechanical models developed by Hwang [25] and Balci [26] are a oneway coupling, which only consider the effects of thermal on mechanical, and Balci's model is suitable for 2D thermomechanical analysis. Comparing to the numerical methods, analytical methods have contributed to understand the essence of thermomechanical coupled problems if the analytical solutions of problems could be obtained. Kennedy [29] reviewed the significant developments in the research of thermomechanical phenomena. Liu and Wang [30] investigated the 3D thermomechanical contact of nonconforming rough surfaces and derived the analytical solutions using the integral forms of known Green's functions; however, the solutions presented are adjusted to small Peclet numbers and steadystate heat transfer. They also illustrated and discussed the thermoelastic stress field inside a translating halfspace with constant velocities by using the discrete convolution and fast Fourier transform method when a parabolic type or an irregularly distributed heat source is applied [31]. Chen and Wang [32] developed a threedimensional thermoelastoplastic contact model for counterformal bodies, which considers steady state heat flux, temperaturedependent strain hardening behavior, and interaction of mechanical and thermal loads. Wang and Ai [33] established an extended precise integration model to analyze the thermomechanical behavior of multilayered transversely isotropic materials, and the actual solutions were obtained by adopting numerical quadrature schemes for the LaplaceFourier transform inversion. Altan et al. [34] studied the elasticplastic thermal stresses of the composite disc under known temperature rise by using the Fortran computer program. Zhang and Xiang [35] proposed a semianalytical method for calculating the displacement and stress induced by heat transfer using Goodier's thermoelastic displacement potential and Laplace transform, but the thermal response was not influenced by the mechanical response. KulchytskyZhyhailo et al. [36] studied the thermal stresses for two extreme cases of an axially symmetric contact problem with frictional heat, i.e. the unloaded bounding surfaces of bodies are thermally insulated and the temperature of outside of contact area is zero. It was found that the thermal stresses in both cases differ considerably.
Although the previous works have contained very useful information and insight into the temperature and stress fields in solids due to the stationary/moving heat source or the sliding friction, the comprehensive study of the semianalytical transient solutions of fully coupled thermomechanical problem for the sliding contact are inadequate. Therefore, the main objective of present research was to deduce the systematic semianalytical transient solutions of fully coupled temperature and stress fields for the sliding contact based on the heat source method and the thermalelastic displacement potential. Results of temperature rise because of different variables and the effects of thermal behavior on stress fields are discussed and analyzed. The phenomenon of sliding friction exists widely in the field of machinery, such as gears, bearings and so on. The work done in this paper can provide some theoretical support for the research of thermomechanical coupled problem in the mechanical field.
2 Hypotheses and model description
The following necessary hypotheses are made in the contact model development for obtaining the semianalytical transient solutions of temperature rise and stress:

the material of contact bodies is same and they are perfectly elastic solids, the deformation is small;

the contact is dry friction or boundary lubrication and the sliding friction coefficient is constant. Without considering the effects of the surface curvature microchange and roughness;

the area outside of the contact patches is assumed to be adiabatic, all of the dissipative mechanical energy due to the frictional work is converted into heat, and it is completely absorbed by two contact bodies with no loss of heat to the surroundings. The bodies have an initial temperature of uniform distribution.
In order to illustrate the universality of the problem, the sliding contact of nonconforming surfaces is the object of the fully coupled thermomechanical investigation, as shown in Figure 1a. The dimensions of the interface Ω between contact bodies, which are subject to the normal load F_{n}, are generally small when compared to the characteristic size of the bodies. Hence, each contact body can be treated as a halfspace elastic solid, as shown in Figure 1b. In the figure, the symbol V is the sliding velocity of body 2 relative to body 1. The symbol V_{p} is the moving velocity of contact patch, V_{p} = V for the body 1, and V_{p} = 0 for the body 2. The coordinate system oxyz is fixed with the contact body, the initial contact point is the coordinate origin and the contact surface is the xyplane; the coordinate system o'x'y'z' is connected with the contact patch Ω and the center of patch is the coordinate origin. The positive of xaxis (or x'axis) is regarded as the moving direction of contact patch. The two coordinate systems coincide under the initial state.
The interface Ω has to be subjected to the distributed normal pressure p(x',y'), tangential traction s(x',y') and heat flux q(x',y') due to the friction of the contact bodies. If the two bodies are same, the tangential traction does not disturb the distribution of normal pressure, and contact characteristics are determined by the profiles of the two surfaces and the normal force [1], so the law between them can be expressed as (1a) (1b) (1c) where p_{0} is the centeral pressure of the Hertzian contact ellipse; a and b are the semiaxes of the ellipse in the x' and y' direction respectively; f is the sliding friction coefficient on the interface.
Fig. 1 The sliding contact between nonconforming surfaces. 
3 Decouple of governing equations
Based on the above assumptions and simplified model, ignoring the influence of the volume force and the inertia term, the thermomechanical coupling governing equations can be denoted as (2) where λ and G are Lame's constants. β is thermal stress coefficient, β = Eα/(1−2µ). E, µ and α are the Young's modulus, Poisson ratio and thermal expansion coefficient, respectively. ρ is the density of material. k_{d} is the thermal diffusivity, k_{d} = k/ρc. k and c are the thermal conductivity and specific heat capacity respectively. u_{ζ} (ζ = x, y, z) is the displacement components. ε is the volumetric strain. θ is the temperature rise, which equals to the difference of the temperature at time t and the initial temperature, θ = T_{t}−T_{0}.
The temperature field and the stress field (or strain field) are fully coupled, they interact with each other. Therefore, the governing equations should be decoupled firstly for obtaining the transient solutions of the coupled thermomechanical contact if we want to use the analytical method. The approach that introduces the thermoelastic displacement potential Φ is used to do it, and then the decoupled governing equations can be obtained (3) (4)
4 Semianalytical transient solutions for contact bodies
4.1 Temperature rise
The basis for the heat source method is the solution for the instantaneous point heat source [3], i.e. (5) where Q is the amount of heat liberated by the instantaneous point heat source; R^{2} = (x−x_{0})^{2} + (y−y_{0})^{2} + (z−z_{0})^{2}, (x_{0}, y_{0}, z_{0}) is the coordinate of the point heat source. The following classical Fourier heat conduction equation can be satisfied with the equation (5) (6)
It is obvious that the equations (4) and (6) are identical in form, and the only difference is the coefficients in each equation. Therefore, the temperature rise considering the thermomechanical coupling due to the plane heat source of halfspace elastic solid is deduced based on the solution of the instantaneous point heat source.
For the body 1, it can be equivalent to a halfspace elastic solid which is subjected to an elliptical moving plane heat source on the surface, as shown in Figure 2. In order to analyse the temperature rise of the arbitratry point in contact body at the observation time t, a new time parameter τ (known as the historical time) is required. The elliptic heat source is moved from the initial position to the dottedline position in the graph at the τ moment. The infinitesimal Aarea on the moving elliptic heat source is taken, and the heat acting on the Aarea in the dτ time is expressed as (7) where γ is the heating partition coefficient (the amount of heat absorbed by the body 1 accounts for the proportion of the total frictional heating). It is inconvenient to determine the complete heating partition coefficient by doing a specific numerical analysis for each sliding cantact case, most investigators have used an approximation, such as the one first made by Blok [17]. Blok's criterion determined the heating partition coefficient can be expressed as (8)
The temperature rise of the arbitrary point (x, y, z) in contact body 1 at the observation time t, which is induced by the heat dQ_{A}, can be written as (9)
On the level of theory, the integration of the equation (9) over the entire contact ellipse and the entire time obtains the temperature rise at any point as well as at any observation time in the body 1. The integral operation may be easier if the shape of the heat source is rectangular, but for the general shape of the heat source, even with the use of numerical integration, this integral operation may also be very difficult. The method of direct integration is not adopted in this paper. The semianalytical transient solution of the temperature rise for the contact body is deduced by using an equivalent moving point heat source instead of the moving plane heat source. As shown in Figure 2, an equivalent moving point heat source dQ_{e} is assumed at the center of the contact ellipse, and the associated temperature rise is equal to the result generating by the dQ_{A} at arbitrary point and arbitrary time in the contact body. Then the problem of moving plane heat source can be equivalent to the issue of moving point heat source, and the works may become relatively simple.
Due to the hypothetical equivalent moving point heat source dQ_{e}, the temperature rise at arbitrary point and arbitrary time in the contact body 1 is expressed as (10)
As shown in Figure 3, the threedimensional region near the interface is selected as the computational target zone for the correlation fields. The region is discretized into M × N × L elements. The pressure on the area outside the contact ellipse is zero. It is assumed that the pressure acting on each element on the surface is uniform distribution [37]. The assumed equivalent moving point heat source Q_{e} corresponding to the plane heat source at the τ moment can be expressed as the following semianalytical form by uniting the equations (7), (9)–(10) and integrating over the interface (11)
To make according to the solution of instantaneous point heat source, the temperature rise of the arbitrary point in the contact body 1 at the t moment is expressed as (12)
The equation (12) can be regarded as the operation of discrete convolution with considering the transformation relation between coordinate systems oxyz and o'x'y'z', therefore the discrete convolution and fast Fourier transform technique can be applied to improve the computational efficiency. To facilitate the presentation, the p_{i}_{,j} is used to denote the pressure p(x'_{i}, y'_{j}) of the element (i, j) on the surface, and the is used to denote the influence of the pressure p_{i}_{,j} on the temperature rise of the specified point (x_{m}, y_{n}, z_{l}), i.e. p_{i}_{,j} = p(x_{i}, y_{j}), . The discrete convolution form of the equation (12) is expressed as (13) where m = 1, 2,…, M; n = 1, 2,…, N; l = 1, 2,…, L; The sign “**” denotes the twodimensional discrete convolution operation; “≈” denotes the twodimensional discrete fast Fourier transform; “IFFT” represents the inverse fast Fourier transform. For Matlab, the twodimensional discrete convolution operation is only a directive, which is very convenient to use. Of course, we can also calculate the temperature rise using the fast Fourier transform based on the discrete convolution theorem.
For the body 2, it can be treated as a halfspace elastic solid which is subjected to the elliptic stationary plane heat source on the surface due to V_{p} = 0. Then the temperature rise of the arbitrary point and arbitrary time in the body 2 is expressed as (14)
And the discrete convolution form (15)
Although the above semianalytical transient solutions are developed based on the elliptical heat source and the Hertzian pressure distribution, they are equally effective for other shaped heat sources, (e.g. rectangular, circular or even irregular shapes) and pressure distribution patterns (parabolic, homogeneous, etc.). For various heat sources, the corresponding temperature rise can be obtained by altering the pressure distribution in the abovementioned formulas.
Fig. 2 The elliptic moving heat source on the surface of the body 1. 
Fig. 3 The computational target region of the temperature and stress. 
4.2 Stress
Total stress of contact body can be decomposed into two parts: elastic stress produced by the surface distributed tractions p(x', y') and s(x', y'), and the thermal stress caused by temperature variation θ, (16)
The elastic stress near the contact interface is given by [37–39] (17) where and are the influence coefficients for the normal pressure and the tangential traction, respectively, the detailed expressions are listed in Appendix.
Because of the similar expressions of the temperature rise for contact bodies, the identical approach and principle is employed to analyze the thermal stress field of the contact body 1 or body 2 due to the temperature variation. For the body 1, the alternative expression of the temperature rise generated by the moving plane heat source is given as (18)
Substituting the equation (18) into the (3), then (19) where , ▽^{2} is Laplace operator.
The following transformations are conducted for the integral variables
And then the particular solution of the equation (19) is expressed as (20) here, or, (21)
The thermal stress associated with the special solution of thermoelastic displacement potential can be calculated using the constitutive equation of elasticity, (22) the corresponding secondorder numerical partial derivatives are calculated by using the central differences, for example (23a) (23b)
The rest of secondorder partial derivatives can be calculated in a similar fashion.
The area inside the contact interface can not expand freely under the action of frictional heat and the pressure on the contact patch will increase due to the mutual constraints; while the area outside the interface is a free surface where the distribution pressure is zero. Therefore, the boundary conditions on the surface of the contact body can be expressed as (24)
The stress formed by the superposition of and can only satisfy the first boundary condition of equation (24), the following assumed normal distribution pressure should be applied to the surface of the contact body in order to satisfy the second boundary condition (25)
The stress produced by the assumed pressure in the contact body can also be calculated by the influence coefficient method. The superposition of and is the required thermal stress, i.e. (26)
The total stress considering the thermomechanical coupling for the contact body 1 is obtained by substituting equations (26) and (17) into (16).
For the contact body 2, the relevant mathematical formulas are as follows (27) (28)
5 Results and discussion
Numerical simulation results of transient temperature rise and stress of the contact interface which is simulated by the sliding contact of an ellipsoid and a halfspace solid is presented in this section. The materials of two contact bodies are assumed to be bearing steel (GCr15), and the mechanical and thermal properties of GCr15 [40] and the simulation parameters are summarized in Table 1. The Cartesian coordinate equation of the ellipsoid is , and the coordinates of the contact point with the halfspace body on the surface of the elliposid is (0,0,4), then the corresponding main curvature radiuses at the contact point are R_{1} = 9 mm, R_{2} = 4 mm. The Hertzian solutions [1] of contact elliptical semiaxes and peak pressure are a = 0.1918 mm, b = 0.1279 mm and p_{0} = 1946.8 MPa when the load F_{n} = 100 N.
Parameters and material (GCr15) properties in simulations.
5.1 Verification
In order to verify the current method, the steadystate maximum temperature rise (Fig. 4) and the pressure distribution (Fig. 5) on the contact interface are calculated based on different approaches [41–44], and the results obtained are compared. The comparison results indicate that the solutions obtained by present method have an approximate agreement with those calculated by other investigators' methods.
Fig. 4 Comparison of the steadystate maximum temperature rise. 
Fig. 5 Comparison of surface pressure along the xaxis for body 1. 
5.2 Temperature rise
Owing to the heat allocation scheme from Blok's work is adopted in this paper, the maximum surface temperature rises of the two contacting bodies can only be guaranteed to be equal, as shown in Figure 6, and it is impossible to ensure that the temperatures at all points within the contact interface are equal. It can be seen from the data in the figure that although the maximum surface temperature rises of two contacting bodies are equal, the gradient of surface temperature rise for the body 1 is greater than the temperature rise gradient of the body 2, this indicats that the movement of heat source will aggravate the uneven degree of temperature distribution on the contact surface.
For the sliding contact problems with elliptical contact patch, the dimensionless Péclet number P_{e} is usually defined as P_{e} = Va/2k_{d}, it is an extremely useful parameter for the relative velocity of motion of the heat source considering the thermal properties of the conduction medium which determines the speed of dissipation of heat in the medium. According to Blok's heat allocation scheme, Figure 7a presents that the variation of heating partition coefficient with sliding time under various Péclet numbers, and the variation of corresponding maximum contact temperature rise is shown in Figure 7b. As shown in Figure 7a, it is known that the heating partition coefficient is approximately equal to 0.5 when the Péclet number P_{e} ≤ 0.1; for P_{e} ≥ 0.1, the greater the relative sliding velocity and the longer the sliding time, the higher the proportion of heat distributed to the contact body 1. Figure 7b manifests that the contact interface temperature rise can reach the quasisteady state in a very short time during the sliding process under the current constant velocity and constant force, the current quasisteady state time is about 0.1 s. However, there are few cases of constant force and constant velocity in practical engineering. In most cases, the force and speed are timevarying, such as alternating load in bearing and frequent variations of velocity in drive system, and it will take a long time to reach quasisteady state in these cases. The transient models developed here can solve the thermal problems of these cases and contribute to better understand the changing process of interface temperature. After the quasisteady state time, the change of the contact temperature rise is very small and the increase of corresponding heating partition coefficient has become smaller.
The difference of surface temperature rises between the coupled and uncoupled models for body 1 is depicted in Figure 8, which can reflect the effect of thermomechanical coupled term on the temperature rise. When P_{e} ≤ 1, it means that the relative sliding velocity between contact bodies is small, the difference between coupled and uncoupled models on the temperature rise is small, hence the coupled term can be neglected for this situation, but there are few cases of small sliding velocity in actual sliding contact process, the sliding velocity is larger in most cases of sliding contact; for the case with the P_{e} ≥ 5, the effect of coupled term must be considered, it is obvious that the effect of the thermomechanical coupled term is more significant with the increase of the sliding velocity. According to the trend shown in the figure, for the cases with large sliding velocity (e.g. automotive braking, grinding process), the coupled model can more accurately calculate the temperature rises. Besides, the figure also shows that the higher the temperature rise on the contact patch, the greater the effect of thermomechanical coupled term.
The effects of sliding velocity, friction coefficient, initial temperature on the temperature rise are described in Figure 9. The curves in the figure demonstrate that the variations of sliding velocity, frictional coefficient can significantly affect the temperature rise, and the temperature rise is not sensitive to the variation of the initial temperature. Also, the surface temperature rise of the contact body is rapidly reduced from the peak position to the surrounding area, this phenomenon indicates that the heat affected zone of the contact bodies is concentrated only in a small region near the contact patches.
For the case with a sliding velocity which belongs to the Péclet number P_{e} ≤ 1, as shown in Figure 9a_{1}, the curves reveal that the characteristics of surface temperature rise distributions of the body 1 for the low sliding velocity are still similar to the symmetrical distribution of the contact center, and the peak temperature rises are located on the contact center. However, the maximum temperature rises markedly shift from the center to the rear edge of the contact region when P_{e} ≥ 5, and the trend of temperature with the sliding velocity is enhanced.
Fig. 6 Surface temperature rise distribution. (a): body 1; (b): body 2. 
Fig. 7 Variation of the heating partition coefficient and the maximum temperature rise. 
Fig. 8 Difference of surface temperature rise along the xaxis between coupled and uncoupled models for body 1. 
Fig. 9 The effects of sliding velocity, friction coefficient, initial temperature on temperature rise, (a_{1})∼(a_{3}): body 1; (b_{1})∼(b_{3}): body 2. 
5.3 Stress
The surface pressure distribution of the body 1 calculated by the different models is shown in Figure 10. Figure 10a is obtained by the Hertzian model and Figure 10b is obtained by the model discribed herein. Figure 10 manifests that the boundary conditions on the surface of the contact body can be satisfied by applying the assumed normal distribution pressure listed in equation (24), i.e. the surface outside the contact zone is free, the pressure distribution on it is zero, and the pressure on the contact patch is increased due to the surface inside the contact patch can not expand freely. Figure 11 shows that the surface pressure considering the thermal effect of the body 1 under various sliding velocities. For the case of P_{e} ≤ 1, due to the contact temperature rise caused by friction heat is small, the effect of thermal behavior on the surface pressure is weak; but for P_{e} ≥ 5, the friction heat has a significant effect on the surface contact pressure.
On the surface of the contact body, the effect of the thermal stress produced by the uneven temperature distribution on the each stress component along the xaxis is illustrated in Figure 12. Negative values in the figures mean compressive stress and positive values represent the tensile stress. For the three normal stress components, as shown in Figure 12a–c, the thermal stress increases the compressive stress on the contact patch, and the effect of thermal stress on the normal stresses σ_{1x} and σ_{1y} is greater than that on the normal stress σ_{1z}; as the free expansion of the surface outside the contact patch, the thermal stress only increases the value of the normal stresses σ_{1x} and σ_{1y}, and the normal stress σ_{1z} remains at zero.
Figure 12d–f show the effect of thermal stress on the three shear stress components on the surface of contact body, the shear stress component τ_{1xy} has both tensile shear stress and compressive shear stress on the contact patch, and the thermal stress reduces the tensile shear stress and increases the compressive shear stress, the value of τ_{1xy} is small by comparing with the normal stress; the surface pressure without considering the thermal effect is symmetric about the xaxis, so the elastic shear stress along the xaxis is zero, but the slight compressive shear stress along the xaxis will be generated when the thermal effect is considered; for the shear stress component τ_{1zx}, the thermal stress increases the compressive shear stress on most contact areas and slightly decreases the value of τ_{1zx} on the small part of the contact area, the thermal stress causes the compressive shear stress on the area ahead of the patch and causes the tensile shear stress on the rear area of the patch.
The abovementioned analyses indicate that the thermal stress has different effects on the various stress components, and it can not be ascertained whether the stress state of a specified point in the elastic body is enhanced or weakened by the thermal stress, also the stress state of a specified point can not be determined only by a certain stress component. The Mises equivalent stress is defined as it takes into account the influence of each stress component, and its value is independent of the selection of the coordinate system. Therefore, the best choice to measure the stress state of the elastic body is Mises equivalent stress. The effect of thermal stress on the Mises equivalent stress on the surface of contact body is shown in Figure 13. By comparing three curves in the figure, it can be seen that the thermal stress will decrease the Mises equivalent stress on the contact patch and increase the value of equivalent stress on the surface outside the patch. The reason why Mises equivalent stress shows this feature can be explained by the definition of Mises equivalent stress and Figure 14. The total Mises equivalent stress is not a simple superposition of the elastic Mises equivalent stress and the thermal Mises equivalent stress when considering the thermal effect, it should be determined by the value of the difference between the normal stress components and the value of the shear stress components. As shown in Figure 12, the shear stress is small compared to normal stress and thus the normal stress is considered as a major factor affecting Mises equivalent stress; according to the information reflected in Figure 14, it is known that the difference value between the normal stress components on the contact patch is decreased and it is increased for the surface outside the contact patch when considering the thermal effect, therefore the Mises equivalent stress on the surface of contact body shows the feature shown in Figure 13.
Fig. 10 Distribution of surface pressure for body 1. 
Fig. 11 Surface pressure along the xaxis under various sliding velocities for body 1. 
Fig. 12 The effect of thermal stress on each stress component along the xaxis on the surface of body 1. 
Fig. 13 The effect of thermal stress on the Mises equivalent stress on the surface of body 1 along the xaxis. 
Fig. 14 The difference between the normal stress components on the surface of body 1 along the xaxis. 
6 Conclusions
In this paper, the planepoint heat source equivalent method, the thermoelastic displacement potential and the central difference were used to derive the semianalytical solutions of the transient temperature rise and stress of the fully thermomechanical coupling for the sliding contact. The presented method had be verified by using comparative analysis approach, the effects of thermomechanical coupled term, sliding velocity, frictional coefficient and initial temperature on temperature rise and the effect of thermal stress on stress fields were discussed, the results indicate that
during the sliding contact process, for the case of the Péclet number P_{e} ≤ 1, the effect of thermomechanical coupled term can not be considered when calculating the contact temperature rise, but it will have a significant effect on the temperature rise for the P_{e} ≥ 5;

the contact temperature rise of bodies is greatly affected by the sliding velocity and the friction coefficient, the effect of the initial temperature on the temperature rise is small;

when P_{e} ≤ 1, the effect of thermal behavior on the surface pressure is weak, and increasing sliding velocity leads to a noticeable increase in surface contact pressure for the case of P_{e} ≥ 5;

in general, the thermal stress increases the various stress components on the surface of contacting bodies. However, it has different effects on the Mises equivalent stress on the surface inside and outside the contact area: the Mises equivalent stress outside the contact area is increased and the value on the contact area is decreased.
Nomenclature
a, b: Semiaxes of the contact ellipse (mm)
c: Specific heat capacity (J kg^{−1} K^{−1})
f: Sliding friction coefficient
(i, j): Element number on the surface, i = 1, 2,…, M; j = 1, 2,…, N
k: Thermal conductivity (W mm^{−1} K^{−1})
k_{d}: Thermal diffusivity, k_{d} = k/ρc (mm^{2} K^{−1})
m, n, l: m = 1, 2,…, M; n = 1, 2,…, N; l = 1, 2,…, L
M, N, L: Amount of elements in the xʹ, yʹ and zʹ directions
p_{0}: The centeral pressure of contact ellipse (MPa)
p^{*t}: Assumed pressure to satisfy the boundary condition (MPa)
Q: Amount of heat liberated by the instantaneous point heat source (J)
Q_{e}: Hypothetical equivalent point heat source (J)
R_{1}, R_{2}: Main curvature radiuses (mm)
T_{t}, T_{0}: Temperature at time t and the initial temperature (K)
V: Sliding velocity of body 2 relative to body 1 (mm s^{−1})
V_{p}: Moving velocity of contact patch (mm s^{−1})
x, y, z, xʹ, yʹ, zʹ: Space coordinates
Δxʹ, Δyʹ, Δzʹ: Intervals in the xʹ, yʹ and zʹ directions
α: Thermal expansion coefficient (K^{−1})
β: Thermal stress coefficient, β = Eα/(12µ)
ρ: Density of material (kg mm^{−3})
γ: Heating partition coefficient
: Total stress, elastic stress, thermal stress (MPa)
σ_{x}, σ_{y}, σ_{z}: Normal stress components (MPa)
τ_{xy}, τ_{yz}, τ_{zx}: Shear stress components (MPa)
σ_{vm}: Mises equivalent stress (MPa)
Φ, Φ^{*}: Thermoelastic displacement potential and its a particular solution
P_{e}: Peclet numbers, P_{e} = Va/2k_{d}
IFFT: Inverse fast Fourier transform
: Influence coefficients for normal pressure and tangential traction
**: Twodimensional discrete convolution operation
≈: Twodimensional discrete fast Fourier transform
Acknowledgements
The authors would like to gratefully acknowledge the financial support from the National Natural Science Foundation of China (No.51275440).
Appendix
The influence coefficients and can be expressed as [39] (A.1) (A.2) where x_{+} = x_{m}−x_{i}+Δx/2, x_{−} = x_{m}−x_{i}−Δx/2, y_{+} = y_{n}−y_{j} + Δy/2, y_{−} = y_{n}−y_{j}−Δy/2. (A.3)and, (A.4)where R = [x^{2} + y^{2} + z^{2}]^{1/2}; .
References
 K.L. Johnson, Contact mechanics, Cambridge University Press, London, 1985 [Google Scholar]
 Y. Zongda, W. Hongli, Thermal stress, Higher Education Press, Beijing, 1993 [Google Scholar]
 J.C. Jaeger, Moving sources of heat and the temperature at sliding contacts, J. Proc. R. Soc. New South Wales, 76 (1942) 203–224 [Google Scholar]
 H. Chen, Y. Hu, H. Wang, et al., Calculation of temperature fields of bodies in sliding contact without lubrication, J. Tsinghua Univ. Sci. Technol. 47 (2007) 1962–1964 [Google Scholar]
 Z.B. Hou, R. Komanduri, General solutions for stationary/moving plane heat source problems in manufacturing and tribology, Int. J. Heat Mass Transfer 43 (2000) 1679–1698 [CrossRef] [Google Scholar]
 P. Levin, A general solution of 3D quasisteadystate problem of a moving heat source on a semiinfinite solid, Mech. Res. Commun. 35 (2008) 151–157 [CrossRef] [Google Scholar]
 R. Bosman, M.B. de Rooij, Transient thermal effects and heat partition in sliding contacts, J. Tribol. 132 (2010) 021401–021409 [CrossRef] [Google Scholar]
 W.L. Kuo, J.F. Lin, General temperature rise solution for a moving plane heat source problem in surface grinding, Int. J. Adv. Manuf. Technol. 31 (2006) 268–277 [CrossRef] [Google Scholar]
 G. Araya, G. Gutierrez, Analytical solution for a transient, threedimensional temperature distribution due to a moving laser beam, Int. J. Heat Mass Transfer 49 (2006) 4124–4131 [CrossRef] [Google Scholar]
 H. Belghazi, M. El Ganaoui, J.C. Labbe, Analytical solution of unsteady heat conduction in a twolayered material in imperfect contact subjected to a moving heat source, Int. J. Therm. Sci. 49 (2010) 311–318 [CrossRef] [Google Scholar]
 M. Hamraoui, T. Osman, A. Boucheffa, et al., Analytical modelling of the threedimensional steadystate temperature in a bearing ring, Mech. Ind. 12 (2011) 1–4 [Google Scholar]
 F.E. Kennedy, X. Tian, Modeling sliding contact temperatures, including effects of surface roughness and convection, J. Tribol. 138 (2016) 042101 [CrossRef] [Google Scholar]
 F.E. Kennedy, Y. Lu, I. Baker, Contact temperatures and their influence on wear during pinon disk tribotesting, Tribol. Int. 82 (2015) 534–542 [CrossRef] [Google Scholar]
 N. Alilat, A. Bairi, N. Laraqi, Threedimensional calculation of temperature in a rotating disk subjected to an eccentric circular heat source and surface cooling, Numer. Heat Transfer, 46 (2004) 167–180 [CrossRef] [Google Scholar]
 N. Laraqi, N. Alilat, J.M. Garcia de Maria, A. Bairi, Temperature and division of heat in a pinondisc frictional derive−exact analytical solution, Wear, 266 (2009) 765–770 [CrossRef] [Google Scholar]
 D.G. Bansal, J.L. Streator, On estimation of maximum and average interfacial temperature rise in sliding elliptical contacts, Wear, 278–279 (2012) 18–27 [CrossRef] [Google Scholar]
 H. Blok, Theoretical study of temperature rise at surfaces of actual contact under oiliness lubricating condition, in: Proceeding of the Institute of Mechanical Engineers General Discussion of Lubrication, Institute of Mechanical Engineers, London, 1937, pp. 222–235 [Google Scholar]
 X. Tian, F.E. Kennedy, Maximum and average flash temperatures in sliding contacts, J. Tribol. 116 (1994) 167–174 [CrossRef] [Google Scholar]
 H.W. Wu, Y.Y. Chen, J.H. Horng, Contact temperature under threebody dry friction conditions, Wear, 330–331 (2015) 85–92 [CrossRef] [Google Scholar]
 A. Belhocine, M. Bouchetara, Transient thermal Ansys analysis of dry contactsApplication to automotive braking, Mech. Ind. 13 (2012) 45–57 [CrossRef] [Google Scholar]
 Y. Wang, W. Tang, Y. Chen, et al., Investigation into the meshing friction heat generation and transient thermal characteristics of spiral beral gears, Appl. Therm. Eng. 119 (2017) 245–253 [CrossRef] [Google Scholar]
 H. Xiao, D. Tang, Z. Deng, et al., Thermal analysis and experimental verification of the transmission in a lunar drilling system, Appl. Therm. Eng. 113 (2017) 765–773 [CrossRef] [Google Scholar]
 W. Li, J. Tian, Unsteadystate temperature field and sensitivity analysis of gear transmission, Tribol. Int. 116 (2017) 229–243 [CrossRef] [Google Scholar]
 J. Wen, M.M. Khonsari, Thermomechanical effects on transient temperature in nonconformal contacts experiencing reciprocating sliding motion, Int. J. Heat Mass Transfer, 52 (2009) 4390–4399 [CrossRef] [Google Scholar]
 P. Hwang, X. Wu, Investigation of temperature and thermal stress in ventilated disc brake based on 3D thermomechanical coupling model, J. Mech. Sci. Technol. 24 (2010) 81–84 [Google Scholar]
 M.N. Balci, S. Dag, B. Yildirim, Subsurface stresses in graded coatings subjected to frictional contact with heat generation, J. Therm. Stresses, 40 (2017) 517–534 [CrossRef] [Google Scholar]
 G. Zhou, L. Hua, D. Qian, et al., Effects of axial rolls motions on radialaxial rolling process for largescale alloy steel ring with 3D coupled thermomechanical FEA, Int. J. Mech. Sci. 59 (2012) 1–7 [CrossRef] [Google Scholar]
 A. Zahedi, M.R. Movahhedy, Thermomechanical modeling of high speed spindles, Sci. Iran. B, 19 (2012) 282–293 [CrossRef] [Google Scholar]
 F.E. Kennedy Jr, Thermal and thermomechanical effects in dry sliding, Wear, 100 (1984) 453–476 [CrossRef] [Google Scholar]
 S. Liu, Q. Wang, A threedimensional thermomechanical model of contact between nonconforming rough surfaces, J. Tribol. 123 (2001) 17–26 [CrossRef] [Google Scholar]
 S. Liu, Q. Wang, Transient thermoelastic stress field in a halfspace, J. Tribol. 125 (2003) 33–43 [CrossRef] [Google Scholar]
 W.W. Chen, Q.J. Wang, Thermomechanical analysis of elastoplastic bodies in a sliding spherical contact and the effects of sliding speed, heat partition, and thermal softening, J. Tribol. 130 (2008) 041402–041410 [CrossRef] [Google Scholar]
 L.J. Wang, Z.Y. Ai, Plane strain and threedimensional analysis for thermomechanical behavior of multilayered transversely isotropic materials, Int. J. Mech. Sci. 103 (2015) 199–221 [CrossRef] [Google Scholar]
 G. Altan, M. Topcu, N.B. Bektas, et al., Elasticplastic thermal stress analysis of an aluminum composite disc under parabolic thermal load distribution, J. Mech. Sci. Technol. 22 (2008) 2318–2327 [CrossRef] [Google Scholar]
 Y. Zhang, Y.Y. Xiang, A semianalytical method and its application for calculating the thermal stress and displacement of sparsely fractured rocks with water flow and heat transfer, J. Zhejiang Univ. Sci. A Appl. Phys. Eng. 16 (2015) 922–934 [CrossRef] [Google Scholar]
 R.D. KulchytskyZhyhailo, Z.S. Olesiak, Stress distribution in rotating solids with frictional heat excited over contact region, J. Therm. Stresses, 29 (2006) 957–972 [CrossRef] [Google Scholar]
 J.J. Kalker, Numerical calculation of the elastic field in a halfspace, Commun. Appl. Numer. Methods, 2 (1986) 401–410 [CrossRef] [Google Scholar]
 S. Liu, Q. Wang, G. Liu, A versatile method of discrete convolution and FFT (DCFFT) for contact analyses, Wear, 243 (2000) 101–111 [CrossRef] [Google Scholar]
 S. Liu, Q. Wang, Studying contact stress fields caused by surface tractions with a discrete convolution and fast Fourier transform algorithm, J. Tribol. 124 (2002) 36–45 [CrossRef] [Google Scholar]
 X. Li, H. Ding, Z. Tang, Study of themophysical properties for GCr15 bearing steel in continuous casting, J. Mater. Metall. 9 (2010) 241–244 [Google Scholar]
 K. Knothe, S. Liebelt, Determination of temperatures for sliding contact with application for wheelrail systems, Wear, 189 (1995) 91–99 [CrossRef] [Google Scholar]
 F.E. Kennedy, Contact temperature of a moving solid surface, Springer US, 2013, pp. 544–548 [Google Scholar]
 M.A. Tanvir, Temperature rise due to slip between wheel and railan analytical solution for Hertzian contact, Wear, 61 (1980) 295–308 [CrossRef] [Google Scholar]
 J.R. Barber, Distortion of the semiinfinite solid due to transient surface heating, Int. J. Mech. Sci. 14 (1972) 377–393 [CrossRef] [Google Scholar]
Cite this article as: Z. An, P. Sun, Semianalytical transient solutions of fully thermomechanical coupled problem for sliding contact, Mechanics & Industry 19, 405 (2018)
All Tables
All Figures
Fig. 1 The sliding contact between nonconforming surfaces. 

In the text 
Fig. 2 The elliptic moving heat source on the surface of the body 1. 

In the text 
Fig. 3 The computational target region of the temperature and stress. 

In the text 
Fig. 4 Comparison of the steadystate maximum temperature rise. 

In the text 
Fig. 5 Comparison of surface pressure along the xaxis for body 1. 

In the text 
Fig. 6 Surface temperature rise distribution. (a): body 1; (b): body 2. 

In the text 
Fig. 7 Variation of the heating partition coefficient and the maximum temperature rise. 

In the text 
Fig. 8 Difference of surface temperature rise along the xaxis between coupled and uncoupled models for body 1. 

In the text 
Fig. 9 The effects of sliding velocity, friction coefficient, initial temperature on temperature rise, (a_{1})∼(a_{3}): body 1; (b_{1})∼(b_{3}): body 2. 

In the text 
Fig. 10 Distribution of surface pressure for body 1. 

In the text 
Fig. 11 Surface pressure along the xaxis under various sliding velocities for body 1. 

In the text 
Fig. 12 The effect of thermal stress on each stress component along the xaxis on the surface of body 1. 

In the text 
Fig. 13 The effect of thermal stress on the Mises equivalent stress on the surface of body 1 along the xaxis. 

In the text 
Fig. 14 The difference between the normal stress components on the surface of body 1 along the xaxis. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.