Tribological behavior of coated spur gear pairs with tooth surface roughness

Coating is an effective way to reduce friction and wear and to improve the contact-fatigue lives of gear components, which further guarantees a longer service life and better reliability of industrial machinery. The fact that the influence coefficient linking the tractions and stress components could not be expressed explicitly increases the difficulty of coated solids contact analysis. The complicated tribological behavior between tooth surfaces influenced by lubrication and surface roughness further adds difficulty to the coated gear pair contact problems. A numerical elastohydrodynamic lubricated (EHL) contact model of a coated gear pair is proposed by considering the coupled effects of gear kinematics, coating properties, lubrication, and surface roughness. The frequency response function and the discrete convolute, fast Fourier transformation (DC-FFT) method are combined to calculate the surface deformation and the subsurface stress fields at each meshing position along the line of action (LOA). The Ree-Eyring fluid is assumed to incorporate the non-Newtonian effect, which is represented in the generalized Reynolds equation. Influences of the ratio between the Young’s modulus of the coating and the substrate on the contact performance, such as pressure, film thickness, tooth friction coefficient, and subsurface stress field, are studied. The effect of the root mean square (RMS) value of the tooth surface roughness is studied by introducing the roughness data, deterministically measured by an optical profiler.


Introduction
The development of the machinery industry is characterized by high speed, high load, and high power trends. Hence, higher requirements have been proposed for the power density and service life of key mechanical components, such as gears and bearings. Mechanical equipment accidents caused by gear failures [1] occur from time to time in fields of wind turbines, helicopters, and warships, etc. In order to get higher power densities and fatigue lives of gears, advanced technologies, such as case hardening [2] and shot peening [3], have been developed to augment the service lives or reduce the friction and wear. Among these technologies, coating is an effective method that has been widely used in many industrial fields [4,5].
Owing to the difference between the mechanical properties of coating and substrate, both the surface tribological behavior and subsurface stress field would change compared with the uncoated case. Furthermore, the fatigue strength would be affected by the application of the coating. Numerically, the elastic deformation or stress components can be calculated by the pressure (or shear traction) distribution and the influence coefficients linked to them. However, in the case where coating exists, the influence coefficients that link the normal pressure (or the shear traction) to the elastic deformation (or the stress components) cannot be expressed explicitly [6]. Hence, the effects of the coating mechanical properties on contact behavior are rarely studied. An exception is that the finite element method is applied in coated contact problems [7,8]. Recently, Wang et al. developed a model for elastohydrodynamic lubrication of multilayered materials under point contact cases, in which the frequency response functions, relating pressure to surface displacements and stress components, are derived from the Papkovich-Neuber potentials [9]. Then, they extended this method for solving three-dimensional fretting contact problems involving functionally graded materials [10]. Liu et al. considered both the thermal and the mechanical properties of the coating on a line contact [11]. The contact analysis of coated gear pairs would be even more difficult when considering the complicated gear kinematics and tooth surface roughness. Liu et al. [12,13], Li et al. [14,15], and Evans et al. [16,17] have studied gear tribology problems without the consideration of coating. Liu et al. [18] provided a discrete convolute, fast Fourier transform (DC-FFT) method for fast calculation of the surface deformation and subsurface stress components. A numerical elastohydrodynamic lubrication model based on the frequency response function and the DC-FFT method is presented in this paper to evaluate the distribution of pressure, film thickness, friction, subsurface stress field, etc., during the meshing process considering gear coating properties. The tooth surface roughness is involved deterministically based upon the measured tooth surface roughness data. Effects of the modulus ratio between the coating and the substrate and the root mean square (RMS) value of the surface roughness are studied.

Model parameters
The sample used in the study is the intermediate parallel stage of a megawatt-level wind turbine gearbox, as showed in Fig. 1 For simplification, the helix angle is not taken into account, and hence, the plane strain condition can be assumed to simulate gear contacts. The contact at any meshing moment is simplified as two circles with different radius of curvature contacting with each other. The variations of the radius, rolling and sliding velocity, and normal load during the meshing process could be calculated based on the gear meshing theory. It is necessary to determine the parameters of the mechanical properties of coating, such as the Young's modulus and Poisson's ratio, when studying the coated gear pair. The diagram of the coated gear contact is shown in Fig. 2, in which the coordinates are indicated. In this paper, the pinion and the wheel are assumed to have the same coating thickness and mechanical properties.
The gear parameters and lubrication properties are listed in Table 1. The Ree-Eyring fluid is assumed to represent the non-Newtonian behavior by employing the generalized Reynolds equation proposed by Yang and Wen [19].
A quasi-steady tooth normal load distribution is assumed, as follows. Around the pitch point, the meshing point belongs to the single teeth-meshing zone, which means that the load is carried by the studied tooth pair, while at the double tooth-meshing zone, the tooth load is shared between two adjacent pairs. At the beginning of the engagement from the engage-in point, the tooth normal load increases linearly from 1/3 of the total load at the engage-in point to 2/3 of the total load at the lowest point of single tooth contact along the line of action (LPSTC). Once the contact passes the highest point of single tooth contact along the line of action (HPSTC), the normal tooth load decreases linearly from 2/3 of the total load to 1/3 of the total load at the recess point. Therefore, sudden changes of tooth normal load occur at the HPSTC and LPSTC points.
The governing equations mainly consist of the Reynolds equation, film thickness equation, force balance equation, and viscosity-pressure equation of the oil film. The dimensionless parameters of this model are as follows: The generalized Reynolds equation, which does not consider the transient squeeze effect, could be expressed as [19]: The meaning of parameters in the Reynolds equation could refer to Ref. [20]. The load balance equation, film viscosity-pressure equation, discretization scheme, and iteration method could be found in Ref. [21]. In this work there are 513 equally-spaced points For the coated elastic contact problem, an explicit Green's function is not available for the normal displacement in the space domain. Its frequency response function, however, is available in the frequency domain [22]. Previous studies have derived the 2D frequency response function of the Green function for a single body of line contact as [22]: For the line contact plane-strain problem, the body forces are assumed to be negligible, and the stress components can be written as where  denotes the Airy stress function. x, y, and z represent the rolling direction, tooth width direction, and direction of depth, respectively. By introducing the Fourier transform x , the stresses and surface displacements can be given [24]. Detailed derivation of the frequency response function of the stress components can be found in Ref. [25]. Once the frequency response functions, the elastic deformation and stresses can be achieved by using the DC-FFT algorithm [18]. To avoid the singularity of the frequency response functions in the origin point, a 16-point Gaussian quadrature integration is applied to evaluate the response function at this point.
From the gear parameters showed in Table 1, the equivalent radius, rolling and sliding velocity, and normal load can be obtained at each engaging point along the line of action (LOA) through the gear meshing theory. The calculation of those parameters can be found in Ref. [21]. Figure 3 shows variations within a whole meshing period of the rolling velocity r u , sliding velocity s u , equivalent radius R , and Hertzian parameters ( h p and h b ) during the meshing process under the given working condition. It can be seen that, during the meshing process, the rolling velocity decreases almost linearly, while the sliding velocity first decreases to zero at the pitch point and then increases gradually, which means that the slide-to-roll ratio at the pitch point is zero. The equivalent radius decreases during the meshing process remarkably. It should be noted that the equivalent radius of a gear pair affects the location of the dangerous zone of the material, thus determining which specific mode of contact failures, such as macro-or micro-pitting, occurs first. The sudden changes of Hertzian parameters at the LPSTC and HPSTC points are observed. Figure 4 shows the distribution of the pressure and oil film shear stress within an entire meshing period of the uncoated gear pair. It can be seen that the contact area changes abruptly at the LPSTC and HPSTC points, the contact pressure is higher at the single teeth-meshing area (near the pitch point), and the maximum pressure occurs when the meshing reaches almost 0.8 GPa, around the pitch point engaging position. Figure 5 shows the comparison of the pressure distribution at the pitch point meshing position, without coating, between the lubrication and dry contact results. As can be seen, under such a working condition, the pressure profile of the lubrication condition only slightly differs from the pressure profile of the dry contact case at the inlet and outlet zones. The numerical curve represents typical characteristics of the inlet zone and the second pressure spike. Figure 6 shows the comparison of the von-Mises stress field at the pitch point meshing position, without coating, between the lubrication and dry contact results. The maximum values of the von-Mises stress under the lubrication and dry contact cases are identical. Only a slight difference occurs at the outlet zone, which is caused by the second pressure spike of the lubrication case.    Fig. 7. They also show a close agreement. The error between the maximum Mises  of the two cases is around 9%, which is acceptable. It is worth to note that the numerical result (left) shows stress concentration at the outlet zone which is caused by the second pressure spike.

Effects of mechanical properties of coatings
The effects of different coating Young's modulus on the gear lubrication contact condition are studied. Figure 8 shows the comparison of the effects of three different coating Young's modulus on the contact pressure when the coating thickness is 100 μm, constantly. It can be seen that the contact width of the soft coating (left) is larger than the contact width of the hard coating (right), and the amplitude of | https://mc03.manuscriptcentral.com/friction pressure in the contact area of hard coating is the largest among these three conditions. The pressure distribution with coatings is different from the Hertzian contact pressure distribution. This issue should be considered in the subsequent calculation of stress field. Figure 9 shows the distribution of subsurface orthogonal shear stress at the pitch point when using three different kinds of coatings. It can be seen that different coating modulus ratios have small impacts on the distribution or the amplitude of the shear stress.
However, the distribution of Mises  shows that the influence of different modulus ratios cannot be neglected. As can be seen in Fig. 10  the big stress appears not only at the depth of approximately 0.4 mm, but also occurs in the coating at the subsurface, especially near the recess point. The position of maximum stress is very important to the specific gear contact fatigue failure modes because a different position of maximum stress may decide which fatigue failure, such as micro-pitting, pitting, or tooth interior fatigue fracture, will occur first.
In the case of relative sliding, the gear tooth surface is subject to the shear traction and normal traction, which results from the action of the oil film shear. Figure 11 shows the oil film shear stress distribution under three different modulus ratio conditions. The oil film shear stress is almost zero at the pitch point, because the sliding velocity is zero at the same position. The sliding velocities on both sides of the pitch point have opposite directions. Hence, the corresponding shear stresses also have opposite directions. The position of maximum shear stress is approximately at the recess point. Because of the high pressure gradient and the high viscosity caused by high pressure in the hard coating, as the Young's modulus of the oil film increases, the amplitude of oil film shear stress increases too.
The condition of this study is based on EHL, and thus, direct contact between asperities does not happen, and the integral of the oil film shear stress is applied to calculate the tooth surface friction. Figure 12 shows how the friction coefficient changes with the Young's modulus of coatings. The mean friction coefficient is defined as: where L is the distance between the engage-in and recess points. This definition is practical and useful in engineering design;  represents the current coefficient of friction: where 1  and 2  are the shear stresses of the film at the interfaces between the oil film and the two solids, respectively. F is the current normal load (N/m), xin and xout are the inlet and the outlet boundaries, respectively. It is noticed that a full elastohydrodynamic lubrication region is assumed and no dry asperity contact occurs. As can be seen in Fig. 12, the mean friction coefficient or the friction coefficient at engage-in point, pitch point, recess point, or any other point along the LOA increases with the Young's modulus of coating, and this kind of change is unaffected by the rollingsliding condition. This means that the soft coating is better when considering the reduction of friction. The effect of friction on the stress field is not significant because the friction coefficient is not larger than 0.8 at any meshing position under the given working condition. Figure 13 shows how the minimum film thickness changes with the Young's modulus of coating at three  typical meshing positions. As can be seen, the minimum film thickness slightly increases with the Young's modulus of coating because of the effect of pressure and viscosity based on the governing equations. Taking the pitch point as an example, the modulus ratio increases from 0.25 to 4, which covers the domain of actual engineering practice. The minimum film thickness rises by only 18% from 0.5 μm to 0.59 μm.
The change in minimum film thickness can be ignored considering the small coating thickness in industrial machinery.

Effect of surface roughness with coatings
The tooth surface has roughness, which is on the same order of magnitude as the oil film thickness in the case of machining or other operations. Roughness has significant influence on the surface tribological behavior. In order to study the effect of surface roughness on the contact properties of the coated gear, a sample surface roughness, taken from a forminggrinding gear tooth surface, is obtained through the Alicona optical profiler, and it is studied instead of the coated gear roughness. As shown in Fig. 14, at one fixed section along the tooth width, data is collected from the tooth tip to the root and there are in total 10,835 measuring points within a 3.5 mm profile, which is enough for the deterministic study. The original root mean square of the surface roughness is q 0.5 μm R  . In order to study the effect of roughness RMS, original roughness data multiplied by a specific scaling factor are taken to obtain a group of roughness data with different RMS.
The surface roughness would lead to marked fluctuations of the surface pressure distribution. Hence, the local pressure may be several times the maximum Hertzian pressure. Figure 15 shows the distribution of contact pressure under three different modulus ratio conditions when roughness RMS is 0.25 μm. Compared with Fig. 8, the local pressure with roughness is much bigger than the local pressure on a smooth surface, and its distribution becomes more complex. When  c s / 2 E E (right), the maximum contact pressure achieves 1.2 GPa, whereas the maximum pressure on a smooth surface is about 0.8 GPa. It means that if the tooth surface gets rougher, the local pressure would be much larger.
The dramatic variation of pressure causes fluctuations of the subsurface stress field. Figure 16 shows the distribution of orthogonal shear stress at the pitch point according to different roughness RMS with coatings during the meshing period. As can be seen, a significant fluctuation of the subsurface orthogonal shear stress occurs as the roughness RMS increases.
When considering roughness, the maximum orthogonal shear stress appears at two positions. One of them coincides with the position of maximum Hertzian pressure, which is approximately at the depth of 0.4 mm; the other is in the near surface. Compared with Fig. 9, a conclusion can be made that the effect of roughness on orthogonal shear stress is much greater than the effect of coatings. Figure 17 shows the effect of roughness RMS on    is about 700 MPa when the roughness RMS takes 0.5 μm, which is 40% bigger than the former in amplitude. Furthermore, the position of the maximum Mises  is closer to the surface, which makes the gear contact fatigue life more sensitive to roughness. Superfinishing technologies have been proven as effective for the improvement of the gear contact fatigue life. Figure 18 shows the effects of the roughness RMS on the subsurface maximum Mises  at the pitch and engage-in points according to two different coating thickness (100 μm and 200 μm). As can be seen, the maximum Mises  at these two meshing positions increases as the roughness RMS increases, and the rise in the maximum Mises  at the engage-in point is much larger than that at the pitch point. However, the effect of coating thickness on the maximum Mises  is not obvious. When the roughness RMS is 0.5 μm, the peak Mises  at the engage-in point reaches almost 1 GPa, which surpasses the yield strength of most steels. Therefore, future works on the gear rough surface contacts should consider the elasto-plastic behavior of materials.

Conclusions
A numerical lubricated-contact model of a gear pair with surface roughness and coatings is developed. Effects of the Young's modulus of the coating and the RMS value of the surface roughness are studied. The following conclusions are drawn: (1) Compared with the effect of pressure, the effect of coating mechanical properties on the orthogonal shearing stress distribution is not notable, but its effect on the Mises  distribution is remarkable. In addition, the existence of coating leads to a discontinuity of the distribution of Mises  at the coating/substrate interface.  (2) Owing to the high-pressure gradient and the high oil viscosity under hard coating conditions, the amplitude of the oil film shear stress and the friction coefficient increase as the Young's modulus of coating increases.
(3) Surface roughness causes a dramatic fluctuation in surface pressure and the subsurface stress field. The effect of roughness on the amplitude and distribution of subsurface stress is more significant than the effect of coatings.