Analysis of the Surface / Interface Damage Evolution Behavior of a Coating–Substrate System under Heavy-Load Elastohydrodynamic Lubrication

: Keeping a coating–substrate system undamaged during heavy-load elastohydrodynamic lubrication (EHL) conditions is challenging. To overcome this problem, an EHL model with a coated gear pair was built. Firstly, based on the full-system ﬁnite element method, the e ﬀ ect of the coating elastic modulus on the oil ﬁlm pressure was obtained. Secondly, the failure mode was predicted after the stress analysis. Finally, the surface / interface damage evolution behavior of the coating–substrate system was analyzed visually by embedding cohesive zone elements. In the numerical calculation, sti ﬀ er coatings tended to increase the ﬁlm pressure and secondary pressure spike, compared with more compliant coatings. As the coating sti ﬀ ness decreased, the maximum equivalent stress in the system reduced, and its location tended to develop close to or at the substrate. The coating cracking and interfacial delamination were individually caused by the shear stress in the coating and shear stress on the interface, and both of them initiated in the region of the secondary pressure peak. The interfacial delamination increased the crack failure probability of coating and vice versa. Therefore, through analyzing the EHL model, the exact damage growth location and its evolution in the coated solids can be determined, and the failure mechanism can be comprehensively revealed.


Introduction
As an effective way to achieve excellent surface load carrying capacity, coating technology has been widely employed in the field of rolling-sliding contact industrial friction pairs which are under heavy-duty elastohydrodynamic lubrication (EHL) conditions [1], such as gears, bearings, cams, etc. [2,3]. It is particularly important to reveal the EHL mechanism of coated couples. On the other hand, there are two main failure modes of thin coatings in service: coating cracking and interfacial delamination. The lack of a method to analyze the failure behavior of the coating-substrate system under heavy-load EHL conditions is an urgent problem that needs to be solved.
In recent years, many scholars have studied the EHL tribological behavior of coating-substrate systems. Based on fast Fourier inverse transform or the finite element (FE) method, Liu et al. [4] carried out point contact EHL analysis for a coating-substrate system. Wang et al. [5] studied the effect of coating material parameters on the stress distribution of coated parts. Liu et al. [6] and effect of coating material parameters on the stress distribution of coated parts. Liu et al. [6] and Xiao et al. [7] analyzed the EHL problem of line contact in coated gear transmissions. Habchi [8] used the FE method combined with the solution of the Reynolds equation to solve the point contact EHL issue for coated solids. Alakhramsing et al. [9] investigated EHL characteristics of coating-substrate systems with finite line contact taking non-Newtonian fluid into account. Stahl et al. [10] acquired the friction reduction mechanism of diamond-like carbon (DLC)-coated gear pairs coupling it with the temperature field. From the aspect of coating-substrate failure under the normal or tangential load, previous studies mostly focused on nano-indentation or scratch characterization experiments combined with cohesive zone model (CZM) [11][12][13] or virtual crack closure technology (VCCT) [12] to study the law of coating cracking and interface debonding damages [14]. Lin et al. [15] used the FE method to predict the interfacial delamination of a multi-layered system under indentation loading at the macroscale with the cohesive zone being considered around the delamination or crack tips. The failure properties of a DLC/steel system were evaluated using a combined nano-indentation and CZM coupled FE approach [16]. Early in the development lifecycle, the authors performed the running-in test of the DLC-coated gears under the heavy-load EHL conditions. As with other hard coatings [14], DLC deposited via physical vapor deposition usually results in a columnar grain structure. To provide additional insight into the cracked and delaminated area induced by the gear meshing force, the corresponding gear sample was milled by a focused ion beam (FIB) to obtain a cross-sectional view. The microstructural image ( Figure 1) was obtained by scanning electron microscopy (SEM). It can be seen that the cracks vertically extend through the thickness of the entire coating (intergranular sliding of columnar grains), as well as a typical feature of delaminated area on the interface. It means that a coating's cracking often is often followed by interfacial delamination under meshing load. Previous studies on the damage evolution behavior of coating-substrate systems in heavy-load EHL service ( Figure 2) are rare, and there are few studies that consider the interaction between coating cracking and interfacial delamination. Thus, these gaps were the focus of the this paper, and the above methods were the foundation for the present research. Among the above scholars, Habchi's full-system FE method not only has an outstanding advantage in regards to solving speed, but can also establish a good interface for subsequent modelling of coating-substrate failure. Moreover, the cohesive zone model (CZM) was adopted because it has superiority in the following two aspects: 1) stress is a function of crack displacement in the failure process, which avoids the singularity of crack tip stress; and 2) there is no need for a pre-crack before calculation.
The main problem now is that when the CZM is coupled as fracture/damage behavior in the coating-substrate system, the non-linearity becomes stronger when solving the heavy-load EHL model. In the current paper, the full-system FE method combined with CZM was employed to analyze the complex mechanical problems involved in the surface contact and failure behavior of the coatings. The inverse analysis [17] of the coating-substrate can make up for the defect of the over-reliance on trial-and-error methods in the design of coatings, and can predict failures in the system in advance, which can reduce the selection time of coatings as well as save costs related to testing. Therefore, accurately obtained internal and interfacial stressors of the coated system, timely Previous studies on the damage evolution behavior of coating-substrate systems in heavy-load EHL service ( Figure 2) are rare, and there are few studies that consider the interaction between coating cracking and interfacial delamination. Thus, these gaps were the focus of the this paper, and the above methods were the foundation for the present research. Among the above scholars, Habchi's full-system FE method not only has an outstanding advantage in regards to solving speed, but can also establish a good interface for subsequent modelling of coating-substrate failure. Moreover, the cohesive zone model (CZM) was adopted because it has superiority in the following two aspects: 1) stress is a function of crack displacement in the failure process, which avoids the singularity of crack tip stress; and 2) there is no need for a pre-crack before calculation.
The main problem now is that when the CZM is coupled as fracture/damage behavior in the coating-substrate system, the non-linearity becomes stronger when solving the heavy-load EHL model. In the current paper, the full-system FE method combined with CZM was employed to analyze the complex mechanical problems involved in the surface contact and failure behavior of the coatings. The inverse analysis [17] of the coating-substrate can make up for the defect of the over-reliance on trial-and-error methods in the design of coatings, and can predict failures in the system in advance, which can reduce the selection time of coatings as well as save costs related to testing. Therefore, accurately obtained internal and interfacial stressors of the coated system, timely prediction of the failure of the coated components, and improvement of the surface load carrying capacity of the coated friction pairs creates a foundation for the application of the coating in the field of heavy-duty transmission. Furthermore, inverse analysis emphasizes the usefulness and significance of the FE model. prediction of the failure of the coated components, and improvement of the surface load carrying capacity of the coated friction pairs creates a foundation for the application of the coating in the field of heavy-duty transmission. Furthermore, inverse analysis emphasizes the usefulness and significance of the FE model.

EHL Model
The computational domain of the numerical model is shown in Figure 3. The EHL problem is solved under the assumed plane strain condition. It is appropriate to define the range of two-dimensional calculation domain Ω as ∈ X [- 30,30] and h ∈ c Y a [-/ -60, 0] for solving elastic deformation from References [8,18]. On the other hand, to ensure that the oil film pressure at the inlet and outlet is 0, the calculation of the lubricant domain can be calculated on the one-dimensional contact area c Ω ( ∈ X [-4.5, 3] ). The fixed-boundary condition is applied to the substrate base ( ). In order to solve the generalized Reynolds and elastic equations in the EHL problem conveniently, the parameters of these equations are firstly dimensionless and treated as follows:

EHL Model
The computational domain of the numerical model is shown in Figure 3. The EHL problem is solved under the assumed plane strain condition. It is appropriate to define the range of two-dimensional calculation domain Ω as X ∈ [− 30,30] and Y ∈ [−h c /a − 60, 0] for solving elastic deformation from References [8,18]. On the other hand, to ensure that the oil film pressure at the inlet and outlet is 0, the calculation of the lubricant domain can be calculated on the one-dimensional contact area Ω c (X ∈ [−4.5, 3]). The fixed-boundary condition is applied to the substrate base (∂Ω b = 0). prediction of the failure of the coated components, and improvement of the surface load carrying capacity of the coated friction pairs creates a foundation for the application of the coating in the field of heavy-duty transmission. Furthermore, inverse analysis emphasizes the usefulness and significance of the FE model.

EHL Model
The computational domain of the numerical model is shown in Figure 3. The EHL problem is solved under the assumed plane strain condition. It is appropriate to define the range of two-dimensional calculation domain Ω as ∈ X [- 30,30] and h ∈ c Y a [-/ -60, 0] for solving elastic deformation from References [8,18]. On the other hand, to ensure that the oil film pressure at the inlet and outlet is 0, the calculation of the lubricant domain can be calculated on the one-dimensional contact area c Ω ( ∈ X [-4.5, 3] ). The fixed-boundary condition is applied to the substrate base ( ). In order to solve the generalized Reynolds and elastic equations in the EHL problem conveniently, the parameters of these equations are firstly dimensionless and treated as follows:  In order to solve the generalized Reynolds and elastic equations in the EHL problem conveniently, the parameters of these equations are firstly dimensionless and treated as follows: where P, τ, η, ρ, X, Y, U, V, and H are dimensionless forms of the oil film pressure p, oil film shear stress τ, oil viscosity η, lubricant density ρ, x and y components of the coordinates, solid elastic deformation along x and y directions, and oil film thickness h, respectively. η 0 and ρ 0 individually correspond to the ambient viscosity and density of the lubricant. p h and a denote the maximum Hertzian contact pressure and half-width of uncoated solids, and are defined as below: where F is the external load. Two substrates are made of a same material, so the reduced elastic modulus of the contacting bodies E can be described by the following equation: where E s and ν s are the elastic modulus and Poisson's ratio of the substrate, respectively. Constitutive equation of non-Newtonian fluid lubricants is constructed by the Ree-Eyring model. The equivalent oil film viscosity η * is shown as: where τ 0 is the Eyring characteristic shear stress. The relationship between viscosity and pressure is described by the Roelands [19] viscosity-pressure equation: where α is the Barus viscosity-pressure coefficient.
The density-pressure equation of Dowson-Higginson [20] is adopted: The dimensionless generalized Reynolds equation for calculating oil film pressure on Ω c is: The boundary conditions of the Reynolds equation are shown as follows [21]: on Ω c ∂P/∂X = 0 on the cavitation boundary (8) Dimensionless oil film thickness H is expressed as: where, H 0 is the approach distance of the rigid body, which is determined by the load balance equation: Coatings 2019, 9, 642 5 of 14

Cohesive Zone Model
In this paper, the CZM (bilinear traction-separation law) [22] was adopted to simulate the surface/interface failure behavior. Although an advanced CZM with more accurate behavior in the shear direction was mentioned by previous researchers [23,24], the currently selected model is a simple and effective way to solve the problem. Figure 4 shows the CZM coupled with tensile and shear loads, and its constitutive equation can be expressed as: where, t N,T is the separation traction stress. δ N,T is the separation traction displacement. δ N,T 0 is the critical displacement when damage occurs. δ N,T f is the displacement when the cohesive element fails completely. K N,T is the cohesive initial stiffness. D N,T is the cohesion cumulative damage stiffness (SDEG). When SDEG = 1, the cohesive element is completely invalid. The subscripts N and T in the parameters represent the normal and tangential directions, respectively.

Cohesive Zone Model
In this paper, the CZM (bilinear traction-separation law) [22] was adopted to simulate the surface/interface failure behavior. Although an advanced CZM with more accurate behavior in the shear direction was mentioned by previous researchers [23,24], the currently selected model is a simple and effective way to solve the problem. Figure 4 shows the CZM coupled with tensile and shear loads, and its constitutive equation can be expressed as: is the displacement when the cohesive element fails completely. N T K , is the cohesive initial stiffness. N T D , is the cohesion cumulative damage stiffness (SDEG). When SDEG = 1, the cohesive element is completely invalid. The subscripts N and T in the parameters represent the normal and tangential directions, respectively. Under mixed load, damage may occur when any stress component in normal or tangential directions reaches the limit value. Therefore, this paper selects the maximum nominal stress criterion to judge the failure: where < > N t stands for Macaulay parentheses.
The damage criterion of power function is used to evolve the damage process of coating and interface under mixed load, which can be expressed as: Where N Γ and T Γ are individually the work done by the normal and tangential forces in the process of damage. α is a power exponent. In this study, the power exponential α is defined as 1 [25][26][27], which is the most commonly used value and corresponds to a linear fracture criterion.  Under mixed load, damage may occur when any stress component in normal or tangential directions reaches the limit value. Therefore, this paper selects the maximum nominal stress criterion to judge the failure: where < t N > stands for Macaulay parentheses. The damage criterion of power function is used to evolve the damage process of coating and interface under mixed load, which can be expressed as: Where Γ N and Γ T are individually the work done by the normal and tangential forces in the process of damage. α is a power exponent. In this study, the power exponential α is defined as 1 [25][26][27], which is the most commonly used value and corresponds to a linear fracture criterion.

Numerical Calculation
The current paper takes the Forschungsstelle für Zahnräder und Getriebebau (FZG) standard gear pair (20CrMnTi) as the research object in the light of ISO 14635-1. Compared to the light load, the coating-substrate system is more likely to fail under heavy load. In addition, the solution to the Reynolds' equation is known to be unstable in the central contact area for highly loaded contacts. To make the solution method universal, a load stage of 12th was selected, which stands for heavy load. The DLC coating commonly used in the gear industry was chosen. The coating was obtained by low temperature vapor deposition and its thickness was 2 µm [28]. For a DLC coating, the great variety of its structures and sp 3 content leads to a wide range of mechanical property. Detailed geometric parameters, material properties, operating conditions, and lubricant properties are shown in Table 1. The global numerical solution scheme is shown in Figure 5. Based on the full-system FE method, all the equations to be solved are discretized under the steady-state calculation condition. The numerical calculation was carried out under given initial values H 0 , P, U, and V, and then the EHL equations of each time step were solved by direct FE and Newton-Raphson iterations to solve them. An implicit relationship was established between the fluid and solid domains and H 0 was continuously updated through the form of full coupling. Lagrangian fifth-order and second-order elements were applied to the Reynolds and elasticity equations, respectively. The specific EHL solution method is described in Reference [21]. Then, the EHL pressure and shear stress were coupled into the coated system to obtain the surface/interface damage evolution. The residual stress in the thin coating is usually approximately −0.02 to −2 GPa. Although it may be a bit large in the coating, it has only a little effect on failure due to the compressive residual stress. On the other hand, the residual stress is arduous to simulate using the FE method, and it does not affect the fracture/damage law explored in this work. This means that it is reasonable to consider the ideal situation of the system (residual stress-free) when it comes to studying the problem of the failure mechanism. However, it should be emphasized that residual stress will be taken into account in further works for improvement, which can then obtain more accurate and realistic results.

EHL Response
The full-system FE method was applied to analyze the line contact EHL model of coated solids, and the effect of coating elastic modulus on oil film pressure was studied. As shown in Figure 6, the calculated oil film pressure distribution curve conformed to the typical EHL characteristics. The lubricant was acting on the gear surface under the rolling-sliding contact. At a high contact pressure of 1841 MPa, the material deformed locally, and a parallel clearance was formed among the contact surfaces. Unlike Hertzian contact, there was a secondary peak pressure at the oil film outlet at the edge of the contact surface. On the whole, the five oil pressure curves basically coincided because of the thin coating, which means that its coating elastic modulus had little effect on the oil film pressure distribution. However, from the local magnification observation, with the increase in the coating elastic modulus (stiffness), the central pressure of the oil film increased, and the secondary pressure peak became more prominent. This conclusion is consistent with the results of the analysis of the uncoupled CZM on the interface. In order to predict the specific failure mode and location of the coating-substrate system, a stress analysis will be carried out during follow-up. Two typical cases were selected in order to analyze the subsequent damage evolution behavior of the system: 1) soft coating and hard substrate ( c E = 50 GPa ); and 2) hard coating and soft substrate ( c E = 500 GPa ).

EHL Response
The full-system FE method was applied to analyze the line contact EHL model of coated solids, and the effect of coating elastic modulus on oil film pressure was studied. As shown in Figure 6, the calculated oil film pressure distribution curve conformed to the typical EHL characteristics. The lubricant was acting on the gear surface under the rolling-sliding contact. At a high contact pressure of 1841 MPa, the material deformed locally, and a parallel clearance was formed among the contact surfaces. Unlike Hertzian contact, there was a secondary peak pressure at the oil film outlet at the edge of the contact surface. On the whole, the five oil pressure curves basically coincided because of the thin coating, which means that its coating elastic modulus had little effect on the oil film pressure distribution. However, from the local magnification observation, with the increase in the coating elastic modulus (stiffness), the central pressure of the oil film increased, and the secondary pressure peak became more prominent. This conclusion is consistent with the results of the analysis of the uncoupled CZM on the interface. In order to predict the specific failure mode and location of the coating-substrate system, a stress analysis will be carried out during follow-up. Two typical cases were selected in order to analyze the subsequent damage evolution behavior of the system: 1) soft coating and hard substrate (E c = 50 GPa); and 2) hard coating and soft substrate (E c = 500 GPa).

Equivalent Stress Analysis
Coupling the above EHL pressure and shear stress into the failure model of the coated solids, the equivalent stress of the system was calculated first. Figure 7 shows the von Mises equivalent

Equivalent Stress Analysis
Coupling the above EHL pressure and shear stress into the failure model of the coated solids, the equivalent stress of the system was calculated first. Figure 7 shows the von Mises equivalent stress distribution nephogram. The smaller the coating elastic modulus, the smaller the stress that will cause certain elastic deformation in the coating. Obviously, under the same external load, the maximum equivalent stress of the soft coating (1029 MPa) was smaller than that of the hard coating (1575 MPa). On the other hand, for soft coatings, the maximum equivalent stress mainly distributed in the substrate (sub-surface). For hard coatings, the maximum equivalent stress mainly concentrated on the surface of the coatings and the interface between coatings and substrates, thus increasing the yield probability of these areas. In view of this reason, the subsequent failure analysis of the coating-substrate system was performed with a hard coating and a soft substrate as the focus of this paper. It can also be seen from the nephogram that the equivalent stress of the coated solid was larger at the contact edge (i.e. near the oil film outlet) due to the influence of the secondary pressure peak of the EHL. In addition, Wang et al. [5] showed that the harder coating would cause greater von Mises stress as well as stress differences between the upper and lower interfaces under the EHL point contact, which proved the correctness of the current numerical calculations.
Coatings 2019, 9, x FOR PEER REVIEW 9 of 14 stress distribution nephogram. The smaller the coating elastic modulus, the smaller the stress that will cause certain elastic deformation in the coating. Obviously, under the same external load, the maximum equivalent stress of the soft coating (1029 MPa) was smaller than that of the hard coating (1575 MPa). On the other hand, for soft coatings, the maximum equivalent stress mainly distributed in the substrate (sub-surface). For hard coatings, the maximum equivalent stress mainly concentrated on the surface of the coatings and the interface between coatings and substrates, thus increasing the yield probability of these areas. In view of this reason, the subsequent failure analysis of the coating-substrate system was performed with a hard coating and a soft substrate as the focus of this paper. It can also be seen from the nephogram that the equivalent stress of the coated solid was larger at the contact edge (i.e. near the oil film outlet) due to the influence of the secondary pressure peak of the EHL. In addition, Wang et al. [5] showed that the harder coating would cause greater von Mises stress as well as stress differences between the upper and lower interfaces under the EHL point contact, which proved the correctness of the current numerical calculations.

Failure Prediction
The coating and interfacial stress of ideal non-destructive coating-substrate system were analyzed, and the failure mode and location were predicted. Figure 8 shows the stress distribution along the x-direction, x σ , on the surface of the coating and the side near the interface. The tensile stress in the contact width region was small, mainly presenting as compressive stress. It was found that the maximum shear stress in the coating was much larger than the tensile stress, which distributed near the secondary pressure peak (x = 0.153 mm). Figure 9 shows the distribution of shear stress along the thickness direction in the coating. It appeared that the maximum shear stress was about 500 MPa at the position of 1.25 μm thickness of the coating. It can be seen that the crack of thin coating was mainly caused by the shear stress at the second pressure peak under the heavy-load EHL pressure, which shows the failure of mode II (sliding type).

Failure Prediction
The coating and interfacial stress of ideal non-destructive coating-substrate system were analyzed, and the failure mode and location were predicted. Figure 8 shows the stress distribution along the x-direction, σ x , on the surface of the coating and the side near the interface. The tensile stress in the contact width region was small, mainly presenting as compressive stress. It was found that the maximum shear stress in the coating was much larger than the tensile stress, which distributed near the secondary pressure peak (x = 0.153 mm). Figure 9 shows the distribution of shear stress along the thickness direction in the coating. It appeared that the maximum shear stress was about 500 MPa at the position of 1.25 µm thickness of the coating. It can be seen that the crack of thin coating was mainly caused by the shear stress at the second pressure peak under the heavy-load EHL pressure, which shows the failure of mode II (sliding type). Coatings 2019, 9, x FOR PEER REVIEW 10 of 14   Figure 10 shows the normal and tangential stress distribution curves of the interface under heavy-load EHL pressure. As mentioned above, the influence of compressive stress (negative value) on the system failure can be neglected. Compared with the normal stress at the interface, shear stress plays a dominant role in the delamination failure of the interface. Near the secondary pressure peak, the shear stress can reach 400 MPa, so the region is where the interfacial delamination initiation occurs. Like the coating cracking, it also shows mode II failure.    Figure 10 shows the normal and tangential stress distribution curves of the interface under heavy-load EHL pressure. As mentioned above, the influence of compressive stress (negative value) on the system failure can be neglected. Compared with the normal stress at the interface, shear stress plays a dominant role in the delamination failure of the interface. Near the secondary pressure peak, the shear stress can reach 400 MPa, so the region is where the interfacial delamination initiation occurs. Like the coating cracking, it also shows mode II failure.   Figure 10 shows the normal and tangential stress distribution curves of the interface under heavy-load EHL pressure. As mentioned above, the influence of compressive stress (negative value) on the system failure can be neglected. Compared with the normal stress at the interface, shear stress plays a dominant role in the delamination failure of the interface. Near the secondary pressure peak, the shear stress can reach 400 MPa, so the region is where the interfacial delamination initiation occurs. Like the coating cracking, it also shows mode II failure.   Figure 10 shows the normal and tangential stress distribution curves of the interface under heavy-load EHL pressure. As mentioned above, the influence of compressive stress (negative value) on the system failure can be neglected. Compared with the normal stress at the interface, shear stress plays a dominant role in the delamination failure of the interface. Near the secondary pressure peak, the shear stress can reach 400 MPa, so the region is where the interfacial delamination initiation occurs. Like the coating cracking, it also shows mode II failure.

Damage Evolution
Firstly, the cohesive elements were implanted into the coating-substrate system to simulate the cracks (perpendicular to the coating, coordinate x = 0.153 mm) and interfacial delamination failure. Then, the EHL pressure and shear stress corresponding to hard coating and soft substrate were coupled into the system. Finally, the surface/interface damage evolution behavior was analyzed visually. Table 2 exhibits the main parameters of the CZM for the DLC coating and steel substrate. The failure parameters for coatings and interfaces are usually determined by an indentation test and pull-out tests, etc. According to the maximum failure stress (~500 MPa and~400 MPa) in Figures 9 and 10, to visualize the damage evolution behavior of the system, it was necessary to select the smaller stress values in Table 2. For coatings, it was assumed that 450 MPa (weak crack resistance) was required to be close to 500 MPa in order to accumulate crack damage. For the interface, 350 MPa (weak bonding capacity) was chosen. Table 2. Main parameters of cohesive zone model.
Coating [29,30] Fracture energy (J/m 2 ) 30 Fracture strength (MPa) 1440~1920 Interface [29,31] Adhesion energy (J/m 2 ) 150 Bonding strength (MPa) 285~1940 The calculated initiation and propagation behaviors of coating cracks and the interface are shown in Figure 11, which was evaluated by the cumulative damage stiffness of the CZM. The location of the maximum SDEG value of the coating and interface was observed, which was consistent with the failure prediction of the coating-substrate system by stress analysis in Section 4.3. Coating cracks and interfacial delamination initiated near the secondary pressure peak, and then propagated from the SDEG maximum region to the vicinity.

Damage Evolution
Firstly, the cohesive elements were implanted into the coating-substrate system to simulate the cracks (perpendicular to the coating, coordinate x = 0.153 mm) and interfacial delamination failure. Then, the EHL pressure and shear stress corresponding to hard coating and soft substrate were coupled into the system. Finally, the surface/interface damage evolution behavior was analyzed visually. Table 2 exhibits the main parameters of the CZM for the DLC coating and steel substrate. The failure parameters for coatings and interfaces are usually determined by an indentation test and pull-out tests, etc. According to the maximum failure stress (~500 MPa and ~400 MPa) in Figures 9 and 10, to visualize the damage evolution behavior of the system, it was necessary to select the smaller stress values in Table 2. For coatings, it was assumed that 450 MPa (weak crack resistance) was required to be close to 500 MPa in order to accumulate crack damage. For the interface, 350 MPa (weak bonding capacity) was chosen. Coating [29,30] Fracture energy (J/m 2 ) 30 Fracture strength (MPa) 1440~1920 Interface [29,31] Adhesion energy (J/m 2 ) 150 Bonding strength (MPa) 285~1940 The calculated initiation and propagation behaviors of coating cracks and the interface are shown in Figure 11, which was evaluated by the cumulative damage stiffness of the CZM. The location of the maximum SDEG value of the coating and interface was observed, which was consistent with the failure prediction of the coating-substrate system by stress analysis in Section 4.3. Coating cracks and interfacial delamination initiated near the secondary pressure peak, and then propagated from the SDEG maximum region to the vicinity. On the other hand, comparing Figure 11c with Figure 11a,b, it was found that the maximum cumulative damage stiffness of the coating crack increased (0.843→0.913) under the EHL pressure compared with the coating-substrate system with the ideal bonded interface. Compared with the ideal system without coating crack damage, the maximum cumulative damage stiffness of the system with crack decreases (0.753→0.653). It also means that if the interfacial delamination damage exists in the system beforehand, the probability of cracking in the coating increases. If the coating cracking occurs, the interfacial delamination damage will be weakened. This conclusion is consistent with the relationship between coating cracking and interfacial delamination damage studied by the indentation method in Reference [32]. It can be explained by the energy transfer theory: if the anti-crack ability of the coating is weak, the coating crack will consume more energy under the EHL pressure, which leads to the interface delamination consuming less energy. This means that the interfacial delamination increased the crack failure probability of coating, otherwise, the reverse. It confirms the reliability of the implanted cohesive elements for simulation evaluation of coating cracking and interfacial delamination.

Conclusions
This research investigated the surface/interface damage evolution behavior of a coatingsubstrate system under heavy-load EHL. The main findings from this study are listed below.

•
On the whole, when the load was constant, the elastic modulus of the thin coating (on the steel substrate) had little effect on the oil film pressure distribution. Locally, with the increase of the coating elastic modulus, the oil center pressure increased, and the secondary pressure peak became more significant due to the change in coating stiffness. This conclusion is consistent with the results of the analysis of the uncoupled CZM on interface. • Under the same external load, the maximum equivalent stress of the soft coating system was smaller than that of the hard coating system, which was located on the subsurface of the system, while the maximum equivalent stress of the hard coating system was distributed on the surface of the coating and the interface. This was because of the stiffness difference between the coating and substrate from the point of view of mechanics.

•
When the thin coating and steel gear substrate system was in service under heavy-load EHL, the shear stress in the coating and shear stress on the interface were the main causes of the coating cracking and interfacial delamination, respectively, which appeared near the secondary pressure peak and all show sliding type failure. The reason is that the geometric dimensions of the contact solids varied greatly, which resulted in shear stress in the dangerous section being much larger than the tensile stress.

•
The surface/interface damage evolution behavior of the coating-substrate system was analyzed visually by embedding cohesive zone elements. In light of the energy transfer theory, the interfacial delamination increased the crack failure probability of the coating and vice versa. On the other hand, comparing Figure 11c with Figure 11a,b, it was found that the maximum cumulative damage stiffness of the coating crack increased (0.843→0.913) under the EHL pressure compared with the coating-substrate system with the ideal bonded interface. Compared with the ideal system without coating crack damage, the maximum cumulative damage stiffness of the system with crack decreases (0.753→0.653). It also means that if the interfacial delamination damage exists in the system beforehand, the probability of cracking in the coating increases. If the coating cracking occurs, the interfacial delamination damage will be weakened. This conclusion is consistent with the relationship between coating cracking and interfacial delamination damage studied by the indentation method in Reference [32]. It can be explained by the energy transfer theory: if the anti-crack ability of the coating is weak, the coating crack will consume more energy under the EHL pressure, which leads to the interface delamination consuming less energy. This means that the interfacial delamination increased the crack failure probability of coating, otherwise, the reverse. It confirms the reliability of the implanted cohesive elements for simulation evaluation of coating cracking and interfacial delamination.

Conclusions
This research investigated the surface/interface damage evolution behavior of a coating-substrate system under heavy-load EHL. The main findings from this study are listed below.

•
On the whole, when the load was constant, the elastic modulus of the thin coating (on the steel substrate) had little effect on the oil film pressure distribution. Locally, with the increase of the coating elastic modulus, the oil center pressure increased, and the secondary pressure peak became more significant due to the change in coating stiffness. This conclusion is consistent with the results of the analysis of the uncoupled CZM on interface. • Under the same external load, the maximum equivalent stress of the soft coating system was smaller than that of the hard coating system, which was located on the subsurface of the system, while the maximum equivalent stress of the hard coating system was distributed on the surface of the coating and the interface. This was because of the stiffness difference between the coating and substrate from the point of view of mechanics.

•
When the thin coating and steel gear substrate system was in service under heavy-load EHL, the shear stress in the coating and shear stress on the interface were the main causes of the coating cracking and interfacial delamination, respectively, which appeared near the secondary pressure peak and all show sliding type failure. The reason is that the geometric dimensions of the contact solids varied greatly, which resulted in shear stress in the dangerous section being much larger than the tensile stress.

•
The surface/interface damage evolution behavior of the coating-substrate system was analyzed visually by embedding cohesive zone elements. In light of the energy transfer theory, the interfacial delamination increased the crack failure probability of the coating and vice versa.