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

Abstract

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.

1. 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₃N₄ 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].

2. Theory

2.1. 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.

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⁻²; then, the friction force distribution, q(x,y) N·m⁻², and the distribution of frictional heat flux, H(x,y) W·m⁻², are generated within the contact region, with f being the friction coefficient. Their relationships are as follows:

q(x,y) = f_p(x,y), H(x,y) = f_vp(x,y)

(1)

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⁻³ 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:

$$\begin{array}{l}p(i,j)>0,g(i,j)=0,(i,j)\in {\mathrm{\Omega }}_c\\ p(i,j)=0,g(i,j)>0,(i,j)\notin {\mathrm{\Omega }}_c\\ {\mathrm{\Delta }}_x{\mathrm{\Delta }}_y{\displaystyle \sum _{j=0}^{N_y-1}{\displaystyle \sum _{i=0}^{N_x-1}{p}_{(i,j)}}}=F_z\\ g_{(i,j)}={g^0}_{(i,j)}+{u^{pz}}_{(i,j)}+{u^{tz}}_{(i,j)}-\delta \end{array}$$

(2)

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⁰(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.

The initial gap g⁰(i,j) can be described as the distance between the central point Q and the outer surface of the ball in the z direction:

$$\begin{array}{l}{x}_{i,j}^2/d^2+y_{i,j}^2/d^2+{(z_{i,j}-d)}^2/d^2=1\\ z_{i,j}=g^0(i,j)\end{array}$$

(3)

Due to z_i,j << d within the contact zone, the relationship between g⁰(i,j) and x_i,j, y_i,j can be described as

$$g^0(i,j)=x_{i,j}^2/2d+y_{i,j}^2/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).

$${u^z}_{(i,j)}={u^{pz}}_{(i,j)}+{u^{tz}}_{(i,j)}={\mathrm{\Delta }}_x{\mathrm{\Delta }}_y{\displaystyle \sum _{k=0}^{N_y-1}{\displaystyle \sum _{l=0}^{N_x-1}{G}_{z-(i-l,j-k)}^p{p}_{(l,k)}}}+{\mathrm{\Delta }}_x{\mathrm{\Delta }}_y{\displaystyle \sum _{k=0}^{N_y-1}{\displaystyle \sum _{l=0}^{N_x-1}{G}_{z-(i-l,j-k)}^t{H}_{(l,k)}}}$$

(5)

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

$$V^{*}=\frac{1}{2}{\mathrm{\Delta }}_x{\mathrm{\Delta }}_y{\displaystyle \sum _{i=0}^{N_y-1}{\displaystyle \sum _{j=0}^{N_x-1}{p}_{(l,k)}}}{u^z}_{(i,j)}+{\mathrm{\Delta }}_x{\mathrm{\Delta }}_y{\displaystyle \sum _{i=0}^{N_y-1}{\displaystyle \sum _{j=0}^{N_x-1}{p}_{(l,k)}}}({g^0}_{(i,j)}-\delta )$$

(6)

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:

$${\sigma }_{AA-(x_0,y_0,z_0)}={\mathrm{\Delta }}_x{\mathrm{\Delta }}_y{\displaystyle \sum _{k=0}^{N_y-1}{\displaystyle \sum _{l=0}^{N_x-1}{G}_{AA}^p{p}_{(l,k)}}}+{\mathrm{\Delta }}_x{\mathrm{\Delta }}_y{\displaystyle \sum _{k=0}^{N_y-1}{\displaystyle \sum _{l=0}^{N_x-1}{G}_{AA}^v{q}_{(l,k)}}}+{\mathrm{\Delta }}_x{\mathrm{\Delta }}_y{\displaystyle \sum _{k=0}^{N_y-1}{\displaystyle \sum _{l=0}^{N_x-1}{G}_{AA}^t{H}_{(l,k)}}}$$

(7)

where σ_AA is an arbitrary stress component of point (x₀,y₀,z₀) 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₀,y₀,z₀), 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.

2.2. 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Ⅱ, 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₁ 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₁. 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₁, 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 those in the substrate space are marked as C_k. Local coordinate systems are established and originate from each image point.

The relationship between the global coordinate and the local coordinates can be written as:

$$\{\begin{array}{l}{R}_k=\sqrt{r^2+z_k^2},z_k=z-(2k-1)h\\ {\mathsf{\rho }}_k=\sqrt{r^2+s_k^2},s_k=z+(2k-1)h\end{array}k=1,2,3....\infty .$$

(8)

Love’s stress and displacement functions can be written as:

$$\{\begin{array}{l}2\mathsf{\mu }{u}_r=\frac{\partial \varphi }{\partial r}+\frac{\partial \mathrm{\Psi }}{\partial r},2\mu u_z=\frac{\partial \varphi }{\partial z}-(3-4\nu )\mathrm{\Psi }+z\frac{\partial \mathrm{\Psi }}{\partial z}\\ {\sigma }_r=-\frac{1}{r}\frac{\partial \varphi }{\partial r}-\frac{{\partial }^2\varphi }{\partial z^2}-\frac{z}{r}\frac{\partial \mathrm{\Psi }}{\partial r}-z\frac{{\partial }^2\mathrm{\Psi }}{\partial r^2}-2\nu \frac{\partial \mathrm{\Psi }}{\partial z}\\ {\sigma }_{\theta }=\frac{1}{r}\frac{\partial \varphi }{\partial r}+\frac{z}{r}\frac{\partial \mathrm{\Psi }}{\partial r}-2\nu \frac{\partial \mathrm{\Psi }}{\partial z}\\ {\sigma }_z=\frac{{\partial }^2\varphi }{\partial z^2}-2(1-\nu )\frac{\partial \mathrm{\Psi }}{\partial z}+z\frac{{\partial }^2\mathrm{\Psi }}{\partial z^2}\\ {\tau }_{rz}=\frac{{\partial }^2\varphi }{\partial r\partial z}-(1-2\nu )\frac{\partial \mathrm{\Psi }}{\partial r}+z\frac{{\partial }^2\mathrm{\Psi }}{\partial r\partial z}\end{array}$$

(9)

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 φ¹ and Ψ¹. For coating materials, the harmonic functions are noted as φ² and Ψ². Due to the singular point O₁ 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

$$\{\begin{array}{l}{\varphi }^1={\displaystyle \sum _{i=1}^{\infty }{\varphi }_i^1(r,z_i)+{\displaystyle \sum _{i=1}^{\infty }{\delta }_i^1(r,s_i),{\varphi }^2={\displaystyle \sum _{i=1}^{\infty }{\varphi }_i^2(r,z_i)}}}\\ {\mathrm{\Psi }}^1={\displaystyle \sum _{i=1}^{\infty }{\mathrm{\Omega }}_i^1(r,z_i)+{\displaystyle \sum _{i=1}^{\infty }{\eta }_i^1(r,s_i)}},{\mathrm{\Psi }}^2={\displaystyle \sum _{i=1}^{\infty }{\mathrm{\Omega }}_i^2(r,z_i)}\end{array}$$

(10)

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.

$$u_r^1=u_r^2,u_z^1=u_z^2,{\sigma }_z^1={\sigma }_z^2,{\tau }_{rz}^1={\tau }_{rz}^2,z=0$$

(11)

Another boundary condition can be established by satisfying the free surface.

$${\sigma }_z^1=0,{\tau }_{rz}^1=0,z=h$$

(12)

Substituting Equation (10) into Equation (9) using boundary condition Equation (11) and geometry condition z = 0, z_i = s_i = −(2i − 1)h, R_i = ρ_i on the interface; recursive relations can be obtained:

$$\begin{array}{l}\{\begin{array}{l}{\eta }_i^1(r,s_i)=\frac{\beta -\alpha }{1-\beta }[2\frac{\partial {\varphi }_i^1(r,z_i)}{\partial z}-{\kappa }_1{\mathrm{\Omega }}_i^1(r,z_i)]\\ {\delta }_i^1(r,s_i)=-\frac{\beta -\alpha }{1-\beta }{\varphi }_i^1(r,z_i)-\frac{1-\alpha }{2\mathrm{\Gamma }}[\frac{{\kappa }_1}{1-\beta }-\frac{{\kappa }_2}{1+\beta }]{\displaystyle \int {\mathrm{\Omega }}_i^1(r,z_i)\mathrm{d}z}\end{array}\\ z_i=-s_i,R_i={\rho }_i\\ \{\begin{array}{l}{\varphi }_i^2(r,z_i)=\frac{1-\alpha }{1-\beta }{\varphi }_i^1(r,z_i)-\frac{1-\alpha }{2}[\frac{{\kappa }_1}{1-\beta }-\frac{{\kappa }_2}{1+\beta }]{\displaystyle \int {\mathrm{\Omega }}_i^1(r,z_i)\mathrm{d}z}\\ {\mathrm{\Omega }}_i^2(r,z_i)=\frac{1-\alpha }{1+\beta }{\mathrm{\Omega }}_i^1(r,z_i)\end{array}\end{array}.$$

(13)

Substituting Equation (10) into Equation (9) using boundary condition Equation (12) and geometry condition z = h, z_i = –(2i − 1)h, s_i = 2ih, R_i+1 = ρ_i on the surface, another recursive relation can be obtained:

$$\begin{array}{l}\{\begin{array}{l}{\mathrm{\Omega }}_{i+1}^1(r,z_{i+1})=-2\frac{\partial {\delta }_{ii-1}^1(r,s_i)}{\partial z}+(3-4{\nu }_1){\eta }_i^1(r,s_i)-2h\frac{\partial {\eta }_i^1(r,s_i)}{\partial z}\\ {\varphi }_{i+1}^1(r,z_{i+1})=(3-4{\nu }_1){\delta }_{ii-1}^1(r,s_i)-4(1-2{\nu }_1)(1-{\nu }_1){\displaystyle \int {\eta }_i^1(r,s_i)}dz+h[2\frac{\partial {\delta }_{ii-1}^1(r,s_i)}{\partial z}+2h\frac{\partial {\eta }_i^1(r,s_i)}{\partial z}]\end{array}\\ s_i=-z_{i+1},{\rho }_i=R_{i+1}\end{array}$$

(14)

where κ is the Kappa parameter; κ₁ = 3 − 4v₁ and κ₂ = 3 − 4v₂ are Kappa parameters for the coating material and substrate material. α and β are Dundurs parameters [24] for the coating–substrate system.

Thus, φ¹₁(r,z₁) and Ω¹₁(r,z₁) 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.

$$\begin{array}{l}\{\begin{array}{l}{\varphi }_1^1(r,z_1)=-(1-2v_1)A_n\ln (R_1-z_1)-\frac{h{A}_n}{R_1}\\ {\mathrm{\Omega }}_1^1(r,z_1)=\frac{A_n}{R_1}\end{array}\\ A_n=-\frac{p}{2\mathsf{\pi }}\end{array}$$

(15)

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.

2.3. 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₁ 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₁. 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:

$$\begin{array}{l}\{\begin{array}{l}{u}_x=\frac{3-4v}{4(1-v)}{B}_1-\frac{1}{4(1-v)}[x\frac{\partial B_1}{\partial x}+y\frac{\partial B_2}{\partial x}+z\frac{\partial B_3}{\partial x}+\frac{\partial \delta }{\partial x}]\\ u_y=\frac{3-4v}{4(1-v)}{B}_2-\frac{1}{4(1-v)}[x\frac{\partial B_1}{\partial y}+y\frac{\partial B_2}{\partial y}+z\frac{\partial B_3}{\partial y}+\frac{\partial \delta }{\partial y}]\\ u_z=\frac{3-4v}{4(1-v)}{B}_3-\frac{1}{4(1-v)}[x\frac{\partial B_1}{\partial z}+y\frac{\partial B_2}{\partial z}+z\frac{\partial B_3}{\partial z}+\frac{\partial \delta }{\partial z}]\end{array}\\ \{\begin{array}{l}{\sigma }_x=\frac{v\mu }{1-v}(\frac{\partial B_1}{\partial x}+\frac{\partial B_2}{\partial y}+\frac{\partial B_3}{\partial z})+\frac{\mu (1-2v)}{1-v}\frac{\partial B_1}{\partial x}-\frac{\mu }{2(1-v)}(x\frac{{\partial }^2{B}_1}{\partial x^2}+y\frac{{\partial }^2{B}_2}{\partial x^2}+z\frac{{\partial }^2{B}_3}{\partial x^2}+\frac{{\partial }^2\delta }{\partial x^2})\\ {\sigma }_y=\frac{v\mu }{1-v}(\frac{\partial B_1}{\partial x}+\frac{\partial B_2}{\partial y}+\frac{\partial B_3}{\partial z})+\frac{\mu (1-2v)}{1-v}\frac{\partial B_2}{\partial y}-\frac{\mu }{2(1-v)}(x\frac{{\partial }^2{B}_1}{\partial y^2}+y\frac{{\partial }^2{B}_2}{\partial y^2}+z\frac{{\partial }^2{B}_3}{\partial y^2}+\frac{{\partial }^2\delta }{\partial y^2})\\ {\sigma }_z=\frac{v\mu }{1-v}(\frac{\partial B_1}{\partial x}+\frac{\partial B_2}{\partial y}+\frac{\partial B_3}{\partial z})+\frac{\mu (1-2v)}{1-v}\frac{\partial B_3}{\partial z}-\frac{\mu }{2(1-v)}(x\frac{{\partial }^2{B}_1}{\partial z^2}+y\frac{{\partial }^2{B}_2}{\partial z^2}+z\frac{{\partial }^2{B}_3}{\partial z^2}+\frac{{\partial }^2\delta }{\partial z^2})\\ {\tau }_{xy}=\frac{(1-2v)\mu }{2(1-v)}(\frac{\partial B_1}{\partial y}+\frac{\partial B_2}{\partial x})-\frac{\mu }{2(1-v)}(x\frac{{\partial }^2{B}_1}{\partial x\partial y}+y\frac{{\partial }^2{B}_2}{\partial x\partial y}+z\frac{{\partial }^2{B}_3}{\partial x\partial y}+\frac{{\partial }^2\delta }{\partial x\partial y})\\ {\tau }_{yz}=\frac{(1-2v)\mu }{2(1-v)}(\frac{\partial B_2}{\partial z}+\frac{\partial B_3}{\partial y})-\frac{\mu }{2(1-v)}(x\frac{{\partial }^2{B}_1}{\partial z\partial y}+y\frac{{\partial }^2{B}_2}{\partial z\partial y}+z\frac{{\partial }^2{B}_3}{\partial z\partial y}+\frac{{\partial }^2\delta }{\partial z\partial y})\\ {\tau }_{xz}=\frac{(1-2v)\mu }{2(1-v)}(\frac{\partial B_3}{\partial x}+\frac{\partial B_1}{\partial z})-\frac{\mu }{2(1-v)}(x\frac{{\partial }^2{B}_1}{\partial x\partial z}+y\frac{{\partial }^2{B}_2}{\partial x\partial z}+z\frac{{\partial }^2{B}_3}{\partial x\partial z}+\frac{{\partial }^2\delta }{\partial x\partial z})\end{array}\end{array}.$$

(16)

One of the recursive relations can be expressed as:

$$\begin{gathered} \{\begin{array}{l}{\varphi }_{\rho i}^1=\frac{1-\mathrm{\Gamma }}{1+\mathrm{\Gamma }}{\varphi }_{Ri}^1\\ {\mathrm{\Omega }}_{\rho i}^1=\frac{{\kappa }_1(1-\mathrm{\Gamma })}{{\kappa }_1\mathrm{\Gamma }+1}{\mathrm{\Omega }}_{Ri}^1+\frac{(1-{\kappa }_1)(1-\mathrm{\Gamma })}{({\kappa }_1\mathrm{\Gamma }+1)(\mathrm{\Gamma }+1)}{\displaystyle \int \frac{\partial {\varphi }_{Ri}^1}{\partial x}dz}-\frac{2(1-\mathrm{\Gamma })}{{\kappa }_1\mathrm{\Gamma }+1}[x\frac{\partial {\varphi }_{Ri}^1}{\partial z}+\frac{\partial {\beta }_{Ri}^1}{\partial z}-(z_i+(2i-1)h)\frac{\partial {\varphi }_{Ri}^1}{\partial x}]\\ {\beta }_{\rho i}^1=-\frac{{\kappa }_1+1}{2}(\frac{{\kappa }_1}{{\kappa }_1\mathrm{\Gamma }+1}-\frac{{\kappa }_2}{\mathrm{\Gamma }+{\kappa }_2}){\displaystyle \int {\mathrm{\Omega }}_{Ri}^{1}dz}+\frac{{\kappa }_1(1-\mathrm{\Gamma })}{{\kappa }_1\mathrm{\Gamma }+1}{\beta }_{Ri}^1+\frac{({\kappa }_1-1)(1-\mathrm{\Gamma })}{({\kappa }_1\mathrm{\Gamma }+1)(1+\mathrm{\Gamma })}[x{\varphi }_{Ri}^1-(z_i+(2i-1)h){\displaystyle \int \frac{\partial {\varphi }_{Ri}^1}{\partial x}dz}]+\\ \frac{1}{2(1+\mathrm{\Gamma })}(\frac{{{\kappa }_1}^2-1+2\mathrm{\Gamma }+2{\kappa }_1}{{\kappa }_1\mathrm{\Gamma }+1}+1-{\kappa }_1-2m+\frac{(1-\mathrm{\Gamma })m(1-{\kappa }_2){\kappa }_2}{\mathrm{\Gamma }+{\kappa }_2}){\displaystyle \int {\displaystyle \int \frac{\partial {\varphi }_{Ri}^1}{\partial x}dz}}dz\end{array} \\ \{\begin{array}{l}{\varphi }_{Ri}^2=\frac{2}{1+\mathrm{\Gamma }}{\varphi }_{Ri}^1\\ {\mathrm{\Omega }}_{Ri}^2=\frac{{\kappa }_2+1}{\mathrm{\Gamma }+{\kappa }_2}{\mathrm{\Omega }}_{Ri}^1+\frac{(1-{\kappa }_2)(1-\mathrm{\Gamma })}{({\kappa }_2\mathrm{\Gamma }+1)(\mathrm{\Gamma }+1)}{\displaystyle \int \frac{\partial {\varphi }_{Ri}^1}{\partial x}dz}\\ {\beta }_{Ri}^2=-\frac{{\kappa }_2+1}{2}(\frac{{\kappa }_1}{{\kappa }_1\mathrm{\Gamma }+1}-\frac{{\kappa }_2}{\mathrm{\Gamma }+{\kappa }_2}){\displaystyle \int {\mathrm{\Omega }}_{Ri}^{1}dz}+\frac{{\kappa }_2+1}{{\kappa }_1\mathrm{\Gamma }+1}{\beta }_{Ri}^1+(\frac{({\kappa }_2+1)}{({\kappa }_1\mathrm{\Gamma }+1)}-\frac{2}{1+\mathrm{\Gamma }})(x{\varphi }_{Ri}^1-z{\displaystyle \int \frac{\partial {\varphi }_{Ri}^1}{\partial x}dz})+\\ \frac{1}{2m(1+\mathrm{\Gamma })}(\frac{{{\kappa }_1}^2-1+2\mathrm{\Gamma }+2{\kappa }_1}{{\kappa }_1\mathrm{\Gamma }+1}+1-{\kappa }_1-2m+\frac{(1-\mathrm{\Gamma })m(1-{\kappa }_2){\kappa }_2}{\mathrm{\Gamma }+{\kappa }_2}){\displaystyle \int {\displaystyle \int \frac{\partial {\varphi }_{Ri}^1}{\partial x}dz}}dz\end{array} \end{gathered}$$

(17)

Another recursive relation can be expressed as:

$$\{\begin{array}{l}{\varphi }_{Ri+1}^1={\varphi }_{\rho i}^1\\ {\mathrm{\Omega }}_{Ri+1}^1=-2(1-2v_1){\displaystyle \int \frac{\partial {\varphi }_{\rho i}^1}{\partial x}dz}+{\kappa }_1{\mathrm{\Omega }}_{\rho i}^1-2(x\frac{\partial {\varphi }_{\rho i}^1}{\partial z}-(s_i-2ih)\frac{\partial {\varphi }_{\rho i}^1}{\partial x}+h\frac{\partial {\mathrm{\Omega }}_{\rho i}^1}{\partial z}+\frac{\partial {\beta }_{\rho i}^1}{\partial z})\\ {\beta }_{Ri+1}^1=\frac{{{\kappa }_1}^2-1}{2}({\displaystyle \int {\displaystyle \int \frac{\partial {\varphi }_{\rho i}^1}{\partial x}dz}}dz-{\displaystyle \int {\mathrm{\Omega }}_{\rho i}^{1}dz})+{\kappa }_1{\beta }_{\rho i}+({\kappa }_1-1)[x{\varphi }_{\rho i}^1-(s_i-(2i+1)h){\displaystyle \int \frac{\partial {\varphi }_{\rho i}^1}{\partial x}dz})+\\ 2h(x\frac{\partial {\varphi }_{\rho i}^1}{\partial z}-(s_i-2ih)\frac{\partial {\varphi }_{\rho i}^1}{\partial x}+h\frac{\partial {\mathrm{\Omega }}_{\rho i}^1}{\partial z}+\frac{\partial {\beta }_{\rho i}^1}{\partial z})\end{array}$$

(18)

where μ₂/μ₁ is noted as Γ and (1 − v₁)/(1 − v₂) is noted as m.

The Midlin solutions for a point tangential load, acting on the surface of a semi-infinite space, can be written as:

$$\begin{array}{l}\{\begin{array}{l}{\varphi }_{R1}^1(x,y,z_1)=\frac{A_f}{R_1}\\ {\mathrm{\Omega }}_{R1}^1(x,y,z_1)=-\frac{(1-2v_1)A_{f}x}{R_1(R_1-z_1)}\\ {\beta }_{R1}^1(x,y,z_1)=-\frac{(1-2v_1)A_{f}x}{R_1-z_1}(1-2v_1-\frac{h}{R_1})\end{array}\\ A_f=-\frac{p_x}{2\pi {\mu }_1}\end{array}$$

(19)

where p_x is a tangential load. Now, B₁, B₂, B₃ and δ can be obtained by using Equation (17)–(19). Green’s functions can be obtained by substituting B₁, B₂, B₃, and δ into Equation (16) and then separating variable p_x.

2.4. 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₁ 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₁. 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.

The general solution for displacement and stress within a semi-infinite space for a thermal load can be written as:

$$\begin{array}{l}\{\begin{array}{l}2\mu u_r=\frac{\partial {\psi }_1}{\partial r}+z\frac{\partial {\psi }_2}{\partial r}\\ 2\mu u_z=\frac{\partial {\psi }_1}{\partial z}-(3-4v){\psi }_2+z\frac{\partial {\psi }_2}{\partial r}+4(1-v){\psi }_3\\ \frac{2\mu t}{E_t}=\frac{\partial {\psi }_3}{\partial z}\end{array}\\ \{\begin{array}{l}{\sigma }_r=\frac{{\partial }^2{\psi }_1}{\partial r^2}-2v\frac{\partial {\psi }_2}{\partial z}+z\frac{{\partial }^2{\psi }_2}{\partial r^2}-2(1-v)\frac{\partial {\psi }_3}{\partial z}\\ {\sigma }_{\theta }=\frac{1}{r}\frac{\partial {\psi }_1}{\partial r}-2v\frac{\partial {\psi }_2}{\partial z}+\frac{r}{z}\frac{\partial {\psi }_2}{\partial r}-2(1-v)\frac{\partial {\psi }_3}{\partial z}\\ {\sigma }_z=\frac{{\partial }^2{\psi }_1}{\partial z^2}-2(1-v)\frac{\partial {\psi }_2}{\partial z}+z\frac{{\partial }^2{\psi }_2}{\partial r^2}+2(1-v)\frac{\partial {\psi }_3}{\partial z}\\ {\tau }_{rz}=\frac{{\partial }^2{\psi }_1}{\partial z\partial r}-(1-2v)\frac{\partial {\psi }_2}{\partial r}+z\frac{{\partial }^2{\psi }_2}{\partial r\partial z}+2(1-v)\frac{\partial {\psi }_3}{\partial r}\end{array}\end{array}$$

(20)

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:

$$\begin{array}{l}\{\begin{array}{l}{\psi }_j^1={\displaystyle \sum _{i=1}^{\infty }({\psi }_{ij}^1+{\psi }_{ij}^2),(j=1,2,3)}\\ {\psi }_{ij}^1={\displaystyle \sum _{k=1}^i{\psi }_{ikj}^1,{\psi }_{ij}^2={\displaystyle \sum _{k=1}^i{\psi }_{ikj}^2,(j=1,2,3)}}\end{array}\\ \{\begin{array}{l}{\psi }_{i11}^1=D_{i111}^1(z_t\ln R_i^{*}+R_i)+D_{i112}^1\ln R_i^{*}\\ {\psi }_{i1j}^1=D_{i1j}^1\ln R_i^{*},(j=2,3)\\ {\psi }_{ikj}^1=D_{ijk1}^1\frac{{\partial }^{2(k-2)}}{\partial z^{2(k-2)}}(\frac{1}{R_i})+D_{ijk2}^1\frac{{\partial }^{(2k-3)}}{\partial z^{(2k-3)}}(\frac{1}{R_i}),(k\ge 2,j=1,2,3)\\ R_i^{*}=R_i-z_i\\ {\psi }_{i11}^2=D_{i111}^2(s_t\ln {\rho }_i^{*}+{\rho }_i)+D_{i112}^2\ln {\rho }_i^{*}\\ {\psi }_{i1j}^2=D_{i1j1}^2\ln {\rho }_i^{*}+D_{i1j2}^2\frac{1}{{\rho }_i},(j=2,3)\\ {\psi }_{ik1}^2=D_{ik11}^2\frac{{\partial }^{2(k-2)}}{\partial z^{2(k-2)}}(\frac{1}{{\rho }_i})+D_{ik12}^2\frac{{\partial }^{(2k-3)}}{\partial z^{(2k-3)}}(\frac{1}{{\rho }_i}),(k\ge 2,j=1)\\ {\psi }_{ikj}^2=D_{ikj1}^2\frac{{\partial }^{(2k-3)}}{\partial z^{(2k-3)}}(\frac{1}{{\rho }_i})+D_{ikj2}^2\frac{{\partial }^{2(k-1)}}{\partial z^{2(k-1)}}(\frac{1}{{\rho }_i}),(k\ge 2,j=2,3)\end{array}\\ {\rho }_i^{*}={\rho }_i-s_i\end{array}$$

(21)

Harmonic functions in the substrate material can be written as

$$\begin{array}{l}{\mathrm{\Psi }}_j^1={\displaystyle \sum _{i=1}^{\infty }{\mathrm{\Psi }}_{ij}^1}=={\displaystyle \sum _{i=1}^{\infty }{\displaystyle \sum _{k=1}^i{\mathrm{\Psi }}_{ikj}^1}},(j=1,2,3)\\ \{\begin{array}{l}{\mathrm{\Psi }}_{i11}^1=L_{i111}^1(z_t\ln R_i^{*}+R_i)+L_{i112}^1\ln R_i^{*}\\ {\mathrm{\Psi }}_{i1j}^1=L_{i1j1}^1\ln R_i^{*}+L_{i1j2}^1\frac{1}{R_i},(j=2,3)\\ {\mathrm{\Psi }}_{ik1}^1=L_{ik11}^1\frac{{\partial }^{2(k-2)}}{\partial z^{2(k-2)}}(\frac{1}{R_i})+L_{ik12}^1\frac{{\partial }^{(2k-3)}}{\partial z^{(2k-3)}}(\frac{1}{R_i}),(k\ge 2,j=1)\\ {\mathrm{\Psi }}_{ikj}^1=L_{ikj1}^1\frac{{\partial }^{(2k-3)}}{\partial z^{(2k-3)}}(\frac{1}{R_i})+L_{ikj2}^1\frac{{\partial }^{2(k-1)}}{\partial z^{2(k-1)}}(\frac{1}{R_i}),(k\ge 2,j=2,3)\end{array}\end{array}$$

(22)

One of the boundary conditions can be established by satisfying the continuous conditions of displacement, stress, temperature, and thermal flux on the interface.

$$u_r^1=u_r^2,u_z^1=u_z^2,t^1=t^2,{\sigma }_z^1={\sigma }_z^2,{\tau }_{rz}^1={\tau }_{rz}^2,{\zeta }_1\frac{\partial t^1}{\partial z}={\zeta }_2\frac{\partial t^2}{\partial z}$$

(23)

Another boundary conditions can be established by satisfying the stress free and thermally insult conditions on the surface.

$${\sigma }_z^1=0,{\tau }_{rz}^1=0,{\zeta }_1\frac{\partial t^1}{\partial z}=0$$

(24)

Additionally, mechanical equilibrium and thermal equilibrium of an infinite plane, ε < z < h, should be satisfied.

$$\{\begin{array}{l}-2\pi {\zeta }_1{\displaystyle {\int }_0^{+\infty }\frac{\partial t^1}{\partial z}dz=\frac{qfP{e}_2{\lambda }_2}{a_h}}\\ -2\pi {\sigma }_z{\displaystyle {\int }_0^{+\infty }{\sigma }_{z}rdr=0}\\ P{e}_2=\frac{V{a}_h}{{\lambda }_2}\end{array}$$

(25)

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. λ₂ is the thermal diffusion coefficient for the substrate material.

Substituting Equation (21) into Equation (20), using boundary condition (24) and equilibrium condition (25), one can obtain:

$$\{\begin{array}{l}-D_{1111}^1+2(1-v)D_{112}^1-2(1-v)D_{113}^1=0\\ D_{1111}^1-2(1-v)D_{112}^1+2(1-v)D_{113}^1=0\\ D_{1112}^1+h{D}_{112}^1=0\\ D_{113}^1=-\frac{H{\mu }_1}{\pi {\zeta }_1{E}_{t1}}\end{array}$$

(26)

Thus, D¹₁₁₁₁, D¹₁₁₂, D¹₁₁₁₂, and D¹₁₁₃ can be obtained from Equation (26). Substituting Equation (21) into Equation (20) using boundary condition (23), one can obtain

$$\{\begin{array}{l}{D}_{1111}^1-D_{1111}^2=L_{1111}^1(\frac{{\mu }_1}{{\mu }_2})\\ D_{1112}^1+D_{1112}^2=L_{1112}^1(\frac{{\mu }_1}{{\mu }_2})\\ D_{1111}^1-(3-4v_1)D_{112}^1+4(1-v_1)D_{113}^1+D_{1111}^2-(3-4v_1)D_{1121}^2+4(1-v_1)D_{1131}^1=\\ [L_{1111}^1-(3-4v_2)L_{1121}^1+4(1-v_2)D_{1131}^1](\frac{{\mu }_1}{{\mu }_2})\\ -D_{1112}^1+D_{1112}^2-(3-4v_1)D_{1122}^2+4(1-v_1)D_{1132}^1=[-L_{1112}^1-(3-4v_2)L_{1122}^1+4(1-v_2)D_{1132}^1](\frac{{\mu }_1}{{\mu }_2})\\ -D_{1111}^1+2(1-v_1)D_{112}^1-2(1-v_1)D_{113}^1+D_{1111}^2-2(1-v_1)D_{1121}^2+2(1-v_1)D_{1131}^2=\\ -L_{1111}^1+2(1-v_2)L_{1121}^1-2(1-v_2)L_{1131}^1\\ -D_{1112}^1-D_{1112}^2+2(1-v_1)D_{1122}^2-2(1-v_1)D_{1132}^2=-L_{1112}^1-2(1-v_2)L_{1122}^1+2(1-v_2)L_{1132}^1\\ D_{1111}^1-(1-2v_1)D_{112}^1+2(1-v_1)D_{113}^1+D_{1111}^2-(1-2v_1)D_{1121}^2+2(1-v_1)D_{1131}^2=\\ L_{1111}^1-(1-2v_2)L_{1121}^1+2(1-v_2)L_{1131}^1\\ D_{1112}^1-D_{1112}^2+(1-2v_1)D_{1122}^2-2(1-v_1)D_{1132}^2=L_{1112}^1+2(1-v_2)L_{1122}^1-2(1-v_2)L_{1132}^1\\ -D_{113}^1+D_{1131}^2=-L_{1131}^1(\frac{{\mu }_1{E}_{t2}}{{\mu }_2{E}_{t1}})\\ D_{1132}^2=L_{1132}^1(\frac{{\mu }_1{E}_{t2}}{{\mu }_2{E}_{t1}})\\ -h{D}_{113}^1-h{D}_{1131}^2-D_{1132}^2=(h{L}_{1131}^1-L_{1132}^1)(\frac{{\mu }_1{E}_{t2}{\zeta }_2}{{\mu }_2{E}_{t1}{\zeta }_1})\\ D_{1132}^2=L_{1132}^1(\frac{{\mu }_1{E}_{t2}{\zeta }_2}{{\mu }_2{E}_{t1}{\zeta }_1})\end{array}.$$

(27)

Substituting D¹₁₁₁₁, D¹₁₁₂, D¹₁₁₁₂, and D¹₁₁₃ into Equation (27), twelve coefficients D²_11jl and L¹_11jl (j = 1, 2, 3; l = 1, 2) can be obtained. Using boundary conditions Equations (23) and (24), successively, equations such as Equations (26) and (27) can be established and each coefficient of the harmonic functions can be obtained.

Green’s functions for the thermal load can be obtained by substituting ψ₁, ψ₂, ψ₃ and Ψ₁, Ψ₂, Ψ₃ into Equation (20) and separating variable H.

2.5. 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 δ⁰, and initial gap g⁰.
• 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⁰, 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₀, 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 − p₀| to judge convergence. If ε^iter < ε₀, 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 G_AA^p, G_AA^v, and G_AA^t.

3. 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:

$$P{e}_2=V{a}_h/{\lambda }_s,p_n=p/c_h,t_n=t{\zeta }_s/f{c}_{h}V,{\sigma }_{n-AA}={\sigma }_{AA}/c_h,x_n=x/a_h,y_n=y/a_h,z_n=z/a_h$$

(28)

where Pe₂ 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⁻⁵ K⁻¹, ζ_s = 50.2 W·m⁻¹·K⁻¹, and λ_s = 1.0 × 10⁻⁵ m²·s⁻¹. 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⁻¹; thus, Pe₂ is 2.4.

3.1. 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.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.

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.

3.2. 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/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/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.

3.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/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.

3.4. 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 indicate that the dependence of interfacial tensile transverse stress on the variations of ζ is negligible.

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⁻¹. 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 load bearing capacity of the coating–substrate system due to its reflecting the real thermal load condition and contact condition well.

4. 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.