Interfacial Mechanics Analysis of a Brittle Coating–ductile Substrate System Involved in Thermoelastic Contact

In this paper, interfacial stress analysis for a brittle coating/ductile substrate system, which is involved in a sliding contact with a rigid ball, is presented. By combining interface mechanics theory and the image point method, stress and displacement responses within a coated material for normal load, tangential load, and thermal load are obtained; further, the Green's functions are established. The effects of coating thickness, friction coefficient, and a coating's thermoelastic properties on the interfacial shear stress, τ xz , and transverse stress, σ xx , distributions are discussed in detail. A phenomenon, where interfacial shear stress tends to be relieved by frictional heating, is found in the case of a coating material's thermal expansion coefficient being less than a substrate material's thermal expansion coefficient. Additionally, numerical results show that distribution of interfacial stress can be altered and, therefore, interfacial damage can be modified by adjusting a coating's structural parameters and thermoelastic properties.


Introduction
Today, the key elements, such as bearings and gears, are under great pressure to meet the legislative demands of long life and high operational speeds [1][2][3].These challenges can be achieved by employing high hardness, anti-corrosion, and wear-resistant brittle coatings, such as metal nitride coatings, metal oxide coatings, diamond-like carbon (DLC), etc., to protect working surfaces, as well as adopting high strength steel [4,5].However, the huge mismatch between coatings and steel substrates in thermoelastic properties always leads to interfacial stress concentrations, which are induced thermally and mechanically, thus, increasing the risk of coating delamination.To avoid blindness in design and application, an interfacial stress analysis of coated solids in dry sliding contacts is essential.
Thus far, many researchers have established Hertzian contact models for coated, elastic semi-infinite space.Linearly, elastic theory uses various forms of transform integration techniques to produce solutions.Gu et al. [6] studied the stress response of coated materials under surface concentrate loading using the image point method.By comparison of calculation results and scratch tests, the notion that the load bearing capacities of DLC-coated Si 3 N 4 was 1.5-3 times larger than that of DLC-coated M50.O'Sullivan and King studied the contact of coated materials using Papkovich-Neuber potentials [7].Wang and Liu [8,9] applied fast Fourier transformation (FFT) to speed up calculations based on the work of O'Sullivan and King.A combination of FFT and Papkovich-Neuber potentials were further extended to study sliding contacts and partial slip contact problems for solids coated with a monolayer or a gradient layer [10,11].The finite element method (FEM) was also used in studying interfacial adhesion strengths between coatings and substrates [11][12][13][14][15].The Hertzian contact problem of coated materials has been extensively studied during past few decades; and a great number of useful conclusions have been drawn.However, the engineering applications of coatings still have a certain degree of blindness due to the small amount of research on the stress responses of coated materials involved in a thermoelastic contact.Ke et al. [16] discussed the effect of the friction coefficient and sliding velocity on the temperature distribution of semi-infinite materials coated with a gradient material.Choi et al. [17] established a two-dimensional thermoelastic contact model for a rigid flat punch sliding over the surface of a graded coating/substrate system.However, little attention has been paid to the interfacial stress distributions of coated solids involved in thermoelastic contacts.As a matter of fact, a significant amount of frictional heat is generated on the surface of coated materials involved in dry sliding contacts.Thermally induced deformation, subsequently, changes the Hertzian stress distribution of coated materials from the surface to several micrometers in depth, especially for hard coatings, which are distinguished from substrates in regard to thermoelastic properties.Thus, to fulfill a thorough optimization of a hard coating-substrate system, it is necessary to establish an interfacial stress analysis model for thermoelastic contacts.
In this paper, a stress analysis model for coated materials under thermoelastic contact conditions is built, based on interface mechanics.The Green's functions of the stress and displacement fields of coated materials for a normal force, tangential force, and thermal loading are established.With the aim of exploring ways to reduce the risks of delamination of hard coatings from steel substrates, the effects of the friction coefficient, coating thickness, and thermoelastic properties on interfacial shear stress τ xz and transverse stress σ xx distributions are investigated, as these two stresses are believed to be strongly related to the propagation of interface cracks in hard coating-substrate systems [18][19][20].

Description of the Thermoelastic Contact Model
As shown in Figure 1a, a rigid insulated ball is subjected to force F z and slides on the surface of a coated, isotropic, elastic half space with a constant speed, V, in the positive x-direction.
Coatings 2017, 7, 21 2 of 21 in studying interfacial adhesion strengths between coatings and substrates [11][12][13][14][15].The Hertzian contact problem of coated materials has been extensively studied during past few decades; and a great number of useful conclusions have been drawn.However, the engineering applications of coatings still have a certain degree of blindness due to the small amount of research on the stress responses of coated materials involved in a thermoelastic contact.Ke et al. [16] discussed the effect of the friction coefficient and sliding velocity on the temperature distribution of semi-infinite materials coated with a gradient material.Choi et al. [17] established a two-dimensional thermoelastic contact model for a rigid flat punch sliding over the surface of a graded coating/substrate system.However, little attention has been paid to the interfacial stress distributions of coated solids involved in thermoelastic contacts.As a matter of fact, a significant amount of frictional heat is generated on the surface of coated materials involved in dry sliding contacts.Thermally induced deformation, subsequently, changes the Hertzian stress distribution of coated materials from the surface to several micrometers in depth, especially for hard coatings, which are distinguished from substrates in regard to thermoelastic properties.Thus, to fulfill a thorough optimization of a hard coating-substrate system, it is necessary to establish an interfacial stress analysis model for thermoelastic contacts.
In this paper, a stress analysis model for coated materials under thermoelastic contact conditions is built, based on interface mechanics.The Green's functions of the stress and displacement fields of coated materials for a normal force, tangential force, and thermal loading are established.With the aim of exploring ways to reduce the risks of delamination of hard coatings from steel substrates, the effects of the friction coefficient, coating thickness, and thermoelastic properties on interfacial shear stress τxz and transverse stress σxx distributions are investigated, as these two stresses are believed to be strongly related to the propagation of interface cracks in hard coating-substrate systems [18][19][20].

Description of the Thermoelastic Contact Model
As shown in Figure 1a, a rigid insulated ball is subjected to force Fz and slides on the surface of a coated, isotropic, elastic half space with a constant speed, V, in the positive x-direction.The moving coordinate system (x,y,z) is attached to the ball with plane xOy coinciding with the surface of the coated material; that is, the contact area will remain stationary with respect to the movement of the ball.The normal pressure distribution is p(x,y) N•m −2 ; then, the friction force distribution, q(x,y) N•m −2 , and the distribution of frictional heat flux, H(x,y) W•m −2 , are generated within the contact region, with f being the friction coefficient.Their relationships are as follows: The moving coordinate system (x,y,z) is attached to the ball with plane xOy coinciding with the surface of the coated material; that is, the contact area will remain stationary with respect to the movement of the ball.The normal pressure distribution is p(x,y) N•m −2 ; then, the friction force distribution, q(x,y) N•m −2 , and the distribution of frictional heat flux, H(x,y) W•m −2 , are generated within the contact region, with f being the friction coefficient.Their relationships are as follows: Coatings 2017, 7, 21

of 20
Consider this thermoelastic contact problem at room temperature in ambient humidity.Because of the thermal conductive coefficient of the air being always less than 10 −3 times of that parameter for the coating material according to Reference [21], the free surface is assumed to be thermally insult.Due to punch sliding with a constant velocity, the temperature will very quickly approach a steady-state value, thus, the convection term in the heat equation is neglected.The coating and substrate are both isotropic materials, with E c , v c , ζ c , α c , and γ c being the elastic modulus, Poisson ratio, thermal conductive coefficient, thermal expansion coefficient, and thermal diffusion coefficient for the coating material, respectively, and E s , v s , ζ s , α s , and γ s being the same parameters for the substrate material.
According to Reference [22], as shown in Figure 1b, in order to obtain contact pressure p(x,y), an area with [-5a h , 5a h ] × [-5a h , 5a h ] in the x and y directions is uniformly divided into N x × N y surface elements, centered on the grid nodes.Here, a h is the Hertz contact radius without coating and thermal effects, and N x and N y are the number of elements in the x and y directions.As being proved in Reference [23], the contact pressure distribution can be approximated by a piecewise constant function that is uniform within each surface element with errors being less than 0.01% in the case of N x and N y being equal to 32.A more finely meshed surface network, with N x × N y being equal to 128 × 128, is adopted in this paper.With the displacements of the element being represented by its central grid; the contact problem can be expressed as a discretized form: For arbitrarily surface element (i,j), p(i,j) is the normal pressure, g(i,j) is the gap between the surfaces of two bodies in contact after external contact load and thermal load was applied, g 0 (i,j) is the initial gap before external mechanical load and heat flux applied on the ball, u pz (i,j) is the surface displacement of the coated material, which was induced mechanically, u tz (i,j) is the surface displacement of the coated material, which was induced thermally, ∆ x and ∆ y are the discretization lengths in the x and y directions, δ is the normal approach, and Ω c is the contact area.
As illustrated in Figure 2, the ball radius is noted as d.The central point of elements (i,j) is noted as Q, with x i,j and y i,j being the positions of point Q in the xOy plane.Consider this thermoelastic contact problem at room temperature in ambient humidity.Because of the thermal conductive coefficient of the air being always less than 10 −3 times of that parameter for the coating material according to Reference [21], the free surface is assumed to be thermally insult.Due to punch sliding with a constant velocity, the temperature will very quickly approach a steadystate value, thus, the convection term in the heat equation is neglected.The coating and substrate are both isotropic materials, with Ec, vc, ζc, αc, and γc being the elastic modulus, Poisson ratio, thermal conductive coefficient, thermal expansion coefficient, and thermal diffusion coefficient for the coating material, respectively, and Es, vs, ζs, αs, and γs being the same parameters for the substrate material.
According to Reference [22], as shown in Figure 1b, in order to obtain contact pressure p(x,y), an area with [-5ah, 5ah] × [-5ah, 5ah] in the x and y directions is uniformly divided into Nx × Ny surface elements, centered on the grid nodes.Here, ah is the Hertz contact radius without coating and thermal effects, and Nx and Ny are the number of elements in the x and y directions.As being proved in Reference [23], the contact pressure distribution can be approximated by a piecewise constant function that is uniform within each surface element with errors being less than 0.01% in the case of Nx and Ny being equal to 32.A more finely meshed surface network, with Nx × Ny being equal to 128 × 128, is adopted in this paper.With the displacements of the element being represented by its central grid; the contact problem can be expressed as a discretized form: ( , ) 0, ( , ) 0,( , ) ( , ) 0, ( , ) 0,( , ) For arbitrarily surface element (i,j), p(i,j) is the normal pressure, g(i,j) is the gap between the surfaces of two bodies in contact after external contact load and thermal load was applied, g 0 (i,j) is the initial gap before external mechanical load and heat flux applied on the ball, u pz (i,j) is the surface displacement of the coated material, which was induced mechanically, u tz (i,j) is the surface displacement of the coated material, which was induced thermally, Δx and Δy are the discretization lengths in the x and y directions, δ is the normal approach, and Ωc is the contact area.
As illustrated in Figure 2, the ball radius is noted as d.The central point of elements (i,j) is noted as Q, with xi,j and yi,j being the positions of point Q in the xOy plane.The initial gap g 0 (i,j) can be described as the distance between the central point Q and the outer surface of the ball in the z direction: Due to zi,j << d within the contact zone, the relationship between g 0 (i,j) and xi,j,yi,j can be described as The initial gap g 0 (i,j) can be described as the distance between the central point Q and the outer surface of the ball in the z direction: Coatings 2017, 7, 21 4 of 20 Due to z i,j << d within the contact zone, the relationship between g 0 (i,j) and x i,j , y i,j can be described as g 0 (i, j) = x 2 i,j /2d + y 2 i,j /2d (4) Total normal displacement u z (i,j) can be obtained by the accumulation effects of the normal contact pressure and thermal loading of each surface element on the z-displacement of element (i,j).
where G z p and G z t are Green's functions of displacement in the z-direction for a normal load and a thermal load.By applying the complementary energy principle to the thermal contact problem, we have With the aid of the quadratic programming method, the true contact pressure, p (i,j) , can be obtained, at which point contact complementary energy V* is minimal and p (i,j) > 0 is within the contact area.
After gaining contact pressure distribution and heat flux, the stress field can be obtained as follows: where σ AA is an arbitrary stress component of point (x 0 ,y 0 ,z 0 ) within a semi-infinite space (AA may refer to xx, yy, xy, and so on).p (l,k) , q (l,k) , and H (l,k) are contact pressure distribution, tangential distribution, and heat flux distribution on surface elements (l,k), respectively.G AA p , G AA v , and G AA t are the values of Green's functions for components AA for a normal load, tangential load, and thermal load at point (x 0 ,y 0 ,z 0 ), respectively.For an elastic solid with a coating bonded to its surface, the explicit expressions for G AA p , G AA v , and G AA t are not available.However, based on interface mechanics, they can be deduced and expressed in the accumulated form by employing the image point method [24,25].The detailed deducing is presented in Section 2.2, Section 2.3, and Section 2.4.In this paper, three orders of image points are used to achieve a high accuracy as is proven in Section 3.1.

Green's Function for a Point Normal Load Acting on the Surface of a Coated Isotropic Thermoelastic Material
In the analytical model of interface mechanics (Figure 3), a straight interface is formed along the border of the coating and substrate.From Figure 3, we can see that coating I, with thickness h, is used to cover substrate II, and is connected only by the interface.Above the coating, there is a free surface, and a force, marked p z , is applied at point O 1 on said surface and is along the z-axis.Considering the symmetries of this problem, a cylindrical coordinate (r,θ,z) is chosen with the r-axis being along the interface, and the z-axis being perpendicular to the r-axis and passing through O 1 .The two axes intersect at O, namely the global point of origin.Coordinate plane rOθ coincides with the interface.
The interface stress, shown in Figure 3, can be calculated using the images method from the complex variable function.The boundary conditions and the constraints of the model are equivalent; furthermore, equations were composed using the relationship between stress and strain.As a method of images, the interface and the surface are imagined as mirrors that reflect points O or O 1 , and will generate infinite points of mirror images.These images influence the interface stresses as a superposition forms, in order to fulfill the boundary conditions and constraints.Namely, the conditions of the interface continuum and free surface are satisfied and interface stresses are gained by superposing the stress solutions of these image points.The image points in the coating space are marked as O k and The relationship between the global coordinate and the local coordinates can be written as: Love's stress and displacement functions can be written as: (1 2 ) where μ is the shear elasticity modulus, and its relationships with E and v is in the form of μ = 0.5E/(1 + v).u is the displacement component and σ is the stress component in each directions.φ and Ψ are harmonic functions.For coating materials, the harmonic functions are noted as φ 1 and Ψ 1 .For coating materials, the harmonic functions are noted as φ 2 and Ψ 2 .Due to the singular point O1 on the surface of the coating, the harmonic functions for coating material have a relationship with zi and si at the same time.However, because no singular point exists in the substrate material, the harmonic functions for the substrate material only have a relationship with zi.Thus, the harmonic functions can be written as where i is the sequence number of the image point.
One of the boundary conditions can be established by satisfying the conditions of the interface continuum.
Another boundary condition can be established by satisfying the free surface.The relationship between the global coordinate and the local coordinates can be written as: Love's stress and displacement functions can be written as: where µ is the shear elasticity modulus, and its relationships with E and v is in the form of µ = 0.5E/(1 + v).u is the displacement component and σ is the stress component in each directions.ϕ and Ψ are harmonic functions.For coating materials, the harmonic functions are noted as ϕ 1 and Ψ 1 .For coating materials, the harmonic functions are noted as ϕ 2 and Ψ 2 .Due to the singular point O 1 on the surface of the coating, the harmonic functions for coating material have a relationship with z i and s i at the same time.However, because no singular point exists in the substrate material, the harmonic functions for the substrate material only have a relationship with z i .Thus, the harmonic functions can be written as where i is the sequence number of the image point.
One of the boundary conditions can be established by satisfying the conditions of the interface continuum.
Another boundary condition can be established by satisfying the free surface.
Thus, ϕ 1 1 (r,z 1 ) and Ω 1 1 (r,z 1 ) are the only factors needed to finish the cyclic recursion.Obviously, the Midlin solutions for a point normal load, acting on the surface of a semi-infinite space, satisfy the requirements.
Now, ϕ and Ψ can be obtained using Equations ( 13)-( 15).Green's functions can be obtained by substituting Φ and Ψ into Equation ( 9), and then separating variable p z .

Green's Function for a Point Tangential Load Acting on the Surface of a Coated Isotropic Thermoelastic Material
As shown in Figure 4, we can see that coating I, with thickness h, is used to cover substrate II, and is connected only by the interface.Above the coating, there is a free surface, and a tangential force, marked p x , is applied at point O 1 on said surface and is parallel to the x-axis.A Cartesian coordinate (x,y,z) is chosen with the x-axis and y-axis being along the interface, and the z-axis being perpendicular to the interface and passing through O 1 .The x-axis and y-axis intersect at O, namely the global point of origin.Coordinate plane xOy coincides with the interface.The image points in the coating space are marked as O k , and those in the substrate space are marked as C k .Green's function for a point tangential load, acting on the surface of a coated isotropic thermoelastic material, can be obtained in the similar way (shown in Section 2.2).Papokovitch's stress and displacement functions can be written as: 4(1 ) Papokovitch's stress and displacement functions can be written as: One of the recursive relations can be expressed as: ) ) Another recursive relation can be expressed as: where µ 2 /µ 1 is noted as Γ and (1 The Midlin solutions for a point tangential load, acting on the surface of a semi-infinite space, can be written as: (19) where p x is a tangential load.Now, B 1 , B 2 , B 3 and δ can be obtained by using Equation ( 17)-( 19).Green's functions can be obtained by substituting B 1 , B 2 , B 3 , and δ into Equation ( 16) and then separating variable p x .

Green's Function for a Moving Point Heat Resource Acting on the Surface of a Coated Isotropic Thermoelastic Material
As shown in Figure 5, we can see that coating I, with thickness h, is used to cover substrate II, and is connected only by the interface.Above the coating, there is a free surface, and a moving point heat resource, marked H, is applied at point O 1 on said surface.A cylindrical coordinate (r,θ,z) is chosen with the r-axis being along the interface, and the z-axis being perpendicular to the r-axis and passing through O 1 .The r-axis and z-axis intersect at O, namely, the global point of origin.Coordinate plane rOθ coincides with the interface.The image points in the coating space are marked as O k and those in the substrate space are marked as C k .The cyclic recursions and stress and displacement functions for point heat resource are given here.

Green's Function for a Moving Point Heat Resource Acting on the Surface of a Coated Isotropic Thermoelastic Material
As shown in Figure 5, we can see that coating I, with thickness h, is used to cover substrate II, and is connected only by the interface.Above the coating, there is a free surface, and a moving point heat resource, marked H, is applied at point O1 on said surface.A cylindrical coordinate (r,θ,z) is chosen with the r-axis being along the interface, and the z-axis being perpendicular to the r-axis and passing through O1.The r-axis and z-axis intersect at O, namely, the global point of origin.Coordinate plane rOθ coincides with the interface.The image points in the coating space are marked as Ok and those in the substrate space are marked as Ck.The cyclic recursions and stress and displacement functions for point heat resource are given here.The general solution for displacement and stress within a semi-infinite space for a thermal load can be written as: The general solution for displacement and stress within a semi-infinite space for a thermal load can be written as: where α is thermal expansion coefficient.λ is the Lame constant [26], its relationships with E and v are in the form of λ = vE/(1 + v)(1 − 2v).E t is the thermal modulus, its relationships with E, v, and α are in the form of E t = 2(λ + 2µ)/α.Harmonic functions in the coating material can be written as: Coatings 2017, 7, 21 Harmonic functions in the substrate material can be written as , (j = 2, 3) One of the boundary conditions can be established by satisfying the continuous conditions of displacement, stress, temperature, and thermal flux on the interface.
Another boundary conditions can be established by satisfying the stress free and thermally insult conditions on the surface.
Additionally, mechanical equilibrium and thermal equilibrium of an infinite plane, ε < z < h, should be satisfied.
where H is the heat flux.As in Grylitsky and Paul [27], we also introduce the Pe number into the above-mentioned equation for the thermoelastic problem.λ 2 is the thermal diffusion coefficient for the substrate material.

Solution Procedure
Figure 6 is a flowchart for solving thermoelastic contact problem of coated solids.The solving procedure can be conducted as follows: • Initialize material parameters and structural parameters.Calculate Hertz contact radius a h , Hertz contact pressure c h , Hertz contact approach δ 0 , and initial gap g 0 .

•
According to the material parameters and the coating thickness, construct a matrix of Green's functions for normal load G p AA using Equation ( 9) and Equations ( 13)-( 15); construct matrix of Green's function for thermal load G t AA using Equations ( 20)-( 27); construct matrix of Green's function for tangential load G v AA using Equations ( 16)-( 19).

•
Employ G p AA and the quadratic programming method, setting g iter0 = g 0 , calculate Hertzian contact pressure distribution p c using Equations ( 2)-( 6).Set p iter0 = p c .

•
Obtain surface heat flux H iter , according to Equation ( 1) and p iter .Employing G t AA and H iter , calculate thermally induced surface displacement u t iter .Update gap g iter = g iter + u t iter .

•
Save p iter as p 0 , then update surface contact pressure p iter by employing g iter , G p AA , and the quadratic programming method.
Green's functions for the thermal load can be obtained by substituting ψ1, ψ2, ψ3 and Ψ1, Ψ2, Ψ3 into Equation (20) and separating variable H.

Solution Procedure
Figure 6 is a flowchart for solving thermoelastic contact problem of coated solids.The solving procedure can be conducted as follows:


Initialize material parameters and structural parameters.Calculate Hertz contact radius ah, Hertz contact pressure ch, Hertz contact approach δ 0 , and initial gap g 0 . According to the material parameters and the coating thickness, construct a matrix of Green's functions for normal load G p AA using Equation ( 9) and Equations ( 13)-( 15); construct matrix of Green's function for thermal load G t AA using Equations ( 20)-( 27); construct matrix of Green's function for tangential load G v AA using Equations ( 16)-( 19). Employ G p AA and the quadratic programming method, setting g iter0 = g 0 , calculate Hertzian contact pressure distribution p c using Equations ( 2)-( 6).Set p iter0 = p c . Obtain surface heat flux H iter , according to Equation ( 1) and p iter .Employing G t AA and H iter , calculate thermally induced surface displacement u t iter.Update gap g iter = g iter + u t iter. Save p iter as p0, then update surface contact pressure p iter by employing g iter , G p AA, and the quadratic programming method. Calculate ε iter = ΔxΔy| p iter − p0| to judge convergence.If ε iter < ε0, go to step 7, otherwise go to step 4 for the next iteration. Obtain contact pressure p and heat flux H. Calculate σAA using Equation ( 7) by employing GAA p , GAA v , and GAA t .

Results and Discussion
The ratio of shear modulus between the coating and substrate is µ = µ c /µ s , the ratio of thermal expansion coefficient between the coating and substrate is α = α c /α s , and the ratio of the thermal conductive coefficient between the coating and substrate is ζ = ζ c /ζ s .The following non-dimensional temperature rise, stress components, and coordinates are used: where Pe 2 is the Pecelet Number, λ s is the thermal diffusion coefficient for the substrate material, and p n is the non-dimensional contact pressure.t n is the non-dimensional temperature rise.σ n-AA is the non-dimensional stress components in an arbitrary direction.a h and c h are the Hertz contact radius and max contact pressure.Coating and substrate material's Poisson ratios are set to be 0.3.In this study, the substrate material is 52100 steel with E s = 210 GPa, α s = 1.17 × 10 and λ s = 1.0 × 10 −5 m 2 •s −1 .Max Hertz contact pressure c h is 1.09 GPa and the Hertz contact radius, a h , is 100 µm.Sliding velocity, V, is 2.4 m•s −1 ; thus, Pe 2 is 2.4.

Model Verification
In order to verify the validity of the present model, the numerical results, obtained by setting µ = 1, α = 1, and ζ = 1, are compared with the exact analytical solution given by Johnson and Midlin [21] for the thermoelastic contact problem of solo material semi-space at f = 0.2 and Pe 2 = 2.4.The contour plot of the Von Mises stress and temperature on the xz plane are shown in Figure 7a,b, respectively, with numerical results being implied by the dotted line and the analytical results being implied by the solid line.All numerical simulations agree with those from the analytical solutions.

Results and Discussion
The ratio of shear modulus between the coating and substrate is μ = μc/μs, the ratio of thermal expansion coefficient between the coating and substrate is α = αc/αs, and the ratio of the thermal conductive coefficient between the coating and substrate is ζ = ζc/ζs.
where Pe2 is the Pecelet Number, λs is the thermal diffusion coefficient for the substrate material, and pn is the non-dimensional contact pressure.tn is the non-dimensional temperature rise.σn-AA is the non-dimensional stress components in an arbitrary direction.ah and ch are the Hertz contact radius and max contact pressure.Coating and substrate material's Poisson ratios are set to be 0.3.In this study, the substrate material is 52100 steel with Es = 210 GPa, αs = 1.17 × 10 and λs = 1.0 × 10 −5 m 2 •s −1 .Max Hertz contact pressure ch is 1.09 GPa and the Hertz contact radius, ah, is 100 μm.Sliding velocity, V, is 2.4 m•s −1 ; thus, Pe2 is 2.4.

Model Verification
In order to verify the validity of the present model, the numerical results, obtained by setting μ = 1, α = 1, and ζ = 1, are compared with the exact analytical solution given by Johnson and Midlin [21] for the thermoelastic contact problem of solo material semi-space at f = 0.  Note that, for verifying the validity of the present model further, the work by Wang [23] for the contact problem of a coated solids without thermal effect can be recovered by assuming V = 0 in the present model.Here, we choose the same material parameters as in Reference [9] and calculate the contact pressure distribution.Figure 8 shows the dimensionless contact stress at h/ah = 0.5, f = 0.2, V = 0 with μc/μs varying from 2 to 4. It is seen that the present results agree very well with Liu's results.Note that, for verifying the validity of the present model further, the work by Wang [23] for the contact problem of a coated solids without thermal effect can be recovered by assuming V = 0 in the present model.Here, we choose the same material parameters as in Reference [9] and calculate the contact pressure distribution.Figure 8 shows the dimensionless contact stress at h/a h = 0.5, f = 0.2, V = 0 with µ c /µ s varying from 2 to 4. It is seen that the present results agree very well with Liu's results.

Effects of Friction Coefficient and Coating Thickness on Contact Pressure and Interface Stress Distribution
The effects of the friction coefficient on contact pressure and surface temperature are shown in Figure 9.With the increase of f, more frictional heats are generated within the contact zone, and the thermal expansion of the contact surface in the z-direction is enlarged.As shown in Figure 9b, this leads to a decrease in the contact area by 23% and an increase in maximum contact pressure by 307% when f increases from 0.1 to 0.5.A comment to be made regarding Figure 9a is that the magnitude of the surface temperature in proportion to f is offset by greater concentrations in the tail of the contact zone due to the moving of heat resources.Because of the higher thermal expansion in the tail of the contact zone, contact pressure is more concentrated in this area with the increase of f, as shown in Figure 9b.

Effects of Friction Coefficient and Coating Thickness on Contact Pressure and Interface Stress Distribution
The effects of the friction coefficient on contact pressure and surface temperature are shown in Figure 9.With the increase of f, more frictional heats are generated within the contact zone, and the thermal expansion of the contact surface in the z-direction is enlarged.As shown in Figure 9b, this leads to a decrease in the contact area by 23% and an increase in maximum contact pressure by 307% when f increases from 0.1 to 0.5.A comment to be made regarding Figure 9a is that the magnitude of the surface temperature in proportion to f is offset by greater concentrations in the tail of the contact zone due to the moving of heat resources.Because of the higher thermal expansion in the tail of the contact zone, contact pressure is more concentrated in this area with the increase of f, as shown in Figure 9b.

Effects of Friction Coefficient and Coating Thickness on Contact Pressure and Interface Stress Distribution
The effects of the friction coefficient on contact pressure and surface temperature are shown in Figure 9.With the increase of f, more frictional heats are generated within the contact zone, and the thermal expansion of the contact surface in the z-direction is enlarged.As shown in Figure 9b, this leads to a decrease in the contact area by 23% and an increase in maximum contact pressure by 307% when f increases from 0.1 to 0.5.A comment to be made regarding Figure 9a is that the magnitude of the surface temperature in proportion to f is offset by greater concentrations in the tail of the contact zone due to the moving of heat resources.Because of the higher thermal expansion in the tail of the contact zone, contact pressure is more concentrated in this area with the increase of f, as shown in Figure 9b.The effect of the friction coefficient on shear stress τn−xz along the intersecting line of plane xoz and the interface is shown in Figure 10, with variations of h/ah at μ = 2, α = 0.3, ζ = 0.5.With the increase of f, the maximum interfacial shear stress increases sharply for coatings with h/ah = 0.02, as shown in Figure 10a, and moderately for coatings with h/ah = 0.25 and h/ah = 0.5, as shown in Figure 10c,d.For coatings with h/ah = 0.1, the maximum interfacial shear stress decreases slightly, and then increases, as shown in Figure 10b.The maximum interfacial shear stress of coated solids with thin coatings (e.g., h/ah = 0.02) is more sensitive to the change of friction coefficient than those that have The effect of the friction coefficient on shear stress τ n−xz along the intersecting line of plane xoz and the interface is shown in Figure 10, with variations of h/a h at µ = 2, α = 0.3, ζ = 0.5.With the increase of f, the maximum interfacial shear stress increases sharply for coatings with h/a h = 0.02, as shown in Figure 10a, and moderately for coatings with h/a h = 0.25 and h/a h = 0.5, as shown in Figure 10c,d.For coatings with h/a h = 0.1, the maximum interfacial shear stress decreases slightly, and then increases, as shown in Figure 10b.The maximum interfacial shear stress of coated solids with thin coatings (e.g., h/a h = 0.02) is more sensitive to the change of friction coefficient than those that have thick coatings.Traditionally, employing thin coatings is preferable, as lower interfacial shear stress can be obtained as soon as f < 0.3.However, due to being affected by frictional heat, thin coatings (e.g., h/a h = 0.02) are only favored when f is below 0.1.The effect of the friction coefficient on transverse stress σn−xx along the intersecting line of plane xoz and the coating bottom can be inferred from Figure 11, which shows various h/ah at μ = 2, α = 0.3, ζ = 0.5.As shown in Figure 11c,d, coatings with a thickness of h/ah > 0.25 may render interfacial transverse stress σn−xx compressive within the contact zone for various f, counteracting the brittle failure of the coating in general.Meanwhile, thin coatings with a thickness of h/ah = 0.02 may encounter brittle failure at the end of the tail of the contact zone because σn−xx is always tensile in that zone, in all cases of f.For coatings with h/ah = 0.1, brittle failure at the end of the tail of the contact zone may happen in the case where f > 0.2.
By compromising between interfacial shear stress and tensile transverse stress, coatings with thicknesses of about h/ah = 0.02 are preferred in cases where f < 0.1; coatings with a thickness of about h/ah = 0.1 are preferred in the case of 0.1 < f < 0.3, and coatings with a thickness of about h/ah = 0.25 are proposed in cases where f > 0.3.The effect of the friction coefficient on transverse stress σ n−xx along the intersecting line of plane xoz and the coating bottom can be inferred from Figure 11, which shows various h/a h at µ = 2, α = 0.3, ζ = 0.5.As shown in Figure 11c,d, coatings with a thickness of h/a h > 0.25 may render interfacial transverse stress σ n−xx compressive within the contact zone for various f, counteracting the brittle failure of the coating in general.Meanwhile, thin coatings with a thickness of h/a h = 0.02 may encounter brittle failure at the end of the tail of the contact zone because σ n−xx is always tensile in that zone, in all cases of f.For coatings with h/a h = 0.1, brittle failure at the end of the tail of the contact zone may happen in the case where f > 0.2.
By compromising between interfacial shear stress and tensile transverse stress, coatings with thicknesses of about h/a h = 0.02 are preferred in cases where f < 0.1; coatings with a thickness of about h/a h = 0.1 are preferred in the case of 0.1 < f < 0.3, and coatings with a thickness of about h/a h = 0.25 are proposed in cases where f > 0.3.

Effect of a Coating's Elastic Modulus and Thermal Expansion Coefficient on Interface Stress Distribution
The effects of a coatings' elastic modulus ratio, μ, on shear stress τn−xz along the intersecting line of plane xoz and the interface are shown in Figure 12, with variations of α at h/ah = 0.02, ζ = 0.5, f = 0.1.As can be inferred from Figure 12c,d, the increase in stiffness of the coating results in a higher interfacial shear stress when α is greater than 1.However, when α is less than 1, the interfacial shear stress decreases dramatically with the increase in the stiffness of the coatings, as can be inferred from Figure 12a,b.With the increase of μ, thermally induced interfacial shear stress τt−xz always increases, as shown in Figure 13.However, the direction of the thermally induced interfacial shear stress component is opposite to that caused by tangential traction when α < 1, as is shown in Figure 13a.Thus, the increase in the thermally induced interfacial shear stress results in the decrease of total interfacial shear stress.While α > 1, the direction of thermally induced interfacial shear stress is the same as that caused by tangential traction, as shown in Figure 13b.Thus, the increase in the interfacial shear stress, which is induced by fractional heat, results in the increase of total interfacial shear stress under the conditions of α > 1.
Another revealed feature is the effects of a coatings' elastic modulus and thermal expansion coefficient on the tensile transverse stress along the intersecting line of plane xoz and the interface, as shown in Figure 14, with variations of α and μ at h/ah = 0.02, ζ = 0.5, and f = 0.1.Figure 14b-d indicate that coatings with a thermal expansion coefficient of α > 0.6 may render the interfacial transverse stress σn−xx compressive within the contact area on the interface for variations of μ, counteracting, in general, the brittle failure of the coating.

Effect of a Coating's Elastic Modulus and Thermal Expansion Coefficient on Interface Stress Distribution
The effects of a coatings' elastic modulus ratio, µ, on shear stress τ n−xz along the intersecting line of plane xoz and the interface are shown in Figure 12, with variations of α at h/a h = 0.02, ζ = 0.5, f = 0.1.As can be inferred from Figure 12c,d, the increase in stiffness of the coating results in a higher interfacial shear stress when α is greater than 1.However, when α is less than 1, the interfacial shear stress decreases dramatically with the increase in the stiffness of the coatings, as can be inferred from Figure 12a,b.With the increase of µ, thermally induced interfacial shear stress τ t−xz always increases, as shown in Figure 13.However, the direction of the thermally induced interfacial shear stress component is opposite to that caused by tangential traction when α < 1, as is shown in Figure 13a.Thus, the increase in the thermally induced interfacial shear stress results in the decrease of total interfacial shear stress.While α > 1, the direction of thermally induced interfacial shear stress is the same as that caused by tangential traction, as shown in Figure 13b.Thus, the increase in the interfacial shear stress, which is induced by fractional heat, results in the increase of total interfacial shear stress under the conditions of α > 1.
Another revealed feature is the effects of a coatings' elastic modulus and thermal expansion coefficient on the tensile transverse stress along the intersecting line of plane xoz and the interface, as shown in Figure 14, with variations of α and µ at h/a h = 0.02, ζ = 0.5, and f = 0.1.Figure 14b-d indicate that coatings with a thermal expansion coefficient of α > 0.6 may render the interfacial transverse stress σ n−xx compressive within the contact area on the interface for variations of µ, counteracting, in general, the brittle failure of the coating.Thus, in order to obtain a lower interfacial shear stress and compressive transverse stress at the same time, 0.6 < α c /α s < 1 can be recommended as the goal in order to optimize the thermoelastic properties of thin hard coatings.Thus, in order to obtain a lower interfacial shear stress and compressive transverse stress at the same time, 0.6 < αc/αs < 1 can be recommended as the goal in order to optimize the thermoelastic properties of thin hard coatings.Thus, in order to obtain a lower interfacial shear stress and compressive transverse stress at the same time, 0.6 < αc/αs < 1 can be recommended as the goal in order to optimize the thermoelastic properties of thin hard coatings.

Effect of a Coating's Thermal Conductiveness and Thermal Expansion Properties on Interface Stress Distribution
The effects of a coatings' thermal conductive ratio, ζ, on shear stress τn−xz along the intersecting line of plane xoz and the interface are shown in Figure 15, with variations of α at h/ah = 0.02, μ = 2, and f = 0.1.The results shown in Figure 15 confirm that the increase in the thermal conductive coefficient of a coating is always beneficial in reducing the interfacial shear stress.As mentioned above, contact pressure distribution is co-determined by contact load and surface temperature rise in thermoelastic contact conditions.Figure 16 illustrates that the increase in ζ results in decreases in contact pressure and interfacial shear stress.With respect to Figure 16a, the magnitude of the contact pressure grows in an inverse proportion to ζ, and the magnitude of the contact area is proportional to ζ.The fact that coatings with a higher thermal conductivity can facilitate heat transfer within a coated system is understood to be accountable for contact pressure relaxation.With the dispersing of surface shear traction, interfacial shear stress decreases substantially, as shown in Figure 16b.
The effects of a coatings' thermal conductive coefficient and thermal expansion coefficient on the tensile transverse stress along the intersecting line of plane xoz and the interface are provided in Figure 17, with variations of α and ζ at h/ah = 0.02, μ = 2, and f = 0.1.Figure 17a-d

Effect of a Coating's Thermal Conductiveness and Thermal Expansion Properties on Interface Stress Distribution
The effects of a coatings' thermal conductive ratio, ζ, on shear stress τ n−xz along the intersecting line of plane xoz and the interface are shown in Figure 15, with variations of α at h/a h = 0.02, µ = 2, and f = 0.1.The results shown in Figure 15 confirm that the increase in the thermal conductive coefficient of a coating is always beneficial in reducing the interfacial shear stress.As mentioned above, contact pressure distribution is co-determined by contact load and surface temperature rise in thermoelastic contact conditions.Figure 16 illustrates that the increase in ζ results in decreases in contact pressure and interfacial shear stress.With respect to Figure 16a, the magnitude of the contact pressure grows in an inverse proportion to ζ, and the magnitude of the contact area is proportional to ζ.The fact that coatings with a higher thermal conductivity can facilitate heat transfer within a coated system is understood to be accountable for contact pressure relaxation.With the dispersing of surface shear traction, interfacial shear stress decreases substantially, as shown in Figure 16b.
The effects of a coatings' thermal conductive coefficient and thermal expansion coefficient on the tensile transverse stress along the intersecting line of plane xoz and the interface are provided in Figure 17, with variations of α and ζ at h/a h = 0.02, µ = 2, and f = 0.1.Figure 17a-d     Thus, in order to obtain a lower interfacial shear stress and compressive transverse stress at the same time, a higher value of ζc/ζs can be recommended in order to improve the resistance to interfacial delamination of the coating material, when subjected to sliding contact with friction heat generation.
The quasi-static thermoelastic contact problem between a ball and a coated solid is discussed in the present study.The contact mode between the two bodies is assumed to be conformal, with coating and substrate all being isotropic and without defects on the interface.In considering those assumptions, the present model can be applied to the brittle coating-ductile substrate system, with the coating's elastic modulus being greater than the elastic modulus of the substrate material.The diameter of the ball should be 20 times greater than the width of the contact area, and the velocity of the ball should be less than 100 m•s −1 .The work present herein offers room for further viable and promising extensions.For instance, the assumptions made above could be relaxed to cover a broader range of applications.Moreover, in the sense that the frictional contact on the surface of a brittle coating may eventually give rise to damage patterns in the form of interface cracking, the feasible coupled crack/contact analysis in the brittle coating-ductile substrate system incorporating the contact-induced frictional heating effect would pose another interesting research topic of technological significance.Such issues need to be addressed and will be reported in forthcoming papers.The analytical work in this paper can help in the design of the coating-substrate system with a proper coating thickness to obtain a relative low interface stress.In addition, the analytical work in this paper can propose targets on the modification of the coating material in elastic and thermoelastic properties.In the sense that the thermoelastic properties of the coating material have effects on the loading bearing capacity of the coating-substrate system.Tribology tests, such as the ball-on-disk test, with relative velocity being similar to real working conditions, could be applied to study the Thus, in order to obtain a lower interfacial shear stress and compressive transverse stress at the same time, a higher value of ζ c /ζ s can be recommended in order to improve the resistance to interfacial delamination of the coating material, when subjected to sliding contact with friction heat generation.
The quasi-static thermoelastic contact problem between a ball and a coated solid is discussed in the present study.The contact mode between the two bodies is assumed to be conformal, with coating and substrate all being isotropic and without defects on the interface.In considering those assumptions, the present model can be applied to the brittle coating-ductile substrate system, with the coating's elastic modulus being greater than the elastic modulus of the substrate material.The diameter of the ball should be 20 times greater than the width of the contact area, and the velocity of the ball should be less than 100 m•s −1 .The work present herein offers room for further viable and promising extensions.For instance, the assumptions made above could be relaxed to cover a broader range of applications.Moreover, in the sense that the frictional contact on the surface of a brittle coating may eventually give rise to damage patterns in the form of interface cracking, the feasible coupled crack/contact analysis in the brittle coating-ductile substrate system incorporating the contact-induced frictional heating effect would pose another interesting research topic of technological significance.Such issues need to be addressed and will be reported in forthcoming papers.The analytical work in this paper can help in the design of the coating-substrate system with a proper coating thickness to obtain a relative low interface stress.In addition, the analytical work in this paper can propose targets on the modification of the coating material in elastic and thermoelastic properties.In the sense that the thermoelastic properties of the coating material have effects on the loading bearing capacity of the coating-substrate system.Tribology tests, such as the ball-on-disk test, with relative velocity being similar to real working Coatings 2017, 7, 21 19 of 20 conditions, could be applied to study the load bearing capacity of the coating-substrate system due to its reflecting the real thermal load condition and contact condition well.

Conclusions
This paper investigated the interfacial mechanical characteristics of a hard coating/substrate system, which was involved in a thermoelastic contact.The image point method and interface mechanics were employed in order to obtain Green's functions for a normal load, tangential load, and thermal load.Based on Green's functions and the contact mechanics, the temperature field and stress field were obtained with a high degree of accuracy.The effects of coating thickness, friction coefficient, and thermoelastic properties on the interfacial stress characteristics were discussed.It was found that:

•
Adequate design of coating thickness for different friction coefficients can be help in obtaining relatively low interfacial shear stress and tensile transverse stress at the same time.For brittle coatings (E c > E s ), in the case where the friction coefficient is less than 0.1, a smaller coating thickness, with h/a h < 0.02, is recommended.With a friction coefficient between 0.1 and 0.3, a coating with a thickness h/a h of about 0.1 is proposed.When the friction coefficient is greater than 0.3, coating with a thickness h/a h above 0.25 is preferable.

•
Interfacial shear stress tends to be relieved by frictional heating in the case of α c /α s < 1 for brittle coatings.Meanwhile, tensile transverse stress on the interface tends to be eliminated in the case of α c /α s > 0.6.A thermal expansion coefficient ratio of 0.6 < α c /α s < 1 can be recommended as a goal for the optimization of thin hard coatings-substrate systems.

•
Coating materials with enhanced thermal conductivity are shown to lower surface contact pressure concentration, and, thus, decrease interfacial shear stress.Although the interfacial tensile transverse stress is affected, to a negligible extent, by the variations of the thermal conductive coefficient of coating materials, it is predictable that the increased thermal conductivity of a coating could improve the resistance to interfacial delamination of the coating material.
Finally, we believe that this paper provide an understanding of the mechanism of modifying interfacial failure for hard coatings by employing adequate analysis.

Figure 1 .
Figure 1.Schematic of the thermoelastic contact problem.(a) A rigid insulated ball sliding over the surface of a coated solid; (b) Contact area meshing.

Figure 1 .
Figure 1.Schematic of the thermoelastic contact problem.(a) A rigid insulated ball sliding over the surface of a coated solid; (b) Contact area meshing.

Figure 2 .
Figure 2. Schematic description of the initial gap between surface elements and ball.

Figure 2 .
Figure 2. Schematic description of the initial gap between surface elements and ball.

Figure 3 .
Figure 3.A point normal force pz acting on the surface of coated material.

Figure 3 .
Figure 3.A point normal force p z acting on the surface of coated material.

Figure 4 .
Figure 4.A point tangential force px acting on the surface of coated material.

Figure 4 .
Figure 4.A point tangential force p x acting on the surface of coated material.

Figure 5 .
Figure 5.A moving point heat resource H acting on the surface of coated material.

Figure 5 .
Figure 5.A moving point heat resource H acting on the surface of coated material.

Figure 6 .
Figure 6.Flowchart for solving the thermoelastic contact problem of a coated solid.

Figure 6 .
Figure 6.Flowchart for solving the thermoelastic contact problem of a coated solid.
2 and Pe2 = 2.4.The contour plot of the Von Mises stress and temperature on the xz plane are shown in Figure 7a,b, respectively, with numerical results being implied by the dotted line and the analytical results being implied by the solid line.All numerical simulations agree with those from the analytical solutions.

Figure 7 .
Figure 7.Comparison between the numerical solution and the analytical solution at μ = 1, α = 1, ζ = 1 and f = 0.2; Pe2 = 2.4.(a) Contour plot of σn−v obtained using the present model and analytical solution.(b) Contour plot of tn obtained by the present model and analytical solution.

Figure 7 .
Figure 7.Comparison between the numerical solution and the analytical solution at µ = 1, α = 1, ζ = 1 and f = 0.2; Pe 2 = 2.4.(a) Contour plot of σ n−v obtained using the present model and analytical solution.(b) Contour plot of t n obtained by the present model and analytical solution.

Figure 8 .
Figure 8.Comparison of contact pressure obtained from Wang's result and the present model at h/ah = 0.5, f = 0.2, V = 0 with μc/μs varying from 2 to 4.

Figure 9 .
Figure 9.Effect of f on contact characteristics.(a) Effect of f on contact pressure distribution; (b) Effect of f on surface temperature distribution.

Figure 8 .
Figure 8.Comparison of contact pressure obtained from Wang's result and the present model at h/a h = 0.5, f = 0.2, V = 0 with µ c /µ s varying from 2 to 4.

Figure 9 .
Figure 9.Effect of f on contact characteristics.(a) Effect of f on contact pressure distribution; (b) Effect of f on surface temperature distribution.

Figure 9 .
Figure 9.Effect of f on contact characteristics.(a) Effect of f on contact pressure distribution; (b) Effect of f on surface temperature distribution.

Figure 10 .
Figure 10.τ n−xz along the intersecting line of plane xoz and interface with various f and h/a h when µ = 2, α = 0.3, ζ = 0.5.(a) With various f at h/a h = 0.02.(b) With various f at h/a h = 0.1.(c) With various f at h/a h = 0.25.(d) With various f at h/a h = 0.5.

Figure 11 .
Figure 11.σ n−xx along the intersecting line of plane xoz and interface with various f and h/a h when µ = 2, α = 0.3, ζ = 0.5.(a) with various f at h/a h = 0.02; (b) with various f at h/a h = 0.1; (c) with various f at h/a h = 0.25; (d) with various f at h/a h = 0.5.

Figure 13 .Figure 14 .
Figure 13.Distributions of thermally induced and tangential traction caused interfacial shear stresses for different α.(a) Direction of thermally induced interfacial shear stress is opposite to that caused by tangential traction when α < 1; (b) Direction of thermally induced interfacial shear stress is same as that caused by tangential traction when α > 1.
indicate that the dependence of interfacial tensile transverse stress on the variations of ζ is negligible.
indicate that the dependence of interfacial tensile transverse stress on the variations of ζ is negligible.

Figure 16 .Figure 17 .
Figure 16.Effect of ζ on the contact pressure and interfacial shear stress component, caused by friction force.(a) Increase in ζ results in decreases in contact pressure; (b) Increase in ζ results in decreases in interfacial shear stress component, caused by friction force.
D 1 112, D 1 1112, and D 1 113 into Equation ( • Obtain contact pressure p and heat flux H. Calculate σ AA using Equation (7) by employing G AA p , G AA v , and G AA t .
The following non-dimensional temperature rise, stress components, and coordinates are used: