Dynamic characteristics analysis of planetary gear system with internal and external excitation under turbulent wind load

Abstract

According to the working characteristics of a 1.5 MW wind turbine planetary gear system under complex and random wind load, a two-parameter Weibull distribution model is used to describe the distribution of random wind speed, and the time-varying load caused by random wind speed is obtained. The nonlinear dynamic model of planetary gear transmission system is established by using the lumped parameter method, and the relative relations among various components are derived by using Lagrange method. Then, the relative relationship between the components is solved by Runge Kutta method. Considering the influence of random load and stiffness ratio on the planetary gear transmission system, the nonlinear dynamic response of cyclic load and random wind load on the transmission system is analyzed. The analysis results show that the variation of the stiffness ratio makes the planetary gear have abundant nonlinear dynamics behavior and the planetary gear can get rid of chaos and enter into stable periodic motion by changing the stiffness ratio properly on the premise of ensuring transmission efficiency. For the variable pitch wind turbine, the random change of external load increases the instability of the system.

Keywords

Turbulent wind load planetary gear nonlinear dynamic characteristics stiffness

Introduction

Planetary gear system is a key component of wind turbines. It is inconvenient for real-time monitoring and detection because the external load is time-varying and the working environment is in the engine room which is set up in high altitude. In addition, due to the time-varying internal and external excitation, the gear transmission system generates vibration and noise during operation, so it is particularly important to study the dynamic characteristics of the gearbox transmission system under random wind load more accurately.

At present, domestic and foreign scholars have conducted a lot of research on the dynamic characteristics of gear transmission systems. Kahraman¹ established a nonlinear torsional vibration model of single-stage planetary gear transmission and analyzed the torsional vibration characteristics of planetary gear transmission. Al-shyyab et al.² established a nonlinear torsional vibration model of a single-stage planetary gear system, solved the dynamic response of the system by using the incremental harmonic balance method, and compared the analytical solution with the numerical solution. Sun et al.³ found that there are abundant strong nonlinear dynamic behaviors in the torsional vibration system of planetary gears while ignoring the meshing stiffness fluctuation and eccentric error of each gear pair. Then, scholars have carried out in-depth research on the influence of internal and external excitation on the gear system. Peeters et al.⁴ established a flexible multi-body dynamics model for the wind turbine transmission system under variable wind loads with the help of multi-body dynamics software, and studied the dynamic characteristics of the system. Liew and Lim⁵ considered the influence of the time-varying stiffness of the rolling bearing on the transmission system, and analyzed the response of the time-varying stiffness of the rolling bearing to the coupled vibration of a single-stage parallel shaft gear-rolling bearing. Qin et al.⁶ used trapezoidal waves to express the time-varying characteristics of meshing stiffness, analyzed the phase relationship in the meshing process of planetary gears, and calculated the dynamic response of the three-stage planetary gear reducer of a shield machine at different meshing positions. Zhou et al.⁷ considered the coupling of gear and rolling bearing, established a nonlinear dynamic model for the transmission system of wind turbine, and obtained the dynamic response of the system. Chen et al.⁸ analyzed the nonlinear dynamic response characteristics of wind turbine drive system from the aspect of backlash theoretically. Sun et al.⁹ obtained the nonlinear frequency response characteristics of the system by adopting the harmonic balance method on the basis of considering time-varying meshing stiffness, comprehensive meshing error, and backlash. Chen et al.¹⁰ took the comprehensive transmission error and external load of gears as random variables and considered the influence of time-varying meshing stiffness, backlash, and other factors, established a single pair of gear torsion dynamic model and analyzed the statistical characteristics of each response quantity and dynamic meshing force of the system. Li et al.¹¹ established a pure torsion model for the transmission system of wind turbines by taking into account nonlinear factors such as time-varying meshing stiffness, comprehensive meshing error of gears, and backlash, and analyzed the dynamic response of the system with backlash, excitation frequency, and damping ratio. Zhou et al.¹² analyzed the inherent characteristics of the system from the combination of bending and torsion coupling, considering the gravity of the gear and supporting bearings. Kahnamouei and Yang¹³ established a planetary gear hybrid power model considering stiffness variation and backlash nonlinearity. Based on the statistical linearization method, the nonlinear backlash was linearized and the system was simulated. Chen et al.¹⁴ established a gear nonlinear model with random wind load. The dynamic characteristics of the system and the influence of the backlash on the response of the system were analyzed.

The linear and nonlinear characteristics of planetary gears have been studied from different perspectives and some important conclusions have been obtained in the above literatures. In this paper, a direct drive single-stage planetary gear transmission system is taken as the research object. Based on the comprehensive consideration of strong nonlinear factors such as time-varying meshing stiffness, backlash, and comprehensive meshing error, the time-varying wind speed in the natural wind field is simulated by combining the average and variance of wind speed provided in a certain area, the dimension-dynamic equation of the planetary gear transmission system in generalized coordinates is derived, the Runge Kutta method is used to solve the dynamic equation. Furthermore, the nonlinear dynamic response characteristics of the system under cyclic load and two-parameter Weibull distribution random wind speed are discussed, and the influence of the above excitation factors on the system characteristics is analyzed.

Torsional dynamic model of planetary gear transmission system

Lumped mass model of planetary gear

This paper takes the 1.5 MW direct-drive wind turbine planetary gear transmission as the research object, adopts the structure form of planet carrier input, sun gear output, fixed ring gear, and uses three planetary gears to split the power and bear the load. The structure diagram of the planetary gear transmission system is shown in Figure 1. The gear model in the text is a straight tooth meshing, so the axial load is very small and can be ignored.

Figure 1.

1.5 MW torsional dynamics model of planetary gear for direct drive wind turbines.

In addition, due to the large bearing stiffness, the comparative analysis of the natural mode calculation results of planetary gear in Kahraman¹ shows that when the ratio of the bearing support stiffness to the gear tooth meshing stiffness is greater than 10, the influence of the bearing support stiffness in planetary gear can be ignored. Since the main task of the transmission system in direct drive wind turbine is to undertake the transmission of torsional power, and the structural load is mainly borne by the casing and bearing. Therefore, only the dynamic characteristics of the torsion direction of the planetary gear are considered here. In order to study the gear drive system of variable pitch variable speed wind turbine, the external and internal excitations caused by the change of wind speed and output torque are considered. In order to simplify the model, the planetary gear system is regarded as a lumped parameter system, and the meshing part of the gear teeth is transformed into an undamped spring, and the inner ring gear is completely fixed without considering the elastic deformation of the inner ring gear. It is considered that the mass, moment of inertia, meshing stiffness and bearing support stiffness of the three planetary gears are identical, and the eccentricity error and roundness error are ignored. The nonlinear factors such as torsional flexibility of the planet carrier, tooth surface friction, and component mass eccentricity are ignored. The force analysis of the sun gear S, planet gear P, and planet carrier C is carried out, and the dynamic equation is listed as shown in equation (1).

{\begin{matrix} J_{s} \frac{d^{2} θ_{s}}{d t^{2}} + \sum_{i = 1}^{N} (F_{sp i} + F_{sp i}^{d}) . r_{bs} = T_{S} \\ J_{pi} \frac{d^{2} θ_{pi}}{d t^{2}} - (F_{sp i} + F_{sp i}^{d}) r_{bp} + (F_{rp i} + F_{rp i}^{d}) . r_{bp} = 0 \\ (J_{c} + \sum_{i = 1}^{N} m_{p i} r_{c}^{2}) \frac{d^{2} θ_{c}}{d t^{2}} - \sum_{i = 1}^{N} [(F_{sp i} + F_{sp i}^{d}) + (F_{rp i} + F_{rp i}^{d})] . r_{bc} = - T_{P} \end{matrix}

(1)

where J_s, J_pi, and J_c are the moment of inertia of the sun gear, planet gear, and planet carrier respectively; θ_s, θ_pi, and θ_c are the angles of the sun gear, planet gear, and planet carrier respectively; $F_{sp i}$ and $F_{sp i}^{d}$ represent the elastic restoring force and damping force of the meshing pair between the sun gear and the ith planetary gear as shown in equation (2); $F_{rp i}$ and $F_{rp i}^{d}$ represent the elastic restoring force and damping force of the meshing pair between the planet carrier and the ith planetary gear as shown in equation (2); r_bs, r_bpi, and r_br are the base circle radius of the sun gear, planetary gear, and internal gear respectively; r_c is the radius of planet carrier; m_s, m_pi, and m_c are the mass of sun gear, planet gear, and planet carrier respectively; N is the number of planet gears.

{\begin{matrix} F_{sp i} = k_{sp i} f ({\bar{x}}_{sp i}, {\bar{b}}_{sp i}) \\ F_{rp i} = k_{rp i} f ({\bar{x}}_{rp i}, {\bar{b}}_{rp i}) \\ F_{sp i}^{d} = c_{sp i} \frac{d {\bar{x}}_{sp i}}{d t} \\ F_{rp i}^{d} = c_{rp i} \frac{d {\bar{x}}_{rp i}}{d t} \end{matrix}

(2)

where k_spi and k_rpi are the stiffness between the sun gear, the inner ring gear, and the ith planetary gear respectively; $f ({\bar{x}}_{sp i}, {\bar{b}}_{sp i})$ and $f ({\bar{x}}_{rp i}, {\bar{b}}_{rp i})$ are the backlash function; c_spi and c_rpi are the meshing damping between the sun gear and the ith planetary gear.

Due to the gap between the gear teeth, the constraint of the planetary gear system is incomplete, so equation (1) represents a positive semidefinite second-order nonlinear system. In order to eliminate the influence of rigid body displacement and indefinite solution on the system, the number of equations is reduced and the dynamic transfer error is introduced as shown in equation (3).

{\begin{matrix} x_{sp i} = x_{s} - x_{p i} - x_{c} - e_{sp i} (t) \\ x_{rp i} = x_{p i} - x_{c} - e_{rp i} (t) \end{matrix}

(3)

where x_spi and x_rpi are the relative displacements on the meshing line between the sun gear, the inner ring gear, and the ith planetary gear respectively; e_spi and e_rpi are the integrated transfer errors of the external and internal meshing pairs respectively; b_sp is the half backlash of sun gear and planetary gear pair; b_rp is the half backlash between the inner ring gear and the planetary gear pair; N is the number of planetary gears; x_c is the projection of the planetary line displacement on the corresponding meshing line, that is, x_c = r_cθ_c cos α, Introducing the equivalent base circle radius of the planet carrier r_bc = r_c cos α, then x_c = r_bcθ_c; Therefore, the linear displacement x generated by the angular displacement of the sun gear and the planetary gear on the corresponding meshing line is shown in equation (4).

{\begin{matrix} x_{s} = r_{bs} θ_{s} \\ x_{pi} = r_{bpi} θ_{pi} \end{matrix}

(4)

Transform equation (1) and obtain $\frac{d^{2} {\bar{x}}_{s}}{d t^{2}}$ , $\frac{d^{2} {\bar{x}}_{pi}}{d t^{2}}$ , and $\frac{d^{2} {\bar{x}}_{c}}{d t^{2}}$ respectively as shown in Appendix B. Then, substituting $\frac{d^{2} {\bar{x}}_{s}}{d t^{2}}$ , $\frac{d^{2} {\bar{x}}_{pi}}{d t^{2}}$ , and $\frac{d^{2} {\bar{x}}_{c}}{d t^{2}}$ into equation (3) to obtain $\frac{d^{2} {\bar{x}}_{sp i}}{d t^{2}}$ and $\frac{d^{2} {\bar{x}}_{rp i}}{d t^{2}}$ , as shown in Appendix B.

where $m_{s} = J_{s} / r_{bs}^{2}$ , $m_{p i} = J_{p i} / r_{bp i}^{2}$ , $m_{c} = J_{c} + \sum_{i = 1}^{N} m_{p i} r_{c}^{2} / r_{bc}^{2}$ ; supposing that T_p and T_s are the input and output torque of the system, T_p = P₀/ω, T_s = T_p/j, ${\bar{F}}_{L} = T_{P} / r_{bs}$ , ${\bar{F}}_{D} = T_{S} / r_{bs}$ , planetary gear ratio j = 5.44; r_bp is the base circle radius of the planetary gear; b is the half backlash of the gear pair of the road, and b_sp or b_rp is taken according to the inside and outside of the gear pair.

The differential equation in the Appendix B (equation (B.2)) is dimension, and the dimensionless program shown in the following equation is obtained. Introducing dimensionless-time τ = ω_pt, where $ω_{p} = \sqrt{k_{ms 1} / m_{1}}$ (k_ms1 is the amplitude of the external meshing stiffness), let ${\bar{x}}_{sp i, rp i} = c x_{sp i, rp i}$ , ${\bar{b}}_{sp i, rp i} = c b_{sp i, rp i}$ , $m_{1} = m_{p i} m_{s} / m_{p i} + m_{s}$ , $m_{2} = m_{p i} m_{r} / m_{p i} + m_{r}$ , $ω_{ss} = \sqrt{k_{ms 1} / m_{s}}$ , $ω_{sc} = \sqrt{k_{ms 1} / m_{c}}$ , $ω_{sp} = \sqrt{k_{ms 1} / m_{p i}}$ (i = 1, 2, 3), $ω_{rc} = \sqrt{k_{mr 1} / m_{c}}$ , $ω_{rp} = \sqrt{k_{mr 1} / m_{p i}}$ , $c_{sp i} = 2 m_{1} ω_{1} ζ_{sp i}$ , $c_{rp i} = 2 m_{2} ω_{2} ζ_{rp i}$ , $ω_{1} = \sqrt{k_{ms 1} / m_{1}}$ , $ω_{2} = \sqrt{k_{mr 1} / m_{2}}$ , $F_{L} = {\bar{F}}_{L} / c m_{c}$ , $F_{D} = {\bar{F}}_{D} / c m_{s}$ . In the following text, b and e are used directly for the assignment of backlash and transmission error. c is 1/2b, then the dimensionless equation is obtained as shown in equation (5).

{\begin{matrix} \frac{d^{2} {\bar{x}}_{sp i}}{d t^{2}} + \frac{ω_{ss}^{2}}{ω_{p}^{2}} \sum_{i = 1}^{N} k_{sp i} (τ) f (x_{sp i}, b_{sp i}) + \frac{2 m_{1} ω_{1}}{m_{s} ω_{p}} \sum_{i = 1}^{N} ξ_{sp i} \frac{d {\bar{x}}_{sp i}}{d t} + \frac{ω_{sc}^{2}}{ω_{p}^{2}} \sum_{i = 1}^{N} k_{sp i} (τ) f (x_{sp i}, b_{sp i}) \\ + \frac{2 m_{1} ω_{1}}{m_{c} ω_{p}} \sum_{i = 1}^{N} ξ_{sp i} \frac{d {\bar{x}}_{sp i}}{d t} + \frac{ω_{rc}^{2}}{ω_{p}^{2}} \sum_{i = 1}^{N} k_{rp i} (τ) f (x_{rp i}, b_{rp i}) + \frac{2 m_{2} ω_{2}}{m_{c} ω_{p}} \sum_{i = 1}^{N} ξ_{rp i} \frac{d {\bar{x}}_{rp i}}{d t} \\ + \frac{ω_{sp}^{2}}{ω_{p}^{2}} k_{sp i} (τ) f (x_{sp i}, b_{sp i}) + \frac{2 m_{1} ω_{1}}{m_{p i} ω_{p}} ξ_{sp i} \frac{d {\bar{x}}_{sp i}}{d t} - \frac{ω_{rp}^{2}}{ω_{p}^{2}} k_{rp i} (τ) f (x_{rp i}, b_{rp i}) - \frac{2 m_{2} ω_{2}}{m_{p i} ω_{p}} ξ_{rp i} \frac{d {\bar{x}}_{rp i}}{d t} \\ = \frac{F_{D}}{ω_{p}^{2}} + \frac{F_{L}}{ω_{p}^{2}} - \frac{d^{2} {\bar{e}}_{sp i} (t)}{d t^{2}} \\ \frac{d^{2} {\bar{x}}_{rp i}}{d t^{2}} + \frac{ω_{sc}^{2}}{ω_{p}^{2}} \sum_{i = 1}^{N} k_{sp i} (τ) f (x_{sp i}, b_{sp i}) + \frac{2 m_{1} ω_{1}}{m_{c} ω_{p}} \sum_{i = 1}^{N} ξ_{sp i} \frac{d {\bar{x}}_{sp i}}{d t} + \frac{ω_{rc}^{2}}{ω_{p}^{2}} \sum_{i = 1}^{N} k_{rp i} (τ) f (x_{rp i}, b_{rp i}) \\ + \frac{2 m_{2} ω_{2}}{m_{c} ω_{p}} \sum_{i = 1}^{N} ξ_{rp i} \frac{d {\bar{x}}_{rp i}}{d t} - \frac{ω_{sp}^{2}}{ω_{p}^{2}} k_{sp i} (τ) f (x_{sp i}, b_{sp i}) - \frac{2 m_{1} ω_{1}}{m_{p i} ω_{p}} ξ_{sp i} \frac{d {\bar{x}}_{sp i}}{d t} + \frac{ω_{rp}^{2}}{ω_{p}^{2}} k_{rp i} (τ) f (x_{rp i}, b_{rp i}) \\ + \frac{2 m_{2} ω_{2}}{m_{p i} ω_{p}} ξ_{rp i} \frac{d {\bar{x}}_{rp i}}{d t} = \frac{F_{L}}{ω_{p}^{2}} - \frac{d^{2} {\bar{e}}_{rp i} (t)}{d t^{2}} \end{matrix}

(5)

Turbulent wind load

At present, there are many models used to describe wind speed distribution, such as Rayleigh distribution, two parameter Weibull distribution, three parameter Weibull distribution, lognormal distribution, Pearson curve fitting, and so on. A large number of studies show that the wind energy curve fitted by the two parameter Weibull density function is closer to the actual wind speed distribution. Therefore, this paper uses the two parameter Weibull distribution to simulate the change of local wind speed.

According to Liu et al.¹⁵ and Ye,¹⁶ the two parameter Weibull distribution curve is expressed as a unimodal curve, the expression of probability density function is shown in equation (6), the calculation equation of mathematical expectation (average wind speed) of wind speed is shown in equation (7), and the variance of wind speed is shown in equation (8).

F (v) = \frac{k}{c} {(\frac{v}{c})}^{k - 1} \exp [- {(\frac{v}{c})}^{k}]

(6)

E (v) = \int_{0}^{+ \infty} v . F (v) dv = c . Γ (1 + \frac{1}{k})

(7)

Inverse transformation of the above equations (6) and (7), one can get the equation (8).

D (v) = \int_{0}^{+ \infty} {[v - E (v)]}^{2} f (v) dv = c^{2} {Γ (\frac{2}{k} + 1) - {[Γ (\frac{1}{k} + 1)]}^{2}}

(8)

where k is the shape parameter; c is the scale parameter; v is the wind speed; Г is the Gamma function.

According to the measured wind speed of a wind farm, we can obtain the average wind speed E(v) = 12m/s and the wind speed variance D(v) = 9m²/s². By substituting the maximum average wind speed and the wind speed variance into the corresponding equation, the shape parameters k = 4.5065, and c = 13.1485 can be obtained.

The shape and scale parameters of Weibull distribution can be obtained by equations (6)–(8), Invoking the Weibull distribution random number generation function rad (1) in Matlab to generate a random number between [0, 1] to get random wind speed data. The random wind speed time series can be obtained by matching the simulated wind speed data with the sampling time points, the two parameter Weibull random wind speed model in 10 min is generated as shown in Figure 2.

Figure 2.

Time history diagram of random wind speed within 600 s.

At present, variable pitch variable speed wind turbines are mainly used for large-scale wind power generation above megawatt level. For variable pitch variable speed wind turbines, the natural wind speed can be divided into four regions according to Vladimír and Tomáš,¹⁷ as shown in Figure 3.

Figure 3.

Variable speed wind turbine operating area.

where v_{cut in} is the cut in wind speed of the wind turbine, v_rate is the rated wind speed of the wind turbine, v_{cut off} is the cut out wind speed of the wind turbine, and the power absorbed by the wind wheel from the wind can be expressed as

P = \frac{1}{2} C_{p} A ρ v^{3}

(9)

where P is the input power of the wind turbine main shaft; C_p is the wind energy utilization coefficient, which is known as C_p = 0.55 from the literature⁹; A is the sweeping area of the wind turbine, A = πR² = 3847 m²; R is the blade radius R = 35 m; ρ is the air density; v is the wind speed which unit is m/s.

The tip speed ratio is a very important parameter used to describe the characteristics of a wind turbine. It refers to the ratio of the tip speed of the wind turbine blade to the wind speed. According to its definition, the expression can be written as follows:

λ = \frac{ω R}{v}

(10)

where λ is the tip speed ratio; ω is the rotor speed.

It can be known from the structure of the wind power generator that the rotational speed, torque, and output power of the wind turbine impeller are the input rotational speed, input torque, and input power of the wind power gear box. By substituting equation (10) into equation (9), the input power of the wind turbine gearbox can be rewritten as equation (11).

P = \frac{1}{2} ρ π R^{5} \frac{C_{p}}{λ^{3}} ω^{3}

(11)

When the wind speed is higher than the cut-in wind speed and lower than the rated wind speed, the input power P of the wind-electric gearbox can also be determined by the following equation:

P = T_{in} ω

(12)

If the wind speed is higher than the fan cut-in wind speed and lower than the rated wind speed, the input torque of the speed increaser as a function of wind speed, input speed, and wind speed can be expressed by equation (13).

{\begin{matrix} T_{in} = \frac{T_{rate}}{v_{rate}^{2}} . v^{2} \\ ω = \frac{ω_{rate}}{v_{rate}} . v \end{matrix}

(13)

where T_rate is the rated torque, N/m; v_rate is the rated wind speed, m/s; ω_rate is the rated speed, r/min.

The random wind speed intercept in 120 s is shown in Figure 4.

Figure 4.

Random wind speed time-history diagram in 120 s.

The data provided by a wind power plant is shown in Table 1. From equations (11)–(13), the curves of input torque and wind speed, input speed, and wind speed of 1.5 MW wind turbine speed increaser are shown in Figure 5. At the same time, the input torque and speed curves of 1.5 MW fan speed increaser planetary gear train corresponding to external time-varying wind speed within 120 s can be obtained as shown in Figure 6.

Table 1.

Technical parameters of wind.

Rated powerP_rate (MW)	Rated speedω_rate(r/min)	Cut in wind speedv_{cut in}(m/s)	Cut out windspeed v_{cut off}(m/s)	Planetary gear ratio j
1.65	12	3	24	5.44

Figure 5.

Input torque and input speed change curve with wind speed: (a) input torque and wind speed and (b) input speed and wind speed.

Figure 6.

Random wind input speed/torque change curve in 120 s: (a) curve of random wind input torque variation in 120 s and (b) curve of random wind input speed change in 120 s.

Internal excitation analysis of planetary gear system

Calculation of time-varying meshing stiffness

There are many methods to calculate the meshing stiffness of gears. Vladimír and Tomáš¹⁷ use the empirical formula and finite element method provided by ISO standard to calculate the meshing stiffness of gears respectively. Al-shyyab and Kahraman,² Hbaieb et al.,¹⁸ and Wang¹⁹ considered the time-varying meshing stiffness as a function of time and assumed that the period of rectangular wave is constant, which is effective in analyzing the natural characteristics of the planetary gear steady-state motion. Due to the constant change of external wind speed, the wind turbine has frequent acceleration and deceleration which causes the speed of the speed increaser to change frequently. If the time-varying meshing stiffness is regarded as a fixed period function of time, it is obviously unable to reflect the change of time-varying meshing stiffness in the process of acceleration and deceleration, so it is different from the real meshing stiffness. In order to be more suitable for engineering application, an improvement is made on the basis of Pan et al.²⁰ The specific method is to express the time-varying meshing stiffness as a mode of adding the mean value and the fluctuation value.

k (θ_{pi}) = \bar{k} + Δ k (θ_{pi}) (i = 1, 2, . . ., n)

(14)

where θ_pi is the rotation angle of the planetary gear around its rotation axis in the relative inertial coordinate system, and the period of the time-varying meshing stiffness is a constant value, which is only related to the number of teeth z_p of the planetary gear.

T_{p} = \frac{p_{bp}}{r_{bp}} = \frac{π m_{p} \cos α_{p}}{\frac{1}{2} m_{p} z_{p} \cos α_{p}} = \frac{2 π}{z_{p}}

(15)

where P_bp is the pitch of planetary gear base circle, r_bp is the radius of planetary gear base circle, m_p is the module of planetary gear, α_p is the pressure angle of planetary gear indexing circle.

The wave component of time-varying meshing stiffness is expanded by Fourier series, and equation (16) is obtained.

Δ k (θ_{pi}) = \sum_{in = 1}^{\infty} [a_{n} \cos (n z_{p} θ_{pi}) + b_{n} \sin (n z_{p} θ_{pi})]

(16)

where

a_{0} = \frac{z_{p}}{2 π} \int_{θ_{0}}^{θ_{0} + \frac{z_{p}}{2 π}} Δ k (θ_{pi}) d θ_{pi} = 0

(17)

a_{n} = \frac{z_{p}}{π} \int_{θ_{0}}^{θ_{0} + \frac{z_{p}}{2 π}} Δ k (θ_{pi}) \cos (n z_{p} θ_{pi}) d θ_{pi}

(18)

\begin{matrix} b_{n} = \frac{z_{p}}{π} \int_{θ_{0}}^{θ_{0} + \frac{z_{p}}{2 π}} Δ k (θ_{pi}) \sin (n z_{p} θ_{pi}) d θ_{pi} \\ n = 1, 2, 3, . . . \end{matrix}

(19)

where θ₀ is the initial phase of time-varying mesh stiffness.

Calculation of meshing error

The meshing error of gear system is generally described by the sum of geometric error and intercept error. Qin et al.⁶ and Vladimír and Tomáš¹⁷ assume that the meshing error is a sine wave of meshing frequency, and then use Fourier series to fit the meshing error. For planetary gear systems, due to the existence of internal and external meshing at the same time, the meshing error includes not only static meshing error, but also eccentricity error of sun gear and planetary gear, collectively called the integrated error of the planetary gear system, which is one of the important internal excitations. Assuming that the meshing error is a sine function of the rotation angle of the planetary gear, the meshing errors between the sun gear and the planetary gear, and the inner ring gear and the planetary gear are shown in equation (20).

{\begin{matrix} e_{sp i}^{'} (t) = A_{sp i} \sin (ω t + φ_{s i}) \\ e_{rp i}^{'} (t) = A_{rp i} \sin (ω t + φ_{ri} + γ_{sr}) \end{matrix}

(20)

where $e_{sp i}^{'} (t)$ and $e_{rp i}^{'} (t)$ are the meshing errors between the ith sun gear and the planetary gear, the planetary gear and the internal gear respectively; A_spi is the transmission error variation amplitude between the ith sun gear and the planetary gear; A_rpi is the transmission error variation amplitude between the ith planetary gear and the internal gear; $φ_{si}$ is the initial phase of the transmission error change between the ith sun gear and the planetary gear; $φ_{ri}$ is the initial phase of the transmission error change between the ith planetary gear and the internal gear; γ_sr is the phase between the internal and external tooth meshing pairs in the same path.

According to the description in Li et al.,¹¹ the phase difference of different external meshing pairs is $φ_{s i} = z_{s} ψ_{i}$ , the phase difference of different internal meshing pairs is $φ_{r i} = z_{r} ψ_{i}$ , where $ψ_{i} = 2 π (i - 1) / N$ is the position angle of the ith planetary gear relative to the first planetary gear. According to the parity of planetary gears, γ_sr is taken as 0 or π respectively.

Calculation of backlash

In actual engineering, due to the existence of the backlash, a certain impact will be generated during gear meshing transmission, which is not conducive to the stability and reliability of the system. The value of this clearance is generally very small. Kahraman and Singh²¹ found that the gear side clearance can be described in the form of a piecewise function as shown in equation (21).

f (X, b) = {\begin{matrix} X - b X > b \\ 0 - b \leq X \leq b \\ X + b X < - b \end{matrix}

(21)

where X is the relative displacement between the sun gear and the ith planetary gear and the ring gear and the ith planetary gear in the direction of the meshing line; b is the half-backlash of the gear pair, b = b₀ + Δb, b₀ is the initial backlash, which is the amount of backlash change caused by tooth surface wear, which is equal to the amount of tooth surface wear in value, and b_sp or b_rp is taken according to the inside and outside of the gear pair.

Simulation calculation and analysis

Basic parameters

In this paper, the planetary gear transmission system of 1.5 MW direct drive wind turbine is taken as the model, and the geometric parameters of the gear transmission system of wind turbine are shown in Table 2.

Table 2.

Geometric parameters of wind turbine gear transmission system.

	Sun gear	Equivalent 1	Planet gear	Equivalent 2	Ring gear
Tooth number	25	—	43	—	111
Helix angle β (°)	0	—	0	—	0
Module (mm)	12	—	12	—	12
Pressure angle α(°)	20°	—	20°	—	20°
Center distance (mm)	—	408	—	924	—
Mass (kg)	211	—	241	—	1363
Contact ratio	—	1.8	—	1.8	—
Average meshstiffness (N/m)	—	1.52e9	—	1.8e9	—
Backlash (mm)	—	0.2	—	0.2	—

System simulation and results

Effect of stiffness ratio on system response

In this paper, the Runge Kutta method is used to solve equation (5). Figures 7 to 12 show the nonlinear dynamic characteristics of the system under different stiffness ratios and external loads. In Vladimír and Tomáš,¹⁷ it can be found that x_spi and x_rpi have the same bifurcation characteristics. In order to save space, this paper only studies the bifurcation characteristics of x_spi.

Figure 7.

Bifurcation diagram and 3D spectrum diagram of the system under cyclic load with the stiffness ratio k_i = 1: (a) bifurcation diagram and (b) three dimensional spectrum.

Figure 8.

Phase diagram, Poincaré cross-section, and spectrum diagram under cyclic load with the stiffness ratio k_i = 1: (a) phase, (b) Poincaré map and (c) FFT spectrum.

Figure 9.

Bifurcation diagram and 3D spectrum diagram of the system under cyclic load with the stiffness ratio k_i = 1.5: (a) bifurcation diagram and (b) three dimensional spectrum.

Figure 10.

Phase diagram, Poincaré cross-section, and spectrum diagram under cyclic load with the stiffness ratio k_i = 1.5: (a) phase, (b) Poincaré map and (c) FFT spectrum.

Figure 11.

Bifurcation diagram and 3D spectrum diagram of the system under turbulent wind load with the stiffness ratio k_i = 1.5: (a) bifurcation diagram and (b) three dimensional spectrum.

Figure 12.

Phase diagram, Poincaré cross-section and spectrum diagram under turbulent wind load with the stiffness ratio k_i = 1.5: (a) phase, (b) Poincaré map and (c) FFT spectrum.

Under the action of cyclic load: when the stiffness ratio k_i = 1, it can be seen from the bifurcation diagram that when the dimensionless excitation speed ω of the system is low, the motion of the system is more complex, with single cycle, quasi cycle, double cycle, and multi cycle state changes. When ω = 3.25, chaotic motion becomes the main form of motion state. When the rotational speed is further increased to 3.5, the system is in stable harmonic 1T-periodic motion; Figure 8 is the phase diagram, Poincaré section, and spectrum diagram of the system when ω = 1.5. From the phase diagram (Figure 8(a)), it can be seen that the system response is shown as a limit cycle; Since the ordinate is relatively small, the Poincaré section (Figure 8(b)) is a set of points. When the coordinates are the same as those in Figures 10 and 12, it should be a point; At the same time, the main frequency in the two-dimensional spectrum of the system (Figure 8(c)) is f_m, which indicates that the system is in the state of single periodic motion.

Under the action of cyclic load: when the stiffness ratio k_i = 1.5, it can be seen from the bifurcation diagram that the quasi-periodic motion of the system is advanced to ω = 2.75, the quasi-periodic motion is obviously advanced, and the system quickly enters into the single periodic motion, the results show that increasing the stiffness ratio properly can make the system enter a stable one-cycle cycle quickly. However, when the speed is in the range of [3.2, 3.5], the system appears quasi-bifurcated motion. As the speed further increases to 3.5, the system is in a stable simple harmonic single periodic motion, and there are some chaotic fluctuations; Figure 10 is the phase diagram, Poincaré cross section, and spectrum diagram of the system when ω = 1.5. From the phase diagram (Figure 10(a)), it can be seen that the system response is multiple closed curves, Poincaré cross section (Figure 10(b)) is represented as a ring composed of multiple concentration points, and the frequency diagram (Figure 10(c)) has obvious peaks, but there are broadband components, which is preliminarily judged as quasi periodic motion.

Above all, when the stiffness ratio k_i = 1, ω = 1, 1.75, 3.25 are the key points for the system to enter chaos. In the actual work, the working speed of the system will have a small fluctuation; when the stiffness ratio k_i = 1.5, ω = 1.5, 1.75, 2.75, and 3.25 are the key points for the system to enter chaos. Thus, increasing the stiffness ratio can shorten the time of chaotic motion and lengthen the interval of periodic motion, quasi periodic motion, and chaotic motion appear alternately. Moreover, it is similar to the conclusion in Li et al.,²² which proves the correctness of the model and the research content in this paper.

Effect of turbulent wind load on system response

It can be seen from the three-dimensional spectrum that there are obvious meshing frequency f_m and composite frequency of sun gear and planet gear in Figures 7(b), 9(b) and 11(b). The peak value at position 4 is larger and the fluctuation of load torque has a great influence on the relative displacement of sun gear and planetary gear. With the increase of rotating speed, the high frequency component f_m first increases, then decreases, and then tends to be stable, which shows that there is no obvious change in the frequency component of the system response in the normal working range.

When the stiffness ratio k_i = 1.5 under turbulent wind load as shown in Figure 12. The phase diagram (Figure 12(a)) is still a number of closed curves, Poincaré map section (Figure 12(b)) is a ring formed by discrete points, the set of points is continuously diffused and fused, the system tends to be unstable, and the frequency diagram (Figure 12(c)) has obvious peaks, but there are also broadband components. Thus, it can be determined that the system is in quasi periodic motion.

When the stiffness ratio remains unchanged at k_i = 1.5, under the action of turbulent wind load, its motion is very similar to that under the action of cyclic load from the bifurcation diagram of the system. However, it can be seen that the motion state gradually changes from quasi-periodic to chaotic, and the turbulent wind load increases the instability of the system. At the same time, it is corresponding to the conclusion that random excitation increases the instability of the system obtained in Chen et al.,²³ which proves the correctness of the study.

Conclusion

In this paper, the wind speed model of random wind field based on Weibull distribution is established. According to the relationship between the actual operation area and wind speed of variable wind speed wind turbine, the actual input load of wind turbine gearbox caused by variable wind speed is obtained. Considering the influence of time-varying meshing stiffness of gears and internal excitation of planetary gear meshing error comprehensively, the nonlinear dynamic model of planetary gear transmission system is established by using lumped parameter method. The Runge-Kutta method is used to solve the dynamic equation of the system, and the meshing characteristics of the system with different stiffness ratio under cyclic load and variable wind load are calculated. The results show the following regularity

For planetary gear system, increasing the stiffness ratio can shorten the time of chaotic motion and lengthen the interval of periodic motion, quasi periodic motion, and chaotic motion appear alternately. From the point of view of suppressing chaos, increasing the stiffness ratio of the system can make the system gradually get rid of chaos and enter stable periodic motion. The transmission system enters chaotic motion from periodic motion. The evolution path from periodic motion to quasi periodic motion and finally to chaos is a conventional path. When the actual equipment works, it should avoid complex motion areas such as quasi periodic motion and chaotic motion, so as to ensure the normal and reliable operation of the system in the stable periodic motion area; the displacement response at the meshing point of the transmission system is closely related to the response of the external load, and has a similar variation law with the time-varying torque. The larger the external excitation, the larger the vibration displacement; when the stiffness ratio is increased properly, the motion of the system under cyclic load is very similar to that under turbulent wind load. However, the motion state gradually changes from quasi-periodic to chaotic, and the turbulent wind load increases the instability of the system.

Footnotes

Appendix A

Appendix B

Some equations of gear dynamics.

(B.1)

{\begin{matrix} \frac{d^{2} {\bar{x}}_{s}}{d t^{2}} = - \frac{1}{m_{s}} \sum_{i = 1}^{N} k_{sp i} (t) f ({\bar{x}}_{sp i}, {\bar{b}}_{sp i}) - \frac{1}{m_{s}} \sum_{i = 1}^{N} c_{sp i} \frac{d {\bar{x}}_{sp i}}{d t} + \frac{F_{D}}{m_{s}} \\ \frac{d^{2} {\bar{x}}_{pi}}{d t^{2}} = \frac{1}{m_{p i}} k_{sp i} (t) f ({\bar{x}}_{sp i}, {\bar{b}}_{sp i}) + \frac{1}{m_{p i}} c_{sp i} \frac{d {\bar{x}}_{sp i}}{d t} - \frac{1}{m_{p i}} k_{rp i} (t) f ({\bar{x}}_{rp i}, {\bar{b}}_{rp i}) - \frac{1}{m_{p i}} c_{rp i} \frac{d {\bar{x}}_{rp i}}{d t} \\ \frac{d^{2} {\bar{x}}_{c}}{d t^{2}} = \frac{1}{m_{c}} \sum_{i = 1}^{N} k_{sp i} (t) f ({\bar{x}}_{sp i}, {\bar{b}}_{sp i}) + \frac{1}{m_{c}} \sum_{i = 1}^{N} c_{sp i} \frac{d {\bar{x}}_{sp i}}{d t} + \frac{1}{m_{c}} \sum_{i = 1}^{N} k_{rp i} (t) f ({\bar{x}}_{rp i}, {\bar{b}}_{rp i}) \\ + \frac{1}{m_{c}} \sum_{i = 1}^{N} c_{rp i} \frac{d {\bar{x}}_{rp i}}{d t} - \frac{F_{L}}{m_{c}} \end{matrix}

(B.2)

{\begin{matrix} \frac{d^{2} {\bar{x}}_{sp i}}{d t^{2}} + \frac{1}{m_{s}} \sum_{i = 1}^{N} k_{sp i} (t) f ({\bar{x}}_{sp i}, {\bar{b}}_{sp i}) + \frac{1}{m_{s}} \sum_{i = 1}^{N} c_{sp i} \frac{d {\bar{x}}_{sp i}}{d t} + \frac{1}{m_{c}} \sum_{i = 1}^{N} k_{sp i} (t) f ({\bar{x}}_{sp i}, {\bar{b}}_{sp i}) \\ + \frac{1}{m_{c}} \sum_{i = 1}^{N} c_{sp i} \frac{d {\bar{x}}_{sp i}}{d t} + \frac{1}{m_{c}} \sum_{i = 1}^{N} k_{rp i} (t) f ({\bar{x}}_{rp i}, {\bar{b}}_{rp i}) + \frac{1}{m_{c}} \sum_{i = 1}^{N} c_{rp i} \frac{d {\bar{x}}_{rp i}}{d t} + \frac{1}{m_{p i}} k_{sp i} (t) f ({\bar{x}}_{sp i}, {\bar{b}}_{sp i}) \\ + \frac{1}{m_{p i}} c_{sp i} \frac{d {\bar{x}}_{sp i}}{d t} - \frac{1}{m_{p i}} k_{rp i} (t) f ({\bar{x}}_{rp i}, {\bar{b}}_{rp i}) - \frac{1}{m_{p i}} c_{rp i} \frac{d {\bar{x}}_{rp i}}{d t} = \frac{F_{D}}{m_{s}} + \frac{F_{L}}{m_{c}} - \frac{d^{2} {\bar{e}}_{sp i} (t)}{d t^{2}} \\ \frac{d^{2} {\bar{x}}_{rp i}}{d t^{2}} + \frac{1}{m_{c}} \sum_{i = 1}^{N} k_{sp i} (t) f ({\bar{x}}_{sp i}, {\bar{b}}_{sp i}) + \frac{1}{m_{c}} \sum_{i = 1}^{N} c_{sp i} \frac{d {\bar{x}}_{sp i}}{d t} + \frac{1}{m_{c}} \sum_{i = 1}^{N} k_{rp i} (t) f ({\bar{x}}_{rp i}, {\bar{b}}_{rp i}) + \frac{1}{m_{c}} \sum_{i = 1}^{N} c_{rp i} \frac{d {\bar{x}}_{rp i}}{d t} \\ - \frac{1}{m_{p i}} k_{sp i} (t) f ({\bar{x}}_{sp i}, {\bar{b}}_{sp i}) - \frac{1}{m_{p i}} c_{sp i} \frac{d {\bar{x}}_{sp i}}{d t} + \frac{1}{m_{p i}} k_{rp i} (t) f ({\bar{x}}_{rp i}, {\bar{b}}_{rp i}) + \frac{1}{m_{p i}} c_{rp i} \frac{d {\bar{x}}_{rp i}}{d t} = \frac{F_{L}}{m_{c}} - \frac{d^{2} {\bar{e}}_{rp i} (t)}{d t^{2}} \end{matrix}

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is supported by Natural Science Foundation of Ningxia Hui Autonomous Region in 2020(2020AAC03279) and Key research and development projects of Ningxia Hui Autonomous Region in 2020(2020BDE03014)

ORCID iD

Hexu Yang

Author biographies

Hexu Yang, born in 1982, received the BE degree from College of Mechanical Engineering and Automation, North West Agriculture and Forestry University, China, in 2005, and master degree from College of Mechanical Engineering and Automation, Northeastern University, China, in 2017. He is PhD candidate at Northeastern University. His research are nonlinear dynamics and wind turbine. Tel: 13895022927; E-mail: yanghxneu@hotmail.com.

Xiaopeng Li, born in 1976, is currently a Professor of School of Mechanical Engineering and Automation, Northeastern University. He received his PhD degree from Northeastern University in 2006. He specializes in machine design theory and methodology, mechanical dynamics, friction-induced vibration. He is the assistant director of the Institute of Mechanical Design and Theory in Northeastern University. The e-mail id of author is xpli@me.neu.edu.cn.

Jinchi Xu, born in 1994, received the BE degree from College of Mechanical Engineering and Automation, Northeastern University, China, in 2016, and master degree from College of Mechanical Engineering and Automation, Northeastern University, China, in 2019. He is PhD candidate at Northeastern University. His research interests are nonlinear dynamics and PI control strategy. Tel: 15541568088; E-mail: xjc19940612@163.com.

Renzhen Chen, born in 1987, received the BE degree from College of Mechanical Engineering and Automation, Northeastern University, China, in 2010, and master degree from College of Mechanical Engineering and Automation, Northeastern University, China, in 2017. He is PhD candidate at Northeastern University. His research interests are nonlinear dynamics and control of robot. Tel: 18804027618; E-mail: crz871025@126.com.

Zemin Yang, born in 1993, received the master degree from College of Mechanical Engineering and Automation, Northeastern University, China, in 2019, is PhD candidate at Northeastern University. His research interests are variable stiffness actuator and mechanical seals. Tel: 13840338295; E-mail: 18689847260@163.com.

References

Kahraman

. Natural modes of planetary Gear trains. J Sound Vib 1994; 173(1): 125–130.

Al-Shyyab

Kahraman

. A non-linear dynamic model for planetary gear sets. Proc IMechE, Part K: J Multi-body Dynamics 2007; 221: 567–576.

Sun

Shen

. Nonlinear dynamics of 2K-H planetary gear transmission. J Tsinghua Univ 2003; 43(5): 636–639.

Peeters

JLM

Vandepitte

Sas

. Analysis of internal drive train dynamics in a wind turbine. Wind Energy 2006; 9(1–2): 141–161.

Liew

Lim

. Analysis of time-varying rolling element bearing characteristics. J Sound Vib 2005; 283(3–5): 1163–1179.

Qin

Xiao

Wang

. Dynamic characteristics of multi-stage planetary gear transmission of shield machine reducer based on meshing phase analysis. J Mech Eng 2011; 47(23): 20–29.

Zhou

Qin

Yang

Chen

. Dynamic characteristics of gear bearing coupling in planetary gear transmission system of wind turbine under variable load. J Chongqing Univ 2012; 35(12): 7–14.

Chen

Zhang

Liu

, et al. Nonlinear dynamic response analysis of planetary transmission system of wind turbine gearbox considering backlash. Acta Solar Energy Sin 2016; 37(6): 1587–1593.

Sun

Shen

sun

, et al. Solution of nonlinear dynamic equation and analysis of dynamic characteristics of planetary gear transmission. J Mech Eng 2002; 3: 11–15.

10.

Chen

Qin

, et al. Influence analysis of random internal and external excitation on dynamic characteristics of gear system. China Mech Eng 2013; 24(4): 533–537.

11.

Zhu

Bao

, et al. Study on torsional nonlinear vibration modeling and motion bifurcation characteristics of planetary gear train. Chin J Mech Eng 2011; 47(21): 76–83.

12.

Zhou

Wang

, et al. Dynamic modeling of wind power gearbox transmission system. J Northeast Univ 2014; 35(9): 1301–1305.

13.

Kahnamouei

Yang

. Development and verification of a computationally efficient stochastically linearized planetary gear train model with ring elasticity. Mech Mach Theory 2021; 155: 104061.

14.

Chen

Wang

Gao

, et al. Dynamic characteristics of wind turbine gear transmission system with random wind and the effect of random backlash on system stability. Proc IMechE, Part C: J Mechanical Engineering Science 2017; 231: 2590–2597.

15.

Liu

Wang

Zhao

, et al. Induction generator for wind power generation. Beijing: China Electric Power Press, 2009, pp.42–76.

16.

. Control technology of wind turbine. Beijing: China Machine Press, 2006, pp.20–25.

17.

Vladimír

Tomáš

. Notes to design of the cylindrical gears with high contact ratio (HCR). J Middle Eur Construct Des Cars 2006; 4(1): 1–5.

18.

Hbaieb

Chaari

Fakhfakh , et al. Dynamic stability of a planetary gear train under the influence of variable meshing stiffnesses. Proc IMechE, Part D: J Automobile Engineering 2006; 220: 1711–1725.

19.

Wang

. Theoretical and experimental study on dynamics of spur gear transmission based on phase tuning. Tianjin University, 2005.

20.

Pan

Sun

, et al. Dynamic modeling and analysis of spur gear. J Dyn Control 2018; 16(2): 121–128.

21.

Kahraman

Singh

. Non-linear dynamics of a geared rotor-bearing system with multiple clearances. J Sound Vib 1991; 144(3): 469–506.

22.

Pan

. Influence of fractal backlash on dynamic behavior of gear-bearing system. J Mech Eng 2018; 54(1): 153.

23.

Chen

Qin

, et al. Dynamic characteristics of gear transmission system subjected to random internal and external excitation. Zhongguo Jixie Gong/China Mech Eng 2013; 24(4): 533–537.