Time Domain Modeling and Analysis of Dynamic Gear Contact Force in a Wind Turbine Gearbox with Respect to Fatigue Assessment

Abstract

The gearbox is one of the most expensive components of the wind turbine system. In order to refine the design and hence increase the long-term reliability, there has been increasing interest in utilizing time domain simulations in the prediction of gearbox design loads. In this study, three problems in time domain based gear contact fatigue analysis under dynamic conditions are discussed: (1) the torque reversal problem under low wind speed conditions, (2) statistical uncertainty effects due to time domain simulations and (3) simplified long term contact fatigue analysis of the gear tooth under dynamic conditions. Several recommendations to deal with these issues are proposed based on analyses of the National Renewable Energy Laboratory’s 750 kW land-based Gearbox Reliability Collaborative wind turbine.

1. Introduction

With the growth of the share of wind energy in the energy market, the design and implementation of large scale wind turbines has become a common occurrence. However, since its inception the wind energy industry has experienced high gearbox failure rates [1]. In order to achieve their stated design life goals of 20 years, most systems require significant repairs or overhauls well before the intended life is reached [2,3,4]. Some firm conclusions about the nature of the failures have been made based on the work of Musial et al. [5]: (i) gearbox failures are not specific to a single gear manufacturer or turbine model, which are general; (ii) poor adherence to accepted gear industry practices, or otherwise poor workmanship, is not the primary source of failures; (iii) most gearbox failures do not begin as gear failures or gear-tooth design deficiencies, and up to 10% of gearbox failures may be manufacturing anomalies and quality issues that are gear related; (iv) the majority of wind turbine gearbox failures appear to initiate in the bearings, and (v) it is believed that the gearbox failures observed in the earlier 500 kW to 1000 kW sizes 5 to 10 years ago may still occur in many of the larger 1 to 2 MW gearboxes being built today with the same architecture. With larger wind turbines, the cost of gearbox rebuilds, as well as the down time associated with these failures, has become a significant portion of the overall cost of wind energy [6]. Presently NREL is performing a long-term project to improve the accuracy of dynamic gearbox testing to assess gearbox and drivetrain options, problems, and solutions under simulated field conditions. In addition, in order to increase the long-term reliability of gearboxes and make their design more reasonable, there is increasing interest in utilizing time domain simulations and physical tests in the prediction of gearbox design loads, with the continual development of computer technology, simulation tools and measurement equipments. In several previous studies Klose et al. [7] performed an integrated analysis of wind turbine behavior and structural dynamics of a jacket support structure under combined wind and wave loads in the time domain. Seidel et al. [8] used the sequential coupling and the full coupling methods to simulate offshore loads on jacket wind turbines, and validated these methods using measurement data from the DOWNVInD project. Gao and Moan [9] performed long-term fatigue analysis of offshore fixed wind turbines using time domain simulations. Dong et al. [10] performed long-term fatigue analysis of multi-planar tubular joints for jacket-type offshore wind turbines using time domain simulations. Peeters et al. [11,12], Xing et al. [13] performed a detailed analysis of internal drive train dynamics in a wind turbine using multi-body simulations. However, there is, at present, limited literature concerning long-term time-domain based analysis of mechanical components, e.g., main shaft, gears and bearings, in the wind turbine drive train system under dynamic conditions. This is mainly due to the complexities involved in modeling and simulating the drivetrain with respect to the computation efforts and scale. Recently, Dong et al. [14] established and applied a long-term time domain based gear contact fatigue analysis of a wind turbine under dynamic conditions. In the present study, several practical problems of time domain based gear contact fatigue analysis encountered in [14] are described and discussed. These are: (1) the rotation reversal problem of gears under low wind speed conditions, (2) the statistical uncertainty effect due to the time domain simulation and (3) simplified long term contact fatigue analysis of gear teeth under dynamic conditions. Several useful suggestions to address these issues are proposed.

2. Time-Domain Simulation of Wind Turbines

In this study, a 750 kW land-based wind turbine from the Gearbox Reliability Collorative (GRC) project coordinated by the National Renewable Energy Laboratory (NREL) in Colorado, USA is used as the case study. This is a three-bladed upwind turbine. The nominal hub height is 55 m, the rotor diameter is 48.2 m; the rated generator speed is 22/15 rpm, which represents a two-speed generator [four-pole (4P) and six-pole (6P) generator] with rated power of 750 kW and 200 kW, respectively; the nominal cut-in wind speed is 3 m/s; the rated wind speed is 16 m/s; the cut-out wind speed is 25 m/s; stall-regulated control is applied; the design wind class is IEC Class II, and the design life is 20 years. The configuration of the drive train system of the wind turbine is as shown in Figure 1. The performance property of the wind turbine is as shown in Figure 2. More details about this wind turbine can be found in [15,16].

Figure 1. Drive train configuration of the wind turbine [16].

Figure 1. Drive train configuration of the wind turbine [16].

Figure 2. Drive train configuration of the wind turbine [17].

Figure 2. Drive train configuration of the wind turbine [17].

The analysis proceeds in two steps. First, global aero-elastic simulations are performed using the FAST code from NREL [18]. FAST is an aero-elastic code that computes the coupled wind turbine structural response under aerodynamic load effects. The time series of the main shaft torques are obtained and used as inputs in a multibody gearbox model in SIMPACK [19]. SIMPACK is a multi-purpose multibody code with special features available to model gearboxes. Figure 3 shows the gearbox internal components and Figure 4 shows the relative gearbox model in SIMPACK. Figure 5 shows an example of the calculation results of the torques in the mainshaft when using the 4P generator.

Figure 3. View of the gearbox internal components [16].

Figure 3. View of the gearbox internal components [16].

Figure 4. Gearbox model primitive in SIMPACK.

Figure 4. Gearbox model primitive in SIMPACK.

Figure 5. Time histories of main shaft torques when using the 4P generator.

Figure 5. Time histories of main shaft torques when using the 4P generator.

In SIMPACK, each component of the gearbox is modeled as a rigid body and interconnected using joints and force elements. The topology diagram of the gearbox model is shown in Figure 6.

Figure 6. Topology diagram of gearbox model in SIMPACK.

Figure 6. Topology diagram of gearbox model in SIMPACK.

The gear pair force element in SIMPACK, FE 225, is used for modeling gear contact. FE225 models gear contact as a series of discrete springs and dampers. The gear stiffness is calculated in accordance with ISO 6336-1 [20]. This stiffness parameter depends on the location of the contact point and the gear geometry. The forces and torques acting on each individual gear are calculated as a function of the gear stiffness and the penetration depth at the gear teeth. FE225 also considers normal damping, coulomb friction, backlash and micro-geometry. The contact forces in meshing gears are obtained using the classical slicing method. Helical gears can be regarded as several very thin cylindrical gear wheels mounted on a mutual axis and individually rotated a small angle around their common axis [21]. In this way, a helical tooth is sliced into several independent cylindrical teeth. More details about the slicing method can be found in [22,23]. In this study, the time series of the gear contact forces at the contact point on a certain slice generated in the meshing gears for different wind speeds are used to perform the studies in Section 3, Section 4 and Section 5

3. The Gear Torque Reversal Problem

The two step analysis procedures as described in Section 2 are used. The time series of the wind turbine main shaft torques and the gear contact forces generated in meshing gears are obtained from FAST and SIMPACK simulations, respectively. In the analysis shown in Figure 5, the torques are negative for short time periods at low wind speeds, i.e., 4 m/s, 6 m/s and 8 m/s. This means that the gears do not mesh on a single side the entire time. Generally speaking, this phenomenon is bad for gearbox reliability. Figure 7 shows an example of the time series of the gear contact forces when using the 4P generator and a wind speed of 4 m/s. The notations, “tooth −1”, “tooth 0” and “tooth +1”, are also defined in Figure 8. They represent the different stages of gear meshings, which are engagement, middle and recess stage, respectively. In the analysis reported in Figure 7, the contact forces experienced by the gear teeth are not always positive, which is in consistent with the cases reported in Figure 5. This means that the gear contact is not merely on one side of the gear teeth, e.g., surface A. Contact occurs on surface B when the contact forces are negative, as shown in Figure 9(a). This is beneficial from a contact fatigue point of view. However, this is very bad from a tooth root bending fatigue point of view. Cracks could be initiated at the root location of gear tooth and propagate into the interior, which leads to failure. This is shown in Figure 9(a). Another important problem caused by this phenomenon is the postprocessing of the contact forces. In order to do the time domain based gear contact fatigue analysis, the time series of gear contact forces at a selected contact position on a certain gear tooth is required. However, the torque reversal problem makes the postprocessing difficult and inaccurate. A simplified method has to be applied, which is described in the following paragraph.

Figure 7. Time histories of gear contact force based on time domain simulation in SIMPACK (U_w = 4 m/s).

Figure 7. Time histories of gear contact force based on time domain simulation in SIMPACK (U_w = 4 m/s).

Figure 8. Scheme of meshing gears.

Figure 8. Scheme of meshing gears.

Figure 9. Scheme of a gear tooth.

Figure 9. Scheme of a gear tooth.

In this study, three different generator controllers are considered and compared at low wind speeds (4 m/s, 6 m/s and 8 m/s). These are the 4P controller, 6P controller and a simple variable speed controller (VS). The 4P controller and the 6P controller are obtained from NREL directly. The simple variable speed controller is designed in-house. The variable speed controller is designed based on the principle of maximizing the power production at the low wind speeds. This means that the torque speed curve fit as shown on the right hand side of Figure 10 passes through maxima of the individual curves shown on the left hand side. Figure 11 shows the comparison of the torques for different generator controllers at U_w = 4 m/s. As observed in the figure, there are no negative torques. The 6P controller is better than the 4P controller, but there are still many negative torques. As the torque reversal problem is severe for the 4P and the 6P controllers at low wind speeds, a continuous segment taken from the entire time series of the contact forces for each wind speed is identified in terms of the maximum value of the mean contact forces, and is used to get the mean value and the standard deviation of the contact forces for each simulation sample from here on, which is a simplified method used in this study. It is noted that the segment lengths for different simulation samples are not uniform, which should be determined case by case.

Figure 10. Design of a simple variable speed controller.

Figure 10. Design of a simple variable speed controller.

Figure 11. Time histories of main shaft torques using different generator controllers.

Figure 11. Time histories of main shaft torques using different generator controllers.

Figure 12 shows the comparison of the power curves for different generator controllers with respect to three different wind speeds. In the analysis reported in this figure, there are no significant differences among the 4P, the 6P and the VS controllers at U_w = 4 m/s and U_w = 6 m/s. The 4P and VS controllers produce more power than the 6P controller at U_w = 8 m/s. Furthermore, there are no negative values for the VS controller at all the wind speeds simulated.

Figure 13 (F_c represents the gear contact forces obtained from time domain simulations) shows the comparisons of the variation of mean contact forces with the increment of simulation samples at the three different contact points at different wind speeds and using different generators. As previously mentioned, p1 is the engagement point, pitch0 is the pitch point and m1 is the recess point. 20 simulations are performed for each wind speed. The simulation length is 700 s, with the first 100 s discarded. In the analysis reported in Figure 13, the variations of the mean contact forces at three different contact points are very similar. The mean contact force of the VS controller is smaller than those of the 4P and the 6P controllers. The simplified treatment of the time series of the contact forces for the 4P and the 6P controllers mentioned above could be conservative. The mean contact force for the 4P controller increases faster than the 6P and the VS controllers, with increasing wind speeds. Generally speaking, from a contact fatigue and bending fatigue point of view, the 6P and the VS controllers could be better than the 4P controller. From a power generation point of view, the 4P controller might be better than the 6P and the VS controllers, especially when the wind speeds are higher than 6 m/s.

Figure 12. Power curves at different wind speeds using different generators.

Figure 12. Power curves at different wind speeds using different generators.

4. Statistical Uncertainty Effect Due to Time-Domain Simulation

For wind turbines, their performances are subject to a number of uncertainties, which are caused by inherent physical randomness of the system or environment (named as aleatory undertainties, e.g., physical wind process) and lack of knowledge of the system or the environment (named as epistemic uncertainties, e.g., model uncertainty and statistical uncertainty). A rational treatment of these uncertainties in a quantitative manner associated with an engineering problem and its physical representation in an analysis are the essence of reliability analysis. In this section, the statistical uncertainty due to time domain simulation is considered, which is mainly due to that different sample data sets usually produce different statistical estimators such as the sample mean, std., etc.

Figure 13. Variation of mean contact forces when using different generators at the different contact points, wind speeds and number of simulation samples.

Figure 13. Variation of mean contact forces when using different generators at the different contact points, wind speeds and number of simulation samples.

For time domain simulation of wind turbines, one of the key components is the simulation of turbulent wind field. In this study, the turbulent wind field is simulated using the TurbSim code [24], which is a stochastic, full-field, turbulent-wind simulator. It uses a statistical model to numerically simulate time series of three-component wind-speed vectors at points in a two-dimensional vertical rectangular grid that is fixed in space. Figure 14 shows an example of a TurbSim wind field. In TurbSim, random phases (one per frequency per grid point per wind component) for the wind velocity time series are created using the random numbers generated by the pseudorandom number generator. For the same mean wind speed, if the random numbers are different, the wind velocity time series will be different. This will affect the calculations of aerodynamic forces applied on the blades and torques in the main shaft of wind turbine. Figure 15 shows an example of the torque time series obtained in FAST simulation for the same wind speed (U_w = 16 m/s) with respect to different random numbers. Generally speaking, the more simulation samples for the same mean wind speed, the better the simulation results. However, time domain simulations are usually very time consuming and the data size is very big, therefore it is necessary to identify a suitable sample size for a certain wind speed. This will cause the so-called statistical uncertainty. In this study the statistical uncertainty of the contact force calculation with respect to the number of simulations employed, at different wind speeds are estimated. The range of 1-h mean wind speed U_w is 4–24 m/s with an increment of 2 m/s. For each wind speed, 20 simulations are performed using different random seed numbers. For each simulation, a time history of the gear contact forces is obtained and used as a sample. Therefore, totally 20 samples for each wind speed are obtained.

Figure 14. Example of wind field simulation in TurbSim [24].

Figure 14. Example of wind field simulation in TurbSim [24].

Figure 15. Time series of main shat torques in FAST simulation at the same wind speed with different simulations (random seeds are different).

Figure 15. Time series of main shat torques in FAST simulation at the same wind speed with different simulations (random seeds are different).

Figure 16 shows the variation of the coefficient of variation (C.O.V.) of the gear contact forces at the three different contact points for different wind speeds and different simulation samples. The quantities (C.O.V.) in this figure are determined from the ensemble of x samples (x = 1, 2, 3, …, 20), which are the same as the quantities in Figure 17 and Figure 18. The 4P generator controller is used. The parameter C.O.V. is defined as follows:

$$C.O.V.=\frac{{\sigma }_F}{{\mu }_F}$$

(1)

where μ_F represents the mean value of the contact forces; σ_F represents the standard deviation of the contact forces.

Figure 16. Variation of C.O.V. of the contact forces when using the 4P controller at different contact points, wind speeds and number of simulation samples.

Figure 16. Variation of C.O.V. of the contact forces when using the 4P controller at different contact points, wind speeds and number of simulation samples.

Figure 17 shows the variation of the parameter ζ at the three different contact points for different wind speeds and different simulation samples. The 4P generator controller is used. The parameter ζ is defined as follows:

$$\zeta =\frac{{C.O.V.}_c^i-{C.O.V.}_{20}^c}{{C.O.V.}_{20}^c}\cdot 100\%$$

(2)

where ${C.O.V.}_c^i$ represents the coefficient of variation of the contact forces based on i simulation samples (i = 1, 2, 3, …, 20).

Figure 17. Variation of ζ [Equation (2)] when using the 4P controller at different contact points, wind speeds and number of simulation samples.

Figure 17. Variation of ζ [Equation (2)] when using the 4P controller at different contact points, wind speeds and number of simulation samples.

In the analysis reported in Figure 16, the variations of the C.O.V. of the contact forces at three different contact points are very similar. At low wind speeds, the C.O.V. value can be up to 0.6. It is decreased significantly with the increment of wind speeds, which represents that the turbulence effect is decreased with the increment of wind speeds. In the analysis reported in Figure 17, the variations of the parameter ζ [Equation (2)] at three different contact points are very similar. The values of ζ at low wind speeds (U_W < 12 m/s) are much higher than those at other wind speeds, especially when the simulation samples are less than 10. The maximum value of ζ can be up to 24%. The findings here suggest that there might be higher levels of uncertainty when using the simplified method discussed in the previous section. This is because the simplified method employs a portion of the original simulation length, i.e., shorter simulation length. It is, therefore, suggested that more samples be used when employing the simplified method. Based on the work in this study, at least 10 simulation samples might be used at low wind speeds, in order to keep the value of ζ within 5%.

In Figure 16, the 4P generator controller is used. In this study, when the 6P and the VS controllers are used, the C.O.V. values for U_W = 4 m/s, 6 m/s and 8 m/s are compared against the cases of the 4P controller, as shown in Figure 18. In the analysis reported in this figure, the variations of the C.O.V. values of the contact forces at three different contact points are similar for each controller. The C.O.V. values for the 6P and the VS controllers are less than that of the 4P controller. From an uncertainty level point of view, the benefits of the 6P and the VS controllers are increased with the increment of wind speeds. One reason is that the time series of the contact forces used for the 6P and the VS controllers are much longer than that of the 4P controller, especially for the VS controller. Figure 18 shows the comparison results at the pitch point. The cases at other points, e.g., the engagement point and the recess point, are similar to that of the pitch point, and not repeated herein. Based on the work in this study, if the 6P or the VS controller is applied, at least 6-8 simulation samples might be used at low wind speeds (U_W < 12 m/s) in order to obtain a relatively stable C.O.V. value of the contact forces.

Figure 18. Variation of C.O.V. of the contact forces when using different controllers at the pitch contact point.

Figure 18. Variation of C.O.V. of the contact forces when using different controllers at the pitch contact point.

In this section, the statistical uncertainty of the gear contact force calculation due to time domain simulation is investigated. It should be noted that a decoupled dynamic response analysis method and a simplified rigid gearbox model are applied as described in Section 2. These simplifications will cause the so-called model uncertainties, which can be reduced by using more refined method (e.g., fully coupled method) and more accurate model (e.g., gearbox model with flexible components). The availability of the decoupled analysis method used in this study is verified in [14]. The detailed analysis of model uncertainties is out of the scope in this paper and might be investigated in further work. In addition, the statistical uncertainty analysis in this study is mainly with respect to time domain based gear contact fatigue assessment in normal operations. It should be noted that the cases of wind turbine extreme loads analysis in normal operations are different. For extreme loads analysis, the peak values obtained from time histories of a certain variable should be used and the exceedance probability theory is applied. More details about wind turbine extreme loads analysis can be found in [25,26,27,28].

5. Simplified Gar Contact Fatigue Analysis

Recently Dong, et al. [14] presented a simplified predictive pitting model for estimating gear service lives and verified its validity by comparing with the published experimental evidence. The model sets up a link between the global dynamic response analysis of wind turbine and the detailed contact fatigue analysis of gears in the drive train, which is given as follows:

$$N_p={\displaystyle \underset{a_0}{\overset{a_c}{\int }}\frac{1}{C\cdot G_{2a}^m\cdot U^m\cdot \Delta {\overline{p}}_{\max }^m}}da$$

(3)

where N_p represents the number of cycles for crack propatation. C and m represent the material constants. G_2a represents the geometry function. U is a factor related to crack closure. $\Delta {\overline{p}}_{\max }$ represents the equivalent maximum contact pressure range, which is given as follows:

$$\Delta {\overline{p}}_{\max }={\left({\displaystyle \underset{0}{\overset{\infty }{\int }}\Delta p_{\max }^m\cdot f_{\Delta p_{\max }}\left(\Delta p_{\max }\right)d\Delta p_{\max }}\right)}^{1/m}$$

(4)

where Δp_max represents the maximum contact pressure range, which can be calculated from time domain simulations or field tests. f_Δpmax represents the probability density distribution of Δp_max.

In order to obtain Δp_max, accurate postprocessing of the tooth contact loads is important. Figure 19 shows the sun gear with the gear teeth numbered. In order to perform time-domain based gear contact fatigue analysis, the time series of the contact forces for each gear tooth should be obtained by postprocessing the MBS simulation results. In reality, this procedure is very time consuming and inconvenient, especially for low wind speeds where the torque reversal problem is severe. A simplified approach is, therefore, desired. In this study the mean values and the standard deviations of the contact forces at three different locations on a representative tooth surface of the sun gear using three different methods are calculated. The three locations are as shown in Figure 9(b). The three different methods are described as follows:

• Maximum damage method: for each wind speed, a most dangerous gear tooth is identified in terms of the maximum mean contact force. A dummy gear tooth is then defined. It is further assumed that the time series of the contact forces of the most dangerous gear tooth for each wind speed are all applied on this dummy gear tooth, which is taken as the representative of all 21 gear teeth.
• Minimum damage method: this method is similar to the first one, for each wind speed, a safest gear tooth with the minimum mean gear tooth contact force is identified. A dummy gear tooth is then also defined. It is also further assumed that the time series of contact forces of the most safe gear tooth for each wind speed are all applied on this dummy gear tooth, which is also taken as the representative of all 21 gear teeth.
• Simple average method: the time series of the gear contact forces from SIMPACK simulations are used directly. It is assumed that the time series of the contact forces at the same contact position on each tooth surface of the sun gear are the same.

Figure 19. Front view of the sun gear in the SIMPACK model.

Figure 19. Front view of the sun gear in the SIMPACK model.

Figure 20 shows the comparisons of the mean contact forces using the three different methods at the three different contact points at different wind speeds.

Figure 20. Comparison of mean contact forces when using different methods at different contact points and wind speeds (4P generator).

Figure 20. Comparison of mean contact forces when using different methods at different contact points and wind speeds (4P generator).

Figure 21 shows the comparisons of the standard deviations of the contact forces. The 4P controller is applied here. In the analysis reported in these figures, at the contact points p1 and m1, the mean contact forces and the standard deviations for different wind speeds are almost the same with respect to the three different methods mentioned above. At the contact point pitch0, the mean contact forces and the standard deviations for different wind speeds using the simple average method are a little smaller than those using the maximum damage method and minimum damage method, which is slightly not conservative. In addition, the cases of the 6P and the VS controllers are also considered in this study, and compared with the results of the 4P controller, as shown in Figure 22 and Figure 23. In the analysis reported in these figures, the cases of the 6P and the VS controllers are very similar to the case of the 4P controller. Figure 20, Figure 21, Figure 22 and Figure 23 show the validity of the simple average method with respect to different wind speeds and different generator controllers.

Figure 21. Comparison of standard deviations of contact forces when using different methods at different contact points and wind speeds (4P generator).

Figure 21. Comparison of standard deviations of contact forces when using different methods at different contact points and wind speeds (4P generator).

Figure 22. Comparison of mean contact forces when using different methods at different contact points, generators and wind speeds.

Figure 22. Comparison of mean contact forces when using different methods at different contact points, generators and wind speeds.

Figure 23. Comparison of standard deviations of contact forces when using different methods at different contact points, generators and wind speeds.

Figure 23. Comparison of standard deviations of contact forces when using different methods at different contact points, generators and wind speeds.

6. Conclusions

In this paper, three practical problems encountered in time domain based gear contact fatigue analyses are discussed. The main conclusions are presented below:

(1) The reverse rotation problem is severe at low wind speeds (e.g., 4 m/s–8 m/s) for the 4P and the 6P generator controllers. In order to avoid it, the VS generator controller could be used. At low wind speeds, a continuous segment taken from the entire time series of the gear contact forces with respect to the 4P and the 6P generator controllers, which has the maximum mean contact force, might be used to do the time-domain based gear contact fatigue analysis, and the results could be on the conservative side.

(2) The variations of the C.O.V. values at different contact points are similar. In order to keep the value of ζ [Equation (2)] within 5%, at least 10 simulation samples might be used at low wind speeds (U_w < 12 m/s) for the 4P generator controller if the simplified method mentioned in Section 3 is employed, and at least 6–8 simulation samples might be used for the 6P and the VS generator controllers. For high wind speeds (U_w > 12 m/s), 5–6 simulation samples might be ok for the 4P, the 6P and VS generator controllers.

(3) The simple average method is equally accurate for different wind speeds and different generator controllers, and is more efficient than the maximum damage method and the minimum damage method.

In this study, only the main shaft torque loads are considered. The effects of non-torque loads could be investigated in future work. In addition, a simple rigid body gearbox model is used in this paper, more refined gearbox models, e.g., with flexible components modeled, could be applied in future work. Furthermore, time-domain simulations of larger megawatt wind turbines (onshore and offshore) could be also performed in future work.