Effect of Turbulent Wind Conditions on the Dynamic Characteristics of a Herringbone Planetary Gear System of a Wind Turbine

Abstract

Due to complex environmental factors, the gear transmission systems of wind turbines are continuously affected by large torque load excitation with periodic and random properties. This paper shares the load-sharing and dynamic characteristics of a herringbone planetary gear system applied in a wind turbine. The gear dynamic model is established using a typical lumped parameter method, in which the nonlinear transmission errors of the gear pairs and left and right-side coupling stiffness of the herringbone gears are included. With the help of the blade element momentum theory, the precise calculation of the hub load of the wind turbine, which is the external excitation of the gear system, is implemented, in which the wind shear, tower shadow, turbulent effect, and tip loss correction are taken into consideration. The nonlinear dynamic characteristics of the system are obtained using the Runge-Kutta method and then discussed. The results show that the turbulent effect plays a major role in the impact on the load-sharing characteristics, and a reasonable set of the support stiffness of rotational components can improve the load-sharing characteristics of the system. The purpose of this research is to provide some useful references in numerical modelling and methods for designers and researchers of wind turbine transmission systems.

1. Introduction

Currently, herringbone planetary gear systems have been developed and widely applied in wind turbine, aero-engine, helicopter propeller, and ship propeller transmission systems because of their excellent carrying capacity. In wind power engineering applications, the stability of the gear transmission system is always the focus of research because of the complex and changeable external environmental load. Therefore, the load analysis and dynamic characteristics of planetary transmission systems have been widely studied in the past thirty years. The research of early scholars focused on the linear modelling of the gear transmission system and the effect of the fixed value excitation on the dynamic characteristics [1,2,3]. At the same time, experimental and measurement techniques of planetary gear dynamics had been developed to some extent [4,5]. Some interesting mechanisms and features of dynamics can be revealed and found using the nonlinear dynamic theory with the discovery of experiments and the development of dynamics. Thus, many experts have begun to pay attention to the nonlinearity of the system and the influence of internal/external excitation nonlinearity.

Gou et al. [6] focused on the nonlinear effect of a planetary gear system of a wind turbine caused by gravity. Through the nonlinear dynamic model, including the gravity, time-varying meshing stiffness, clearance of the bearing, and nonlinear contact of the tooth surface, a source of planetary bearing failures and its solution strategy were discovered. Zhao et al. [7] studied the nonlinear vibration features of a gear system with different pitting faults. Xiang et al. [8] and Zhang et al. [9] revealed the crack and wear failure mechanisms of a multi-stage gear transmission system of a wind turbine, respectively. Zhao et al. [10] analyzed the effect of transmission error on the planetary gear system of a wind turbine using a multi-gap gear finite element model. Because of the heavy load capacity of the herringbone gear, this structure was developed and applied to wind turbines. Currently, there are many outstanding scholars who have made outstanding contributions to the dynamic research of herringbone planetary gear. Mo et al. [11,12,13] conducted a series of investigations on a herringbone planetary gear system to reveal the influence of different structural parameters on the system’s dynamic characteristics. Wang et al. [14] and Hou et al. [15] investigated the effect of a flexible pin and friction excitation on a herringbone planetary gear system using the finite element method. Wang et al. [16] presented a herringbone planetary gearbox failure model in a wind turbine to reveal the pitting failure mechanism. The results of the research indicated that precise predictions could be obtained using the description of precise excitation considered in the dynamic model. Xu et al. [17] investigated the nonlinear behavior of herringbone planetary gear transmission systems with double-sided meshing impact and found that reducing the backlash ensured stable system operation. Wang [18] established the relationship between transmission system characteristics and gear modification. The optimum modification parameters of the axial and tooth shape were determined, which could achieve the purpose of vibration and noise reduction. Wu et al. [19] investigated the dynamic characteristics of a herringbone gear system with uncertain parameters. Although there are plenty of excellent studies on the nonlinear dynamic characteristics of herringbone planetary gear systems, dynamics research on herringbone planetary gear systems were applied to wind turbines. In addition, the external load of a wind turbine gearbox in service is more complex than that of the general gearbox. In actual wind power engineering, the load acting on the hub, which is the external excitation of the gear transmission system, has strong periodicity and randomness because of the wind field atmospheric environment and wind turbine structure itself, especially turbulent wind.

Thus, actual wind load and time-varied load excitation have gradually gained the attention of researchers. Chen et al. [20] presented a multi-source external environmental load modelling method and investigated resonance identification and dynamic responses under electromechanical coupling conditions. Wang et al. [21] studied the dynamic characteristics of both conventional and compact wind turbine gearboxes, including load-sharing and fatigue damage assessment. Tan et al. [22,23] performed a series of developments on an external environmental load acting on a gear transmission system of a wind turbine, revealing the effect of complex load conditions on the system, such as the voltage dip and platform motion and time-varying wind load. Abo et al. [24] presented a turbulent wind speed model under different influence factors, including wind shear, tower shadow, and turbine inertia. Wang et al. [25] articulated the modelling of wind and wave coupling for a floating wind turbine drivetrain. Porté-Agel et al. [26] summarized recent experimental, computational, and theoretical research efforts for wind farm prediction. Doagou et al. [27] established an erosion model of a wind turbine blade, fully considering random liquid impact. Da Silva et al. [28] studied the nonlinear dynamics of a floating offshore wind turbine, considering turbulent wind and irregular waves and focusing on the natural frequency of the system during surges. Bangga et al. [29,30] established an aerodynamic model of a wind turbine load to obtain the exact load of the system using the CFD and BEM methods under different wind conditions. The randomness of wind speed can be precisely described using the above studies. However, the calculation of hub load is mostly simplified using an approximate expression according to the Bates theory, and wind farm environments, blades, the wind shear effect, and the tower shadow effect were also ignored. As a consequence, the periodicity of the hub load was inexact.

Therefore, on the basis of the literature review above, the description of the precise hub load of a wind turbine is the first purpose of this paper. A calculation model of the hub load is established using the blade element momentum theory (BEM), which includes the wind shear effect, tower shadow effect, and turbulent effect. The coordinate transformation method is used to describe the wind speed uniformly at different locations. At the same time, the second purpose is to emphasize the importance of the turbulent wind condition in herringbone planetary gear modelling and dynamic characteristics. Furthermore, it is also one of the purposes of this paper to find out the cause of the fluctuation of the system’s dynamic characteristics in the external excitation and analyze the solution. It is hoped that this research can provide some useful reference for the design and manufacture of the herringbone planetary gear systems of wind turbines.

2. Modeling and Theories

2.1. Research Framework and Assumption

In this study, a herringbone planetary gear transmission system of a megawatt wind turbine was employed as the object of this study. The key elements of this study are shown in Figure 1, which mainly includes different wind conditions, an impeller, and a gear system. In the process of modeling, the key elements can be simplified into three subsystem models, called the wind velocity model, the load calculation model of the hub, and the dynamic model of the herringbone planetary gear system. Firstly, the turbulent wind can be calculated, involving the wind shear effect, tower shadow effect, and turbulent effect. Then, the hub load can be calculated as the turbulent wind load based on the BEM. Finally, the dynamic characteristics can be acquired using the dynamic model of a herringbone planetary gear transmission system considering the exact input load.

The assumptions for a dynamics model of a herringbone planetary gear system are as follows:

The input torque is directly loaded onto the planetary carrier. The herringbone gear system is considered to be a combination of two helical gear systems. The planetary gears are evenly arranged, and each member is rigid. The cross-angle of the herringbone gears due to manufacturing errors is not considered. The support bearing is reduced to a spring damping system. In this study, 4.5 m/s was the cut-in wind velocity, and the cut-out wind velocity was set to 25 m/s. The rated wind speed of the impeller of a wind turbine was 12 m/s.

2.2. A Dynamic Model of the Gear System

2.2.1. Motion Equations

A dynamic model of the planetary gear dynamics model in the OXY plane is shown in Figure 2. There were three kinds of coordinate systems. The motions of all parts of the hybrid wheel system were described in a uniform fixed coordinate system, OXYZ, which is the first kind of coordinate system. The centre of the planetary carrier was the origin of the coordinate system. The X-axis points to the initial position of planetary gear 1. The Z-axis goes out along the axis. The follow-up coordinate system of the carrier was the second, which was Oxyz. Its initial position coincides with OXYZ. The third coordinate system was o_pix_piy_piz_pi of the planetary gear ($i=1,2,\cdots N$). N is the number of the planetary gear. The horizontal and vertical coordinates are the radial direction and tangential direction of the planetary gear. The z-axis sets out along the axis. The position angles of the planetary gears change with the increase in time. $x, y, z, \theta$ stands for the translational and rotational displacement of each component. r, s, c, and pi are used as their subscript, and they respectively denote the ring, sun gear, carrier, and planetary gears. L and R are used as the positional labels of the herringbone gears. Thus, the column vector X can be expressed using the following equation:

$$X=\left\{\begin{array}{c}{x}_s^L,y_s^L,z_s^L,{\theta }_s^L,x_s^R,y_s^R,z_s^R,{\theta }_s^R\\ x_{pi}^L,y_{pi}^L,z_{pi}^L,{\theta }_{pi}^L,x_{pi}^R,y_{pi}^R,z_{pi}^R,{\theta }_{pi}^R\\ x_r^L,y_r^L,z_r^L,{\theta }_r^L,x_r^R,y_r^R,z_r^R,{\theta }_r^R\\ x_c,y_c,z_c,{\theta }_c\end{array}\right\}$$

(1)

In order to intuitively display the meshing relationship, the meshing dynamic model between the sun and planetary gear is shown in Figure 3a. Hypothetically, the direction of the sun gear output torque is clockwise. ${\phi }_{spi}$ can be written using the following equation:

$${\phi }_{spi}={\phi }_{pi}+{\alpha }_{spi},$$

(2)

where ${\phi }_{pi}$ denotes the position angle of the ith planetary gear. ${\alpha }_{spi}$ represents the transverse pressure angle. According to the geometric relationship, the projection of their meshing displacements along the meshing action line can be deduced as follows:

$$\begin{array}{l}{\delta }_{spi}^L={\mathit{V}}_{spi}^L{\mathit{X}}_{spi}^L-e_{spi}^L\left(t\right)\\ {\delta }_{spi}^R={\mathit{V}}_{spi}^R{\mathit{X}}_{spi}^R-e_{spi}^R\left(t\right)\end{array},$$

(3)

where $e_{spi}^L\left(t\right)$ and $e_{spi}^R\left(t\right)$ denote the meshing errors of the external meshing gear. ${\mathit{X}}_{spi}^L$ and ${\mathit{X}}_{spi}^R$ are the column vectors for the meshing element:

$$\begin{array}{l}{\mathit{X}}_{spi}^L=\left(x_s^L,y_s^L,z_s^L,{\theta }_s^L,x_{pi}^L,y_{pi}^L,z_{pi}^L,{\theta }_{pi}^L\right)\\ {\mathit{X}}_{spi}^R=\left(x_s^R,y_s^R,z_s^R,{\theta }_s^R,x_{pi}^R,y_{pi}^R,z_{pi}^R,{\theta }_{pi}^R\right)\end{array}.$$

(4)

The meshing vectors ${\mathit{V}}_{spi}^L$ and ${\mathit{V}}_{spi}^R$ can be written as:

$$\begin{array}{l}{\mathit{V}}_{spi}^L={\left(\sin {\varphi }_{spi}\cos {\beta }_b,\cos {\varphi }_{spi}\cos {\beta }_b,-\sin {\beta }_b,-r_s\cos {\beta }_b,\sin {\alpha }_{spi}\cos {\beta }_b,-\cos {\alpha }_{spi}\cos {\beta }_b,\sin {\beta }_b,r_{pi}\cos {\beta }_b\right)}^T\\ {\mathit{V}}_{spi}^R={\left(\sin {\varphi }_{spi}\cos {\beta }_b,\cos {\varphi }_{spi}\cos {\beta }_b,-\sin {\beta }_b,r_s\cos {\beta }_b,\sin {\alpha }_{spi}\cos {\beta }_b,-\cos {\alpha }_{spi}\cos {\beta }_b,\sin {\beta }_b,-r_{pi}\cos {\beta }_b\right)}^T\end{array}$$

(5)

where $r_s$ and $r_{pi}$ stand for the radii of the pitch circle, which represents the sun and planetary gear. ${\alpha }_{spi}$ denotes the transverse pressure angle. ${\beta }_b$ is the base helix angle, and its expression is provided in the following equation:

$$\begin{array}{c}{\beta }_b=\arccos \left(\sqrt{1-\left(\sin \beta {\cos }^2{\alpha }_{spi}\right)}\right)\end{array},$$

(6)

where $\beta$ is the helix angle.

If the parameters of the left and right sides are the same, superscripts will be not used. Based on Newton’s second law, the motion differential equations of the sun gear can be expressed as the following equations:

$$\begin{array}{c}{m}_s{\ddot{x}}_s^q=-\left(K_{spi}^q{\delta }_{spi}^q+C_{spi}^q{\dot{\delta }}_{spi}^q\right)sin{\varphi }_{spi}cos{\beta }_b-K_s{x}_s^q-C_s{\dot{x}}_s^q\\ m_s{\ddot{y}}_s^q=-\left(K_{spi}^q{\delta }_{spi}^q+C_{spi}^q{\dot{\delta }}_{spi}^q\right)cos{\varphi }_{spi}cos{\beta }_b-K_s{y}_s^q-C_s{\dot{y}}_s^q\\ m_s{\ddot{z}}_s^q=\left(K_{spi}^q{\delta }_{spi}^q+C_{spi}^q{\dot{\delta }}_{spi}^q\right)sin{\beta }_b-K_s{z}_s^q-C_s{\dot{z}}_s^q-K_{sa}^{LR}\left(z_s^q-z_s^p\right)-C_{sa}^{LR}\left({\dot{z}}_s^q-{\dot{z}}_s^p\right)\\ I_s{\ddot{\theta }}_s^q=\left(K_{spi}^q{\delta }_{spi}^q+C_{spi}^q{\dot{\delta }}_{spi}^q\right)r_{s}cos{\beta }_b-K_{sb}^{LR}\left(v_s^q-v_s^p\right)-C_{sb}^{LR}\left({\dot{v}}_s^q-{\dot{v}}_s^p\right)+T_{out}\end{array},$$

(7)

where $q=L,R$ and $p=R,L$. $K_{spi}$ denote the meshing stiffness of the sun gear and planetary gear. $C_{spi}$ represents their meshing damping. $K_s$ is the radial support stiffness of the sun gear. $C_s$ stands for the radial support damping of the sun gear. $K_{sa}$ denotes its axial support stiffness. $C_{sa}$ represents its axial support damping. $K_{sb}$ denotes its torsional stiffness. $C_{sb}$ represents its torsional damping. $K_{sa}^{LR}$ denotes the coupling axial support stiffness between the left-end helix gear and the right helix gear. $C_{sa}^{LR}$ represents their coupling axial support damping. $K_{sb}^{LR}$ denotes the coupling torsional stiffness between the left-end gear and the right gear. $C_{sb}^{LR}$ denotes their coupling torsional damping. $v_s^q$ is the equivalence displacement in the tangent direction of the corresponding base circle of the sun gear:

$$v_s^q=r_{bs}\cdot {\theta }_s^q,(q=L,R).$$

(8)

Relatively, the motion differential equations of the planetary gears can be expressed as the following:

$$\begin{array}{c}{m}_p{\ddot{x}}_{pi}^q=\left(K_{spi}^q{\delta }_{spi}^q+C_{spi}^q{\dot{\delta }}_{spi}^q\right)sin{\alpha }_{pi}cos{\beta }_b-K_p{x}_{pi}^q-C_p{\dot{x}}_{pi}^q\hfill \\ m_p{\ddot{y}}_{pi}^q=\left(K_{spi}^q{\delta }_{spi}^q+C_{spi}^q{\dot{\delta }}_{spi}^q\right)sin{\alpha }_{pi}cos{\beta }_b-K_p\left(y_{pi}^q-r_c{v}_c\right)-C_p\left({\dot{y}}_{pi}^q-r_c{\dot{v}}_c\right)\hfill \\ m_p{\ddot{z}}_{pi}^q=-\left(K_{spi}^q{\delta }_{spi}^q+C_{spi}^q{\dot{\delta }}_{spi}^q\right)sin{\alpha }_{pi}cos{\beta }_b-K_{pa}{z}_{pi}^q-C_{pa}{\dot{z}}_{pi}^q-K_{pa}^{LR}\left(z_{pi}^q-z_{pi}^p\right)-C_{pa}^{LR}\left({\dot{z}}_{pi}^q-{\dot{z}}_{pi}^p\right)\hfill \\ I_{pi}{\ddot{\theta }}_{pi}^q=-\left(K_{spi}^q{\delta }_{spi}^q+C_{spi}^q{\dot{\delta }}_{spi}^q\right)r_{pi}cos{\beta }_b-K_{pb}^{LR}\left(v_{pi}^q-v_{pi}^p\right)-C_{pb}^{LR}\left({\dot{v}}_{pi}^q-{\dot{v}}_{pi}^p\right)+T_{out}\hfill \end{array},$$

(9)

where $r_c$ is the radius of the pitch circle of the carrier. $m_p$ is the mass of the planetary gear.

The dynamic equation of the extreme meshing can be rewritten as:

$${\mathit{M}}_{spi}^q{\ddot{X}}_{spi}^q+\left(C_{spi}^q+C_{sc}^q\right){\dot{X}}_{spi}^q+\left(K_{spi}^q+K_{sk}^q\right)X_{spi}^q=F_{sn}^q+F_{se}^q,$$

(10)

where $C_{sc}^q$ represents the support and torsional damping matrix, $K_{sc}^q$ denotes the support and torsional stiffness matrix. In addition, the mass matrix $M_{spi}^q$, meshing damping matrix $C_{spi}^q$, meshing stiffness matrix $K_{spi}^q$, vector of forces caused by the dynamic transmission errors $F_{se}^q$, and the vector of external forces $F_{sn}^q$ can be expressed as the following equations:

$$\begin{array}{l}{M}_{spi}^q=diag\left(m_s,m_s,m_s,I_s,m_{pi},m_{pi},m_{pi},I_{pi}\right)\\ C_{spi}^q=C_{spi}{\left(V_{spi}^q\right)}^T{V}_{spi}^q\\ K_{spi}^q=K_{spi}{\left(V_{spi}^q\right)}^T{V}_{spi}^q\\ F_{se}^q=C_{spi}{\left(V_{spi}^q\right)}^T{\dot{e}}_{spi}^q\left(t\right)+K_{spi}{\left(V_{spi}^q\right)}^T{e}_{spi}^q\left(t\right)\\ F_{sn}^q={\left[0,0,0,T_{out},0,0,0,T_p\right]}^T\end{array}.$$

(11)

Figure 3b shows the dynamic model of the ring and planetary gear. Similarly, their motion equations can be established based on Newton’s second law, and it can be written as [12]:

$$M_{rpi}^q{\ddot{X}}_{rpi}^q+\left(C_{rpi}^q+C_{rc}^q\right){\dot{X}}_{rpi}^q+\left(K_{rpi}^q+K_{rk}^q\right)X_{rpi}^q=F_{rn}^q+F_{re}^q.$$

(12)

Figure 3c shows the dynamic model of the carrier and planetary gear. In this section, L and R are not used because of the same parameters. The motion equations of the carrier and planetary gears can be deduced as follows:

$$\begin{array}{cc}\hfill m_c{\ddot{x}}_c& =-K_c\left(x_c-x_{pi}-r_c{\theta }_{c}sin{\alpha }_{pi}\right)-C_c\left({\dot{x}}_c-{\dot{x}}_{pi}-r_c{\dot{\theta }}_{c}sin{\alpha }_{pi}\right)\hfill \\ \hfill m_c{\ddot{y}}_c& =-K_c\left(y_c-y_{pi}+r_c{\theta }_{c}cos{\alpha }_{pi}\right)-C_c\left({\dot{y}}_c-{\dot{y}}_{pi}+r_c{\dot{\theta }}_{c}cos{\alpha }_{pi}\right)\hfill \\ \hfill m_c{\ddot{z}}_c& =-K_{ca}\left(z_c-z_{pi}\right)-C_{ca}\left({\dot{z}}_c-{\dot{z}}_{pi}\right)\hfill \\ \hfill I_c{\ddot{\theta }}_c& =r_{c}sin{\alpha }_{pi}\left[K_c\left(x_c-x_{pi}-r_c{\theta }_{c}sin{\alpha }_{pi}\right)+C_c\left({\dot{x}}_c-{\dot{x}}_{pi}-r_c{\dot{\theta }}_{c}sin{\alpha }_{pi}\right)\right]\hfill \\ & -r_{c}cos{\alpha }_{pi}\left[K_c\left(y_c-y_{pi}+r_c{\theta }_{c}cos{\alpha }_{pi}\right)+C_c\left({\dot{y}}_c-{\dot{y}}_{pi}+r_c{\dot{\theta }}_{c}cos{\alpha }_{pi}\right)\right]\hfill \end{array},$$

(13)

where $K_c$ denotes the radial support stiffness of the carrier, $C_c$ represents the radial support damping of the carrier, $K_{ca}$ denotes its axial support stiffness, and $C_{ca}$ represents its axial support stiffness. $m_c$ is its mass.

$$\begin{array}{l}{m}_p{\ddot{x}}_{pi}=K_c\left(x_c-x_{pi}-r_c{\theta }_c\sin {\alpha }_{pi}\right)+C_c\left({\dot{x}}_c-{\dot{x}}_{pi}-r_c{\dot{\theta }}_c\sin {\alpha }_{pi}\right)\\ m_p{\ddot{y}}_{pi}=K_c\left(y_c-y_{pi}+r_c{\theta }_c\cos {\alpha }_{pi}\right)+C_c\left({\dot{y}}_c-{\dot{y}}_{pi}+r_c{\dot{\theta }}_c\cos {\alpha }_{pi}\right)\\ m_p{\ddot{z}}_{pi}=K_{ca}\left(z_c-z_{pi}\right)+C_{ca}\left({\dot{z}}_c-{\dot{z}}_{pi}\right)\\ I_{pi}{\ddot{\theta }}_{pi}=0\end{array}.$$

(14)

According to Equations (13) and (14), the motion equations can be obtained as follows:

$$M_{cpi}^q{\ddot{X}}_{cpi}^q+C_{cc}^q{\dot{X}}_{rpi}^q+K_{ck}^q{X}_{cpi}^q=F_{cn}^q,$$

(15)

where $C_{cc}^q$ represents the support damping matrix, $K_{cc}^q$ denotes the support stiffness matrix. In addition, the mass matrix $M_{cpi}^q$ and vector of the external forces $F_{cn}^q$ can be provided in the following equations:

$$\begin{array}{l}{M}_{cpi}^q=diag\left(m_c,m_c,m_c,I_c,m_{pi},m_{pi},m_{pi},I_{pi}\right)\\ F_{cn}^q=\left[0,0,0,T_{in},0,0,0,0\right]\end{array}.$$

(16)

The motion equations of all the subsystems were coupled into an overall motion equation, and it could be written as:

$$M\ddot{X}+C\dot{X}+KX=F.$$

(17)

2.2.2. Meshing Stiffness

Based on the assumption of dynamic modelling that two helix gears are coupled into a herringbone gear, the meshing stiffness of the herringbone gear can be similarly calculated by coupling the meshing stiffness of the two helix gears, and its expression can be provided by the following equations:

$$K_m=2b{G}_T,$$

(18)

where $b$ is the unilateral helical tooth width, $G_T$ represents the meshing stiffness per helical tooth width, and its calculation expression is:

$$G_T=\frac{0.75{\epsilon }_{\alpha }+0.25}{h},$$

(19)

where ${\epsilon }_{\alpha }$ is the transverse contact ratio and $h$ stands for the flexibility of the unit tooth width, which is related to the tooth number and modification coefficient.

$$C=2\zeta \sqrt{K_m\frac{r_{b1}^2{r}_{b2}^2{m}_1{m}_2}{r_{b1}^2{m}_1+r_{b2}^2{m}_2}},$$

(20)

$\zeta$ is the meshing damping ratio, and its value is in the range of 0.03 to 0.17. $r_{b1}$ and $r_{b2}$ are the radii of the base circles of the driving and driven gears. $m_1$ and $m_2$ are their masses.

2.2.3. Time-Varying Meshing Transmission Errors

The internal and external meshing errors, respectively, are shown as follows:

$$\begin{array}{l}{e}_{spi}^q\left(t\right)={\displaystyle \sum _l^L{e}_{spi}}\cos \left(l{\omega }_{m}t+{\phi }_{spi}-l{Z}_s{\phi }_i+f\left(q\right)\right)\\ e_{rpi}^r\left(t\right)={\displaystyle \sum _l^L{e}_{rpi}}\cos \left(l{\omega }_{m}t+{\phi }_{rpi}+l{\vartheta }_{rs}+l{Z}_r{\phi }_i+f\left(q\right)\right)\end{array},$$

(21)

where $e_{sp}$ denotes the meshing transmission error between the sun gear and planetary gear, ${\phi }_{spi}$ represents their initial phase, $e_{rp}$ denotes the meshing transmission error between the ring gear and planetary gear, ${\phi }_{rpi}$ represents their initial phase, ${\omega }_m$ is their meshing frequency, $Z_s$ represents the tooth number of the sun gear, $Z_r$ represents the tooth number of the ring gear, and ${\vartheta }_{rs}$ is the meshing phase difference between the sun gear and ring gear. $f\left(q\right)$ is the piecewise function, which is used to distinguish the meshing transmission error of the left and right helical gear. Its expression is provided in the following equation:

$$f\left(q\right)=\left\{\begin{array}{cc}0& q=L\\ l{\gamma }_{stg}& q=R\end{array}\right..$$

(22)

2.3. Wind Velocity Model

2.3.1. Turbulent Effect

The three-dimensional turbulent wind was simulated based on the Kaimal spectral model [31]:

$$\frac{f\cdot S{\left(f\right)}_j}{u_f^2}=\frac{{\eta }_a{\kappa }^{{\eta }_{\gamma }}}{{\left({\eta }_c+{\eta }_b{\kappa }^{{\eta }_{\alpha }}\right)}^{{\eta }_{\beta }}}\left(j=u,v,w,uw\right)$$

(23)

where $S{\left(f\right)}_j$ denotes the power spectral ensity (PSD), $u_f$ represents the friction velocity, ${\eta }_a$, ${\eta }_b$, ${\eta }_c$, ${\eta }_{\alpha }$, ${\eta }_{\beta }$ and ${\eta }_{\gamma }$ are six parameters of the wind spectral model, and $\kappa$ stands for the reduced frequency. Its expression is as follows:

$$\kappa =\frac{f\cdot z}{\overline{u}},$$

(24)

where $\overline{u}$ denotes the mean wind velocity, $z$ represents the height in space, and $f$ is the frequency. In addition, the calculational expression of $u_f$ was simulated based on Weber’s research [32,33]:

$$u_f={\left({\overline{uw}}^2+{\overline{uv}}^2\right)}^{0.25}$$

(25)

According to Tieleman’s theory, the so-called blunt model and so-called pointed model are defined by four parameters, which are ${\eta }_c$, ${\eta }_{\alpha }$, ${\eta }_{\beta }$; ${\eta }_{\gamma }$. ${\eta }_c=1$, ${\eta }_{\alpha }=1$; ${\eta }_{\beta }=5/3$; and ${\eta }_{\gamma }=1$ for the blunt model, which can be used to describe topographic factor interference, whereas ${\eta }_{\alpha }=5/3$, ${\eta }_{\beta }=1$ is the pointed model, which is employed for describing flat topography. The Kaimal model is a mixed model because the horizontal wind is established using the blunt model, and the vertical wind is described using the pointed model. Thus, the calculation expression of the PSD from each component can be written as follows:

$$\begin{array}{c}\frac{f\cdot S{\left(f\right)}_u}{u_f^2}=\frac{105\kappa }{{(1+33\kappa )}^{5/3}}\\ \begin{array}{l}\frac{f\cdot S{\left(f\right)}_v}{u_f^2}=\frac{17\kappa }{{(1+9.5\kappa )}^{5/3}}\\ \frac{f\cdot S{\left(f\right)}_w}{u_f^2}=\frac{2.1\kappa }{(1+5.3{\kappa }^{5/3})},\end{array}\\ \frac{f\mathrm{Re}\left(S{\left(f\right)}_{uw}\right)}{u_f^2}=-\frac{14\kappa }{{(1+9.6\kappa )}^{7/3}}\end{array}.$$

(26)

Except for the one-node wind simulation, a co-coherence function was used to describe the correlation of the turbulence from two nodes, for which the distance is $d_j$. Its expression is defined as follows [34]:

$$\gamma \left(f,d_j\right)=\exp \left(-C_{ij}\frac{f{d}_j}{U(z)}\right);\hspace{1em}i=u,v,w,\hspace{1em}j=y,z,$$

(27)

where $\gamma \left(f,d_j\right)$ denotes the co-coherence and $C_{ij}$ stands for the co-coherence exponential decay coefficient. Their values were set to:

$$C_{uy}=7,C_{uz}=10;C_{vy}=7,C_{vz}=10,C_{wy}=6.5;C_{wz}=3.$$

(28)

Thus, the co-coherence power spectrum density matrix could be calculated via the following equation:

$$S^0(\omega )=\left[\begin{array}{cccc}{S}_{11}(\omega )& S_{12}(\omega )& \cdots & S_{1n}(\omega )\\ S_{21}(\omega )& S_{22}(\omega )& \cdots & S_{2n}(\omega )\\ \cdots & \cdots & \cdots & \cdots \\ S_{n1}(\omega )& S_{n2}(\omega )& \cdots & S_{nn}(\omega )\end{array}\right].$$

(29)

Then, the matrix was solved using the Cholesky decomposition method. Finally, the turbulent wind, u, v, w, could be obtained based on the Shinozika theory [35].

2.3.2. Wind Shear Effect

The wind shear effect is defined as the variation in the wind velocity and direction in a position in space. In the study of wind velocity, generally, the change in the wind velocity in the vertical direction is only considered. Based on the IEC standard [36], an exponential model was used to describe the wind shear effect:

$$V_{shear}(h)=V_{\mathrm{H}}{\left(\frac{h}{H}\right)}^{{\alpha }_w}.$$

(30)

where $H$ denotes the hub height, $V_H$ represents the wind velocity at the hub’s position, $V_{shear}(h)$ stands for the wind velocity of the height $h$, and ${\alpha }_w$ is the wind shear coefficient, normally.

Equation (30) can be also rewritten in the polar coordinate form:

$$V_{shear}(r,\theta )=V_{\mathrm{H}}{\left(\frac{H+r\cos \theta }{H}\right)}^{{\alpha }_w}=V_{\mathrm{H}}\left[1+w_{\mathrm{s}}(r,\theta )\right],$$

(31)

where $\theta$ is the position angle of a node in the rotational plane of the impeller, $r$ is the distance between the centre of the impeller hub and this node, and $w_{\mathrm{s}}$ is the wind shear coefficient. Its expression can be expanded using the third-order Taylor series [37]:

$$w_{\mathrm{s}}(r,\theta )\approx {\alpha }_w\frac{r}{H}\cos \theta +\frac{{\alpha }_w({\alpha }_w-1)}{2}{\left(\frac{r}{H}\right)}^2{\cos }^2\theta +\frac{{\alpha }_w({\alpha }_w-1)({\alpha }_w-2)}{6}{\left(\frac{r}{H}\right)}^3{\cos }^3\theta$$

(32)

Clearly, the effect of the height on locational wind velocity is reflected in the size of $w_{\mathrm{s}}(r,\theta )$, as shown in Figure 4.

The mean wind velocity in the rotational plane of the impeller can be calculated using the following equation:

$${\overline{V}}_{shaer}=\frac{1}{\pi R^2}{\int }_0^{2\pi }{\int }_0^R{V}_{shaer}(r,\theta )·r\mathrm{d}r\mathrm{d}\theta .$$

(33)

2.3.3. Tower Shadow Effect

Due to the hidden effect of the tower on the flow, the inflow velocity in the region of the tower’s upstream and downstream reduced. The phenomenon is defined as the tower shadow effect. In this wind turbine study, the tower shadow effect only occurred in the lower haft of the disk surface of the impeller. Thus, the tower shadow effect was expressed as the following equation:

$$V_{ts}(r,\theta )=\left\{\begin{array}{cc}0& 0<\theta <\pi /2\\ {\overline{V}}_{shear}\left[1+{\eta }_{ts}(r,\theta )\right]& \pi /2\le \theta <1.5\pi \\ 0& 1.5\pi \le \theta <2\pi \end{array}\right.,$$

(34)

where ${\eta }_{ts}(r,\theta )$ denotes the tower shadow coefficient and its expression can be provided using [38,39]:

$${\eta }_{ts}(r,\theta )=\left\{\begin{array}{cc}0& 0<\theta <\pi /2\\ r_t{\left(h\right)}^2\frac{r^2{\sin }^2\theta -l^2}{{\left(r^2{\sin }^2\theta +l^2\right)}^2}& \pi /2\le \theta <1.5\pi \\ 0& 1.5\pi \le \theta <2\pi \end{array}\right.,$$

(35)

where $r_t$ donates the tower’s radius, $l$ represents the shaft length between the impeller and hub. In general, the tower of large megawatt-scale wind turbines is designed to be wide at the bottom and narrow at the top. Thus, the tower’s radius changes with increasing height, and its expression is provided in the following equation:

$$r_t(h)=r_{\mathrm{tb}}+\frac{r_{\mathrm{tt}}-r_{\mathrm{tb}}}{H}h,$$

(36)

where $r_t(h)$ denotes the tower’s radius at $h$, $r_{\mathrm{tb}}$ represents the radius of the tower bottom, and $r_{\mathrm{tt}}$ is the radius of the tower top. The tower shadow coefficient curves are shown in Figure 5.

2.3.4. Wind Velocity

Based on the above theory and the sizes of the tower and blades, the sweep plane of the impeller was divided into 32 × 32 (1024) nodes to describe the turbulent wind velocity, as shown in Figure 6a,b. The horizontal coordinate represents the width of space and the vertical coordinate represents the height of space. The size parameters of the wind turbine are shown in

Table 1. Size parameters of the 1.5 MW wind turbine.

Name	Symbol	Value
Hub height	H	83 m
Impeller radius	R	38.75 m
Hub radius	-	1.5 m
Overhang	l	3.3 m
Rotational velocity of the impeller	ωp	11 rmp
Airfoil series	NACA63-4	-

. At 3 s, 6 s, and 9 s, the turbulent wind conditions are shown in Figure 7.

2.4. Aerodynamic Load Calculation

The power of a wind turbine originates from the wind acting on the impeller. Meanwhile, it is also the main load on each component of the wind turbine. Thus, the accurate acquisition of aerodynamic forces on the blades is key to analyzing the loads. Currently, the methods for calculating the aerodynamic loads on the wind turbines are the traditional BEM, the Schmitz theory, and the generalized dynamic wake theory. Since the physical meaning of the BEM is clear and easy to understand, its modified model can be widely used in the wind turbine field. Therefore, the calculation of the aerodynamic load is solved using the BEM.

Before starting with the calculation of the aerodynamic load, the four coordinate systems were defined to display the blade element’s position as shown in Figure 8, and the construction of the coordinate systems could begin with the observation point P of the rotational wind blade 1. The coordinate system $T_1{w}_1{a}_1{h}_1$ is connected to the ground and is a fixed system. The original point is set at the base of the tower. $T_2{w}_2{a}_2{h}_2$, $T_3{w}_3{a}_3{h}_3$, and $P{w}_4{a}_4{h}_4$ represent the tower top position, hub position, and blade element position, respectively. Two inclination angles exist in the $h-a$ plane (${\theta }_{ti}$ and ${\theta }_{co}$), but they are all set to zero for this calculated load task. The position angle, ${\beta }_p$, between the three blades is 120 degrees. ${\omega }_b$ stands for the angular velocity of the impeller.

According to the BEM, the lift and drag forces on each blade element are provided in:

$$\begin{array}{c}\mathrm{d}L=\frac{1}{2}\rho W{\left(r,\theta \right)}^{2}c(r)C_l\mathrm{d}r\\ \mathrm{d}D=\frac{1}{2}\rho W{\left(r,\theta \right)}^{2}c(r)C_d\mathrm{d}r\end{array}$$

(37)

where $\mathrm{d}L$ denotes the lift force unit, $\mathrm{d}D$ represents the drag force unit, $\rho$ is the airflow density, $c(r)$ stands for the chord length on the interface where the blade element is located, and $W\left(r,\theta \right)$ is the actual inflow wind velocity of the blade element; it can be obtained via the coordinate transformation of velocity under the turbulent wind condition. $C_l$ is the lift coefficient and $C_l$ is the drag coefficient.

The torque of the blade element can be written as:

$$F_{torque}=dL\cos \phi -dD\sin \phi .$$

(38)

To calculate the torques acting on the entire blade, it is only necessary to integrate the torque elements spanwise across the blade, as shown in Figure 9, Figure 10 and Figure 11. The inflow velocity of the hub is steady in the results from Figure 9 and Figure 10. The torque of the different blades can be seen in Figure 9a, and the torque variation is mainly affected by the wind shear effect. The torque presents periodic fluctuations. Combined with Figure 5, the value of $w\left(r,\theta \right)$ is at its minimum when the blades rotate to 180 degrees. The blade alternating fluctuation is due to the position angle difference between the blades of 120 degrees.

In addition, another key factor is the tower effect. After considering the combined wind shear and tower effect, the torque of the blades and hub is illustrated in Figure 10. The torque is weakened by this combined effect when the blades rotate to the centre area of the tower shadow, as can be seen in the marked position in Figure 10. This is because the actual wind velocity plunges under the combined effect when the blades’ position is 180 degrees, as shown in Figure 8. The torque drop not only has a huge impact on the transmission system but also has an obvious impact on the power grid. Through carefully comparing Figure 9a and Figure 10a, it can be found that the torque under the sole wind shear effect is smaller than the torque under the combined effect on both sides of the tower shadow centre area. This can make up for the power loss in the centre area of the tower shadow to a certain extent.

When the turbulence, wind shear, and tower effect are all considered, the torque presents a nonlinear fluctuation, as shown in Figure 11. The torque in the centre area of the tower shadow may be enhanced or not diminished because of the turbulent wind effect. Combined with the axial flow wind velocity map in Figure 7, the wind shear is still the main factor among the three. This is because the influence of the wind shear effect lasts, whereas the tower shadow effect works in a small area at its centre.

3. Results and Analysis

The parameters of the herringbone planetary gear of the wind turbine are shown in

Table 2. Wind turbine gearbox parameters [13].

Parameters	Unit	Sun Gear	Planetary Gear	Ring Gear	Carrier
Tooth		22	41	104
Modulus	mm	16	16	16
Mass	kg	149.6	517.8	671.3	275.1
Helix angle	deg	15	15	15
Meshing stiffness	N/m	2.53 × 1010		3.6 × 1010
Moment of inertial	kg·m2	3.7	43.1	527.6
Couling stiffness	N/m	108	108	109
Input power	kW	1500

. Based on the modelling and theories, the hub loads could be precisely calculated under different wind conditions, as shown in Figure 9, Figure 10 and Figure 11. They can be seen as the external excitation acting on the carrier of this herringbone planetary gear system. Equation (17) can be solved using the Runge-Kutta method to obtain the displacements and velocities from each component in different directions. After that, the dynamics and load-sharing characteristics of this system could then be fully and clearly understood.

3.1. Load-Sharing Characteristics

In order to quantify the load-balancing performance of the planetary transmission system, the load-balancing coefficient is defined here as [13]:

$$\begin{array}{l}s{l}_{spi}^q=\frac{3\times {\left(F_{spi}^q\right)}_{\max }}{{\displaystyle \sum _i^3{F}_{spi}^q}}\left(q=L,R\right)\\ s{l}_{rpi}^q=\frac{3\times {\left(F_{rpi}^q\right)}_{\max }}{{\displaystyle \sum _i^3{F}_{rpi}^q}}\left(q=L,R\right)\end{array},$$

(39)

where $s{l}_{spi}^q$ and $s{l}_{rpi}^q$ represent the time-varying load-sharing coefficient. $F_{spi}^q$ stands for the meshing forces between the sun gear and planetary gear and $F_{rpi}^q$ stands for the meshing forces between the ring gear and planetary gear.

In order to evaluate the effect of the support stiffness on the load-sharing coefficient, the global load-sharing coefficient, SL, is present. Its expression can be written using [13]:

$$\begin{array}{l}S{L}_{sp}^q=\max \left(s{l}_{spi}^q\right)\\ S{L}_{rp}^q=\max \left(s{l}_{rpi}^q\right)\end{array}.$$

(40)

3.1.1. Effect of the Support Stiffness

In the section, the wind shear effect and the tower shadow effect are considered in the calculation of the hub load. Because the turbulent wind velocity has significant randomness, it was not considered. In addition, using periodic external excitation could more clearly describe the influence of the support stiffness on the load-sharing characteristics. The support stiffnesses of the main rotational components was set as the controlled parameters, and the damping was 0.2% of the stiffness based on the study of bearings in a wind turbine [40]. The load-sharing coefficient can be calculated using Equation (40). Value orthogonal tests of the sun gear support stiffness and planetary gear support stiffness are shown in Figure 12 and Figure 13. The value orthogonal tests of the ring gear support stiffness and planetary gear support stiffness are shown in Figure 14 and Figure 15.

In summary, the value of $S{L}_{rp}^L$ can be seen in Figure 12a. It can be found that the coefficient decreases as the value of $K_s$ increases and increases as the value of $K_p$ increases. The minimum value of $S{L}_{rp}^L$ can be observed at the position corresponding to the coordinate $\left({10}^6,{10}^8\right)$. Moreover, as the value of the support stiffness increases, the variation degree of the load-sharing coefficient is different. When the $K_s$ is set to 10⁶, the increase in $S{L}_{rp}^L$ is obvious with an increase in the value of $K_p$, whereas this increasing trend is slight at $K_s={10}^8$. Similarly, when the $K_p$ is set to the minimum, the fluctuation decrease in the load-sharing coefficient is obvious with an increase in $K_s$, whereas it exhibits a decline in the form of a parabolic function at $K_p={10}^8$. The variation trend of $S{L}_{rp}^R$ is shown in Figure 12b. Combined with Figure 12a, it can be seen that in the right-side model, the variation trend of the load-sharing coefficient is the same, but the magnitude of the number and the magnitude of the change are different. $S{L}_{rp}^R$ is smaller than $S{L}_{rp}^L$. This is because the torsional coupling stiffness and damping of the left and right ends play a role.

The values of $S{L}_{sp}^L$ and $S{L}_{sp}^R$ can be seen in Figure 13a,b. It can be found that an obvious difference appears, as can be seen in Figure 13, with the variation in $K_s$, whereas with an increase in $K_p$, the decline in $SL$ is slight.

The values of $S{L}_{rp}^L$ and $S{L}_{rp}^R$ with the change in $K_r$ and $K_p$ can be seen in Figure 14a,b. The global load-sharing coefficient $S{L}_{rp}^L$ and $S{L}_{rp}^R$ will decrease if $K_r$ increases. They will increase if $K_p$ increases. On the other hand, the effect of $K_r$ and $K_r$ on the global load-sharing coefficient is present. The values of $S{L}_{sp}^L$ and $S{L}_{sp}^R$ with the change of $K_r$ and $K_p$ can be seen in Figure 15a,b. It can be seen that opposite trends to Figure 14 are shown as the values of $K_r$ and $K_p$ change.

According to the above description and analysis, combined with Figure 12, Figure 13, Figure 14 and Figure 15, the support stiffness of the sun gear, $K_s$, has a great influence on both the load-sharing coefficient between the ring and planetary gear and the load-sharing coefficient between the sun and planetary gear. The support stiffness of the ring gear, $K_r$, only has an obvious effect on the load-sharing coefficient between the ring and planetary gear. Similarly, the support stiffness of the planetary gear, $K_p$, is the same. Therefore, improving the support stiffness of the ring gear is the best scheme for the parameter adjustment. This scheme can not only reduce the load-sharing coefficient between the ring gear and planetary gear but also ensure that it does not have a significant impact on the load-sharing coefficient between the sun gear and planetary gear. Meanwhile, the support stiffness of the planetary gear should be designed to be a low value while ensuring that the design is reasonable. This can make the global load-sharing coefficient closer to 1 than before.

3.1.2. Contribution Rate Analysis

In order to analyze the effect of the different factors on the time-varying load-sharing characteristics, such as the turbulent, tower shadow effect, and wind shear effect, the relative contribution ratio $r_i$ and the real-time rate of change $c_t$ were proposed, and their expressions can be written as [41]:

$$\begin{array}{l}{r}_i(t)=\frac{\left[s_i\left(t\right)-s_t\left(t\right)\right]}{{\displaystyle \sum \left[s_i\left(t\right)-s_t\left(t\right)\right]}}\times 100\%\left(i=Turbulent,tower shadow,wind shear\right)\\ c_t\left(t\right)=\frac{\left[s_t\left(t\right)-s_t\left(V_H\right)\right]}{s_t\left(V_H\right)}\times 100\%\end{array}.$$

(41)

where $s_i\left(t\right)$ stands for the time-varying load-sharing coefficient without considering the i influence factor, $s_t\left(t\right)$ denotes the time-varying load-sharing coefficient with all influence factors, and $s_t\left(V_H\right)$ represents the time-varying load-sharing coefficient under the steady inflow condition ($V_H$ = 11 m/s).

The relative contribution rate of each influence factor is shown in Figure 16 and Figure 17. The relative contribution rate of the turbulent was maintained at the highest level. Strikingly, the mean values of $r_{turbulent}$ in Figure 16a,b and Figure 17a,b are 80.368%, 80.372%, 80.371%, and 80.370%, respectively. Regarding the tower shadow effect, they are 10.314%, 10.312%, 10.312%, and 10.313%. The contribution rate of the wind shear effect were calculated to be 9.318%, 9.316%, 9.317%, and 9.317%. These can describe the whole variation of the relative contribution rate. However, the time-varying properties cannot be reflected.

In Figure 16a, the relative contribution rate ranges of the wind shear, tower shadow, and turbulent are [0.91%, 23.61%], [1.31%, 28.44%], and [47.94%, 93.88%]. It is not difficult to find the law of the contribution rates of the wind shear and tower shadow. The excitation load enters the shadow effect center at t = 11.32 s and t = 13.06 s. When the load excitation enters the region of the tower shadow, their values increase gradually, whereas the load excitation exits the region of the tower shadow, and their values decrease. In particular, the contribution rate of the tower shadow quickly climbed to 28.44%. This is because the tower shadow causes the load excitation to quickly fall. When the blades are not in the region of the tower shadow, the contribution rate of the wind shear and tower is low (mean ≤ 7%). This means that the effect of the wind shear and tower shadow can be ignored at this stage.

In addition, combined with Figure 16 and Figure 17, the real-time rate of change in $s{l}_{sp1}^q$ is smaller than $s{l}_{rp1}^q$. It shows that the load-sharing characteristic between the ring gear and planetary gear is more sensitive to turbulent wind conditions than between the sun gear and planetary gear. This also reminds designers to focus on the design of the gear rings. The cyclic aerodynamic load caused by the wind shear and tower shadow effect always exists, and when the excitation enters the tower shadow region, the wind shear and tower shadow effect will affect the load-sharing characteristics of the system. The vibration and fatigue problems of the wind turbine gearboxes become more complicated. At the same time, due to the periodic change in the load, the output power of the grid-connected wind turbines also has certain fluctuations, and the output power quality will be affected to a certain extent.

3.2. Dynamic Characteristics

To analyze the dynamic characteristics in the herringbone planetary gear system, the different wind conditions, including the steady inflow (V_H = 11 m/s) and turbulent wind, were considered in the motion equations. After solving, the meshing forces, displacements, and centroid trajectories could be obtained. Their description and analysis are shown as follows. In addition, the turbulent wind in this section stands for a wind condition including the wind shear, tower shadow, and turbulent effect.

3.2.1. Meshing Force

In Figure 18 and Figure 19, the internal and external meshing forces are displayed, respectively. It can be seen that the meshing forces acting on the internal meshing line are slightly affected by the turbulent wind condition, and the right side is larger than the left side. This is because the space translational displacement of the right side from the ring gear is larger than the left side, as can be seen in Figure 20 and Figure 21. Only after the blade enters the centre of the tower shadow effect, the meshing force curve will have an obvious impact. In addition, the meshing forces curves under turbulent wind conditions are lower than the meshing forces under steady inflow. This is mainly due to the decline in the wind velocity caused by the wind shear.

In Figure 19, except for the above analysis, the meshing forces acting on the external meshing line are larger than those acting on the internal meshing line. They are strongly affected by the turbulent wind conditions. This is because the torsional displacement fluctuation is large from the sun gear. Notably, the meshing frequency of the sun gear is higher than the other. If the force on the tooth surface is also high, tooth surface damage is unavoidable. Therefore, the gear pair between the sun gear and planetary gear should be handled specifically in manufacturing and design to decrease the meshing stiffness. In the design of the herringbone gearbox, the external meshing gear pair should be specially treated to deal with the problem of meshing impact enhancement under dynamic loads.

3.2.2. Displacement

As can be seen in Figure 20 and Figure 21, the torsional displacements of the ring gear are more obviously influenced by the turbulent wind condition than the axial displacements. This is because the axial excitation is ignored after considering the herringbone gear special structure and the bearing support design of this wind turbine. Therefore, the fluctuations in the axial displacements are generated by the axial–torsion coupling term in the motion equations. This greatly weakens the influence of the turbulent wind condition on the axial vibration in the transmission system. Meanwhile, this means that the axial displacements of the ring gear are slightly affected by the torsional excitation under turbulent wind conditions. In the design of the ring gear, the fluctuation in the torsional excitation can be ignored.

The axial and torsional displacements of the sun gear are shown in Figure 22 and Figure 23. Similarly, the torsional displacement variation trends of the sun gear are the same as the ring, but these are larger than the ring gear. This is because the rotational velocity is larger, the torque is smaller, and the vibration is larger.

In summary, the vibration amplitude of the sun gear in the torsional direction should be concerned. However, both the ring gear and sun gear are sensitive to turbulent wind conditions but the vibration amplitudes of the ring gear are small, whereas, from the sun gear, it is too large. In this way, the sun gear is subjected to both high-frequency loads and low-frequency loads.

4. Discussion

In Section 3, the contribution rate from the different influence factors on the load-sharing coefficient is discussed using statistical methods. It can be found that the most sensitive factor that causes changes in the load-sharing characteristics is the turbulent effect through the time-varying relative contribution rate and rate of change. The contribution rates of the tower shadow and wind shear effect also surge when the blades enter the centre region of the tower shadow. This causes serious fluctuations in the time-varying load-sharing coefficient. This is because the superposition of the two effects makes the external excitation load fluctuate greatly.

The wind shear effect is mainly affected by the ground roughness. Based on the measured data, the mean wind shear coefficient of the sea, plain, hill, mountainous area, and seaside are approximately 0.1, 0.14, 0.17, 0.22, and 0.32, respectively. But these are just the estimates related to ground roughness. For example, the terrain is mostly hilly and plain in the central and eastern parts of China, but the measured wind shear coefficient is above 0.3. Based on the above different terrain, the influence of the wind shear effect on the global load-sharing coefficient is shown in

Table 3. Effect of the wind shear coefficient.

Wind Shear Coefficient	0	0.1	0.14	0.17	0.22	0.32
$S{L}_{rp}^L$	1.2127	1.2133	1.2116	1.2138	1.2141	1.2143
$S{L}_{sp}^L$	1.0360	1.0361	1.0362	1.0363	1.0364	1.0365
Amplitude $S{L}_{rp}^L$/$S{L}_{sp}^L$	0.212/0.036	0.212/0.036	0.213/0.0361	0.214/0.0362	0.214/0.0363	0.2142/0.0364

through the method of this paper. The wind shear effect is only considered.

Flat terrain can improve the load-sharing characteristics, but the improvement is slight. However, the turbulence intensity is different in different terrains, and turbulence is the main cause of uneven load based on above the study.

In addition, the tower shadow coefficient is mainly affected by the tower size. The load-sharing characteristics of the transmission system can be improved by reducing the impact caused by the tower shadow effect. Meanwhile, the theoretical numerical results are provided in

Table 4. Effect of the tower shadow.

$\mathbf{Tower} \mathbf{Size} {\mathit{r}}_{\mathit{t}\mathit{t}}/2$	0.6	0.7	0.8	1
$S{L}_{rp}^L$	1.2181	1..2186	1.2234	1.2350
$S{L}_{sp}^L$	1.0378	1.0378	1.0392	1.0402
Amplitude $S{L}_{rp}^L$/$S{L}_{sp}^L$	0.2176/0.0378	0.2185/0.0378	0.2233/0.0391	0.2338/0.0401

.

Through the table, it can be seen that the tapered tower design can effectively improve the load-sharing characteristics of the transmission system. Three influence factors were considered in this paper, but realistic wind conditions are complex and there are many more factors than were considered in this paper, such as the wake effect, stall effect, voltage drop, extreme wind regime, and so on. If their contribution rates are analyzed, the research mode in this paper can provide a reference. The precise modelling of environmental loads and transmission systems is still the direction of our future research.

5. Conclusions

In this study, the load-sharing and dynamic characteristics of a herringbone planetary gear transmission system in a wind turbine were researched using the simulation method. The main conclusions are as follows:

(1)

The load-sharing coefficients of the left and right sides are similarly sensitive to the support stiffness, and the difference between the two is small. The support stiffness of the planetary gear can alter the effect of the support stiffness of the ring gear and sun gear on the load-sharing coefficient and vice versa. The interaction between the three stiffnesses should be considered when adjusting the support stiffness to control the load-sharing characteristics;

(2)

The turbulent effect is the most critical cause of the time-varying load-sharing characteristic variation, followed by the wind shear effect and tower shadow effect. Although when the blades enter the centre region of the tower shadow, the influence of the wind shear effect and tower shadow effect increased; they still do not exceed the turbulent effect;

(3)

The meshing forces acting on the external meshing line are obviously greater than those acting on the internal meshing line under turbulent wind conditions. In the design of the herringbone gearbox, the external meshing gear pair should be specially treated to deal with the problem of meshing impact enhancement under dynamic loads.

As a whole, the present work investigated the dynamic characteristics of a herringbone planetary gear system applied to a wind turbine under turbulent wind conditions. The study demonstrates that accurately calculating the torque load of the wind turbine gearbox is necessary because the dynamic behavior of the inner and outer meshing pairs is different after the action of a time-varying load. Additionally, the left and right sides of the herringbone gear are not the same. This provides an effective theoretical reference for the design of the herringbone planetary gear of wind turbines. Through the analysis and discussion, the correct geographical position of the wind turbine and the reasonable design of the tower are both good schemes that can effectively ensure the load-sharing performance of the transmission system and reduce the vibration. Further studies on the optimal design of the herringbone planetary gearbox of a wind turbine will be conducted to achieve the purpose of reducing vibration and extending service life.