Evaluation of contact fatigue risk of a carburized gear considering gradients of mechanical properties

Carburized gears are widely used in geared machines such as wind turbines. Contact fatigue problems occur in engineering practice, reducing reliabilities of machines. Contact fatigue failures are related to many factors, such as gradients of mechanical properties of the hardening layer. In this work, an elastic-plastic contact model of a carburized gear is developed based on the finite element method to evaluate contact fatigue failure risk, considering variations in hardness and strength. The Dang Van multiaxial equivalent stress is calculated via Python coding within the Abaqus framework. The gradient of yield strength along the depth from case to core is considered. The concept of local material fatigue failure risk is defined to evaluate the probability of pitting failure. The effects of design factors, such as the case hardening depth (CHD), surface hardness, and contact pressure on fatigue failure risk, are studied. As the CHD increases or the surface hardness decreases, the risk of deep spalling failure reduces. The increase in surface hardness leads to a decreased risk of pitting failure, while the variation in CHD hardly affects the pitting failure risk.


Introduction
Wind turbines are widely used in power generation due to the ever-rising demand for renewable energy. As the requirements of higher power and power density from wind turbines grow, rolling contact fatigue problems in wind turbine gearboxes have become one of vital problems limiting turbine reliability [1]. In order to cope with the harsh operating environment, surface hardening, especially carburization [2], is widely applied in megawatt wind turbine gears. Contact fatigue failures relate to many factors, such as working conditions, mechanical properties, and residual stress. In addition, during gear engagement, the contact undergoes non-proportional multiaxial stress, which means both normal and shear stress components contribute to the gear contact fatigue process. In engineering practice, several failure modes are associated with contact fatigue, such as micropitting [3], pitting [4], and deep spalling [5], as shown in Fig. 1. Furthermore, wind turbines may encounter overloading and other extreme working conditions (such as misalignment-induced stress concentration) during operation, leading to elastic-plastic contact behavior. These factors indicate that the contact fatigue failure mechanism of heavy-duty carburized gears is complex. Therefore, this work focuses on the contact fatigue problem of large-module carburized gears from two main perspectives: the determination of multiaxial stress state, and the gradient characteristics of mechanical properties caused by carburizing and quenching.
The influencing factors of gear contact fatigue behavior have been widely studied. Li et al. [6] and Liu et al. [7] carried out gear tribology studies considering the lubrication state. Evans et al. [8] carried out gear contact fatigue prediction considering tooth surface roughness. Hu et al. [9] considered the effect of dynamic load on gear pitting failure life. Zhou et al. [10] developed an elastoplasticity asperity contact model to investigate gear contact. However, these studies were limited to homogeneous materials, and did not consider the variation in mechanical properties caused by the case hardening process. Some studies have documented the effect of case hardening process on the gradients of mechanical properties of gears along the depth.
Wang et al. [11] predicted the rolling contact fatigue (RCF) life of a carburized wind turbine gear by adopting the Fatemi-Socie multiaxial fatigue criterion, and the effects of compressive residual stress and the tensile residual stress on RCF life were studied. Thomas [12] presented an empirical formula for hardness with depth for carburized gears. Witzig [13] proposed the relationship between hardness and local strength. Zeng and Chiu [14] attempted to link the indentationunloading curve to the elastic modulus and strainhardening properties by using the nanoindentation method. Pavlina and van Tyne [15] proposed the relationship between the yield strength and hardness of local material in steel. Boiadjiev et al. [16] defined the concept of material exposure and evaluated the risk of contact fatigue failure at each material point along the depth direction. However, this work did not consider the impact of elastic-plastic contact and focused only on the action of shear stress; the multiaxial stress state that occurs during gear contact was neglected. In this work, the concept of local fatigue failure risk is defined by considering both the timevarying multiaxial stress state and the hardness gradient along the depth of carburized gears used in wind turbines. The risk of pitting failure and of spalling failure are distinguished based on the critical depth position of contact failure risk value. The influence of design parameters, such as case hardening depth (CHD, defined as the depth at which the hardness is 550 HV) [13,17,18] and surface hardness, on fatigue failure risk values is evaluated.

Methodology
Engineering practice shows that the intermediate parallel stage small cylindrical gear in a 2 MW wind turbine gearbox suffers pitting or spalling failure frequently. Gears then fail within only a few years, causing considerable economic loss. This study focuses on the gear as sample with which to evaluate the risk of contact fatigue failure. Table 1 shows the basic gear parameters required for fatigue evaluation. Because the gear surface is finished, tooth surface roughness is neglected in this study.
Typical reference points are illustrated in Fig. 2(a). The HPSTC represents the highest point of single tooth contact, and the LPSTC is the lowest point of single tooth contact along the line of action. The variations in the tooth radius of curvature and the normal load (F*) of tooth surface can be obtained based on meshing theory [19]. The Hertzian contact halfwidth b and the maximum Hertzian pressure P H can then be calculated via Eq. (1) [20].
The distribution of P H (red solid line) along the line of action is calculated. As illustrated in Fig. 2 (right), a remarkable change is found in the contact pressure at HPSTC and LPSTC. It is true that the gear has varying loads across the involute geometry. However, it should be noted that the contact pressure in the single tooth meshing area is more pronounced than in the double tooth meshing area, and stress concentration is more pronounced when the gear engages at the pitch point. Thus, the gear RCF failure risk is primarily evaluated at this critical position. According to Ref. [20], the gear contact at any meshing moment could be represented as two curved bodies contacting each other. Furthermore, this model could be represented as a rigid circle contacting an infinite deformable half-space with the assumption of plane strain condition. According to Table 1, an equivalent Young's modulus (E) of the infinite halfspace and the equivalent radius of curvature (R) of the rigid circle are calculated by Eq. (2) [21]. A schematic diagram of the gear contact model is shown in Fig. 3.
The proposed fatigue evaluation model mainly consists of three aspects: the mechanical properties of material, a finite element elastic-plastic contact model, and the concept of failure risk due to contact fatigue at local material points, which are explained as follows.

The framework of fatigue failure risk model
Pitting and spalling are the two main fatigue failures in heavy-duty case hardening gears. Evidence shows that pitting mainly occurs at the subsurface (at a depth of several hundred micrometers). On the other hand, spalling mainly occurs at the transition zone between the hardening layer and the core (usually at a depth of several millimeters). The approach compares the local multiaxial equivalent stress with the local strength derived from the hardness curve. The failure risk can be defined as the ratio of equivalent stress to local material strength. Because surface roughness is not considered, this value varies only with depth. The fatigue failure risk value is expressed as [16]: where  local [ ]( ) y represents the strength of local material, and  equi ( ) y is the equivalent stress calculated using stress components combined with the multiaxial stress criterion. They are functions of depth y. The maximum magnitude of material failure risk can be calculated as: where ff ,max A is the position corresponding to the point at which a fatigue crack is most likely to be initiated.
Generally, mechanical properties along depth of a case hardening gear can be divided into three zones, namely, the hardened layer (commonly with a depth of several hundred micrometers), the transitional layer, and the core zone. As shown in Fig. 4, the hardened layer is the area where pitting might occur, while the transitional zone is the area where spalling is prone to initiate. Because the core is far away from the stress-affected zone, there is almost no fatigue failure occurring within the core. The fatigue failure risk can be evaluated by emphasizing the difference between pitting failure and deep spalling failure. A good design should control both the pitting risk and the spalling risk simultaneously.

Gradients of hardness and mechanical properties
To improve contact fatigue resistance, a case hardening process is widely used in gear industry. The hardness gradient along depth depends on the heat treatment process. For carburized gears, CHD and surface hardness are important parameters that determine loading capacity. The hardness of sample is obtained through two methods: calculated via an existing empirical formula, or directly measured via the Vickers indentation test. The testing force for hardness measurement is 4.9 N. The measured results (red hollow circle) are shown in Fig. 5. The hardness gradually decreases from case to core, with the maximum value reaches 650 HV. The core hardness is approximately 400 HV. The CHD is approximately 2.5 mm.
The hardness gradient from surface to core induced by carburizing processes has been studied by many researchers [12,22]. Among them, Thomas presented an empirical formula to depict the hardness curve for carburized case gears [12]: represents the surface hardness, core HV is the core hardness; and HV ,max y is the distance between the depth of maximum hardness and surface, in this study, it equals to 0. core y is the y-coordinate where core ( ) HV y HV  . The solid blue line in Fig. 5 shows the general profile of hardness, calculated by the Thomas empirical formula when the surface hardness, core hardness, and CHD are determined as 650 HV, 400 HV, and 2.5 mm, respectively. As can be seen, the  empirical value is in good agreement with the measured result.
Because multiaxial equivalent stress is applied, the local material strength is defined as the yield strength: Some works have studied graded material properties of hardened case materials [23,24]. Pavlina and van Tyne [15] showed a linear correlation between yield strength and Vickers hardness: To study the effects of CHD and surface hardness on the material contact fatigue failure risk, two sets of models were selected via the Thomas equation, as shown in Fig. 6. Figure 6(a) shows the effect of CHD on the yield strength curve when the surface hardness is fixed at 650 HV and the core hardness is fixed at 400 HV, while Fig. 6(b) shows the effect of surface hardness on yield strength when the core hardness is controlled at 400 HV and the CHD is fixed at 2.5 mm. Because the core hardness is mainly determined by the forging process of original material, the effect of core hardness is not studied in this work.

Stress calculation based on the elastic-plastic finite element model
The FE model used in the work is schematically shown in Fig. 7. The radius of curvature of two circles at the pitch point are calculated using gear meshing theory from Table 1 as R 1 = 0.248 m, R 2 = 0.049 m. The equivalent radius of the rigid body and the equivalent elasticity modulus are then derived as R = 0.0409 m   ,  0 mm  y 20 mm . The x-axis indicates the rolling direction, while the y-axis indicates the depth direction. To introduce gradients of mechanical properties caused by carburizing, the finite element model is geometrically partitioned horizontally. The thickness of a single layer is set to 0.01 mm, which is fine enough to describe gradients along depth. To obtain a high-quality mesh, the mesh size in the purple solid line section of Fig. 7 is set to 0.01 mm × 0.01 mm. The element type is CPE4. During simulation, the contact load moves horizontally from x = -5 mm to x = 5 mm, a range that guarantees that each critical material point undergoes a complete loading cycle. Because material points at the same horizontal level experience the same loading cycle, the present work mainly concerns the 500 equally spaced material points located on the blue dotted line (x = 0 mm, and y  [0, 5] mm).
Due to the highly localized stress under contact surface, plasticity may occur in heavy-duty wind turbine gears. An appropriate material model is required to capture the elastic-plastic response under heavy loads. Figure 8 shows a typical linear kinematic hardening in which E is the Young's modulus, and YS represents the yield strength. M represents the linear hardening modulus, defined as 5% of the Young's modulus E [26]. This stress-strain response is also used in Refs. [25,26] to simulate carburized steel. It is worth noting that, under extremely heavy load conditions, a shakedown state may occur, which means several loading cycles should be applied to obtain a stabilized stress-strain response. In this work, a stabilized stress-strain field after five loading cycles is used for further RCF failure risk determination. Figure 9 shows the equivalent plastic strain after the stressstrain response is stabilized. The high surface strength caused by carburizing leads to plasticity at a deeper location. The stress-strain response over five loading cycles at a depth of approximately 2 mm (which is found to have a relatively higher equivalent plastic strain) is also depicted in Fig. 9(b). As can be seen, a stabilized loop is obtained after the contact moves through five cycles of the half-infinite domain. The loop traced with a dark line is the stabilized loop.
The gear contact process induces typical nonproportional multiaxial stress states. Dang Van et al. [27] constructed a well-accepted multiaxial fatigue failure criterion using the maximum shear stress and the hydrostatic stress. The Dang Van criterion can be expressed as [27,28]: where  max is the maximum shear stress, and  h is the instantaneous hydrostatic component of stress tensor.  and  are both material parameters that can be obtained by [28]: where  1 represents the fatigue limit in reversal torque, and  1 is the fatigue limit of the fully reversed tension. The ratio     1 1 / is set as 0.577 [30].
The fatigue performance of a material can be quantified by the equivalent stress value evaluated by weighted sum of the local maximal shear stress ( max ) and the hydrostatic stress (  h ). The Dang Van equivalent stress can be expressed by [30][31][32]: The Dang Van equivalent stress distribution of target material points is calculated according to Eq. (10). The failure risk considering mechanical property gradients is also depicted. Figure 10 and Table 2 display the results. The maximum equivalent stress is found at a depth of approximately 1.08 mm, while the maximum failure risk value occurs at a depth of 1.25 mm, deeper than the maximum equivalent stress. Simulation results indicate that mechanical property gradients influence the failure depth.

Effect of surface hardness
A fine control of surface hardness is an important step in the manufacturing process. In many situations engineers realize the importance of surface hardness, but they have limited knowledge of how exactly surface hardness influences fatigue performance and service life. In this section, the influence of surface hardness on the distribution of ff A is studied under the following conditions: = 675, 700, 725, and 750 HV experience a purely elastic state. Figure 11 shows the variation in fatigue risk and the equivalent stress as the surface hardness changes. The profiles of equivalent stress and ff A show remarkable discrepancy. This means that, a method without the consideration of the local material strength is not suitable for RCF failure risk evaluation of carburized gears. It is obvious that among the selected   | https://mc03.manuscriptcentral.com/friction cases, the critical position at 2.5 mm (black dotted line) depth is a featured point. In areas with depth less than 2.5 mm, as the surface hardness increases, the risk value decreases. Conversely, in areas with depth less than 2.5 mm, as the surface hardness increases, the risk value increases. Furthermore, when the surface hardness is less than 700 HV, the location of ff ,max A is at the position where pitting failure may happen, but when the surface hardness is higher than 700 HV, ff ,max A occurs at deeper positions, which means the risk of spalling failure increases. The predicted failure pattern under different surface hardness conditions is listed in Table 3. Figure 12 shows the effect of surface hardness from another point of view. The position of ff ,max A becomes deeper as the surface hardness increases. It also clearly shows that the critical value of 700 HV distinguishes the pitting failure risk and the spalling failure risk. From the right section of Fig. 12, it can be seen that the case with surface HV = 700 HV has the lowest predicted value of ff ,max A . The ff ,max A of this carburized gear at a surface hardness of 700 HV is 0.67, which is 9.5% less than the case with a surface hardness of 600 HV.

Effect of CHD
Case thickness is one of the most important parameters during manufacturing. The CHD is often defined as the distance between the surface and the depth at which the hardness is 550 HV. In this section, the effect  Fig. 13. Stress results show that the models with CHD = 2.00, 2.25, and 2.50 mm generate plasticity under the given conditions, while the models with CHD = 2.75, 3.00, 3.25, and 3.50 mm predict a purely elastic response. It is evident that with the increase in CHD the value of ff A decreases at the transition zone between the hardening layer and the core. The resistance to spalling improves as CHD increases, which can be achieved in engineering practice by extending the duration of carburization process. As CHD decreases, the regime with high-risk values enlarges, leading to a higher risk of spalling failure. However, the value of CHD has little impact on pitting risk underneath surface. Figure 14(a) shows the distribution of ff A as CHD changes from 2.0 to 3.5 mm. The blue curve represents the maximum fatigue failure risk value under the selected case using the elastic-plastic model. It can be seen, for the given conditions, when CHD is less than 2.5 mm, plasticity occurs within the substrate, caused by the reduction in the yield strength as CHD decreases. Within the range where plasticity occurs in the elastic-plastic model (CHD ≤ 2.5 mm), the position with the maximum risk value becomes shallower as the CHD decreases. When CHD is large enough (CHD > 2.5 mm in this case), no plasticity occurs within the material. As CHD continues to increase, the position of ff ,max A becomes slightly shallower. Figure 14(b) shows the nonlinear relationship between CHD and ff ,max A . Compared with the result of case with CHD = 2.0 mm, ff ,max A is 4.2% less in the case with CHD = 3.5 mm. Results indicate that as CHD increases, the probability of fatigue failure within the material decreases.

Effect of normal load
In this section, the impact of load on fatigue failure risk is studied under the following conditions: surface HV = 650 HV, core HV = 400 HV, CHD = 2.5 mm. Figure 15 depicts ff A along depth under several normal load cases. Results show that as P 0 increases, ff A increases everywhere, which is easy to understand. As the load increases and plasticity occurs, the ff A curve at the transition zone becomes flat indicating the increasing of deep spalling risk. Figure 16(a) shows how the distribution of ff A varies with P 0 in the range of 1.95-3.45 GPa. As can be seen, in the purely elastic state (in this case, P 0 less than 2.7 GPa), as P 0 increases, ff ,max A increases linearly, but as the load further increases in the areas where plasticity occurs, the increase in ff ,max A slows down. The position of ff ,max A becomes deeper with an almost linear relationship with P 0 , but as the load increases and plasticity enhances, the position of ff ,max A remains at the same depth of approximately 1 mm. Figure 16

Conclusions
In this work, a method for evaluating the contact fatigue failure risk of carburized gears is proposed based on the multiaxial stress state and gradients of mechanical properties. The concept of local material fatigue failure risk is defined to investigate the risk value for pitting and spalling. Conclusions can be made as follows: (1) Keeping the surface hardness surface HV and the core hardness constant, as the effective case hardening depth increases, the resistance to spalling failure improves, while the resistance to pitting is almost unaffected.
(2) Keeping the case hardening depth and core hardness constant, an increase in surface HV leads to a decreased risk of pitting failure, but the risk of deep spalling failure increases. When the surface hardness  exceeds 700 HV, the maximum risk value occurs in the transition area, which means the gear suffers a higher risk of spalling than pitting.
(3) As the load continues to increase leading plasticity, ff A at the transition zone becomes flat, which indicates an increase in spalling risk. In a purely elastic state, as the normal load increases, ff ,max A increases linearly, but as the load further increases where plasticity happens, the increasing trend of ff ,max A slows down.