Random Vibration and Dynamic Analysis of a Planetary Gear Train in a Wind Turbine

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) inwind turbines under the excitation of wind turbulence.The turbulence is representedby theVonKarmon 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.


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][5][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][8][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][11][12][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 loadsharing 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][27][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.

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.

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.

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  ro is the rotation speed of the rotor.

Planet Gears. The rotation of a planet gear is represented as
and   are computed by ( 2) and (3).

Aerodynamic Torque.
A rotor of wind turbines is shown in Figure 2. The wind goes through the rotor and generates the  ro 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  0 and a turbulence   as below: The mean wind speed  0 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  0 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. V and  V are two quantities related to the mean wind speed and parameters on the site.They are computed by where  1 and  2 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.

General Equation of Motion.
Choosing the augmented general coordinate vector as  = (  ,   ,  1 , . . .,   ,   ,   )  , 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.

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  1 and  2 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  1 in the Appendix is assumed as proportional which is a combination of the mass matrix and the stiffness matrix as  1 =    1 +   1 .  and   are two constants, and  1 and  1 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.

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.  =     /(2 cos ); here  = ,  or .  represents the module of the gears and   is the tooth number of the gear.This can be found in any wellwritten 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].

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  0 = 0 and Ẋ0 = 0. 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 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 (⟨ 2  ⟩) 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.
The variances of the rotation of the planet gears and the carrier, ⟨ 2  ⟩ and ⟨ 2  ⟩, 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].

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, ⟨ 2 1 ⟩ and ⟨ 2 1 ⟩, 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 ⟨ 2 1 ⟩ and ⟨ 2 1 ⟩ are not stationary; rather they change with the basic meshing frequency.

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 1/ 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.

Entries of Matrices
The matrices used in the model are as follows: Phase angle between the 1st planet and the 1st rotor blade   ,   : Contact ratio of the sun-planet and ringplanet meshes , : Time and time step   : Time lag between  −  and  −  mesh , : Parameters of Newmark scheme.

Figure 4 :
Figure 4: Variances of the sun gear: (a) variance of   ; (b) variance of   ; (c) variance of   .

Table 2 :
Parameters of the turbine and the wind.
+   (  −   ) cos  Carrier, ring, sun, and the th planet when used as subscript :R a d i io fb a s ec i r c l e s ,u s e dw i t h subscripts   ,   : Meshing forces at the th sun-planet and ring-planet mesh :N u m be ro fp l a n e t s   : [ [  ∑ =1   sin 2   +   −  ∑ =1   sin   cos    1 sin  1 ⋅ ⋅ ⋅   sin    ∑ =1   cos   cos   ∑ =1   cos 2   +   − 1 cos  1 ⋅ ⋅ ⋅ −  cos   −  ∑ =1   cos   cos Nomenclature, , , :