Two-Dimensional Mechanical Behavior Analysis of Multilayered Solids Subjected to Surface Contact Loading Based on a Semi-Analytical Method

In this paper, a two-dimensional semi-analytical method is developed for the mechanical behavior analysis of multilayered solids subjected to surface contact loading, which is indispensable for realizing an optimized tribological performance from the mechanical behavior point of view. Firstly, the explicit analytical frequency response functions of the multilayered solid are derived in a recursive form by analytically solving a system of linear equations established according to the boundary conditions and the interface continuous conditions. Then, the two-dimensional elastic field solution in the subsurface of multilayered solids in the space domain is converted from its corresponding frequency response functions by employing a numerical conversion method based on the inverse fast Fourier transformation. The present method is validated by comparing with the solution given by other methods. Lastly, the stress analysis of multilayered coatings with various structure layouts and various layer number of the multilayers were performed with the present method.


Introduction
With the progress of surface engineering in the past several decades, multilayer coatings are often used in the surface modification of tribo-parts to enhance their tribological performance and service life under extreme service conditions, such as high vacuum, high radiation, high load, and high temperature [1][2][3][4]. Meanwhile, from the mechanical behavior point of view, a deformation and stress analysis of multilayered solids is indispensable to their optimal design and engineering application [5].
Many researches have been carried out theoretically and numerically to understand the mechanical behavior of coated solids by using different methods. Besides the finite element method [6][7][8][9], boundary element method [10,11], and image point method [12][13][14], the integral transform techniques [15] were most popular in producing the elastic filed solution of layered mediums. If the fundamental solutions that satisfy the governing equations are obtained with integral transform techniques or other methods, then a semi-analytical method or model, which always has a higher rate of convergence and consumes less solving time than the other numerical methods, can be developed to analyze the engineering problems [16][17][18]. Therefore, the fundamental solutions that satisfy the governing equations are the key to developing the semi-analytical method or model. In the early research work, the substrate of the coating-substrate system was treated as a rigid body [19,20]. Based on the Airy stress function, a two-dimensional elastic field solution in terms of Fourier integral transform (FT) for solids coated with a monolayer was given and the substrate was treated as an elastic body [21]. However, the tangential traction applied on the surface was not considered. Based on the formulas in Ref. [21], King and O'Sullivan [22] further gave the two-dimensional (2D) frequency response functions (FRF) of monolayer coatings with both the normal pressure and the tangential traction considered, which is the 2D elastic field solution in the frequency domain of the monolayer coating. Based on the 2D FRF of monolayer coatings, Wang et al. built a 2D contact model between a solid coated with a single layer and a rigid cylinder to analyze the influence of the coating thickness on its mechanical behavior under various friction coefficients [23]. Elsharkawy and Hamrock described the basic steps of solving the 2D FRF of the multilayer coating-substrate system by introducing the FT [24]. However, the solving process becomes increasingly complicated with the increase of the layer number, and only the elastic field solution of a bilayer coating-substrate system was presented. Elsharkawy and Hamrock [25] further investigated the effect of the frictional force on the stress field based on their research in Ref. [24]. In order to model the contact problem of a graded layer, Ke et al. [26][27][28] also adopted the FT technique to obtain the 2D FRF of multilayered solids by assuming that the elasticity modulus of each layer meets the exponential variation law along its depth, the elasticity modulus is continuous on the interface and the Poisson ratio is the same to each layer. In the past several decades, considerable efforts have been devoted to the 2D mechanical behavior analysis of solids coated with monolayer or bilayer coatings. However, only a few efforts have been devoted to the 2D mechanical behavior analysis involved solids coated with multilayer coatings. In the research involving multilayer coatings, the 2D FRF was deduced by forcing the elasticity moduli and Possion ratios of two adjacent layers to be continue on their interface. To my best knowledge, too few attentions have been paid to the 2D elastic filed solution multilayer coating-substrate systems without any mandatory requirement on the interface continuity of the elasticity modulus and Possion ratios. Furthermore, the mechanical behavior of multilayered solids with various layouts and layer numbers of the multilayers subjected to the surface contact loading, which may potentially play an important role in the optimal design and performance improvement of multilayered solids, still receives too little attention.
Therefore, efforts are devoted to developing a 2D semi-analytical method for the mechanical behavior analysis of multilayered solids subjected to surface contact loading in this paper. First, based on the 2D general FRF of layered material, a system of linear equations involved the unknown constants in the 2D general FRF of layered materials is established according to the boundary conditions and the interface continuous conditions. Next, the 2D FRF of multilayered solids is obtained by determining the unknown constants in the 2D general FRF of the layered materials through solving a system of linear equations analytically. Then, the 2D elastic field solution of multilayered solids is converted from its 2D FRF by employing a numerical conversion method based on the inverse fast Fourier transform (IFFT). The present method of multilayered solids is validated by comparisons with the finite element method as well as the exact analytical solution given by McEwen. Lastly, several specific conclusions are drawn based on the numerical investigation on the mechanical behavior of multilayered solid with various layouts and various layer number of the multilayered system by utilizing the present method. Figure 1 illustrates a 2D plain strain problem, in which an elastic solid coated with multilayers is subjected to surface contact loading p(x) and q(x). In the figure, N is the layer number of the multilayer coatings, while h k , E k and ν k are the thickness, elasticity modulus and Poisson ratio of the kth layer, respectively. E N+1 and ν N+1 are the elasticity modulus and Poisson ratio of the substrate material, respectively. p(x) and q(x) are the surface normal pressure and the tangential traction applied on the surface of the multilayer coatings, respectively. The normal pressure and the tangential traction applied on the surface are given by Hertz line contact theory as follows [29]:

Statement of the Mechanical Problem
where p H and b H represent the maximum contact pressure and the contact half-width of the substrate material in the Hertz line contact, and µ f represents the friction coefficient.
Coatings 2020, 10, x FOR PEER REVIEW 3 of 17 k th layer, respectively. EN+1 and νN+1 are the elasticity modulus and Poisson ratio of the substrate material, respectively. p(x) and q(x) are the surface normal pressure and the tangential traction applied on the surface of the multilayer coatings, respectively. The normal pressure and the tangential traction applied on the surface are given by Hertz line contact theory as follows [29]: where pH and bH represent the maximum contact pressure and the contact half-width of the substrate material in the Hertz line contact, and μf represents the friction coefficient.

Two-Dimensional Analytical Frequency Response Functions of Multilayered Solids
For a 2D plane strain problem, the stress and displacement components of layered materials in frequency domain are given as Equations (3)-(7) [22,23] based on the Airy stress function by using the FT technique. x k x k where the symbol 'i' represents the imaginary unit, the symbol '~' represents the one-dimensional Fourier integral transform. A1 (k) , A2 (k) , A3 (k) are A4 (k) are the unknown constants in the 2D FRF of the k th layered material. There are 4(N + 1) − 2 unknown constants that need to be determined in total to obtain the 2D FRF of a multilayered solid with N layers, since the unknown constants A3 (k) and A4 (k) are zero because

Two-Dimensional Analytical Frequency Response Functions of Multilayered Solids
For a 2D plane strain problem, the stress and displacement components of layered materials in frequency domain are given as Equations (3)-(7) [22,23] based on the Airy stress function by using the FT technique.
where the symbol 'i' represents the imaginary unit, the symbol '~' represents the one-dimensional Fourier integral transform. A 1 (k) , A 2 (k) , A 3 (k) are A 4 (k) are the unknown constants in the 2D FRF of the kth layered material. There are 4(N + 1) − 2 unknown constants that need to be determined in total to obtain the 2D FRF of a multilayered solid with N layers, since the unknown constants A 3 (k) and A 4 (k) are zero because all the stress and displacement components of the substrate become zero when z N+1 tends to infinity. The 4(N + 1) − 2 unknown constants can be determined with the boundary conditions and the interface continuous conditions of the multilayered solids, which are: where 1 ≤ k ≤ N, p and q are the one-dimensional FT of the normal pressure p(x) and the tangential traction q(x). Using the results given in Ref. [30,31], the one-dimensional FT of the normal pressure p(x) and the tangential traction q(x) can be given as Equations (14) and (15), respectively.
Coatings 2020, 10, 429 The system of linear equations composed of Equations (16)- (25) can be solved simultaneously with a numerical method according to Ref. [32]. However, the numerical method will consume a lot solving time especially when the layer number is large, since an equation group with 4(N + 1) − 2 equations should be solved for λN x times. The meaning of the symols λ and N x can be seen in Section 2.3. Therefore, analogous to the method proposed in Ref. [33], an analytical solution of the unknown constants in Equations (3)-(7) is derived and given in a recursive form in this paper. Equations (22)-(25) based on the continuous conditions of the last interface are used to solved the unknown constants A 1 (N+1) and A 2 (N+1) of the substrate as functions of the unknown constants of the kth layer.
Two more relationships of the unknown constants A 1 (N) , A 2 (N) , A 3 (N) and A 4 (N) are given as: where the S terms in Equations (28) and (29) where the t terms in Equations (30)-(33) are the intermediate constants which are also listed in the Appendix A. Assume that two relationships of the unknown constants A 1 (k+1) , A 2 (k+1) , A 3 (k+1) , and A 4 (k+1) of the (k + 1)th layer have been developed as follows: Coatings 2020, 10, 429 6 of 17 By appropriate algebra operations among the equations shown in Equations (30)- (35), two relationships of the unknown constants A 1 (k) , A 2 (k) , A 3 (k) , and A 4 (k) of the kth layer can be deduced as follows: where the S terms in Equations (36) and (37) are the intermediate constants and listed in the Appendix A. It can be found that the S terms of the kth layer contain those of the (k + 1)th layer and their algebraic expressions are true for every layer of the multilayer system. Therefore, the two relationships of the unknown constants A 1 1 , A 2 1 , A 3 1 , and A 4 1 of the 1st layer are: The S terms of the 1st layer in Equations (36) and (37)  Equations (16), (17), (38), and (39), and given as follows: Then the unknown constants in the 2D general FRF of the other layers and the substrate can be obtained with their recursive formulas shown in Equations (28) and (29) and Equations (22)-(25) layer by layer in a top-down fashion.

Numerical Conversion Method
Having obtained the 2D explicit analytical FRF of multilayered solids, the elastic field solution in the space domain can be converted from its corresponding FRF by using an IFFT based numerical conversion method. The conversion method originally is used to produce the influence coefficients of the stress and displacement components for developing a semi-analytical contact model [34,35]. The detailed steps of using the IFFT based numerical conversion method to produce the stress component σ rs at z depth are as follows: • Select a computing domain {x|x b ≤ x ≤ x e } at z depth and divide it uniformly into N x line elements. The number N x is required to be a positive integral power of 2. x b and x e are usually chosen to be −2b H and 2b H , respectively. • Refine the grid into N ωx uniform elements in the frequency domain {ω x |−π < ω x ≤ π}, where N ωx = λN x and λ should be a nonnegative integral power of 2.
• Construct a discrete seriesσ rs by using the FRF of the stress component σ rs : Coatings 2020, 10, 429 7 of 17 • Reconstruct a new discrete series σ rs by applying a wrap-order operation to σ rs : • Apply the IFFT operation to the discrete series σ rs : • The stress at each line element in the discrete space domain can be obtained as follows: A very high conversion accuracy can be achieved by selecting an enough large λ and κ to eliminate the periodic error and the error caused by aliasing effect. However, the increase of λ and κ will lead to a heavier calculation burden. In this paper, λ and κ are fixed as 4 and 5, respectively because the elastic field solution of the present method agrees well with the exact analytical solution given by McEwen [36] and the solution of the finite element method, when λ = 4 and κ = 5. It can be found in Section 3 of this paper.

Validation of the Present Method
In order to validate the present method, the maximum shear stress in the subsurface of a half plane of the substrate material subjected to surface mechanical loading obtained by the present method is compared with the exact analytical solution given by McEwen seen in Ref. [36]. The input parameters of the present method are listed in Table 1. The elasticity parameters of the layer are the same to those of the substrate when the present method is adopted to simulate the problem of the half plane of the substrate material subjected to surface contact loading. Table 1. Input parameters of the present method for simulating a half plane of the substrate material subjected to subjected to surface contact loading.    In order to further validate the present method, a comparison of the maximum shear stress in the subsurface of a multilayered solid with 40 layers is also conducted between the solution of the present method and that of the finite element method. The parameters of the multilayered solid and the other concerning parameters are shown in Table 2. Table 2. Input parameters of the present method for simulating the multilayered solid with 40 layers subjected to surface contact loading.
The finite element analysis of a multilayered solids subjected to surface mechanical loading with Hertzian distribution is performed with a commercial software ANSYS 15.0. The FEM model of the multilayered solid subjected to surface mechanical loading is shown in Figure 3. In order to make the FEM model better approximate to the half-plane assumption of the semi-analytical method developed in this paper, the geometric model, which is a rectangle with a size of 600b H in x direction and 300b H in z direction, is built and meshed with the finite element named Plane 182, since the two-dimensional mechanical problem of a multilayered solid subjected to surface mechanical loading is a plain strain problem. The displacements of the nodes located at the bottom of the FEM model are set to be 0, and the surface mechanical loading with Hertzian distribution is applied on the surface loading region of the FEM model by using the finite element named Surf 153.   The finite element analysis of a multilayered solids subjected to surface mechanical loading with Hertzian distribution is performed with a commercial software ANSYS 15.0. The FEM model of the multilayered solid subjected to surface mechanical loading is shown in Figure 3. In order to make the FEM model better approximate to the half-plane assumption of the semi-analytical method developed in this paper, the geometric model, which is a rectangle with a size of 600bH in x direction and 300bH in z direction, is built and meshed with the finite element named Plane 182, since the twodimensional mechanical problem of a multilayered solid subjected to surface mechanical loading is a plain strain problem. The displacements of the nodes located at the bottom of the FEM model are set to be 0, and the surface mechanical loading with Hertzian distribution is applied on the surface loading region of the FEM model by using the finite element named Surf 153.    Figure 4a for both µ f = 0 and µ f = 0.2. The relative error of the peak value of the maximum shear stress is not more than 0.15% for both µ f = 0 and µ f = 0.2.
Coatings 2020, 10, x FOR PEER REVIEW 9 of 17 Figure 4 is the contour plots of the maximum shear stress of the multilayered solid with 40 layers subjected to surface contact loading under different friction coefficients. Obviously, the contour plots of the maximum shear stress obtained with the present method shown in Figure 4b are in a good agreement with those obtained with the finite element method shown in Figure 4a for both μf = 0 and μf = 0.2. The relative error of the peak value of the maximum shear stress is not more than 0.15% for both μf = 0 and μf = 0.2. The comparisons presented above not only validate the present method and the IFFT based numerical conversion method adopted in this paper, but also show that the present method adapts to multilayer coatings with various layer numbers and layer thicknesses.

Results and Discussions
In order to get a deep insight into the mechanical behavior of multilayered solids with various structure layouts and various layer numbers of the multilayers subjected to the surface contact loading, numerical simulations are conducted with the present method. In all numerical simulation cases, the maximum contact pressure pH and the contact half-width bH are the same to those given in Table 1. The Poisson ratios of the layers and the substrate are all fixed at 0.3. The elasticity modulus of the top layer is fixed at 800 GPa, the elasticity modulus of the substrate materials is fixed at 200 GPa, and the moduli of the other layers is due to the layout of the multilayers. Figure 5 is the contour plots of the maximum shear stress in the subsurface of multilayered solids with various structure layouts when the friction coefficient is 0 and 0.5, respectively. The elasticity modulus along z axis of various structure layouts of the multilayers is shown in Figure 5a. The layer number N is 10, and all layers from the top layer to the bottom layer have the same thickness of bH/10. The first layout of the multilayer is a monolayer coating since the elasticity modulus remains unchanged from the top layer to the bottom layer and changes abruptly on the interface between the bottom layer and the substrate. The other four layouts of the multilayers afford four alternatives to avoid the sudden change of elasticity modulus as that of the 1st layout through a gradient change layer by layer. For the first layout, a significant discontinuity can be found on the interface between the bottom layer and the substrate. The peak value of the maximum shear stress is 0.567pH for μf = 0 and locates at the bottom of the bottom layer. This is a result of the sudden decrease of elasticity modulus on from 800 GPa of the layers to 200 GPa of the substrate. The peak value of the maximum shear stress increases rapidly to 0.816pH for μf = 0.5 and its location shifts to the surface due to the The comparisons presented above not only validate the present method and the IFFT based numerical conversion method adopted in this paper, but also show that the present method adapts to multilayer coatings with various layer numbers and layer thicknesses.

Results and Discussions
In order to get a deep insight into the mechanical behavior of multilayered solids with various structure layouts and various layer numbers of the multilayers subjected to the surface contact loading, numerical simulations are conducted with the present method. In all numerical simulation cases, the maximum contact pressure p H and the contact half-width b H are the same to those given in Table 1. The Poisson ratios of the layers and the substrate are all fixed at 0.3. The elasticity modulus of the top layer is fixed at 800 GPa, the elasticity modulus of the substrate materials is fixed at 200 GPa, and the moduli of the other layers is due to the layout of the multilayers. Figure 5 is the contour plots of the maximum shear stress in the subsurface of multilayered solids with various structure layouts when the friction coefficient is 0 and 0.5, respectively. The elasticity modulus along z axis of various structure layouts of the multilayers is shown in Figure 5a. The layer number N is 10, and all layers from the top layer to the bottom layer have the same thickness of b H /10. The first layout of the multilayer is a monolayer coating since the elasticity modulus remains unchanged from the top layer to the bottom layer and changes abruptly on the interface between the bottom layer and the substrate. The other four layouts of the multilayers afford four alternatives to avoid the sudden change of elasticity modulus as that of the 1st layout through a gradient change layer by layer. For the first layout, a significant discontinuity can be found on the interface between the bottom layer and the substrate. The peak value of the maximum shear stress is 0.567p H for µ f = 0 and locates at the bottom of the bottom layer. This is a result of the sudden decrease of elasticity modulus on from 800 GPa of the layers to 200 GPa of the substrate. The peak value of the maximum shear stress increases rapidly to 0.816p H for µ f = 0.5 and its location shifts to the surface due to the friction traction applied on the surface. For the second layout, the significant discontinuity of the maximum shear stress as that of the 1st layout does not exist anymore and a relatively small discontinuity is observed on the interface of each two adjacent layers, and the discontinuity on the interface of the top five layers is more visible than that on the interface of the bottom five layers. The maximum shear stress decreases in the layers close to the substrate while increases dramatically in the surface region, so its peak value locates on the surface and reaches 0.582p H for µ f = 0 and 0.960p H for µ f = 0.5 respectively. For the third layout, the contour plots of the maximum shear stress are very similar to those of the second layout. The peak value of the maximum shear stress is 0.548p H for µ f = 0 and 0.886p H for µ f = 0.5 respectively and both locate on the surface of the multilayered solids. For the 4th layout, no discontinuity of the maximum shear stress is observed from the first layer to the fifth layer due to the same elasticity modulus of the top five layers. The maximum shear stress also reduces in the layers close to the substrate and increases in the surface region. The peak value of the maximum shear stress is 0.526p H for µ f = 0 and 0.836p H for µ f = 0.5 respectively and both locate on the surface too. For the fifth layout, the decrease of the maximum shear stress in the layers close to the substrate is not as significant as those of the second layout to fourth layout, and the increase of the maximum shear stress in the surface region is also not as significant as those of the second layout to fourth layout. The peak value of the maximum shear stress is 0.487p H for µ f = 0 and 0.764p H for µ f = 0.5 respectively, which are the smallest among the five layouts of the multilayers investigated numerically in this research.

Influence of the Structure Layout of the Multilayers on the Mechanical Behavior
Coatings 2020, 10, x FOR PEER REVIEW 10 of 17 friction traction applied on the surface. For the second layout, the significant discontinuity of the maximum shear stress as that of the 1st layout does not exist anymore and a relatively small discontinuity is observed on the interface of each two adjacent layers, and the discontinuity on the interface of the top five layers is more visible than that on the interface of the bottom five layers. The maximum shear stress decreases in the layers close to the substrate while increases dramatically in the surface region, so its peak value locates on the surface and reaches 0.582pH for μf = 0 and 0.960pH for μf = 0.5 respectively. For the third layout, the contour plots of the maximum shear stress are very similar to those of the second layout. The peak value of the maximum shear stress is 0.548pH for μf = 0 and 0.886pH for μf = 0.5 respectively and both locate on the surface of the multilayered solids. For the 4th layout, no discontinuity of the maximum shear stress is observed from the first layer to the fifth layer due to the same elasticity modulus of the top five layers. The maximum shear stress also reduces in the layers close to the substrate and increases in the surface region. The peak value of the maximum shear stress is 0.526pH for μf = 0 and 0.836pH for μf = 0.5 respectively and both locate on the surface too. For the fifth layout, the decrease of the maximum shear stress in the layers close to the substrate is not as significant as those of the second layout to fourth layout, and the increase of the maximum shear stress in the surface region is also not as significant as those of the second layout to fourth layout. The peak value of the maximum shear stress is 0.487pH for μf = 0 and 0.764pH for μf = 0.5 respectively, which are the smallest among the five layouts of the multilayers investigated numerically in this research.   Figure 6 shows the influence of the structure layouts of the multilayered solids on the peak value of the maximum shear stress under various friction coefficients. With the increase of the friction coefficient, the peak value of the maximum shear stress increases slightly when µ f ranges from 0 to 0.2 for various layouts of multilayered solids, while increases relatively significantly when µ f ranges from 0.2 to 0.5. For the second layout to fifth layout, the peak value of the maximum shear stress decreases in the order of the second layout, the third layout, the fourth layout and the fifth layout under various friction coefficients. As indicated in the figure, the peak value of the maximum shear stress of the second layout is larger than that of the first layout under various friction coefficients. The peak value of the maximum shear stress of the third layout and the fourth layout is smaller than that of the first layout when the friction coefficient is smaller than 0.1, and it is opposite when the friction coefficient is larger than about 0.2. It can be further observed that only the peak value of the maximum shear stress of the fifth layout is smaller than that of the first layout under various friction coefficients. Further, the peak value of the maximum shear stress of the fifth layout is about 16.3% smaller than that of the second layout when the friction coefficient is 0. Therefore, the peak value of the maximum shear stress, which is closely related to the yield and fatigue damage of materials [37,38], can be reduced effectively by optimizing the structure layout of multilayered solids and minimized by employing a structure layout with elasticity modulus increasing first in the top layers and then decreasing in the bottom layers just as the fifth layout shown in Figure 5a.
Coatings 2020, 10, x FOR PEER REVIEW 11 of 17 Figure 6 shows the influence of the structure layouts of the multilayered solids on the peak value of the maximum shear stress under various friction coefficients. With the increase of the friction coefficient, the peak value of the maximum shear stress increases slightly when μf ranges from 0 to 0.2 for various layouts of multilayered solids, while increases relatively significantly when μf ranges from 0.2 to 0.5. For the second layout to fifth layout, the peak value of the maximum shear stress decreases in the order of the second layout, the third layout, the fourth layout and the fifth layout under various friction coefficients. As indicated in the figure, the peak value of the maximum shear stress of the second layout is larger than that of the first layout under various friction coefficients. The peak value of the maximum shear stress of the third layout and the fourth layout is smaller than that of the first layout when the friction coefficient is smaller than 0.1, and it is opposite when the friction coefficient is larger than about 0.2. It can be further observed that only the peak value of the maximum shear stress of the fifth layout is smaller than that of the first layout under various friction coefficients. Further, the peak value of the maximum shear stress of the fifth layout is about 16.3% smaller than that of the second layout when the friction coefficient is 0. Therefore, the peak value of the maximum shear stress, which is closely related to the yield and fatigue damage of materials [37,38], can be reduced effectively by optimizing the structure layout of multilayered solids and minimized by employing a structure layout with elasticity modulus increasing first in the top layers and then decreasing in the bottom layers just as the fifth layout shown in Figure 5a.  Figure 7 is the contour plots of the maximum shear stress in the subsurface of multilayered solids with a similar structure layout of the multilayers but a different layer number. It can be seen from Figure 7b,c that the discontinuity of the maximum shear stress on the interface of each two adjacent layers becomes more invisible with the increase of the layer number for both μf = 0 and μf = 0.5, which is a result of that the difference in the elasticity moduli of each two neighboring layers becomes smaller with the increase of the layer number as shown in Figure 7a. The peak value of the maximum shear stress of the solid coated with three layers locates at the bottom of the third layer. However, it moves to the surface when the layer number increases to nine layers and 36 layers. A decrease of the peak value of the maximum shear stress can be observed when the layer number increases from three to nine, while almost no difference of the peak value of the maximum shear stress is observed when the layer number increases from 9 to 36.  Figure 7 is the contour plots of the maximum shear stress in the subsurface of multilayered solids with a similar structure layout of the multilayers but a different layer number. It can be seen from Figure 7b,c that the discontinuity of the maximum shear stress on the interface of each two adjacent layers becomes more invisible with the increase of the layer number for both µ f = 0 and µ f = 0.5, which is a result of that the difference in the elasticity moduli of each two neighboring layers becomes smaller with the increase of the layer number as shown in Figure 7a. The peak value of the maximum shear stress of the solid coated with three layers locates at the bottom of the third layer. However, it moves to the surface when the layer number increases to nine layers and 36 layers. A decrease of the peak value of the maximum shear stress can be observed when the layer number increases from three to nine, while almost no difference of the peak value of the maximum shear stress is observed when the layer number increases from 9 to 36. Coatings 2020, 10, x FOR PEER REVIEW 12 of 17 The stress σxx along z axis of multilayered solids with a similar layout of the multilayers but a different layer number is shown in Figure 8 when the friction coefficient is 0. It can be observed that the stress σxx is discontinuous on the interface of each two neighboring layers, and the discontinuity of the stress σxx also becomes more invisible with the increase of the layer number. As indicated in Figure 7, the stress σxx is tensile in the region close to the substrate, and the tensile stress σxx is so large that its maximum value reaches up to 0.5pH when the layer number is three. It also can be found that the maximum value of the tensile stress σxx in the region close to the substrate becomes smaller with the increase of the layer number, although the tensile stress σxx could not be completely eliminated by increasing the layer number. The maximum tensile stress σxx decreases to 0.28pH when the layer number increases to 36. Obviously, the maximum tensile stress σxx, which is believed to be responsible for the initiation and propagation of cracks at the base of hard coatings [39,40], can be reduced effectively by increasing the layer number.

Conclusions
A semi-analytical method for the 2D mechanical behavior analysis of multilayered solids has been developed and validated in this paper. The influence of the structure layout and the layer number on the mechanical behavior of the multilayered solid has also been numerically investigated based on the present method. Several conclusions are drawn as follows: The stress σ xx along z axis of multilayered solids with a similar layout of the multilayers but a different layer number is shown in Figure 8 when the friction coefficient is 0. It can be observed that the stress σ xx is discontinuous on the interface of each two neighboring layers, and the discontinuity of the stress σ xx also becomes more invisible with the increase of the layer number. As indicated in Figure 7, the stress σ xx is tensile in the region close to the substrate, and the tensile stress σ xx is so large that its maximum value reaches up to 0.5p H when the layer number is three. It also can be found that the maximum value of the tensile stress σ xx in the region close to the substrate becomes smaller with the increase of the layer number, although the tensile stress σ xx could not be completely eliminated by increasing the layer number. The maximum tensile stress σ xx decreases to 0.28p H when the layer number increases to 36. Obviously, the maximum tensile stress σ xx , which is believed to be responsible for the initiation and propagation of cracks at the base of hard coatings [39,40], can be reduced effectively by increasing the layer number. The stress σxx along z axis of multilayered solids with a similar layout of the multilayers but a different layer number is shown in Figure 8 when the friction coefficient is 0. It can be observed that the stress σxx is discontinuous on the interface of each two neighboring layers, and the discontinuity of the stress σxx also becomes more invisible with the increase of the layer number. As indicated in Figure 7, the stress σxx is tensile in the region close to the substrate, and the tensile stress σxx is so large that its maximum value reaches up to 0.5pH when the layer number is three. It also can be found that the maximum value of the tensile stress σxx in the region close to the substrate becomes smaller with the increase of the layer number, although the tensile stress σxx could not be completely eliminated by increasing the layer number. The maximum tensile stress σxx decreases to 0.28pH when the layer number increases to 36. Obviously, the maximum tensile stress σxx, which is believed to be responsible for the initiation and propagation of cracks at the base of hard coatings [39,40], can be reduced effectively by increasing the layer number.

Conclusions
A semi-analytical method for the 2D mechanical behavior analysis of multilayered solids has been developed and validated in this paper. The influence of the structure layout and the layer number on the mechanical behavior of the multilayered solid has also been numerically investigated based on the present method. Several conclusions are drawn as follows:

Conclusions
A semi-analytical method for the 2D mechanical behavior analysis of multilayered solids has been developed and validated in this paper. The influence of the structure layout and the layer number on the mechanical behavior of the multilayered solid has also been numerically investigated based on the present method. Several conclusions are drawn as follows: • A 2D analytical frequency response functions of multilayered solids subjected to surface contact loading is derived and given in a recursive form. • A 2D semi-analytical method for the mechanical behavior analysis of multilayered solids is developed by converting from its frequency response functions with an IFFT based numerical conversion method.

•
Among various structure layouts numerically investigated in this paper, the layout with elasticity modulus increasing first in the top layers and then decreasing in the bottom layers can minimize the peak value of the maximum shear stress under various friction coefficients.

•
The maximum tensile stress σ xx in subsurface can be reduced effectively by increasing the layer number of multilayered solids.

Conflicts of Interest:
The authors declare no conflict of interest.