Three-Dimensional Thermoelastic Contact Model of Coated Solids with Frictional Heat Partition Considered

: In this paper, a three-dimensional thermoelastic contact model of coated solids with the frictional heat partition considered is developed by introducing a frictional heat partition model. The inﬂuence coefﬁcients of the temperature rise, normal displacement and stress components in the three-dimensional thermoelastic contact model are converted from their corresponding frequency response functions (FRFs) with a conversion method based on the fast Fourier transform (FFT), and the FRFs of solids coated with a homogeneous coating subjected to a coupled action of the mechanical loading and the frictional heat ﬂux on its surface are deduced in the frequency domain by introducing a two-dimensional Fourier integral transform. The contact pressure and the frictional heat partition between the two bodies are solved by employing a fast numerical algorithm based on the conjugate gradient method (CGM) and a discrete convolution fast Fourier transformation (DC-FFT). Comparison between the solutions of the present model and those of a thermoelastic contact model in literature is conducted in order to validate the present model. Several speciﬁc conclusions on the effect of the sliding speed, thermoelastic properties and thickness of the coating are drawn based on the result of numerical investigation by utilizing the present model.


Introduction
Solid coatings, such as TiN, diamond-like-carbon (DLC) and MoS 2 , are widely employed to improve the tribological performance and service life of tribo-parts [1][2][3].More and more tribo-parts, such as gears, bearings, are coated for anti-scuffing, anti-friction and anti-wear, especially those working under severe performance conditions.A three-dimensional thermoelastic contact model of coated solids is essential for a more accurate analysis of the contact behavior by considering the partition of frictional heat between the bodies in contact, since the thermal effect of the frictional heat can exert a significant impact on the contact behavior [4].
In the theoretical analysis of the mechanics of solids, linear elastic theory uses various integral transform techniques to produce solutions [5].The general theory of the stress and displacement of layered systems was structured by Burmister for analyzing the contact under prescribed axially symmetric surface normal loading [6,7].Chen [8,9] extended the applications to both axisymmetric and non-axisymmetric normal surface loading.O'Sullivian and King [10] studied the contact of layered materials using Papkovich-Neuber potentials.Nogi and Kato [11] improved the computing speed by employing the fast Fourier transform (FFT) technique based on the work of O'Sullivian and King.The Papkovich-Neuber potentials were further employed to study the sliding contact or partial slip contact problem for solids coated with a monolayer or multilayers [12][13][14].Wang et al. [15] also adopted Papkovich-Neuber potentials to conduct a systematic investigation of the effect of the thickness and stiffness of a homogeneous coating on the contact behavior under various friction coefficients.It was also reported that the equivalent inclusion method, which was often employed to analyze the contact concerning inhomogeneities, also can be used to analyze the contact of coated solids [16].Kral and Komvopoulos [17,18] employed the finite element method (FEM) to study the contact of elastic-plastic layered mediums subjected to repeated indentations by a rigid sphere.Kang et al. [19] also adopted FEM to analyze the failure mechanism of plasma sprayed AT40 coatings under different rolling-sliding contact conditions.Besides the research work mentioned above, numerous investigations have also been conducted on the problem of solids with a homogeneous layer subjected to the mechanical loading coupling with the thermal heat flux, as well as that of uncoated solids [20][21][22][23][24]. Leroy et al. [25,26] applied the Fourier integral transform to a finite thickness layered medium subjected to a moving line heat source and obtained a system of equations which links the transformed quantities, and the solution in the space domain by applying the inverse FFT.Ju and Farris [27] deduced an FFT thermoelastic solution for a half plane subjected to a moving heat flux.Shodja and Ghahremaninejad [28], Ke and Wang [29] studied the problem of a half plane coated with an functionally gradient material (FGM) layer subjected to the coupled action of the mechanical loading and the thermal heat flux.Shi et al. [30] investigated the thermal elastic field of a half space with multilayers using the Fourier integral transform.However, in those studies, the heat flux was given by hypothesis rather than a partition model, and the effect of the thermal displacement on the contact pressure was not considered.The FEM was also used to analyze the contact problem of coated solids subjected to the coupled action of the mechanical loading and the thermal heat flux.Kulkarni et al. [31] established a two-dimensional FEM model for the thermoelastoplastic contact under repeated translation loading.Komvopoulos et al. [32,33] developed a three-dimensional FEM model of elasto-plastic layered mediums under the coupled action of the mechanical loading and the thermal heat flux.In the last several decades, meaningful studies have been done on the contact problem of coated solids and efforts have been made to develop a design and optimization technology of coating for improving the performance of anti-wear, anti-friction and fatigue strength etc. [34][35][36][37].However, little attention has been paid to the investigation of the contact problem of coated solids considering the frictional heat partition, which should be one of the essential factors in the contact problem.
In this paper, a three-dimensional thermoelastic contact model of coated solids considering the partition of frictional heat is developed by introducing a frictional heat partition model.The influence coefficients (ICs) of the temperature rise, surface normal displacement and stress components are converted from their corresponding frequency response functions (FRFs) by employing a numerical conversion method based on the FFT, and the FRFs of the temperature rise, displacement and stress components of a solid with a homogenous coating subjected to the mechanical loading coupling with the frictional heat flux are derived and given in explicit formulas.A fast iteration algorithm based on the conjugate gradient method (CGM) and a discrete convolution fast Fourier transformation (DC-FFT) are adopted to obtain the contact pressure, heat partition and stress in the subsurface.The present model is validated by a comparison between the solution of the present model and that of a model in literature.Lastly, a numerical investigation of the effect of the sliding speed, thermoelastic properties and thickness of the coating on the contact behavior is further conducted.

Contact Model
According to the authors of [38], the thermoelastic contact model of an elastic ball brought into sliding contact with a half space with a homogeneous coating subjected to an external normal load W, as shown in Figure 1, can be described as follows: where g is the gap between the surfaces of the two bodies in contact after the external normal load W is applied, g 0 is the initial gap before the external normal load W is applied, p is the contact pressure between the two bodies, u e is the surface normal displacement due to the contact pressure, u t is the surface normal displacement due to the frictional heat flux, δ is the normal approach and Ω c represents the real contact area.

Contact Model
According to the authors of [38], the thermoelastic contact model of an elastic ball brought into sliding contact with a half space with a homogeneous coating subjected to an external normal load W, as shown in Figure 1, can be described as follows: where g is the gap between the surfaces of the two bodies in contact after the external normal load W is applied, g0 is the initial gap before the external normal load W is applied, p is the contact pressure between the two bodies, ue is the surface normal displacement due to the contact pressure, ut is the surface normal displacement due to the frictional heat flux, δ is the normal approach and Ωc represents the real contact area.In order to employ numerical methods to solve the contact model, a square calculation domain is selected and uniformly divided into Nx × Ny surface square elements centered on the grid nodes.Here aH is the contact radius of the substrate material in Hertz point contact, Nx and Ny are the numbers of elements in x and y direction respectively.The contact pressure distribution is approximated by a piecewise constant function that is uniform within each surface element.Then the Equation ( 1) can be converted into a discretized form: where Δx and Δy are the discretization sizes in x and y direction respectively.According to the linear superposition principle, the displacements ue and ut are obtained by: In order to employ numerical methods to solve the contact model, a square calculation domain and uniformly divided into N x × N y surface square elements centered on the grid nodes.Here a H is the contact radius of the substrate material in Hertz point contact, N x and N y are the numbers of elements in x and y direction respectively.The contact pressure distribution is approximated by a piecewise constant function that is uniform within each surface element.Then the Equation (1) can be converted into a discretized form: where ∆x and ∆y are the discretization sizes in x and y direction respectively.According to the linear superposition principle, the displacements u e and u t are obtained by: Coatings 2018, 8, 470 where KN k [i − r, j − s] and KQ k [i − r, j − s] are the ICs of the surface normal displacement due to the contact pressure and the surface heat flux respectively, k = 1 and k = 2 represent the elastic ball and the coated solid, respectively.The overall flow chart of the numerical algorithm is shown in Figure 2a.One solver named NCS in Figure 2a is a normal contact solver for determining the normal contact pressure with a known u t in Equation ( 2), the other solver named HPS is a heat partition solver for determining the heat flux of Q 1 and Q 2 within the calculation domain based on a heat partition model seen below in the current section.Several different iteration schemes have been proposed for the solver NCS, and the iteration scheme based on the CGM [39,40] is employed here for its comparatively high rate of convergence and explicit iteration format, which makes it compatible with the fast multi-summation method, namely, the DC-FFT [41].
where KNk[i − r, j − s] and KQk[i − r, j − s] are the ICs of the surface normal displacement due to the contact pressure and the surface heat flux respectively, k = 1 and k = 2 represent the elastic ball and the coated solid, respectively.The overall flow chart of the numerical algorithm is shown in Figure 2a.One solver named NCS in Figure 2a is a normal contact solver for determining the normal contact pressure with a known ut in Equation ( 2), the other solver named HPS is a heat partition solver for determining the heat flux of Q1 and Q2 within the calculation domain based on a heat partition model seen below in the current section.Several different iteration schemes have been proposed for the solver NCS, and the iteration scheme based on the CGM [39,40] is employed here for its comparatively high rate of convergence and explicit iteration format, which makes it compatible with the fast multi-summation method, namely, the DC-FFT [41].After the contact pressure and the heat flux applied on the surface of coated solids have been obtained, the stress components and the temperature rise in the subsurface of coated solids can be obtained by:  After the contact pressure and the heat flux applied on the surface of coated solids have been obtained, the stress components and the temperature rise in the subsurface of coated solids can be obtained by: where SN mn [i − r, j − s, z], ST mn [i − r, j − s, z], SQ mn [i − r, j − s, z] and TR 2 [i − r, j − s, z] are the stress components and temperature rise influence coefficients of point at z depth lying directly below the grid node [i, j] for the normal force p, tangential traction q, and heat flux Q 2 applied on the surface element [r, s] of coated solids, respectively.According to the Coulomb's friction law, the tangential traction q[r, s] is given as follows: where µ f is the frictional coefficient.
In Equations ( 3)-( 6), the surface heat fluxes Q 1 and Q 2 should be determined with a consideration of the frictional heat partition between the two bodies in contact.For simplicity, the previous work [42,43] assumed that the frictional heat flux is solely absorbed by one body or evenly partitioned between the two bodies in contact regardless of the real surface temperature distribution.Blok [44] proposed a frictional heat partition model by matching the maximum flash temperature of the two bodies, which has been adopted by Tian and Kennedy [45].Jaeger [46] proposed a heat partition model by matching the average temperature within the contact zone instead of matching the maximum temperature, which was adopted by Archard and Rowntree [47].Francis [48] also proposed a heat partition model by holding the interfacial temperature as the harmonic mean of the two surface temperatures when each body receives all of the frictional heat flux.Gao et al. [49] determined the heat partition between the two bodies by matching the temperature in the whole simulation domain.In the present model, the heat partition model shown as Equation ( 8), which is established by matching the temperature within the real contact area and has been employed by Bos and Moes [50], Chen and Wang [4], is adopted to determine the frictional heat partition, in which the heat convection and thermal radiation are ignored: where T 1 and T 2 are the bulk temperatures of the elastic ball and the coated solid respectively, T 1 and T 2 are the surface temperature rises of the two bodies due to the frictional heat flux, Q 1 and Q 2 are the surface heat fluxes of the two bodies, and Q is the frictional heat flux.The frictional heat flux is as follows: where V s is the sliding speed.Here the heat partition coefficient H pc is defined to describe the ratio between the heat flux flowing into the coated solid and the total frictional heat flux, shown as Equation (10).
The surface temperature rise T k [i, j] can be obtained by: where By substituting Equation ( 11) into the Equation ( 8), the heat partition model can be rewritten as follows: The bulk temperatures T 1 and T 2 in Equation ( 12) can be input according to practical conditions and it is assumed that T 1 equals to T 2 in this paper for simplicity.The unknown Q 2 [i, j] in Equation ( 12) can be solved by a modified iterative algorithm based on a standard CGM [51].The detail of the modified iterative algorithm can be found in Figure 2b, its difference from the iterative algorithm of a standard CGM can be summarized as follows:

•
In the modified iterative algorithm, the calculations of the iteration parameter δ n and the update length of heat flux α in the conjugate gradient direction are conducted over the real contact area Ω c rather than the whole calculation domain Ω as follows: where r[i, j] and b[i, j] are intermediate variables used by the flow chart of heat partition solver as shown in Figure 2b.In the standard CGM, the calculations of δ n and α are conducted over the whole calculation domain Ω.

•
The update of the heat flux Q 2 [i, j] is performed only within the real contact area Ω c in the conjugate gradient direction with the update length α.Outside the real contact area Ω c , the heat flux Q 2 [i, j] is enforced to be zero.
For the elastic ball, the ICs of its displacement and stress components due to the mechanical loading are obtained by applying Boussinesq and Cerruti solutions [52].However, the explicit formulas for all ICs of the coated solid are unavailable as well as the ICs of the elastic ball due to the surface heat flux.In this paper, these ICs are converted from their corresponding FRFs by a numerical conversion method based on FFT [53,54].The detailed steps of the numerical conversion method adopted in this paper can be found in [15].

Frequency Response Functions
For solids with a homogeneous coating subjected to the coupled action of the mechanical loading and the frictional heat flux, including the normal pressure, tangential traction and frictional heat flux moving on its surface, as shown in Figure 3, the governing differential equations of the temperature rise and thermoelastic field can be given as Equations ( 15)-( 18) for quasi-static states: Coatings 2018, 8, 470 where T is the temperature rise, u is the displacement, V is the moving speed of the mechanical loading and the heat flux, and ν, γ and α represent the Poisson ratio, thermal diffusivity and thermal expansion coefficient respectively.The thermal diffusivity is the thermal conductivity κ divided by the volumetric heat capacity c, and the volumetric heat capacity c is a product of the density and the specific heat capacity.where T is the temperature rise, u is the displacement, V is the moving speed of the mechanical loading and the heat flux, and ν, γ and α represent the Poisson ratio, thermal diffusivity and thermal expansion coefficient respectively.The thermal diffusivity is the thermal conductivity κ divided by the volumetric heat capacity c, and the volumetric heat capacity c is a product of the density and the specific heat capacity.
then Equations ( 15)-( 18) can be rewritten as follows: where and   k F can be deduced as follows: e e e e e e e e e e e e where Analogous to the method employed in [55], a two-dimensional Fourier integral transform Φe −i(ω x x+ω y y) dxdy and two intermediate functions y /∂y are introduced, then Equations ( 15)-( 18) can be rewritten as follows: (1 where The general solutions of T (k) , u z (k) , G (k) and F (k) can be deduced as follows: Coatings 2018, 8, 470 8 of 20 where The other two displacement components u x (k) , u y (k) and the six stress components σ xx (k) , σ yy (k) , k) and F (k) , which can be found in Appendix A.
The elastic ball is subjected to a stationary normal force, tangential traction and heat flux, its FRFs of temperature rise T (k) and surface normal displacement u z (k) in frequency domain can be found in Appendix C.

Verification of the Present Model
In order to validate the present model, the solutions obtained with it for the thermoelastic contact between a homogeneous solid of substrate material and an elastic ball are compared with those obtained with Chen's model in [4].The input parameters for the present model are shown in Table 1.
Table 1.Input parameters of present model for simulating the thermoelastic contact between a homogeneous solid and an elastic ball.

Name and Unit of the Input Parameter Value
Elasticity modulus of the ball, E (GPa) 2.1 × 10 2 Poisson ratio of the ball, ν 3 × 10 −1 Thermal conductivity of the ball, κ (W/m•K) 5.02 × 10 1 Volumetric heat capacity of the ball, c (J/m 3 The coating has the same thermoelastic properties as those of the substrate, therefore The symbol a H is the contact radius of the substrate material in Hertz point contact with the ball, which has been mentioned in Section 2.1.In this simulation, the coating thickness equals to the contact radius of the substrate material in Hertz point contact with the ball.
Coatings 2018, 8, 470 9 of 20 The thermoelastic properties of the coating are the same to those of the substrate when the present model is adopted to simulate the contact between a homogeneous solid of substrate material and an elastic ball.Figure 4 compares the contact pressure and the surface temperature rise along x axis under various sliding speeds, where p H and a H represent the maximum contact pressure and the contact radius of the substrate material in Hertz point contact.As indicated in Figure 4, the contact pressure and the surface temperature rise obtained with the present model are in an excellent match with those obtained with Chen's model for various sliding speeds.The comparisons of the temperature rise T 2 and the von Mises stress σ V are shown in Figures 5 and 6, respectively.It can be seen that the temperature rise and the von Mises stress obtained with the present model are also in an excellent agreement with those obtained with Chen's model for different sliding speeds.
Coatings 2018, 8, x FOR PEER REVIEW 9 of 20 The thermoelastic properties of the coating are the same to those of the substrate when the present model is adopted to simulate the contact between a homogeneous solid of substrate material and an elastic ball.Figure 4 compares the contact pressure and the surface temperature rise along x axis under various sliding speeds, where pH and aH represent the maximum contact pressure and the contact radius of the substrate material in Hertz point contact.As indicated in Figure 4, the contact pressure and the surface temperature rise obtained with the present model are in an excellent match with those obtained with Chen's model for various sliding speeds.The comparisons of the temperature rise T2 and the von Mises stress σV are shown in Figures 5 and 6, respectively.It can be seen that the temperature rise and the von Mises stress obtained with the present model are also in an excellent agreement with those obtained with Chen's model for different sliding speeds.The thermoelastic properties of the coating are the same to those of the substrate when the present model is adopted to simulate the contact between a homogeneous solid of substrate material and an elastic ball.Figure 4 compares the contact pressure and the surface temperature rise along x axis under various sliding speeds, where pH and aH represent the maximum contact pressure and the contact radius of the substrate material in Hertz point contact.As indicated in Figure 4, the contact pressure and the surface temperature rise obtained with the present model are in an excellent match with those obtained with Chen's model for various sliding speeds.The comparisons of the temperature rise T2 and the von Mises stress σV are shown in Figures 5 and 6, respectively.It can be seen that the temperature rise and the von Mises stress obtained with the present model are also in an excellent agreement with those obtained with Chen's model for different sliding speeds.

Numerical Results and Discussion
In this section, the effect of the sliding speed, thermoelastic properties and thickness of the coating on the contact behavior is numerically investigated with the present model.The thermoelastic properties of the coated solid and the elastic ball and the other concerning parameters are listed in Table 2.
Table 2. Thermoelastic properties of the coated solid and the elastic ball and the other concerning parameters.

Name and Unit of the Input Parameter Value
Elasticity modulus of the ball, E (GPa) 2.1 × 10 2 Poisson ratio of the ball, ν 3 × 10 −1 Thermal conductivity of the ball, κ (W/mK) 5.02 × 10 1 Volumetric heat capacity of the ball, c (J/m 3

Effect of the Sliding Velocity
The sliding velocity V s has a significant effect on the contact behavior of coated solids as shown in Figure 7.When V s increases from 0.5 to 6 m/s, the maximum contact pressure increases to 1.8p H from 1.2p H and the maximum temperature rise also increases to 546 K from about 10 K.The main reason for the increase of the maximum temperature rise is that the increase of the sliding speed produces more frictional heat on the contact interface, as shown in Equation ( 6), which also leads to an increase of the surface thermal displacement as well as the maximum contact pressure.The heat partition coefficient H pc also increases with the increase of the sliding speed, since more cool surfaces of the coated solid pass through the contact zone in a given time and it absorbs and takes away more frictional heat flux when passing through the contact zone.The transverse stress σ xx at the tail side of the contact zone changes to compressive stress from tensile stress with the increase of the sliding speed as shown in Figure 7d, since the thermal stress σ xx caused by the frictional heat flux is compressive and increases with the increase of the sliding speed.The transverse stress σ xx at the base of the coating is tensile around the contact center and its maximum basically increases with the increase of V s .The main reason should be the concentration of the contact pressure distribution as shown in Figure 7a.The maximum shear stress σ zx on the interface also increases with the increase of V s .In addition to the concentration of the contact pressure distribution, another reason is that the thermal shear stress caused by the frictional heat flux has the same direction as that caused by the surface mechanical loading at the position where the maximum σ zx locates.

Effect of the Coating Thermal Conductivity
The effect of the coating thermal conductivity κ1 on the contact behavior is shown in Figure 8.

Effect of the Coating Thermal Conductivity
The effect of the coating thermal conductivity κ 1 on the contact behavior is shown in Figure 8.
When the coating thermal conductivity κ 1 increases from 0.5κ 2 to 2κ 2 , the maximum contact pressure decreases to 1.5p H from 2.2p H and the maximum temperature rise of the coated solid decreases to 260 K from 716 K.With the increase of κ 1 , the frictional heat flux transfers more quickly into the coated solid from the contact interface, therefore the maximum temperature rise of the coated solid decreases.That is also the reason why H pc increases with the increase of the coating thermal conductivity, as shown in Figure 8c.The surface thermal displacement decreases as a consequence of the decrease of the temperature rise caused by the increase of the coating thermal conductivity.That is the reason why the maximum contact pressure decreases with the increase of the coating thermal conductivity.The transverse stress σ xx on the surface is compressive and its maximum decreases with the increase of κ 1 , the transverse stress σ xx at the base of the coating around the contact center is tensile and its maximum decreases with the increase of κ 1 .The maximum of the shear stress σ zx also decreases with the increase of κ 1 .One reason is the deconcentration of the contact pressure distribution caused by the increase of κ 1 , the other reason may be that the temperature gradient from the surface to the interface becomes relatively small with the increase of κ 1 .
rise along x axis, (c) frictional heat partition coefficient, (d) σxx on the surface, (e) σxx at the base of the coating, (f) σzx on the interface.

Effect of the Coating Thermal Conductivity
The effect of the coating thermal conductivity κ1 on the contact behavior is shown in Figure 8.When the coating thermal conductivity κ₁ increases from 0.5κ2 to 2κ2, the maximum contact pressure decreases to 1.5pH from 2.2pH and the maximum temperature rise of the coated solid decreases to 260 K from 716 K.With the increase of κ1, the frictional heat flux transfers more quickly into the coated solid from the contact interface, therefore the maximum temperature rise of the coated solid decreases.That is also the reason why Hpc increases with the increase of the coating thermal conductivity, as shown in Figure 8c.The surface thermal displacement decreases as a consequence of the decrease of the temperature rise caused by the increase of the coating thermal conductivity.That is the reason why the maximum contact pressure decreases with the increase of the coating thermal conductivity.The transverse stress σxx on the surface is compressive and its maximum decreases with the increase of κ1, the transverse stress σxx at the base of the coating around the contact center is tensile and its maximum decreases with the increase of κ1.The maximum of the shear stress σzx also

Effect of the Coating Volumetric Heat Capacity
The effect of the coating volumetric heat capacity c 1 on the contact behavior is shown in Figure 9.When the coating volumetric heat capacity c 1 increases from 0.5c 2 to 2c 2 , the maximum contact pressure decreases from 1.9p H to 1.59p H and the maximum temperature rise of coated solid decreases to 375 from 519 K.The frictional heat partition coefficient increases with the increase of the coating volumetric heat capacity since a unit volume material of the coating absorbs more heat for 1 K temperature rise with the increase of the coating volumetric heat capacity.That is also the explanation for the influence of the coating volumetric heat capacity on the contact pressure and the temperature rise.Around the contact center, the maximum of the compressive stress σ xx on the surface decreases with the increase of c 1 , while the maximum tensile stress σ xx at the base of the coating increases.The maximum of the shear stress σ zx on the interface also decreases with the increase of c 1 , as shown in Figure 9f.decreases with the increase of κ1.One reason is the deconcentration of the contact pressure distribution caused by the increase of κ1, the other reason may be that the temperature gradient from the surface to the interface becomes relatively small with the increase of κ1.

Effect of the Coating Volumetric Heat Capacity
The effect of the coating volumetric heat capacity c1 on the contact behavior is shown in Figure 9.When the coating volumetric heat capacity c1 increases from 0.5c2 to 2c2, the maximum contact pressure decreases from 1.9pH to 1.59pH and the maximum temperature rise of coated solid decreases to 375 from 519 K.The frictional heat partition coefficient increases with the increase of the coating volumetric heat capacity since a unit volume material of the coating absorbs more heat for 1 K temperature rise with the increase of the coating volumetric heat capacity.That is also the explanation for the influence of the coating volumetric heat capacity on the contact pressure and the temperature rise.Around the contact center, the maximum of the compressive stress σxx on the surface decreases with the increase of c1, while the maximum tensile stress σxx at the base of the coating increases.The maximum of the shear stress σzx on the interface also decreases with the increase of c1, as shown in Figure 9f.

Effect of the Coating Thermal Expansion Coefficient
The effect of the coating thermal expansion coefficient α1 is shown in Figure 10.When the coating thermal expansion coefficient α1 increases from 0.5α2 to 2α2, both the maximum contact pressure and the maximum temperature rise increase, as indicated in Figure 10a,b.However, its effect on the heat

Effect of the Coating Thermal Expansion Coefficient
The effect of the coating thermal expansion coefficient α 1 is shown in Figure 10.When the coating thermal expansion coefficient α 1 increases from 0.5α 2 to 2α 2 , both the maximum contact pressure and the maximum temperature rise increase, as indicated in Figure 10a,b.However, its effect on the heat partition coefficient is very small, which is shown in Figure 10c.Actually, the thermal expansion coefficient itself has no effect on the heat partition coefficient H pc , and the slight decrease of H pc may be due to the concentration of the contact pressure, since nothing of the heat partition model has a direct relationship with the thermal expansion coefficient except that the frictional heat flux distribution is related to the contact pressure distribution.The effect of the coating thermal expansion coefficient on the contact pressure can be easily understood, since the increase of α 1 leads to an increase of the surface thermal displacement in the heated zone.With the increase of α 1 , the increase of the maximum temperature rise is mainly caused by the concentration of the frictional heat flux.This result is consistent with that in [27], although the research in [27] is on a two-dimensional thermoelastic problem.The transverse stress σ xx at the base of the coating is tensile around the contact center, however its maximum is almost the same for various α 1 .This result may be explained by a mutual cancellation of the compressive thermal stress σ xx and the tensile mechanical stress σ xx , since not only the tensile mechanical stress σ xx caused by the mechanical loading increases with the increase of α 1 , but also the compressive thermal stress σ xx caused by the frictional heat flux increases.The maximum shear stress σ zx on the contact interface increases with the increase of α 1 .One reason is the concentration of contact pressure caused by the increase of α 1 , the other reason could be that the thermal shear stress α 1 on the interface caused by the frictional heat flux has the same direction with that caused by the mechanical loading at the position where the maximum shear stress σ zx locates and increases with the increase of α 1 .
Coatings 2018, 8, x FOR PEER REVIEW 13 of 20 mutual cancellation of the compressive thermal stress σxx and the tensile mechanical stress σxx, since not only the tensile mechanical stress σxx caused by the mechanical loading increases with the increase of α1, but also the compressive thermal stress σxx caused by the frictional heat flux increases.The maximum shear stress σzx on the contact interface increases with the increase of α₁.One reason is the concentration of contact pressure caused by the increase of α1, the other reason could be that the thermal shear stress α1 on the interface caused by the frictional heat flux has the same direction with that caused by the mechanical loading at the position where the maximum shear stress σzx locates and increases with the increase of α1.

Effect of the Coating Thickness
For various κ1/κ2, the effect of the coating thickness on the maximum contact pressure, the maximum temperature rise and the heat partition coefficient is show in Figure 11.Among various κ1/κ2, the maximum contact pressure increases fastest with the increase of coating thickness for κ1/κ2 = 0.5.One reason is that the coating-substrate system becomes stiffer with the increase of coating thickness since the elasticity moduli of the coating is two times of the substrate, the other reason is κ1/κ2 = 0.5, which leads to the maximum temperature rise increases with the increase of the coating thickness shown as Figure 11b since the frictional heat flux transfers more slowly into the coated solid with the increase of the coating thickness and less frictional heat flows into the coated solid as indicated in Figure 11c.As a result of the increase of the temperature rise, the surface thermal displacement in the contact zone increases, which will lead to the concentration of the contact

Effect of the Coating Thickness
For various κ 1/ κ 2 , the effect of the coating thickness on the maximum contact pressure, the maximum temperature rise and the heat partition coefficient is show in Figure 11.Among various κ 1/ κ 2 , the maximum contact pressure increases fastest with the increase of coating thickness for κ 1/ κ 2 = 0.5.One reason is that the coating-substrate system becomes stiffer with the increase of coating thickness since the elasticity moduli of the coating is two times of the substrate, the other reason is κ 1/ κ 2 = 0.5, which leads to the maximum temperature rise increases with the increase of the coating Coatings 2018, 8, 470 14 of 20 thickness shown as Figure 11b since the frictional heat flux transfers more slowly into the coated solid with the increase of the coating thickness and less frictional heat flows into the coated solid as indicated in Figure 11c.As a result of the increase of the temperature rise, the surface thermal displacement in the contact zone increases, which will lead to the concentration of the contact pressure and the increase of the maximum contact pressure.On the contrary, the maximum temperature rise of the coated solid with κ 1/ κ 2 = 2 decreases with the increase of the coating thickness as shown in Figure 11b, since the frictional heat transfers more quickly into the coated solid with the increase of the coating thickness and more frictional heat flows into the coated solid as indicated in Figure 11c.As a result of the decrease of the temperature rise, the surface thermal displacement also decreases in the contact zone of the coating-substrate system with κ 1/ κ 2 = 2, therefore a low point of the maximum contact pressure is observed for κ 1/ κ 2 = 2 when the coating thickness is about 0.25a H , as shown in Figure 11a.For coated solids with κ 1/ κ 2 = 1, the increase of the coating thickness leads to a concentration of the contact pressure distribution and an increase of the maximum contact pressure as indicated in Figure 11a, since the coating-substrate system becomes stiffer as a whole.As a result of the concentration of the contact pressure distribution, the maximum temperature rise increases as indicated in Figure 11b, and the heat partition coefficient decreases slightly as indicated in Figure 11c.This result is similar to the effect of the coating thermal expansion coefficient on the temperature rise and the heat partition coefficient as shown in Section 4.4.It also can be observed that the increase of the coating thickness ranging from 0 to a H has a significant effect on the contact behavior of coated solids, however a further increase of coating thickness only brings a marginal effect on the contact behavior when the coating thickness is larger than a H . and the heat partition coefficient as shown in Section 4.4.It also can be observed that the increase of the coating thickness ranging from 0 to aH has a significant effect on the contact behavior of coated solids, however a further increase of coating thickness only brings a marginal effect on the contact behavior when the coating thickness is larger than aH.

Conclusions
Based on the FRFs of coated solids subjected to the mechanical loading coupling with the heat flux on its surface, a three-dimensional thermoelastic contact model of coated solids with the frictional heat partition considered has been developed by introducing a heat partition model, and validated through a comparison of the contact pressure, temperature rise and von Mises stress in the subsurface between the present model and a model in literature.A numerical investigation of the effect of the sliding speed, thermoelastic properties and thickness of the coating on the contact behavior has been conducted by utilizing the present model.The numerical results show that the contact behavior of coated solids can be severely affected by the thermoelastic properties and thickness of the coating as well as the sliding speed.Several specific conclusions can be drawn:


With the increase of the sliding speed, the maximum contact pressure, the maximum temperature rise, the heat partition coefficient, the maximum tensile stress of σxx at the base of the coating and the maximum shear stress σzx on the interface all increase. With the increase of the coating thermal conductivity, the maximum contact pressure, the maximum temperature rise, the maximum tensile stress of σxx at the base of the coating and maximum shear stress σzx on the interface decrease, the heat partition coefficient increases.


With the increase of the coating volumetric heat capacity, the maximum contact pressure, the maximum temperature rise and the maximum shear stress σzx on the interface decrease, the heat partition coefficient and the maximum tensile stress of σxx at the base of the coating increase.


With the increase of the coating thermal expansion coefficient, the maximum contact pressure, the maximum temperature rise and the maximum shear stress σzx on the interface increase, the heat partition coefficient and the maximum tensile stress of σxx at the base of the coating vary

Conclusions
Based on the FRFs of coated solids subjected to the mechanical loading coupling with the heat flux on its surface, a three-dimensional thermoelastic contact model of coated solids with the frictional heat partition considered has been developed by introducing a heat partition model, and validated through a comparison of the contact pressure, temperature rise and von Mises stress in the subsurface between the present model and a model in literature.A numerical investigation of the effect of the sliding speed, thermoelastic properties and thickness of the coating on the contact behavior has been conducted by utilizing the present model.The numerical results show that the contact behavior of coated solids can be severely affected by the thermoelastic properties and thickness of the coating as well as the sliding speed.Several specific conclusions can be drawn:

•
With the increase of the sliding speed, the maximum contact pressure, the maximum temperature rise, the heat partition coefficient, the maximum tensile stress of σ xx at the base of the coating and the maximum shear stress σ zx on the interface all increase.

•
With the increase of the coating thermal conductivity, the maximum contact pressure, the maximum temperature rise, the maximum tensile stress of σ xx at the base of the coating and maximum shear stress σ zx on the interface decrease, the heat partition coefficient increases.

•
With the increase of the coating volumetric heat capacity, the maximum contact pressure, the maximum temperature rise and the maximum shear stress σ zx on the interface decrease, the heat partition coefficient and the maximum tensile stress of σ xx at the base of the coating increase.

•
With the increase of the coating thermal expansion coefficient, the maximum contact pressure, the maximum temperature rise and the maximum shear stress σ zx on the interface increase, the heat partition coefficient and the maximum tensile stress of σ xx at the base of the coating vary slightly.

•
The increase of the coating thickness h ranging from 0 to a H exerts a significant effect on the contact behavior, however a further increase of the coating thickness h only brings a marginal effect.

Figure 1 .
Figure 1.Description of the thermoelastic contact between an elastic ball and a half space coated with a homogeneous coating considering the frictional heat partition.

Figure 1 .
Figure 1.Description of the thermoelastic contact between an elastic ball and a half space coated with a homogeneous coating considering the frictional heat partition.

Figure 2 .
Figure 2. The algorithm flow chart: (a) overall flow chart of the thermal elastic contact model, (b) flow chart of the heat partition model.
SN i r j s z p r s ST i r j s z q r s SQ i r j s z Q r s mn zz yy xx xz yz xy

Figure 2 .
Figure 2. The algorithm flow chart: (a) overall flow chart of the thermal elastic contact model, (b) flow chart of the heat partition model.

(
mn = zz, yy, xx, xz, yz, xy) is the ICs of the surface temperature rise at the grid node [i, j] due to the surface heat flux Q k [r, s].

Figure 3 .
Figure 3. Half space coated with a homogeneous coating subjected to the coupled action of the mechanical loading and the frictional heat flux.

Figure 3 .
Figure 3. Half space coated with a homogeneous coating subjected to the coupled action of the mechanical loading and the frictional heat flux.

Figure 4 .
Figure 4. Comparisons on the contact pressure and the surface temperature rise between the present model and Chen's model under various sliding speeds: (a) contact pressure along x axis, (b) surface temperature rise along x axis.

Figure 5 .
Figure 5.Comparison on the temperature rise in xoz plane between the present model and Chen's model under different sliding speeds: (a) solution of Chen's model, (b) solution of the present model.

Figure 6 .
Figure 6.Comparison on the von Mises stress in xoz plane between the present model and Chen's

Figure 4 .
Figure 4. Comparisons on the contact pressure and the surface temperature rise between the present model and Chen's model under various sliding speeds: (a) contact pressure along x axis, (b) surface temperature rise along x axis.

Figure 4 .
Figure 4. Comparisons on the contact pressure and the surface temperature rise between the present model and Chen's model under various sliding speeds: (a) contact pressure along x axis, (b) surface temperature rise along x axis.

Figure 5 .
Figure 5.Comparison on the temperature rise in xoz plane between the present model and Chen's model under different sliding speeds: (a) solution of Chen's model, (b) solution of the present model.

Figure 6 .
Figure 6.Comparison on the von Mises stress in xoz plane between the present model and Chen's

Figure 5 .
Figure 5.Comparison on the temperature rise in xoz plane between the present model and Chen's model under different sliding speeds: (a) solution of Chen's model, (b) solution of the present model.The comparisons conducted above not only validate the present contact model, but also show that the present model can be utilized to analyze the thermoelastic contact problem of homogeneous solids.

Figure 5 .
Figure 5.Comparison on the temperature rise in xoz plane between the present model and Chen's model under different sliding speeds: (a) solution of Chen's model, (b) solution of the present model.

Figure 6 .
Figure 6.Comparison on the von Mises stress in xoz plane between the present model and Chen's model under different sliding speeds: (a) solution of Chen's model, (b) solution of the present model.

Figure 6 .
Figure 6.Comparison on the von Mises stress in xoz plane between the present model and Chen's model under different sliding speeds: (a) solution of Chen's model, (b) solution of the present model.

Coatings 2018, 8 ,
x FOR PEER REVIEW 11 of 20 thermal shear stress caused by the frictional heat flux has the same direction as that caused by the surface mechanical loading at the position where the maximum σzx locates.

Figure 7 .
Figure 7. Effect of the sliding speed Vs on: (a) contact pressure along x axis, (b) surface temperature rise along x axis, (c) frictional heat partition coefficient, (d) σxx on the surface, (e) σxx at the base of the coating, (f) σzx on the interface.

Figure 8 .
Figure 8.Effect of the coating thermal conductivity κ1 on: (a) contact pressure along x axis, (b) surface temperature rise along x axis, (c) frictional heat partition coefficient, (d) σxx on the surface, (e) σxx at

Figure 7 .
Figure 7. Effect of the sliding speed V s on: (a) contact pressure along x axis, (b) surface temperature rise along x axis, (c) frictional heat partition coefficient, (d) σ xx on the surface, (e) σ xx at the base of the coating, (f) σ zx on the interface.

Figure 8 .
Figure 8.Effect of the coating thermal conductivity κ1 on: (a) contact pressure along x axis, (b) surface temperature rise along x axis, (c) frictional heat partition coefficient, (d) σxx on the surface, (e) σxx at the base of the coating, (f) σzx on the interface.

Figure 8 .
Figure 8.Effect of the coating thermal conductivity κ 1 on: (a) contact pressure along x axis, (b) surface temperature rise along x axis, (c) frictional heat partition coefficient, (d) σ xx on the surface, (e) σ xx at the base of the coating, (f) σ zx on the interface.

Figure 9 .
Figure 9.Effect of the coating volumetric heat capacity c1 on: (a) contact pressure along x axis, (b) surface temperature rise along x axis, (c) frictional heat partition coefficient, (d) σxx on the surface, (e) σxx at the base of the coating, (f) σzx on the interface.

Figure 9 .
Figure 9.Effect of the coating volumetric heat capacity c 1 on: (a) contact pressure along x axis, (b) surface temperature rise along x axis, (c) frictional heat partition coefficient, (d) σ xx on the surface, (e) σ xx at the base of the coating, (f) σ zx on the interface.

Figure 10 .
Figure 10.Effect of the coating thermal expansion coefficient α1 on: (a) contact pressure along x axis, (b) surface temperature rise along x axis, (c) frictional heat partition coefficient, (d) σxx on the surface, (e) σxx at the base of the coating, (f) σzx on the interface.

Figure 10 .
Figure 10.Effect of the coating thermal expansion coefficient α 1 on: (a) contact pressure along x axis, (b) surface temperature rise along x axis, (c) frictional heat partition coefficient, (d) σ xx on the surface, (e) σ xx at the base of the coating, (f) σ zx on the interface.

Figure 11 .
Figure 11.Effect of the coating thickness h on: (a) maximum contact pressure, (b) maximum surface temperature, (c) frictional heat partition coefficient.

Figure 11 .
Figure 11.Effect of the coating thickness h on: (a) maximum contact pressure, (b) maximum surface temperature, (c) frictional heat partition coefficient.