Stress Analysis of Multilayered Coatings Subjected to Surface Point Contact Loading Based on Its Three-Dimensional Elastic Field Solution

: In order to investigate the effect of the structural layout of multilayered coatings on its mechanical behavior, a three-dimensional elastic field solution is developed for multilayered solids subjected to surface point contact loading, which is converted from the elastic field solution in frequency domain by using a numerical conversion algorithm. The elastic field solution in frequency domain is obtained by numerically solving a group of linear equations involving the unknown constants in the general elastic field solution of layered material that is obtained by using Fourier integral transform technique. The present solution is validated by comparing with the exact analytical solution for uncoated solids and finite element solution for solids coated with 30 layers. Lastly, the effect of structural layout of multilayered coatings is further investigated with present solution. The result shows that the gradient structural layout with elasticity modulus decreasing gradually from the top layer to the substrate, which is preferable to a larger friction coefficient for multilayered solids subjected to surface line contact loading, is preferable for a smaller friction coefficient <0.1 for multilayered solids subjected to surface point contact loading, and the gradient structural layout with elasticity modulus increasing first in the top layers and then decreasing in the bottom layers, which is preferable to a smaller friction coefficient for multilayered solids subjected to surface line contact loading, is preferable for a friction coefficient >0.2.


Introduction
Advanced coating structures, such as multilayered coatings, gradient coatings etc., are employed to improve the tribological performance and service life of tribo-parts because of its advantage of anti-friction, anti-wear, anti-scuffing, and anti-corrosion, especially for those performing under severe conditions [1][2][3]. No doubt, an insight into the mechanical behavior of a multilayered coatings subjected to surface contact loading is essential to its optimal design and application [4,5].
Theoretical analysis of the mechanics of layered solids, linearly elastic theory used various integral transform techniques to produce the elastic field solutions for both two-dimensional (2D) and three-dimensional (3D) problems [6], since the independent variables in the elastic field governing differential equations can be reduced by using integral transform techniques. For 2D problems, the Airy stress function is always used to produce the 2D elastic field solution in frequency domain by adopting Fourier integral transform technique for developing a 2D contact model [7][8][9]. For 3D problems, the Papkovich-Neuber potential functions are always used to produce the 3D elastic field solution in frequency domain of coated solids by using Fourier integral transform technique for developing various 3D contact model. Based on the 3D elastic field solution in frequency domain, O'Sullivian and King [10] developed a sliding contact model between a spherical indenter and a monolayered solid. Nogi and Kato [11] further developed a rough contact model and introduced fast Fourier transform (FFT) technique to speed up the calculation speed. The 3D elastic field solution in frequency domain were further employed to study the sliding contact or partial slip contact problem for solids coated with monolayer or bilayers [12][13][14][15]. Besides the integral transform techniques, the finite element method (FEM)-based commercial software [16][17][18] and image point method [19][20][21] were also used to explain the failure mechanism of coatings subjected to contact loading based on the result of stress analysis. In the past several decades, extensive meaningful researches have been carried out theoretically or numerically to understand the mechanical behavior of coated solid mediums subjected to surface contact loading by using different methods. However, most of them were focused on the problems of a monolayer or a bilayers coating-substrate system, since the solving process for multilayered coatings becomes more and more complicated with the increase of the layer number. An expression of the 3D elastic field solution in frequency domain was given in term of recursion by using Papkovich-Neuber potential functions [22], while the derivation process is very tedious, and the expression of the solution in frequency domain is so lengthy and complicated that it is very inconvenient for researchers to implement. The effect of the structural layout on the mechanical behavior of multilayered solids subjected to surface line contact loading has been investigated and the preferable structural layout is recommended for various friction coefficient [4,5]. However, the effect of structural layout on the mechanical behavior of multilayered solids subjected to surface point contact loading, which may be different from that of the multilayered solids subjected to surface line contact loading, still lacks a systematic research.
In this paper, a 3D elastic field solution has been developed for numerical investigation on the mechanical behavior of multilayered coatings with various structural layout, which is converted from the 3D elastic field solution in frequency domain with a 2D Inverse Fast Fourier Transform (IFFT) based conversion algorithm. The 3D elastic field solution in frequency domain for multilayered coatings subjected to surface point contact loading is obtained by numerically solving a group of linear equations involving the unknown constants in the general elastic field solution of layered material that is obtained by using Fourier integral transform technique, and the group of linear equations is established according to the boundary conditions and interface continuous conditions. The present solution was validated by comparisons on the maximum shear stress with other approaches. Lastly, the effect of the structural layout of multilayered coatings is further studied with the present solution, and we found that the preferable structural layout of multilayered solids subjected to surface point contact loading for various friction coefficient is different from that of multilayered solids subjected to surface line contact loading.

Problem Description
The problem that a multilayered coating is subjected to the action of surface point contact loading, as shown in Figure 1, is considered in this work. The normal and tangential traction applied on the surface are denoted with ( ) , p x y and ( ) , q x y respectively, the number of the layers is denoted with N , the elastic modulus and Poisson ratio are denoted with k E and k ν respectively for the k th layer, the layer thickness is denoted with k h for the k th layer, and the k z that starts from the top surface of the k th layer is the z axis for the k th layer. The substate and the layers are homogeneous and perfectively bounded to each other, and the strain of the multilayered solids due to the surface point contact loading is so small that the linear elastic mechanic theory is applicable to analyze current problem.
According to Hertz contact theory [23], the normal pressure, and tangential traction applied on the surface are: where H p and H a are the maximum contact pressure and contact radius of Hertz point contact respectively, and f μ is the friction coefficient.

General Solution of the Elastic Field Governing Equation
For 3D problems, the elastic field governing equations of a homogeneous layered material in the multilayered system are [24]: where 2 2 The elastic field governing equations become: where the symbol " i " is the imaginary unit, the symbol "~" represents Fourier integral transform, The six stress components ( )

General Solution of the Elastic Field Governing Equation
Solving Equation (9), the general solution of the intermediate variable F  is: where ( )    (11). Using the elimination method adopted in Ref. [25], the general solution of the displacement ( ) where ( )    (21) and (22) into Equation (10) or Equation (11) as follows:  is the shear modulus of the k th layer.
For the elastic substrate, tends to infinity. Therefore, there are only 6 +3 N unknown constants needed to be determined for a multilayered system with N layers.

Determination of the Undetermined Constants
The stress and displacement components ( ) are continuous on all interfaces of the multilayered system since the layers deposited on the elastic substrate are perfectively bonded to each other. Therefore, the 6 +3 N unknown constants of the multilayered system with N layers can be determined with the boundary conditions and interface continuous conditions listed as follows: Substituting the Equations (12)- (19) into Equations (27)-(35), we have: Substituting Equation (9) into Equations (40), (43), and (44), a system of linear equations with 2 +1 N equations involving unknown constants marked with symbol A can be established and given in a matrix form as follows: Similarly, substituting Equations (9)-(11) into Equations (38), (39), (41), (42), (45), and (46), a system of linear equations with 4 +2 N equations involving unknown constants marked with symbol B can be established and given in a matrix form as follows: The detail of the matrices with a complicated and tedious derivation process, the algorithm Gaussian full pivoting elimination method seen in Ref. [4,28] is adopted to solve the Equations (47) and (48) to keep the accuracy of its solution in this paper.

2D IFFT Based Numerical Conversion Algorithm
After the 3D elastic field solution in frequency domain is determined, the elastic field solution in space domain can be obtained with a 2D IFFT-based numerical conversion algorithm [10,29], which is proposed originally to produce the influence coefficient matrix of stresses and displacements, which is essential to establishing a contact model based on semi-analytical method. Take the calculation of the stress component σ rs at any depth z for example, the basic procedures of the numerical conversion algorithm based on 2D IFFT are: • The grid in frequency domain ( ) and λ should be a non-negative integer to the power 2.
• A two-dimensional matrix ˆr s σ is constructed by using the solution of σ rs  in frequency domain: and κ is an aliasing phenomenon control parameter.
• Apply 2D IFFT to the matrix rs

Validation of the Present Method
In order to validate the present solution, the maximum shear stress in subsurface of the present solution is compared with that obtained by other approaches. Figure 2 shows the comparison of the maximum shear stress τ in x-o-z plane between the result of present solution and that of the exact analytical solution for an uncoated half-space. For an uncoated half-space subjected to Hertz-type contact pressure, an exact analytical solution of its stress field can be found in Ref. [30]. Table 1 shows the input parameters of the present method to analyze an uncoated half-space subjected to Hertz-type contact pressure, in which the elastic parameters of the layer and the substrate are the same to those of the half-space. An excellent match of the maximum shear stress distribution between the solution of present method and the exact analytical solution can be observed as well as a relative error less than 0.05% of the peak value of the maximum shear stress.   In order to verify that the present solution has no limit on the layer number of the multilayered systems, a comparison of the maximum shear stress τ in x-o-z plane between the present solution and the finite element solution is also conducted for multilayered system with 30 layers. The input parameters of the present solution are shown in Table 2. The comparison of the maximum shear stress is shown in Figure 3. It can be seen that the maximum shear stress in x-o-z plane obtained with present solution shown in Figure 3b matches with that of the finite element solution shown in Figure 3a, and the relative error of the peak value of the maximum shear stress is not more than 0.14%.
The comparisons conducted above not only validate the present solution and its high accuracy, but also show that the present solution is applicable to the stress analysis of multilayered systems with various layer number.

Numerical Investigation and Discussion
For an insight to the effect of the structural layout of the multilayered system on its mechanical behavior, stress analysis on the multilayered systems with different elasticity modulus changing modes from the top layer to the substrate is conducted based on the present solution. The maximum contact pressure H p and the contact radius H a of the surface mechanical loading are the same to those listed in Tables 1 and 2.
The effect of the structural layout of the multilayered coatings on the maximum shear stress H τ / p in x-o-z plane under various friction coefficients is shown in Figure 4, and the location of the symbol "*" in Figure 4 is the position of the peak value of the maximum shear stress. The multilayered system with elasticity modulus changing mode 1 shown in Figure 4a is a homogeneous hard monolayer, since the elastic modulus of each layer is the same and differs from the substrate. Similar to the results of Ref. [10,11], an significant discontinuity and stress concentration in the region close to the substrate is observed in the contour plot, and the peak value of the maximum shear stress, which is located at the bottom of the hard coating, reaches up to about 0.381pH when the friction coefficient μf is 0, which is recognized as a result of the sharp decrease of the elasticity modulus from the hard coating to the substrate. When the friction coefficient μf increases to 0.5, the peak value of the maximum shear stress increases rapidly in the surface region because of the action of frictional traction, which leads to the location of the peak value of the maximum shear stress shifts to the surface from the bottom of the hard coating, and its value increases to 0.519pH. The elasticity modulus changing modes of 2 to 4 shown in Figure 4a are designed as alternatives to avoid the sudden change of elasticity modulus of mode 1 through a gradient change. The elasticity modulus of mode 2 decreases uniformly from top to bottom by 80 GPa. The elasticity modulus of mode 3 decreases rapidly first by 225 GPa in the top three layers, and then decreases slowly from the 3rd layer to the bottom layer by about 19 GPa. The elasticity modulus of mode 4 increases first from the top layer to the 5th layer by 100 GPa, and then decreases gradually from the 5th layer to the substrate by about 167 GPa. The significant discontinuity observed in the maximum shear stress distribution of mode 1 does not exist in the maximum shear stress distribution of mode 2 to 4, and slight discontinuity appears on the interface of multilayers with structural layout of mode 2 to 4. It can be seen from Figure 4b that the peak value of the maximum shear stress of mode 2 to 4 are 0.328, 0.322, and 0.359pH respectively for μf = 0, which are all smaller than that of the mode 1. However, when the friction coefficient increases to 0.5, the peak value of the maximum shear stress of mode 2 and 3 are 0.568 and 0.677pH respectively, which are all larger than that of mode 1. The peak value of the maximum shear stress of mode 4 is 0.518pH, which is slightly smaller than that of mode 1. Figure 5 is the peak value of the maximum shear stress H τ / p and its depth of multilayered coatings with different elasticity modulus changing modes under various friction coefficient. As shown in Figure 5a, the peak value of the maximum shear stress of mode 2 to 4 are all smaller than that of mode 1 when the friction coefficient μf ranges from 0 to 0.1. It implies that the maximum shear stress concentration on the interface of the homogeneous hard monolayer can be avoided by using a multilayered system with an elasticity modulus changing gradually, as shown in Figure 4. As a result, an apparent decrease of the maximum shear stress of mode 2 to 3 in the region near the substrate can be observed in Figure 4b, and the depth in which the peak value of the maximum shear stress of mode 2 and 3 located is closer to the surface than that of the mode 1 and mode 4, as shown in Figure  5b. Among the four changing modes of elasticity modulus, the peak value of the maximum shear stress of mode 3 is the smallest when the friction coefficient tends to 0 and the peak value of the maximum shear stress of mode 2 is the smallest when the friction coefficient tends to 0.1. That means the maximum shear stress can be reduced by using a multilayered structural layout with elasticity modulus decreasing gradually from the top layer to the substrate when the friction coefficient is not smaller than 0.1. For example, the peak value of the maximum shear stress of mode 2 and 3 are about 14% and 15.5% smaller than that of mode 1 respectively when the friction coefficient is 0. It has been commonly recognized in previous research [4,5,8,11] that the maximum shear stress or von Mises equivalent stress increases with the increases of friction coefficient and its location moves toward to the surface, since the maximum shear stress in the surface region increases due to the increase of the surface frictional traction. The increase of the maximum shear stress in the surface region, which is caused by the increase of tangential traction applied on the surface, also strongly depends on the structural layout of multilayered system. For μf = 0, the peak value of mode 2 and 3 is about 9.4% and 30.4% larger than that of mode 1 respectively, however the peak value of mode 4 is still smaller than that of mode 1. The shifting of the maximum shear stress to the surface during the increase of friction coefficient also strongly depends on the elasticity modulus changing mode of multilayered systems. It can be observed in Figure 5b that the friction coefficient value to make the maximum shear stress of mode 2 and 3 shift to surface is smaller than that of mode 1 and 4.   Figure 6 is the stress σxx along z axis of multilayered coatings with different elasticity modulus changing modes when the friction coefficient μf is 0. Actually, the stress σxx along z axis of each elasticity modulus changing mode is the same under various friction coefficient, since the stress σxx in the x-o-z plane caused by tangential traction applied on the surface is asymmetrical at the two sides of z axis. It can be observed that the stress σxx along z axis is markedly different for multilayered coatings with different elasticity modulus changing mode. For mode 1, which represents a homogeneous hard coating, the stress component σxx is tensile in the region close to the substrate and the maximum tensile stress σxx located at the bottom of the coating reaches up to about 0.48pH, which is closely related to the cracks initiated at the bottom of the coating. For mode 2 to mode 4, which represent various multilayered hard coatings with gradient structural layout, the stress σxx is tensile too in the region close to the substrate, however the maximum tensile stress σxx is much smaller than that of mode 1, and the maximum tensile stress σxx decreases successively in the order of mode 4, mode 2 to mode 3.  Figure 7 is the stress σxx along x axis of multilayered coatings with different elasticity modulus changing modes under various friction coefficients. The stress σxx along x axis is tensile at the tailedge of loading area under various friction coefficients, which is believed to be responsible to the cracks initiated on the surface. The tensile stress σxx at the tail-edge of the loading area increases with the increase of friction coefficient [8,10,11], while the maximum tensile stress σxx at the tail-edge not only depends on the friction coefficient but also strongly depends on the elasticity modulus changing mode of multilayered system. As indicated in Figure 7, the maximum tensile stress σxx at the tail-edge of loading area decreases in the order of mode 3, mode 2, mode 4 to mode 1 under various friction coefficients, and the maximum tensile stress σxx at the tail-edge of mode 4 is almost the same to that of mode 1. Under relatively small friction coefficient such as 0 and 0.1, the maximum tensile stress at the tail-edge of mode 2 and mode 3 as shown in Figure 7a,b are smaller than the maximum tensile stress σxx along z axis of mode 1 as shown in Figure 6, although they are slightly larger than the maximum tensile stress σxx at the tail-edge of mode 1 and mode 4. When the friction coefficient increases to a value larger than 0.3, the maximum tensile stress σxx at the tail-edge of mode 2 and mode 3 becomes remarkably larger than the maximum tensile stress σxx at the tail-edge of mode 1 and mode 4 as well as the maximum tensile stress σxx in the region close to the substrate of mode 1.
Based on the numerical results presented above, it is confirmed that the structural layout has a strong influence on the mechanical of multilayered coatings subjected to surface point contact loading. The concentration of the maximum shear stress and the tensile stress σxx in the region close to the substrate of homogenous hard coatings can be effectively eliminated by using hard coatings with gradient structural layout. The structural layout of hard coatings with elasticity modulus decreasing gradually from the top layer to the substrate such as mode 2 and 3 is better than that with elasticity increasing first in the top layers and then decreasing in the bottom layers such as mode 4 in terms of relieving the concentration of the maximum shear stress and the tensile stress σxx in the region close to the substrate. However, the latter is much better than the former in terms of depressing the increase of the maximum shear stress and the increase of the tensile stress σxx caused by the increase of frictional traction applied on the surface. By synthetically considering the influence of the structural layout of multilayered hard coatings on the tensile stress σxx at the tail-edge of loading area and in the region near substrate as well as the peak value of the maximum shear stress, the gradient structural layout with elasticity modulus decreasing gradually from the top layer to the substrate, which is preferable to a larger friction coefficient for multilayered solids subjected to surface line contact loading [4,5], is preferable for a smaller friction coefficient <0.1 for multilayered solids subjected to surface point contact loading, and the gradient structural layout with elasticity modulus increasing first in the top layers and then decreasing in the bottom layers, which is preferable to a smaller friction coefficient for multilayered solids subjected to surface line contact loading [4,5], is preferable for a friction coefficient >0.2 for multilayered solids subjected to surface point contact loading.

Conclusions
A 3D elastic field solution has been promoted and validated for multilayered solids subjected to the surface point contact loading in this paper. A numerical investigation on the mechanical behavior of multilayered solids with various structural layout has also been performed with the present method. The following conclusions are drawn as follows: • A 3D elastic field solution of multilayered solids which has no limit on the layer number of multilayered coatings is developed and validated. • The disadvantages of homogeneous hard coatings can be effectively overcome by adopting a proper gradient structural layout of multilayered coatings. • The preferable structural layout of multilayered coatings subjected to surface point contact loading for specific friction coefficient is different from that of multilayered coatings subjected to surface line contact loading. For multilayered solids subjected to surface point contact loading, the gradient structural layout with elasticity modulus decreasing gradually from the top layer to the substrate is preferable for a friction coefficient smaller than 0.1 and the gradient structural layout with elasticity modulus increasing first in the top layers and then decreasing in the bottom layers is preferable for a friction coefficient larger than 0.2. u  can be obtained as follows: According to the relationship between stresses and displacements, shown in Equations (14)

Appendix B
For clarity and concise of description, the matrices in Equation (47) are composed with sub-matrices as follows: The sub-matrices of matrices