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

Authors: Jianming Yang; Ping Yang
Journal: Shock and Vibration
 (2016)

---

## 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.


---

## Body

## 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: Section1 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 Figure1, 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 y and z 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.Figure 1
Components of a PGT.


### 2.1. Carrier

The carrier receives the aerodynamic torque from the rotor. Its equation of motion is represented by(1)Jc+Nmprc2rc2uc¨=Trorc+∑i=1Nfspi+frpicos⁡α,where the subscripts r, p, s, and c indicate the ring, the planet, the sun, and the carrier, respectively, while i denotes the ith gear mesh.Following the convention that stretching of the mesh spring is taken as positive, the meshing forcesfspi and frpi are calculated by (2)fspi=kspi-yssin⁡ϕi+zscos⁡ϕi-upi-uccos⁡α(3)frpi=krpiupi-uccos⁡α,where upi=θpirp is the linear tooth displacement of the ith planet gear along the line of action. uc=θcrc 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: (4)msys¨=∑i=1Nfspisin⁡ϕi-ksysmszs¨=-∑i=1Nfspicos⁡ϕi-kszs.The following two relations are derived with reference to Figure1 and the PGT kinematics: (5)ϕi=ψi+αψi=ωrot+β;here β represents a possible phase angle between the 1st planet and the 1st rotor blade and ωro is the rotation speed of the rotor.


### 2.3. Planet Gears

The rotation of a planet gear is represented as(6)Jpirp2upi¨=fspi-frpi.fspi and frpi are computed by (2) and (3).


### 2.4. Aerodynamic Torque

A rotor of wind turbines is shown in Figure2. The wind goes through the rotor and generates the Tro which can be computed by the following equation using an equivalent wind speed at the hub height of the turbine [29]: (7)Tro=12ρAVh2RCpλ,where Vh represents the equivalent wind speed at the hub height, Cp is the so-called power coefficient which is a function of λ and pitch angle βr, and λ is called tip speed ratio defined as(8)λ=ωroRVh.Figure 2
Wind field.The wind speedVh is split into a mean wind speed V0 and a turbulence Vt as below: (9)Vh=V0+Vt.The mean wind speedV0 changes slowly; therefore, it can be treated as a deterministic quantity. While the turbulence Vt represents a fast changing component and is regarded as a random quantity.Substituting (9) into (7) and neglecting the term of 2nd order Vt would give (10)Tro≈12ρARCpλV02+ρARCpλV0Vt.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σt is expressed as a product of the mean wind speed V0 and a turbulence intensity It as (11)σt=ItV0.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]: (12)Sω=kv1+ωTv25/6,where ω is the frequency. kv and Tv are two quantities related to the mean wind speed and parameters on the site. They are computed by (13)Tv=ltV0kv=2Tv1-c22c12c2-c2+1-c12-1,where c1 and c2 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: (14)Hiω=kviωTvc1+1iωTv+1iωTvc2+1.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 turbulenceVt: (15)Vt=σtkvxa+kvTvc1x˙ac2Tv2x¨a+c2Tv+Tvx˙a+xa=ζ˙t,where ζ˙(t) is white noise with unit spectrum and xa is an added variable to the general coordinate vector.


### 2.5. General Equation of Motion

Choosing the augmented general coordinate vector asX=(ys,zs,up1,…,upN,uc,xa)T, 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: (16)MX¨+CX˙+KX=Ft.The entries of all the matrices and vectors can be found in the Appendix.


## 2.1. Carrier

The carrier receives the aerodynamic torque from the rotor. Its equation of motion is represented by(1)Jc+Nmprc2rc2uc¨=Trorc+∑i=1Nfspi+frpicos⁡α,where the subscripts r, p, s, and c indicate the ring, the planet, the sun, and the carrier, respectively, while i denotes the ith gear mesh.Following the convention that stretching of the mesh spring is taken as positive, the meshing forcesfspi and frpi are calculated by (2)fspi=kspi-yssin⁡ϕi+zscos⁡ϕi-upi-uccos⁡α(3)frpi=krpiupi-uccos⁡α,where upi=θpirp is the linear tooth displacement of the ith planet gear along the line of action. uc=θcrc 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: (4)msys¨=∑i=1Nfspisin⁡ϕi-ksysmszs¨=-∑i=1Nfspicos⁡ϕi-kszs.The following two relations are derived with reference to Figure1 and the PGT kinematics: (5)ϕi=ψi+αψi=ωrot+β;here β represents a possible phase angle between the 1st planet and the 1st rotor blade and ωro is the rotation speed of the rotor.


## 2.3. Planet Gears

The rotation of a planet gear is represented as(6)Jpirp2upi¨=fspi-frpi.fspi and frpi are computed by (2) and (3).


## 2.4. Aerodynamic Torque

A rotor of wind turbines is shown in Figure2. The wind goes through the rotor and generates the Tro which can be computed by the following equation using an equivalent wind speed at the hub height of the turbine [29]: (7)Tro=12ρAVh2RCpλ,where Vh represents the equivalent wind speed at the hub height, Cp is the so-called power coefficient which is a function of λ and pitch angle βr, and λ is called tip speed ratio defined as(8)λ=ωroRVh.Figure 2
Wind field.The wind speedVh is split into a mean wind speed V0 and a turbulence Vt as below: (9)Vh=V0+Vt.The mean wind speedV0 changes slowly; therefore, it can be treated as a deterministic quantity. While the turbulence Vt represents a fast changing component and is regarded as a random quantity.Substituting (9) into (7) and neglecting the term of 2nd order Vt would give (10)Tro≈12ρARCpλV02+ρARCpλV0Vt.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σt is expressed as a product of the mean wind speed V0 and a turbulence intensity It as (11)σt=ItV0.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]: (12)Sω=kv1+ωTv25/6,where ω is the frequency. kv and Tv are two quantities related to the mean wind speed and parameters on the site. They are computed by (13)Tv=ltV0kv=2Tv1-c22c12c2-c2+1-c12-1,where c1 and c2 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: (14)Hiω=kviωTvc1+1iωTv+1iωTvc2+1.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 turbulenceVt: (15)Vt=σtkvxa+kvTvc1x˙ac2Tv2x¨a+c2Tv+Tvx˙a+xa=ζ˙t,where ζ˙(t) is white noise with unit spectrum and xa is an added variable to the general coordinate vector.


## 2.5. General Equation of Motion

Choosing the augmented general coordinate vector asX=(ys,zs,up1,…,upN,uc,xa)T, 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: (16)MX¨+CX˙+KX=Ft.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]: (17)Xn+1X˙n+1=A1XnX˙n+A2FnFn+1,where A1 and A2 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]: (18)Dn+1=A1DnA1T+A2En1A2T+A1En2A2T+A1En2A2TT.In (18), (19)Dn=XnXnTXnX˙nTX˙nXnTX˙nX˙nTEn1=FnFnTFnFn+1TFn+1FnTFn+1Fn+1TEn2=XnFnTXnFn+1TX˙nFnTX˙nFn+1T.〈〉 here denotes the ensemble average. The diagonal elements in Dn contain the information we are looking for.The meshing stiffness (kspi and krpi) in the stiffness matrix K periodically change with time. The forms of kspi and krpi used in this paper are simplified as a rectangular form as shown in Figure 3. There is a time lag ts between the s-p and r-p mesh. For standard gears the value of ts 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:(20)kspi=k0+∑i=1Nkkiccos⁡iωmt+∑i=1Nkkissin⁡iωmtk0=cspkspi1kic=1π∫0csp-12π2kspi1cos⁡ixdx+∫csp-12π2πkspi1cos⁡ixdx=kspi1iπsin⁡icsp-12πkis=1π∫0csp-12π2kspi1sin⁡ixdx+∫csp-12π2πkspi1sin⁡ixdx=kspi1nπ1-cos⁡icsp-12π.Figure 3
Mesh stiffness.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 matrixC1 in the Appendix is assumed as proportional which is a combination of the mass matrix and the stiffness matrix as C1=αMM1+αKK1. αM and αK are two constants, and M1 and K1 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 αM and αK 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 Table1 [11], while parameters relating to the turbine and wind are given in Table 2.Table 1
Parameters of PGT.
Parameter Value m s 181.6 kg J c 314.7 kgm2 r s 0.102 m r c 0.293 m k s p i 1 3.55 × 109 N/m N k 10 α 24.6 ° c r p 1.5 N p 26 m g 14 mm m p 114 kg J p 51.9 kgm2 r p 0.165 m k s 1.2 × 107 N/m k r p i 1 4.56 × 109 N/m N 3 c s p 1.6 N s 16 N r 68Table 2
Parameters of the turbine and the wind.
Parameter Value V 0 15 m/s ω r o 20 rpm I t 0.16 C p 0.36 ρ 1.2041 kg/m3 R 23 m l t 600 mThe tooth number of the gears is related to the radius by the following equations for SI gears.ri=mgNi/(2cos⁡α); here i=s,p or r. mg represents the module of the gears and Ni 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 ωm=NsNr/Ns+Nrωs [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 X0=0 and X˙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 (ys) and the vertical directions (zs). Combining the two together makes the overall displacement ds as below: (21)ds2=ys2+zs2,where ds denotes the overall transverse displacement of the sun gear. Taking ensemble average of the above equation leads to (22)ds2=ys2+zs2.This indicates that the variance of the transverse displacement of the sun equals the sum of the variances of the displacements alongy and z 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 ϕi contained in these two equations is clearly related with the rotor period. However, ds 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 (〈ds2〉) 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 s-p and r-p meshes changes with time.Figure 4
Variances of the sun gear: (a) variance ofzs; (b) variance of ys; (c) variance of ds.
(a)
Variance ofys, sun gear coordinate (b)
Variance ofzs, sun gear coordinate (c)
Variance ofds, sun gear overall displacementThe variances of the rotation of the planet gears and the carrier,〈τpi2〉 and 〈τc2〉, 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].Figure 5
Variances of the planet gear rotation: (a) variance ofup1; (b) variance of up2; (c) variance of up3.
(a)
Variance ofup1, displacement of 1st planet (b)
Variance ofup2, displacement of 2nd planet (c)
Variance ofup3, displacement of 3rd planetFigure 6
Variances ofuc, carrier rotation displacement.


### 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 matrixD. 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 s-p and r-p meshing forces can be obtained as below: (23)fspi2=kspi2sin2ϕiys2+cos2ϕizs2+upi2+cos2αuc2-2sin⁡ϕicos⁡ϕiyszs+2sin⁡ϕiysupi+2sin⁡ϕicos⁡αysuc+2cos⁡ϕizsxpi-2cos⁡ϕicos⁡αzsuc-2cos⁡αupiucfrpi2=krpi2upi2+cos2αuc2-2cos⁡αupiuc.All the terms of ensemble averages at the right hand side are obtained from the matrixD, including diagonal and nondiagonal elements. While the other parameters are available from the kinematic analysis or the parameter tables. Two quantities, 〈fsp12〉 and 〈frp12〉, 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 〈fsp12〉 and 〈frp12〉 are not stationary; rather they change with the basic meshing frequency.Figure 7
Variances offsp1, meshing force at 1st s-p mesh.Figure 8
Variances offrp1, meshing force at 1st r-p mesh.


## 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]: (17)Xn+1X˙n+1=A1XnX˙n+A2FnFn+1,where A1 and A2 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]: (18)Dn+1=A1DnA1T+A2En1A2T+A1En2A2T+A1En2A2TT.In (18), (19)Dn=XnXnTXnX˙nTX˙nXnTX˙nX˙nTEn1=FnFnTFnFn+1TFn+1FnTFn+1Fn+1TEn2=XnFnTXnFn+1TX˙nFnTX˙nFn+1T.〈〉 here denotes the ensemble average. The diagonal elements in Dn contain the information we are looking for.The meshing stiffness (kspi and krpi) in the stiffness matrix K periodically change with time. The forms of kspi and krpi used in this paper are simplified as a rectangular form as shown in Figure 3. There is a time lag ts between the s-p and r-p mesh. For standard gears the value of ts 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:(20)kspi=k0+∑i=1Nkkiccos⁡iωmt+∑i=1Nkkissin⁡iωmtk0=cspkspi1kic=1π∫0csp-12π2kspi1cos⁡ixdx+∫csp-12π2πkspi1cos⁡ixdx=kspi1iπsin⁡icsp-12πkis=1π∫0csp-12π2kspi1sin⁡ixdx+∫csp-12π2πkspi1sin⁡ixdx=kspi1nπ1-cos⁡icsp-12π.Figure 3
Mesh stiffness.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 matrixC1 in the Appendix is assumed as proportional which is a combination of the mass matrix and the stiffness matrix as C1=αMM1+αKK1. αM and αK are two constants, and M1 and K1 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 αM and αK 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 Table1 [11], while parameters relating to the turbine and wind are given in Table 2.Table 1
Parameters of PGT.
Parameter Value m s 181.6 kg J c 314.7 kgm2 r s 0.102 m r c 0.293 m k s p i 1 3.55 × 109 N/m N k 10 α 24.6 ° c r p 1.5 N p 26 m g 14 mm m p 114 kg J p 51.9 kgm2 r p 0.165 m k s 1.2 × 107 N/m k r p i 1 4.56 × 109 N/m N 3 c s p 1.6 N s 16 N r 68Table 2
Parameters of the turbine and the wind.
Parameter Value V 0 15 m/s ω r o 20 rpm I t 0.16 C p 0.36 ρ 1.2041 kg/m3 R 23 m l t 600 mThe tooth number of the gears is related to the radius by the following equations for SI gears.ri=mgNi/(2cos⁡α); here i=s,p or r. mg represents the module of the gears and Ni 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 ωm=NsNr/Ns+Nrωs [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 X0=0 and X˙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 (ys) and the vertical directions (zs). Combining the two together makes the overall displacement ds as below: (21)ds2=ys2+zs2,where ds denotes the overall transverse displacement of the sun gear. Taking ensemble average of the above equation leads to (22)ds2=ys2+zs2.This indicates that the variance of the transverse displacement of the sun equals the sum of the variances of the displacements alongy and z 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 ϕi contained in these two equations is clearly related with the rotor period. However, ds 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 (〈ds2〉) 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 s-p and r-p meshes changes with time.Figure 4
Variances of the sun gear: (a) variance ofzs; (b) variance of ys; (c) variance of ds.
(a)
Variance ofys, sun gear coordinate (b)
Variance ofzs, sun gear coordinate (c)
Variance ofds, sun gear overall displacementThe variances of the rotation of the planet gears and the carrier,〈τpi2〉 and 〈τc2〉, 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].Figure 5
Variances of the planet gear rotation: (a) variance ofup1; (b) variance of up2; (c) variance of up3.
(a)
Variance ofup1, displacement of 1st planet (b)
Variance ofup2, displacement of 2nd planet (c)
Variance ofup3, displacement of 3rd planetFigure 6
Variances ofuc, carrier rotation displacement.


## 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 matrixD. 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 s-p and r-p meshing forces can be obtained as below: (23)fspi2=kspi2sin2ϕiys2+cos2ϕizs2+upi2+cos2αuc2-2sin⁡ϕicos⁡ϕiyszs+2sin⁡ϕiysupi+2sin⁡ϕicos⁡αysuc+2cos⁡ϕizsxpi-2cos⁡ϕicos⁡αzsuc-2cos⁡αupiucfrpi2=krpi2upi2+cos2αuc2-2cos⁡αupiuc.All the terms of ensemble averages at the right hand side are obtained from the matrixD, including diagonal and nondiagonal elements. While the other parameters are available from the kinematic analysis or the parameter tables. Two quantities, 〈fsp12〉 and 〈frp12〉, 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 〈fsp12〉 and 〈frp12〉 are not stationary; rather they change with the basic meshing frequency.Figure 7
Variances offsp1, meshing force at 1st s-p mesh.Figure 8
Variances offrp1, meshing force at 1st r-p mesh.


## 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, eithers-p or r-p, demonstrate an exact phase difference of 1/N 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.



---

*Source: 6292953-2016-11-24.xml*