Abstract
Premature failure of gearboxes is a big challenge facing the wind power industry. It highly depends on fully understanding the embedded dynamics to solve this problem. To this end, this paper investigates the random vibration and dynamics of planetary gear trains (PGTs) in wind turbines under the excitation of wind turbulence. The turbulence is represented by the Von Karmon spectrum and implemented by passing white noise through a 2nd-order shaping filter. Then, extra equations are formed and added to the original governing equations of motion. With this augmented equation set, a recursive numerical algorithm based on stochastic Newmark scheme is applied to solve for the statistics of the responses starting from initial conditions. After simulation, the variances of the vibration responses and the dynamic meshing forces at gear meshes are obtained.
1. Introduction
With increasing public awareness of the unsustainability of traditional fossil energy and its severe effects on the environment, renewable energy gains more and more share in the total energy consumption. Wind power, as one of the several types of renewable energy, has experienced fast growing in the past decade. More wind power systems are expected to be deployed in the future with the pushing from governments and the public for environment protection. In a wind turbine, a drive train connects the rotor shaft to the generator. It changes the relatively lower rotation speed of the rotor to a much higher one on which the generator works efficiently. There are various types of drive trains in use, such as direct drives and torque splitting transmissions; however, traditional gearboxes are still predominant in application because of economical and/or technological reasons. Due to the high speed ratio required, a gearbox typically consists of two or three stages of planetary gear trains (PGTs) or one stage of PGT combined with parallel shaft gear trains. One of the big challenges faced by the wind power industry and researchers at present is the premature failure of the gearbox [1]. A gearbox which is typically designed with a life expectancy of 20 years has parts failed commonly after only 3-to-5-year operation [2]. Consequently, replacement of the failed gearbox, or main maintenance on it, is required, causing significant downtime and high cost. It is reported that gearbox failure causes the highest downtime in wind turbines and leads to the most lengthy and costly maintenance procedures [3]. The economical consequence of gearbox failure in wind power industry is so significant that insurers are changing their policies toward wind power operators [1]. Given the significant impact of gearbox failures on the industry, great effort has been made in recent years toward understanding the dynamics of wind turbine gearbox and relavent relibility and fatigue issues. Peeters and his collaborators [4–6] investigated the modal characteristics of PGTs in a wind turbine, and the responses to abruptly applied impacts, such as braking. Methods used in modeling general PGTs were followed in these papers [7–9]. Guo et al. [10] investigated the dynamics with extended harmonic balance approach and compared the analytic results against available experimental data. Also possible tooth wedging and the unique effect of gravity on a wind turbine PGT were examined by Guo and Parker [11] and Guo et al. [12, 13]. In these studies [10–13], focus was placed on the complex nonlinear phenomena in the responses. In modeling the system, Qin et al. [14] used multibody method; the dynamic forces in gear meshing and bearings under a deterministic torque are simulated. In recent years, Zhu and his associates looked at the effect of flexible pins [15] and modeled the dynamics based on measured spectrum [16] and the rigid-flexible coupling dynamics of the drive train in a megawatt wind turbine [17]. Nejad et al. [18] investigated the possibility of floating the sun gear to balance geometric imperfection. Guo et al. [19] studied the load-sharing effect under nontorque loads. Xing and Moan [20] examined the effect of flexibility of substructures within the drive train with Simpack. In spite of the plentiful publications, it is still far from fully understanding the essential dynamics of the gear trains. As pointed out in [1], the lack of full accounting of the critical loads and the transfer of the loads between subsystems constitutes a barrier to understanding the underlying mechanisms of gearboxes, preventing the problem from being solved. One critical missing piece to thorough understanding the dynamics, to the author’s understanding, is proper account of the randomness of the wind loading in the wind turbine and its effect on the dynamics. Without a proper consideration of the randomness in the wind load, it is hard to say the key feature of the dynamics is captured in a wind turbine system. In fact, some researchers have noted this point already and conducted some research works correspondingly. Guo [21] developed a data-driven stochastic analytic model for the PGT in a wind turbine. A low-dimensional yet realistic stochastic wind flow model was used. Qin and Tang also developed a stochastic wind model using the Welbull distribution and studied the dynamic effect with Monte Carlo simulation [22]. Torsional vibration of a drive train under the combined effect of a deterministic mean wind speed and a random turbulence was investigated in [23]. The random excitation was modeled as a nonstationary random process following a predefined spectrum. The effect of random parameters on the dynamic responses was investigated in [24], which simplified the whole drive train as a two-stage gear system with 4 degrees of freedom. J. Yang and P. Yang [25] modeled the dynamics of a PGT under an random external excitation of Gaussian white noise, but no consideration was given to the wind turbine application.
There are several types of excitation to a PGT in wind turbines; the main three include(a)the periodic rotor torque generated by the mean wind speed,(b)the displacement excitation caused by the transmission errors of the gear meshes,(c)the random rotor torque caused by the wind turbulence.
This paper focuses on modeling the dynamics of a PGT to the wind turbulence, the excitation of type (c). Justification for this treatment is twofold. First, the vibrations caused by excitation (a) and (b) are basically deterministic; and there are already plenty of publications dealing with these two parts (see the many cited papers above and [26–28]). Second, knowing the responses related to each individual excitation: (a), (b), and (c), the total response can be easily obtained if needed. The random turbulence in this paper is modeled with the widely used Von Karmon spectrum which is different from the wide band white noise used in [25]. The nonwhite spectrum is implemented in the simulation by passing white noise through a 2nd-order filter; therefore, the overall dynamic system is augmented by combining the original system and the filter.
The paper is organized as below: Section 1 gives a brief literature review. In Section 2, the dynamic model for the PGT is developed. Following that, the vibration responses are simulated and analyzed in Section 3. In Section 4 conclusions are drawn, and future works are discussed as well.
2. Modeling PGT Vibration
In a wind turbine, the rotor receives the wind and generates an aerodynamic torque and a thrust which are transmitted to the carrier of the PGT through a low speed shaft. The sun gear connects a generator directly or through other gear pairs. The PGT under investigation is schematically shown in Figure 1, excluding possible parallel shaft gears. In this model, all the parts are regarded as lumped masses interconnected by springs which represent bearings and/or gear meshes. The ring is fixed without motion. Given the fact that lateral and axial motions of the PGT are coupled with the motion of the tower which is hard to determine at this point, only rotational motions of the carrier and the planet gears are considered. On the other hand, the sun gear is often floated in wind turbines for load-sharing purpose; thus, the transverse motions in both and directions are accounted for. If the PGT is broken down into a carrier, a sun gear, and planet gears, the motion of parts could be represented as below.
2.1. Carrier
The carrier receives the aerodynamic torque from the rotor. Its equation of motion is represented by where the subscripts , , , and indicate the ring, the planet, the sun, and the carrier, respectively, while denotes the th gear mesh.
Following the convention that stretching of the mesh spring is taken as positive, the meshing forces and are calculated by where is the linear tooth displacement of the th planet gear along the line of action. is the linear displacement of carrier at the planet center along the tangential direction. Notations of the variables and parameters are given in Figure 1.
2.2. Sun Gear
In wind turbines, complicated control systems are often used to keep the high speed end of the transmission operating at constant speed for the sake of power quality [29]. Thus, in this paper the sun is assumed rotating at a constant speed, and only the transverse motions are considered. The equations of motion are expressed as below:
The following two relations are derived with reference to Figure 1 and the PGT kinematics: here represents a possible phase angle between the 1st planet and the 1st rotor blade and is the rotation speed of the rotor.
2.3. Planet Gears
The rotation of a planet gear is represented as and are computed by (2) and (3).
2.4. Aerodynamic Torque
A rotor of wind turbines is shown in Figure 2. The wind goes through the rotor and generates the which can be computed by the following equation using an equivalent wind speed at the hub height of the turbine [29]: where represents the equivalent wind speed at the hub height, is the so-called power coefficient which is a function of and pitch angle , and is called tip speed ratio defined as
The wind speed is split into a mean wind speed and a turbulence as below:
The mean wind speed changes slowly; therefore, it can be treated as a deterministic quantity. While the turbulence represents a fast changing component and is regarded as a random quantity.
Substituting (9) into (7) and neglecting the term of 2nd order would give The first term in (10), a deterministic component, is the excitation (a), while the second, a random component, is the excitation (c). As pointed out in the introduction section, we focus on the responses under the random excitation (the second term) here.
The standard deviation of the turbulence is expressed as a product of the mean wind speed and a turbulence intensity as
The turbulence component is generally represented by its spectrum. There are several spectra available in the literature. One of the most commonly used spectra in the wind power industry is the Von Karmon spectrum as below [30]: where is the frequency. and are two quantities related to the mean wind speed and parameters on the site. They are computed by where and are two constants with values of 0.4 and 0.25, respectively [29].
This spectrum can be approximated by passing white noise through a shaping filter [29] with the following transfer function:
The variance of the output of the above shaping filter is unity; thus, in order to represent the turbulence completely, the following two equations can be used to simulate the turbulence : where is white noise with unit spectrum and is an added variable to the general coordinate vector.
2.5. General Equation of Motion
Choosing the augmented general coordinate vector as , and combining (1)–(6), (15), an augmented equation set in the following form is obtained. It is noted that a damping matrix has been inserted into the equation: The entries of all the matrices and vectors can be found in the Appendix.
3. Numerical Simulation and Analysis
3.1. Simulation Scheme
Many algorithms exist for numerically solving (16) to get a specific response sample. However, in random vibration, we are more interested in the statistics of the responses, rather than individual response which is often sought in the deterministic case as in many existing works [14, 31, 32]. In this paper, we look specifically for the variances or standard deviations of the responses, so that we can construct the probabilistic distribution of the responses. This can provide a much more solid base for fatigue and reliability study. To this end, we first recast (16) into the following state space form with one of the stochastic Newmark schemes [33]: where and are two matrices. Their entries and components are given in the Appendix.
Multiplying (17) by its transpose and taking ensemble average give the following recursive relation [34]: In (18), here denotes the ensemble average. The diagonal elements in contain the information we are looking for.
The meshing stiffness ( and ) in the stiffness matrix periodically change with time. The forms of and used in this paper are simplified as a rectangular form as shown in Figure 3. There is a time lag between the and mesh. For standard gears the value of is determined by the method proposed by Parker and Lin [35]. In the simulation, the rectangular form is approximated by expanding it into Fourier series as below:
The number of harmonics retained depends on the accuracy required. In this paper, 10 harmonic terms are retained. Once expansion for one mesh is made, the other meshes can be easily obtained by shifting the series to a certain phase difference.
The damping matrix in the Appendix is assumed as proportional which is a combination of the mass matrix and the stiffness matrix as . and are two constants, and and are the mass and stiffness matrices before augmentation. Their entries are given in the Appendix. In the simulation, values of 0.23 and 0.005 are assumed for and respectively. Generally the value of damping has little effect on the stable response, rather it mainly affects the time for the response to reach stable. Through simulation, this is proven that this combination works well giving relatively accurate simulation results yet not too long a time to reach the stable response.
3.2. Simulation Parameters
The PGT simulated in this paper is of spur gear type with three planet gears. The parameters of the PGT are listed in Table 1 [11], while parameters relating to the turbine and wind are given in Table 2.
The tooth number of the gears is related to the radius by the following equations for SI gears. ; here or . represents the module of the gears and is the tooth number of the gear. This can be found in any well-written text on gears, such as the one by Jelaska [36]. Also, it is straight forward to obtain the meshing frequency through kinematic analysis as [26].
3.3. Statistics of Displacement Response
With (18) the variances and covariances of the state space variables at the descritized time instants can be computed recursively starting from initial conditions. In this paper, we take the zero initial values as and . Simulation shows that stable (not stationary) status is quickly reached at the assumed damping level.
For the sun gear, the transverse displacement is represented by two components along the horizontal () and the vertical directions (). Combining the two together makes the overall displacement as below: where denotes the overall transverse displacement of the sun gear. Taking ensemble average of the above equation leads to
This indicates that the variance of the transverse displacement of the sun equals the sum of the variances of the displacements along and directions. Figure 4 shows the variances of the sun gear’s displacements after the stable status is reached. In Figures 4(a) and 4(b), the meshing period is recognizable, but the responses are not exactly periodic in meshing frequency. This is understandable if one takes a look at (4). The parameter contained in these two equations is clearly related with the rotor period. However, in Figure 4(c) shows almost perfect meshing frequency without noticeable effect of the rotor frequency. Another observation is that the displacement variance of the sun gear () is not stationary even though the excitation is stationary. Rather it changes with the meshing frequency. This is reasonable given that the meshing stiffness in the and meshes changes with time.
(a) Variance of , sun gear coordinate
(b) Variance of , sun gear coordinate
(c) Variance of , sun gear overall displacement
The variances of the rotation of the planet gears and the carrier, and , are presented in Figures 5 and 6, respectively. Similar to the sun gear, meshing frequency is obviously seen in these figures. Among the three planets, any adjacent two have a phase difference of 1/3 meshing period (Figure 5). If we compare Figures 5 and 6, it can be seen that the fluctuation in Figure 6 is much smaller than that in Figure 5. This makes sense because the carrier is under simultaneous action of all the three planets which, in general, are not in phase. Thus, the planets are more likely to cancel out each other in excitation to the carrier. This is similar to the mesh phasing phenomenon in the deterministic case [35, 37].
(a) Variance of , displacement of 1st planet
(b) Variance of , displacement of 2nd planet
(c) Variance of , displacement of 3rd planet
3.4. Statistics of Meshing Loads
The statistics of the meshing forces are important for analyzing the reliability of the system and calculating the accumulated fatigue of the parts. These tasks can not be carried out with deterministic response samples. However, the variances of the forces are not directly contained in the matrix . They can be obtained after some straightforward manipulation. Multiplying (2) and (3) with their transpose, respectively, and then taking the ensemble average, the variances of the and meshing forces can be obtained as below:
All the terms of ensemble averages at the right hand side are obtained from the matrix , including diagonal and nondiagonal elements. While the other parameters are available from the kinematic analysis or the parameter tables. Two quantities, and , are shown in Figures 7 and 8, respectively. For the sake of saving space, variances of other meshing forces are not given here, but it is expected that a phase shift of 1/3 the meshing period exists in any two adjacent phases. Obviously, the variances of meshing forces and are not stationary; rather they change with the basic meshing frequency.
4. Conclusions and Future Work
This paper, based on a lumped parameter dynamic model, investigates the vibration and dynamics of a PGT in a wind turbine to the excitation of wind turbulence. The wind turbulence is approximated by passing white noise through a shaping filter. Then, an augmented set of governing equation of motion is obtained, in which the excitation is white noise. The equations are solved by a direct numerical recursive algorithm based on the stochastic Newmark scheme. The variances of all the responses and the meshing forces are obtained through numerical simulation. The following conclusions can be drawn through the simulation.(a)The variances of the displacements and meshing forces are not stationary under the wind turbulence. Instead, they change with time in the meshing frequency. This fact indicates that neglecting the meshing frequency, as done in most wind power simulation, may not give correct results to the dynamics.(b)The variances of meshing forces, either or , demonstrate an exact phase difference of meshing period between two adjacent meshes. This result is obtained without considering possible manufacturing errors or misalignment. If these factors are considered, the phase difference would be shifted.(c)Relatively, the fluctuation in the variance of the rotation of the carrier is much smaller than that in the rotation of the planets, even though they change with the same meshing frequency.
Dynamics of PGT in deterministic framework is well developed; publications are abundant in journals and conference proceedings. In contrast, research works from random viewpoint are very rare. This paper, as an initial effort, intends to reveal the very basic features of PGT random vibration and dynamics for wind turbine application. Future works to improve include the following.(a)In real wind turbines, PGTs are more likely to be used in combination with one or two-stage parallel shaft gears. As such, one of the future works is to account for the interaction between the PGT and possible parallel shaft gears, as well as the interaction between the whole transmission unit and the generator.(b)The rotation of the rotor brings rotation frequency to the excitation; therefore, the effect of this rotation frequency component on the response statistics needs to be investigated and clarified in the future.(c)The model used in this paper is relatively rough, with focus placed on the torsional vibration. In the future, refinement of the model to better represent the wind turbine application is needed.
Appendix
Entries of Matrices
The matrices used in the model are as follows: is a unit diagonal matrix.
Nomenclature
| , , , : | Carrier, ring, sun, and the th planet when used as subscript |
| : | Radii of base circles, used with subscripts |
| , : | Meshing forces at the th sun-planet and ring-planet mesh |
| : | Number of planets |
| : | Tooth number of sun gear |
| : | Tooth number of planet gear |
| : | Tooth number of ring gear |
| , , : | Linear displacements and angular displacement |
| : | Pressure angle of gears |
| : | Torque of rotor |
| : | Stiffness of individual part or mesh |
| : | Density of air |
| : | Swept area of the rotor |
| , : | Moment inertia and mass |
| : | Gear module |
| : | Equivalent wind speed at hub height |
| , : | Mean and turbulent wind speed |
| : | Radius of the rotor |
| : | Power coefficient |
| : | Tip speed ratio |
| : | Standard deviation of the turbulence |
| : | Turbulence intensity |
| : | Turbulence length |
| , , : | Mass, damping, and stiffness matrices |
| , : | Displacements of th planet gear tooth and the carrier along the tangential directions |
| : | General coordinate vector |
| : | Added general coordinate |
| : | Excitation vector |
| : | White noise |
| : | Meshing angular velocity |
| : | Angular velocity of rotor |
| : | Phase angle between the 1st planet and the 1st rotor blade |
| , : | Contact ratio of the sun-planet and ring-planet meshes |
| , : | Time and time step |
| : | Time lag between and mesh |
| , : | Parameters of Newmark scheme. |
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.