Study on contact fatigue of a wind turbine gear pair considering surface roughness

Contact fatigue issues become more and more crucial in gear industry as they significantly affect the reliability and service life of associated mechanical systems such as wind turbine gearboxes. The contact fatigue behavior is mostly determined by the mechanical properties of materials and stress fields near the contact area, which is further influenced by the lubrication and surface roughness due to pressure fluctuations. In this study, a numerical model incorporating the lubrication state, tooth surface roughness, residual stress, and mechanical properties of the material is developed to determine the contact fatigue behavior of a megawatt level wind turbine carburized gear. The variations of the hardness and residual stress along the depth were characterized by the Vickers hardness measurement and X-ray diffraction test, respectively. The elastohydrodynamic lubrication theory was applied to predict the contact pressure distribution, highlighting the influence of the surface roughness that stemed from the original measurement through an optical profiler. The stress histories of the studied material points during a complete contact loading cycle were fast calculated using the discreteconcrete fast Fourier transformation (DC-FFT) method. Modified Dang Van diagrams under different working conditions were determined to estimate the contact fatigue failure risk. The effect of the root mean square (RMS) value of the surface roughness on the failure risk at critical material points were discussed in detail. Results revealed that the surface roughness significantly increases the contact fatigue failure risk within a shallow area, and the maximum risk appears near the surface.


Introduction
Due to the increasing demand for improved performance, reduced weight, and increased service lives of advanced geared machines such as wind turbines, which present significant contact fatigue issues leading to economic repercussions, diverse engineering techniques have been applied such as case carburizing [1], nitriding [2], and coating [3] to improve the resistance of gear contact fatigue failure. Gears failure due to rolling contact fatigue (RCF) occur in a variety of modes such as spalling [4], micropitting [5], and tooth flank fracture [6]. Gear contact fatigue is affected by various factors, including the lubrication state, tooth surface roughness, hardness gradient, and residual stress distribution, which cause the mechanism of gear contact fatigue failure to be unclear.
The impact of the lubrication and surface roughness on the tribological performance and fatigue behavior of contacting elements have been widely investigated in recent years, and significant improvements have been developed [7][8][9]. Qiao et al. [10] compared fatigue model results for elastohydrodynamic lubrication (EHL) with surface roughness. This study reported several multiaxial fatigue criteria such as Findley and Dang Van criteria which were applied to represent the multiaxial stress state and subsequently, to predict the fatigue failure probability. Moreover, according to Liu et al. [11], a thermal starved EHL line contact model for a spur gear pair was derived to study the effect of the working conditions under the starved lubrication, addressing the importance of good lubrication on a gear surface initiated failure resistance. Pu et al. [12] evaluated the contact fatigue performance under combined rolling and sliding motions by considering the sinusoidal roughness. The results imply that sliding could reduce the contact fatigue life dramatically, and rough surface contacts would lead to accelerated pitting failures. Yan et al. [13] numerically analyzed the lubrication characteristics and fatigue life under full film and surface topography. Zhu et al. [14] studied the effect of rough surface topography on mixed lubrication characteristics. Results reveal that under the mixed EHL state, maximum pressure, contact load ratio, friction, flash temperature, and maximum subsurface shear stress increased dramatically as the roughness increased. Zhou et al. [15] discussed the relationship between microstructure characterization parameters and contact characteristic parameters of two cylinders with surface roughness. Due to the development of the industrial manufacturing technologies, the root mean square (RMS) of the gear surface roughness could be controlled at a relatively low value. For instance, the RMS of a wind turbine gear can be controlled at approximately 0.50 μm by applying manufacturing processes. Although the gear lubrication is sometimes guaranteed, micropitting can still be observed due to the significant influence of the surface roughness. Sheng and Kahraman [16] proposed a numerical model considering the lubrication and surface topography to predict the gear micropitting life, while the residual stress, the material properties, etc., were not considered. Cardoso et al. [17] estimated the gear micropitting performances of high pressure nitriding steel gears which were lubricated with a standard mineral lubricant and two biodegradable esters. Experimental results showed that the gears with the ester lubricants presented better micropitting performance than the gears with the mineral oil. Despite numerous theoretical and experimental investigations, the micropitting mechanism has still not been fully determined and requires further study. As revealed by the experimental observations in Ref. [18], the fatigue crack mostly initiates within the contact stressed zone, where the mechanical properties and residual stress highly vary. This implies that the fatigue damage is sensitive to those variations [19]. In addition, a large amount of residual stresses can also be introduced through the machining, heat treatment, or surface treatment processes, which need to be considered. Walvekar et al. [20] studied the RCF lives of case carburized steels and calculated the damage parameters. The effect of the residual stress distribution on the RCF lives was also investigated. Rudenko and Val'ko [21] studied the contact fatigue resistance of carburized gears experimentally, and an expression for the ultimate deep contact fatigue resistance of carburized gears was derived. Lv et al. [22] studied the effect of the microshot peened treatment on the fatigue behavior of gears as well as the effects of the relaxation mechanism of the compressive residual stress on the gear fatigue performance. Inoue and Kato [23] proposed an experimental formula for the estimation of the bending fatigue strength of carburized gears. This formula also presents a way of enhancing gear fatigue strength by increasing the surface hardness and residual stress. MackAldener and Olsson [24] analyzed a gear failure mode called tooth interior fatigue fracture by numerical simulations using Finite element model (FEM) and the critical plane fatigue initiation criterion according to Findley, highlighting the effect of the residual stress. Brandão et al. [25] proposed a micropitting test performed with spur gears considering the residual stress on a Forschungsstelle für Zahnräder und Getriebebau (FZG) test machine to predict the surface damage.
However, these studies do not simultaneously consider the hardness gradient, residual stress gradient, lubrication state, and surface roughness in the prediction of the gear contact fatigue. A more comprehensive numerical model needs to be developed to include complicated contact situations. Herein, a numerical model incorporating the variations of the hardness and the residual stress is proposed to describe the stress responses and to estimate the contact fatigue performance of a carburized gear under the lubrication condition. Based on the numerical results, the influence of the RMS value of the surface roughness on the gear contact fatigue is discussed in detail.

Simulation methodology
A numerical model for the prediction of the contact fatigue performance of a carburized gear is described in detail as follows. Under a given loading and lubrication condition and different surface topographies, the gear contact fatigue risk can be estimated based on the fast-calculated stress histories and modified Dang Van diagram. Moreover, the material properties and residual stress distribution are considered for a more comprehensive analysis. The technical process of this study is displayed in Fig. 1. The model primarily consists of four aspects: gear contact parameters, lubrication model, Dang Van diagram, and explanation of failure risk index.

Gear contact parameters
The gear sample studied in this work arises from the intermediate parallel stage of a 2 MW wind turbine gearbox. Such a gear stress concentration becomes more pronounced when the gear engages at the pitch point compared with the other meshing positions [26], thus the contact fatigue failure is primarily evaluated at this critical position. The gear contact at the pitch point can be equivalent to two deformable circles with different radii of curvature ( 1 2 , R R ) contacting  with each other. Subsequently, this model can be further simplified as a rigid circle with an equivalent radius of curvature ( e R ) contacting with an infinite elastic half-space, with an equivalent Young's modulus ( e E ) [27,28]. The schematic diagram of the gear contact model based on the plane strain assumption is illustrated in Fig. 2, as well as the information on the coordinate system. The x coordinate denotes the rolling direction, the z coordinate represents the direction of the depth, and the y coordinate represents for the gear width direction. The model geometric as well as the operating parameters are listed in Table 1. The gear manufacturing process is also illustrated in Table 1 due to its crucial role on the final fatigue performance of gears. In fact, the heating treatment process was performed at a carburizing temperature of 930 °C , intensive carburization potential (26-30 h) of 1.03%-1.13%, diffusion potential (14-24 h) of The tooth surface roughness, due to its notable influence on the pressure fluctuation, is considered in the proposed model. The deterministic nature of the model allows the explicit description of the surface roughness. The gear surface roughness, processed by the generating grinding method, is measured using the optical device Alicona G4 equipped with highresolution lenses. Figure 3 shows the gear surface topography and processed data used in the calculation. The RMS value of the measured roughness is 0.25 μm, and the gap along the rolling direction between two adjacent points is 0.18 μm.

Lubrication model
Following the developments of the EHL and micro-EHL theories, gear lubrication has been extensively studied in recent years and relevant numerical models have been developed [29,30]. For clarity, essential equations are explained here. For the plane strain assumption, the fluid flow in the contact area is expressed by a one-dimensional generalized Reynolds equation [31]: where p , e h , and  represent the fluid pressure, film thickness, and lubricant density, respectively. The rolling velocity is represented by is the current rolling velocity of the driving and driven gears, respectively. The variation of the lubricant density depends on the pressure which can be derived through the Dowson-Higginson equation [32]: where  0 is the ambient oil density. An effective viscosity  * is introduced to characterize the non-Newtonian behavior, which is expressed as the Ree-Eyring fluid. The equivalent viscosity of the Ree-Eyring fluid is expressed as follows: where  0 is the Eyring characteristic stress and is set to 10 MPa in the present work.  is the current oil shear stress, and  is the viscosity as a function of the pressure [33]: where 0 p is the pressure index   9 0 = 5.1 10 p ,  0 is the oil viscosity at ambient pressure, and e z is fixed as a constant of 0.60.
The oil film thickness h, considering the term of the surface roughness, can be calculated as follows: where 0 h is the initial separation between two interacting surfaces, R is the equivalent radius of curvature, ED is the elastic deformation, and SR denotes the term of the surface roughness. Moreover, the lubrication parameters are listed in Table 2. The previous equations are dimensionless based upon the Hertzian parameters, namely the half-Hertzian contact width H b and maximum Hertzian pressure H p [27]. The discrete convolute, i.e., the fast Fourier transformation (DC-FFT) method, first proposed by Liu et al. [34], is adopted for the fast calculation of the elastic deformation and stress components. Herein, the friction coefficient is assumed to be 0.10. A more detailed explanation on the solving strategy is detailed in Ref. [11].
As shown in Fig. 4, the calculation domain of the material points covers a range of    x z N N equally spaced grids. The time of the gear contact should be guaranteed to keep each individual point completing an entire loading cycle and is estimated as follows [35]:  where H b is 0.86 mm, and r u is 2.04 m/s at the pitch point in the calculation. Therefore, the gear contact time can simply be calculated as 0.0017 s. The distance between two adjacent measured data points is approximately 0.18 μm along the rolling direction. Compared with the distance between two adjacent interested material points within the calculation domain (approximately 54 μm), the amount of the measured surface roughness data points is large enough to obtain high accuracy and efficiency.

Dang Van diagram
According to the complicated multiaxial stress states during the gear contact period, various multiaxial fatigue criteria have been applied in gear contact fatigue studies. Among them, the Dang Van multiaxial fatigue criterion is widely applied to evaluate the contact fatigue performance [36,37]. The significant hydrostatic stress distribution caused by cyclic contact, which is absent in classical fatigue situations, could be explicitly considered in this criterion. The Dang Van equivalent stress  DangVan is defined as a linear function of the maximum shear stress  max and the hydrostatic stress  H : where  D is a material parameter which could be evaluated according to the following equation [38]: where  1 and  1 are the fully reversed bending fatigue limit and fully reversed torsion fatigue limit of the material points, respectively, which depend on the depth location. Figure 5 presents the original Dang Van diagram representing the contact fatigue limit [39]. In this diagram, the x coordinate denotes the hydrostatic stress, and the y coordinate represents the maximum shear stress  max . The original Dang Van curve is characterized as an inclined line with a negative slope  D . For any load path data on or above this Dang Van curve, the fatigue crack initiation occurrence is expected [39,40], and thus the corresponding material points should be highly noticed.
The original Dang Van criterion predicts a detrimental effect of tensile hydrostatic, and an overestimated positive influence is expected from compressive values [41]. Desimone et al. [40] reported an over-optimistic prediction of the effect of the compressive residual stress on the RCF issues, especially when the compressive hydrostatic stress exists during contact. Therefore, when the value of the hydrostatic stress is negative (the negative part of the x coordinate), the original inclined Dang Van curve should be modified by evaluating the RCF failure.
The residual stress plays a non-negligible role, together with the applied stress, on the formation of the final stress field. The residual stress gradient of the rolling direction was characterized by applying a X-ray diffraction test [26]. The electro-polishing is used for the measurement of the residual stress gradient distribution. The profile is assumed to remain negative along the direction of depth, representing the compressive residual stress state. Several groups of gear samples are tested, and based on the measured data and the deviation, the residual stress is comprised in a range of approximately 40 MPa, which may be caused by the stress relaxation or the uneven measured tooth surface. As the residual stress along the z coordinate is much smaller than that in the other two directions, only the residual stresses along the rolling direction  r, x and width direction  r, y are superposed to the corresponding stress components  x and  y linearly for the further calculation, and  r, x ,  r, y profiles along the depth are assumed to have same profiles along the depth [42].
Considering residual stress into account, the Dang Van equivalent stress could be formulated as follows [43]: where the    D H , l o a d is modified by the superimposition of the residual stress distribution through the hydrostatic part of the residual stress tensor  H, residual [41]. Hence a different Dang Van safe locus is proposed and can be identified in two straight line segments. One is horizontal, representing a conservative fatigue limit for the region of compressive hydrostatic stress values, and the other presents a negative slope  D , as illustrated in Fig. 6 [41]. The modified multiaxial fatigue limit does not depend on the variation of the amplitude of the compressive hydrostatic stress, indicating that the over-optimistic positive effect of the compressive hydrostatic stress on the contact fatigue problems is modified. Moreover, based on the RCF experiment in Ref. [44], the conservative horizontal locus shows a more accurate fit than the original Dang Van curve for the compressive hydrostatic stress.

Index of failure risk
According to the characteristics of heat treatment and carburizing, the hardness value decreases gradually from the surface to the core substrate within the hardened layer of a carburized gear. The hardness gradient variations would cause changes in the local material strength, fully-reversed bending limit, and torsion fatigue limit. Therefore, the gradient properties must be considered when analyzing carburized gear contact fatigue issues. Figure 7 shows the measured hardness data using the Vickers indentation test machine, and the measured data is curve-fitted by the third order Fourier series analysis. In order to avoid the possible excessive indentation, a test force of 4.9 N is chosen. Each measurement point along the tooth flank is measured five times and finally, the average hardness value is determined. As observed in Fig. 7, the hardness varies from 650 to 450 HV, from the case to the core. Moreover, the effective case depth (ECD), defined as the hardness value is 550 HV, is approximately 2.20 mm. According to the measurement, the deviation of the measured hardness data is within 10 HV.
The values of the fully reversed bending fatigue limit  1 and fully reversed torsion fatigue limit  1 used in the Dang Van criterion depend on the material depth and can be determined based on the hardness  | https://mc03.manuscriptcentral.com/friction gradient by using the following equations [45]: where z represents the depth from the surface. Figure 8 shows the profiles of hardness, fatigue limits under fully reversed bending (  1 ), and torsion ( 1 ) along the direction of depth. The value of  1 varies from 0.64 GPa on the surface to 0.49 GPa in the core.
In addition,  1 has a similar profile as  1 and varies from 0.37 to 0.28 GPa. The limit gradually decreases with the increase of the depth, and the load path of each material point should be compared with its corresponding limit in order to capture the severity of each instant in the contact period.
For the evaluation of contact fatigue of carburized gears, a concept called "the material exposure" (a ratio of the equivalent shear stress considering residual stress and local material strength derived from the hardness gradient) is described in Refs. [46,47]. Accordingly, an index R is defined to estimate the risk of the gear contact fatigue along the depth as follows: (11) where R is the risk of gear contact fatigue,  Dang Van is the Dang Van equivalent stress,  1 is the torsion fatigue limit mentioned previously, and z represents the depth. Figure 9 shows the schematic diagram of the Dang Van equivalent stress  Dang Van with and without the measured residual stress and the resultant distribution of the fatigue failure risk. It can be observed that the equivalent stress is affected by the residual stress, which can further influence the fatigue failure risk.

Dang Van diagram with smooth surface
The Dang Van equivalent stress history and the modified Dang Van diagrams can be evaluated to demonstrate the effect of the measured compressive residual stress on the stress distribution and gear contact fatigue performance. Figure 10 illustrates  stress and the y coordinate represents the maximum shear stress  max . The torsion fatigue limit line  1 at the gear surface and core, and the corresponding limit line of the maximum value of  max are illustrated.
The Dang Van equation is expressed as a line with a slope  D based on the depth mentioned previously.
In this section, the data of 401 material points along a depth down to H 4.0b at the middle of the contact length is studied.
Overall, the load path data of each material point does not reach its corresponding torsion fatigue limit, when the smooth surface is considered. This indicates that almost no contact fatigue issue is expected to occur under this circumstance. Considering the superposition of the residual stress distribution, the maximum shear stress  max at the subsurface decreases, and the location of the maximum becomes closer to the surface. However, the maximum value of  max almost remains unchanged. Such compressive residual stress is beneficial to the drop of the gear contact fatigue failure possibilities at the subsurface.  than its corresponding limit value  1 . None of the data point at the surface exceeds the corresponding limit line, but the variation of the hydrostatic stress  H distribution at the surface fluctuates within a wide range because of the fluctuation of the surface pressure due to the roughness. It can be summarized that the influence of the surface roughness is more significant as the RMS increases, especially at the near surface where the stress distributions highly fluctuate. The near surface could be regarded as the expected crack initiation positions rather than the surface, although the important fluctuation of  H could be found at the surface. However, the effects of the surface roughness on the stress distribution drop sharply as the depth becomes deeper.

Effect of surface roughness on the risk of contact fatigue
In order to investigate the gear contact failure probability influence range of the surface roughness, the risk of the gear contact fatigue along the depth with different surface asperities is evaluated according to Section 2.4. Figure 13 schematically illustrates the risk of the contact fatigue along the depth for a maximum value of  max during the contact period with different asperity conditions. When the gear surface is smooth, the maximum risk of fatigue is approximately 0.8, which appears at a depth of 0.1 mm. The maximum risk of fatigue is approximately 1.0, 1.2, and 1.4 when the RMS of the surface roughness is 0.2, 0.3, and 0.4 μm, respectively. As the RMS increases to 0.5 μm, the maximum value of the fatigue risk reaches approximately 1.6, and its occurrence location is H 0.01b (0.0086 mm in this work), which is very close to the surface.
Moreover, the range of the influence of the surface roughness is close to the surface. As it can be observed in Fig. 14, the risk of fatigue is highly affected by the surface roughness within a small range. When the RMS is 0.5 μm, the range of obvious influence is approximately 0.05 mm from the surface. Although discrepancies in the risk of fatigue could be found on a deeper location compared with a smooth surface, the effect can be ignored when describing a gear contact fatigue as the maximum risk of fatigue occurs at a very shallow location below the gear tooth surface, which is the area dominated by gear contact fatigue modes such as micropitting.
The maximum fatigue failure risk appears at a depth of approximately H 0.01b , which is more than twice the value of that at the surface. This is mainly attributed to the interaction of the peaks of the  | https://mc03.manuscriptcentral.com/friction roughness being considered as an extremely small contact area which is subjected to the Hertzian contact theory. Consequently, the maximum fatigue risk occurs near the surface rather than in the surface in spite of the effect of surface roughness.
Consistent results are displayed in Fig. 15 in which 26 groups of RMS of surface roughness are chosen from 0 to 0.5 μm in order to investigate the effect of the surface roughness on the risk of fatigue and to generate the failure risk map. The highest value of the risk of contact fatigue is approximately 1.6 when the RMS is 0.5 μm . More specifically, the maximum fatigue risk appears at a closer position to the surface and its value exceeds 1.0 when the roughness RMS increases to 0.28 μm . The maximum risk increases significantly, but the influence of the surface roughness remain within a small range which is fairly close to the gear surface. Consequently, the micropitting is more likely to occur as the RMS value increases during the gear contact process.

Conclusions
Herein, a numerical model was proposed to study the effect of the gear tooth surface roughness on the contact fatigue behavior of a carburized wind turbine gear pair under lubrication condition. The variations of the hardness and the residual stress were both incorporated in the numerical model. The modified Dang Van diagram and the risk of contact fatigue at critical material points were performed for different roughness cases. The main conclusions are summarized as follows: (1) No contact fatigue failure is expected to appear with the smooth surface assumption. When the compressive residual stress is considered, a decrease of the maximum shear stress at the subsurface can be observed. The location of the maximum is found near the surface, increasing the failure risk near the surface. However, the change of the peak value of the maximum shear stress cannot be observed clearly. The compressive residual stress is beneficial to the resistant ability of the gear contact fatigue failures which occur at the subsurface.
(2) The effect of the surface roughness on the variations of the load path data is more important as the value of RMS increases. More important variations of  H at the surface are observed because of the greater fluctuations of the pressure distribution caused by the increasing surface roughness, and all investigated load path data at the depth of H 0.01b are above the modified Dang Van limit. The influence of the surface roughness on the stress distributions sharply decreases as the position becomes deeper.
(3) The increase of the surface roughness RMS increases the maximum risk of fatigue significantly, which occurs in a shallow range. As the roughness RMS increases to 0.5 μm , the effect of the roughness RMS on the failure risk becomes limited within a depth of 0.05 mm. Therefore, the possibility of the appearance of the micropitting increases during the gear meshing period. Moreover, fatigue experiments based on the gear samples or other equivalent components will be conducted in future work to determine several important materials and tribological parameters. Subsequently, the current numerical results could be modified and verified based on the experimental data.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.
The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
To view a copy of this licence, visit http://creativecommons.