Parametric Study of Drivetrain Dynamics of a Wind Turbine Using the Multibody Dynamics

Wind energy has become one of the most cost-effective and environmental friendly renewable energy resources among all exploited renewable energy. However, the failure of gearbox contributes most of the downtime for wind turbine system. Dynamic properties of drivetrain, including gearbox should be investigated in detail further. In present paper, a mathematical model for a horizontal axis wind turbine drivetrain was developed using the torsional multibody dynamic model. The drivetrain in this study consisted of a low-speed planetary gear stage (three identical planets with spur teeth, sun and fixed ring gears) and two high-speed spur gear stages. This typical arrangement has been commonly used in the wind turbine industry. Based on this model, this paper aims to investigate the influence of drivetrain parameters on the dynamic response of the wind turbine. The dynamic responses of the turbine with different rotor inertia and generator inertia are compared. Then the difference due to changing of the shaft stiffness is also investigated during and after the transient condition. This parametric study shows that lower rotor inertia and generator inertia could lead to more oscillations.


Introduction
The growing awareness of the threat of climate change caused by greenhouse gas has brought wind energy to the top of the global consensus, and wind energy has become the most cost-effective of all currently exploited renewable energy resources and has attracted extensive research and business interest [1][2]. Besides the industrial application success, extensive research is also carried out outside of the commercial companies [3][4][5][6].
The increased sizes of rotor and turbine have led to a complex design of the drivetrain mechanism. As one of the most critical components on the operation of wind turbines, the drivetrain is responsible for most of the downtime and results in the increment in operational costs [7]. Various studies have addressed the dynamic modeling of wind turbines [8]. However, few studies have focused on the advanced dynamic modeling of drivetrain [9][10][11][12][13].
As extension of our previous research [9], this study uses the multibody dynamics simulation techniques to extend the two-lumped-masses model into the torsional model for the drivetrain of a wind turbine, where exactly one degree of freedom (DOF) of each drivetrain component is used to investigate the torsional vibration of the drivetrain. Based on this model, dynamic analysis is carried out to study the transient response of the drivetrain structures. Numerical investigation is carried out to study the effects of the drivetrain parameters on the performance of the wind turbine.

Mathematical Model of the Wind Turbine Drivetrain
The lumped parameter dynamic model of the drivetrain investigated in this study is shown in Figure 1. With one DOF for each component, the dynamic model has ten DOFs. Each gear, i, has a polar moment of inertia of J i , i = rotor, c, p1, p2, p3, s, g1, g2g3, g4, gen, for the rotor, the carrier, the planets 1, 2, 3, the sun, the gears 1, 2 and 3, 4 and the generator, respectively; k j , j = 1, 2, 3, is the torsional stiffness for the Low Speed Shaft (LSS), Intermediate Shaft (IMS) and High Speed Shaft (HSS), respectively. The shaft connecting gears 2 and 3 ( Figure 1) was assumed to be rigid due to immediately adjacent gears. The parameter r i is the base circle radius for each gear component and r c is the radius of the circle passing through the planet centers for the carrier.
The angular displacements of the central members, θ rotor , θ c , θ s , θ g1 , θ g2g3 , θ g4 , θ gen , were defined as absolute displacements about a fixed coordinate frame, whereas the angular displacements of planets, θ cp1 , θ cp2 , θ cp3 , were defined relative to the carrier c as:  The linear springs with time-varying stiffness k rp (t), k sp (t), k g12 (t) and k g34 (t) were used to model the gear mesh [14] ( Figure 2). Each mesh stiffness is expressed in the Fourier series as: where k rp and k sp are mean values. The Fourier coefficients are given as: where C sp , C rp are sun-planet contact ratio and ring-planet contact ratio; γ sp , γ rp are phase angle of sun-planet gear pair and ring-planet gear pair. Similar mesh stiffness is defined for the parallel gear stage, except that there are different gear mesh frequencies (GMFs) and contact ratios. The mesh stiffness for the ring-planet, the sun-planet and the gears 1-2 and 3-4 are presented in Figure  3.
The GMF for each gear mesh pair can be determined from kinematic relationships as: The relative gear mesh displacements of the sun-planet and the ring-planet gear mesh are defined as: where r r is the radius of the ring gear. A fixed ring gear is assumed in this study, which results in θ r = 0. Lagrange's equation was used to derive the EOM [15]. The Lagrangian of the undamped dynamic models shown in Figure 2 with three planets is given as (Kahraman 2001): The Lagrange EOMs of the ten-DOF dynamic model are given as: where a dot over the parameter denotes differentiation with respect to time and ∂ denotes a partial derivative. Eq. (14) together with Eq. (13) leads to eight coupled ordinary differential equations, which can be written in the matrix form as: Mu+Ku=Q ɺɺ (15) where u, M, K, and Q are the displacement vector, the mass matrix the stiffness matrix, and the external force.
The aerodynamic torque and electromagnetic torque are considered as external forces.

Numerical Analysis Method
Equations of motion were solved using the direct integration method to obtain the transient response of the wind turbine drivetrain. Solution at the next time step was calculated from the solutions at the previous time steps. The unconditional stable Newmark-β method (Figure 4) was used for the numerical analysis by implementing it into the Matlab based on the Eqs. 13~15. The basic integration equations that relate the positions, the velocities and the accelerations from the time step t to t+∆t are given as follows: Typical parameters of a wind turbine drivetrain from Todorov (2009) are used in this study ( Table 1). The rotor is exited with the angular velocity of 18 rpm.

Effects of Rotor and Generator Inertia Parameters
To investigate the effects of the drivetrain parameters, seven cases are analyzed with the proposed multibody dynamic model. Table 2 presents the cases studied with different rotor inertia. The generator inertia and shaft stiffness are kept as constant.  Figure 5 shows the rotational speed of the rotor, planet 1, sun gear and generator, respectively. It is observed that the rotor inertia has great influence on the rotor and planet 1 and negligible influence on sun gear and generator. Not only the magnitude but also the oscillation period of rotor speed change due to changing of the rotor inertia. The effect of the rotor inertia on the gear contact forces are investigated ( Figure 6) and the results show that there is a large increase of the contact force of the ring-planet and sun-planet gear pair but no effect on the parallel gear contact force.   Table 3 summarizes the simulation cases for the generator inertia parameters in which the rotor inertia and shaft stiffness remain constant. Figure 7 shows the response of three cases using different generator inertias. No effect is found on the rotor speed. But the higher generator inertia reduces the planet 1, sun gear and generator speed magnitude but increases the oscillation periods. Figure 8 shows the effect of the generator inertia on different gear pair contact force. Oscillation periods of all the contact forces are affected by generator inertia but the magnitude remain unchanged.    Table 4~6 summarizes considered cases for different LSS, IMS, and HSS stiffness parameters, respectively. LSS stiffness only has a small influence on the rotor speed while negligible effects can be found for other components ( Figure  9). With increasing of the LSS stiffness, the contact force between ring gear and planet gear and the contact force between sun gear and planet gear increase but no effect on the parallel gear stage ( Figure 10). Due to the different IMS stiffness, the oscillation period of the planet 1, sun gear and generator speed change significantly ( Figure 11). Similar effect on the oscillation period can be identified for all the gear contact forces ( Figure 12). The oscillation period of the planet 1, sun gear and generator change a little because of the increasing of the HSS stiffness ( Figure 13). The reduction of the oscillation period for all the gear stages can be found with increasing of HSS stiffness ( Figure 14).

Conclusions
A mathematical model of the drivetrain in a horizontal wind turbine was proposed. The model used multibody dynamics to investigate the torsional vibration in the drivetrain. The governing equation was derived using Lagrange's equation. The model took into account the flexibility of the gear mesh by using linear springs with timevarying stiffness and shaft torsion.
The equations of motion of the drivetrain were solved numerically by the direct integration method to obtain the transient response of the system. The parameter studies were carried out based on the proposed model. The following conclusions can be drawn from this study: 1) The generator inertia has most significant effect on the response of the components. 2) Higher generator inertia can reduce the oscillation range and delay the oscillation period. 3) LSS stiffness has negligible effects while IMS and HSS stiffness can reduce the period of the component speed response.

4) Rotor inertia and LSS stiffness have a large effect on
the contact force magnitude of the ring-planet and sunplanet gear pair but limited effect on the parallel gear stage contact force.