Nonlinear characteristics of a multi-degree-of-freedom wind turbine’s gear transmission system involving friction

In this study, a 42-degree-of-freedom (42-DOF) translation–torsion coupling dynamic model of the wind turbine’s compound gear transmission system considering time-varying meshing friction, time-varying meshing stiffness, meshing damping, meshing error and backlash is proposed. Considering the different meshing between internal and external teeth of planetary gear, the time-varying meshing stiffness is calculated by using the cantilever beam theory. An improved meshing friction model takes into account the mixed elastohydrodynamic lubrication to calculate the time-varying meshing friction. The bifurcation diagram of the excitation frequency is utilized to analyze the bifurcation and chaos characteristics of the system. Meanwhile, the nonlinear dynamic characteristics of the gear system are identified by using the time domain diagrams, phase diagrams, Poincare maps and amplitude–frequency spectrums of the gear system. The results show that the system has complex bifurcation and chaotic behaviors including periodic, quasi-periodic, chaotic motion. The bifurcation characteristics of the system become complicated and the chaotic region increases considering the effects of friction in the high-frequency region.


Introduction
Compound gear transmission system is one of the most important transmission devices in wind turbine. Due to harsh operation environment, multi-frequencies excitations and complicated structure, wind turbine's compound gear transmission system is considered vulnerable in most cases. Thus, it is necessary to study the nonlinear dynamic characteristics of wind turbine's compound gear transmission system and evolution mechanism by establishing reliable dynamic model of compound gear transmission system to ensure the safe and stable operation for wind turbine.
A number of studies have been conducted on the dynamic modeling of planetary gear systems [1]. Lin and Parker [2] and Ambarisha and Parker [3] adopted analytical model and finite element model to study the nonlinear dynamic characteristics of planetary gear, and proposed a lumped parameter analytical model. By comparing with the finite element results in the complex nonlinear dynamic range in the system, the validity of the model was verified. Zhou et al. [4] established multi-degree-of-freedom dynamic model considering the coupling translation-torsion vibration. The coupled multi-body dynamics of the spur gear rotor bearing system was studied containing backlash, transmission error, eccentricity, gravity. The dynamic behaviors were investigated in different system parameters. Wang et al. [5] established the planetary gear transmission system considering the coupling effects of multiple nonlinear factors. The Runge-Kutta numerical method was used to solve the system vibration equation. Zhu et al. [6] established the lumped parameter dynamic model of compound planetary gear sets. The nonlinear dynamic characteristics of the gear sets are researched by changing the value of the backlash, meshing stiffness and amplitude of error excitations. Wang et al. [7] established 8-DOF gear system and applied improved meshing stiffness model to calculate the time-varying meshing stiffness. The dynamic chaos and bifurcation behavior of translation-torsion coupling system are analyzed. Zhu et al. [8] studied the influence of backlash on the encased differential planetary gear system. The results show that the backlash will cause the bilateral impact of the meshing tooth, resulting in chaotic motion of the system. The change of backlash at any transmission stage has a significant effect on gear meshing force and dynamic characteristics of the system. It can be seen from the above studies that the research on the nonlinear dynamic characteristics of planetary gear transmission system mainly focuses on the dynamic response analysis of singlestage planetary gear transmission system. Few studies have been done on two-stage planetary gear system in published studies, especially the two-stage planetary gear coupling translation-torsion system of wind turbine's compound gear transmission system.
During operation of compound planetary gear transmission system, friction is a strongly nonlinear internal excitation, which is positively correlated with the dynamic characteristics of compound gear transmission system. He et al. [9,10] proposed a new method considering the sliding friction and time-varying meshing stiffness, and studied the influence of its dynamic characteristics. Li and Kahraman [11,12], Xu et al. [13], Cioc et al. [14] proposed a tribo-dynamic model for spur gear pairs. The iterative calculation scheme of lubrication and dynamic model was proposed to make the tribological and dynamic behavior of spur gear pair fully coupled. Liu et al. [15] and Li et al. [16] proposed a single-stage spur gear model considering starved lubrication. By studying the dynamic coupling relationship between friction and dynamic characteristics of gear system, the dynamic meshing force of the system was predicted. A calculation method of meshing stiffness under starved lubrication was provided. The research of gear lubrication dynamics mainly focuses on single pair of meshing gears. Hou et al. [17] proposed a method to solve the different friction tooth surfaces of singlestage gear by subsection, and then, the friction force and friction torque are analyzed by summation. Han et al. [18] investigated the dynamic response of a pair of helical gears considering the interaction between friction and time-varying meshing stiffness, and analyzed the influence of friction on the dynamic characteristics of the system. Zhao et al. [19] established a fractal contact model for gear pair contact considering the influence of roughness on the normal contact stiffness of gear, and analyzed the influence of roughness on the dynamic characteristics of the system. Lin and Parker [20] established a nonlinear dynamic model of gear considering time-varying meshing stiffness, backlash, sliding friction and random external load. The results show that under the condition of invariable load, friction makes the duration of the system in the stable state longer and the periodic stability of the system worse. Li et al. [21] obtained time-varying meshing friction coefficient of tooth surface by fitting the experimental results of double disks. Chen and Lin [22] studied the calculation method of friction coefficient and meshing efficiency of plastic gear pair under dry friction condition, the dynamic characteristics of the system under dry friction condition were obtained. Bae et al. [23] adopted finite element simulation method to calculate the contact stress between coated polymer gears considering different friction coefficients. Keller et al. [24] established a thermal elastohydrodynamic model for the single-tooth gear. Mixed friction is considered based on measured rough surfaces via indirect coupling. Park [25] investigated friction and transmission error of spur gear due to sliding friction under quasi-static condition. Ma et al. [26] proposed an improved analytical model for time-varying meshing stiffness calculation considering friction between teeth surfaces. Ghosh and Chakraborty [27] studied a 6-DOF translation-rotation model considering friction. The influence of different friction parameters on the system was analyzed. It can be seen from the above studies that the research model for multi-degree-of-freedom planetary gear is mainly the purely rotational model of gear system. The study of gear friction is mostly based on single pair gear. There exists compact coupling relationship between the translational and torsional vibrations of gear system. Friction has a great influence on the dynamic characteristics of wind turbine system. The motivation of this paper is to established multi-degreeof-freedom translation-torsion coupling dynamic planetary gear model to accurately describe the dynamic behavior of wind turbine's compound gear transmission system. The consideration of TVMS is implemented by the cantilever beam theory. Another main research content is the calculation of TVMF. In this paper, a 42-DOF translation-torsion coupling dynamic model of the wind turbine's gear transmission system considering TVMF, TVMS, meshing damping, meshing error and backlash is established. The dimensionless system equation is solved by Runge-Kutta numerical method. The dynamic characteristics of the gear system are identified from the bifurcation diagram, time domain diagrams, phase diagrams, Poincare maps and amplitude-frequency spectrums of the gear system. The research can provide some useful information for the operation of wind turbine's gear transmission system.

Dynamic modeling of gear system
Taking the time-varying meshing friction and timevarying meshing stiffness into account, a 42-DOF translation-torsion coupling lumped parameter dynamic model is established, as shown in Fig. 1. Fig-ure 1a presents wind turbine gear's transmission system consisting of two planetary gear transmission operating at the low-speed stage and intermediate stage and one fixed-axis gear transmissions operating at the highspeed stage, respectively, which have large transmission ratio in relatively compact structures. The parameters of the planetary gear transmission system (I and II) are shown in Fig. 1b. In this model, the involute spur gear pairs are simplified as spring-damping system. The planetary gear transmission system consists of components such as carrier (c), ring (r ), sun (s) and planet ( p n (n = 1, 2, ..., N )). The stiffness of springs k st , k rt , k ct can be defined according to the power distribution of the system. Each spur gears i(i = s 1 , r 1 , c 1 , p 1,n (n = 1, 2, ..., N ), g 1 , s 2 , r 2 , c 2 , p 2,n (n = 1, 2, ..., N )) is simplified as a rigid disk of base radius r i , mass moment of inertia I i and angular displacement θ i . The torsional displacements are represented by u c , u r , u s , u p n , respectively. Figure 1c shows the translation-torsion coupling dynamic model of the fixed-axis gear system including the gear(g 1 ) and pinion(g 2 ). The carrier (c 1 ) and the pinion(g 2 ) are identified as input(T c 1 ) and output(T g 2 ) components of the system. Two coordinates (O XY, ox y) are created to establish the dynamic equations of the system easily. The center of the fixed coordinate system (O XY ) is located in the center of the sun gear. The moving coordinate system (ox y) is fixed with the planetary carrier and rotates with the theoretical angular velocity (ω c ). The relationship between the two coordinates system is shown in Fig. 1b. In terms of the analysis above, the relationship could be obtained as follows, Then, the absolute velocity can be represented in the rotational coordinate system as follows,

Equivalent displacements and transmission errors
Define θ i (i = s 1 , r 1 , c 1 , p 1,n , s 2 , r 2 , c 2 , p 2,n , g 1 , g 2 (n = 1, 2, ..., N )) as the angular displacement of each gear body, which is measured in the fixed coordinate system O XY . So the equivalent gear meshing displacements u i in the meshing line direction are defined as, In order to ensure the contact of teeth surface on meshing performance, it is assumed that the relative deformation of gears is completely changed into elastic deformation on teeth surface along with the meshing line direction. According to the geometrical relationships, the relative gear mesh displacements δ j ( j = s 1 p 1,n , r 1 p 1,n , s 2 p 2,n , r 2 p 2,n , g 1 g 2 (n = 1, 2, ..., N )) the moving coordinate system, and the line displacements δ j ( j = c 1 p 1,n , c 2 p 2,n ) of planet gears relative to the carrier on the x-axis and y-axis of the planet carrier can be represented as, sin(ϕ s 2 p 2,n ) + (y s 2 − y p 2,n ) cos(ϕ s 2 p 2,n ) + u s 2 + u p 2,n − u c r s 2 +r p 2,n r c 2 − E s 2 p 2,n δ r 2 p 2,n = (x p 2,n − x r 2 ) sin(ϕ r 2 p 2,n ) + (y r 2 − y p 2,n ) cos(ϕ r 2 p 2,n ) + u r 2 − u p 2,n − u c 2 r s 2 +r p 2,n r c 2 − E r 2 p 2,n δ p 2,n c 2 x = x p 2,n − x c 2 δ p 2,n c 2 y = y p 2,n − y c 2 δ g 1 g 2 = (x g 2 − x g 1 ) sin(ϕ g 1 g 2 ) + (y g 1 − y g 2 ) cos(ϕ g 1 g 2 ) + u g 1 + u g 2 − E g 1 g 2 (4) where E j is the transmission errors, and ϕ j is the relative engagement angle. The vibration of the system caused by errors mainly exists in the manufacturing and assembly processes. The position of the planet gear is changed as a function of gear position error. O P n is the ideal location of center of the planet gear. The gear position error of magnitude A P n makes the ideal location point O P n move toŌ P n . The equivalent displacement of the position error of planet gearĀ P n can be expressed as, where α s , α r , ϕ p n are the pressure angle of the sun gear, the pressure angle of the ring gear and the initial angle of position error and the planet spacing angles, respectively.
Similarly, the run-out error on the reference circle of the gear can be expressed as, where ω i , ϕ ei are the theoretical angular velocity of the component and the initial angle of the run-out error on the reference circle, respectively. where ϕ p n = 2π n−1 N . The transmission errors on the line of action can be represented as, n −s 1 p 1,n +ē s 1 −s 1 p 1,n +ē p 1,n −s 1 p 1,n E r 1 p 1,n =Ā p 1,n −r 1 p 1,n +ē r 1 −r 1 p 1,n +ē p 1,n −r 1 p 1,n E s 2 p 2,n =Ā p 2,n −s 2 p 2,n +ē s 2 −s 2 p 2,n +ē p 2,n −s 2 p 2,n E r 2 p 2,n =Ā p 2,n −r 2 p 2,n +ē r 2 −r 2 p 2,n +ē p 2,n −r 2 p 2,n The gear backlash exists in the meshing gear teeth due to the need of lubrication, the errors caused by manufacturing and assembly processes, and the wear during operation. And the piecewise-linear backlash functions f (δ j ) can be defined as, where b j is the gear backlash, and δ j is the relative gear mesh displacement. The dimension tolerance must be considered that the translational vibration is taken into consideration for the components. The radial support forces f s 1 x and f s 1 y are changed because of the radial clearance ( b) caused by dimension tolerance of the shaft and the hole for sun gear(s 1 ). The principle of the relationship between the force and the clearance can be expressed as, where k s 1 , x s 1 , y s 1 are the radial support stiffness and radial vibrational displacements. γ s 1 determines whether the shaft and the hole are in contact, according to,

Equation of motion
The general coordinate q i is introduced to establish the equations of motion of the gear system.
Lagrange function of planetary transmission system (L 1 ), fixed-axis transmissions(L 2 ) and coupling relations between levels(L 3 ) are expressed as, According to lagrange equation, where Q i is generalized force of the system.
Taking into account the meshing damping, gear backlash, TVMS and TVMF, the translation-torsion coupling equation of motion involving 42-DOF for wind turbine's compound gear transmission dimensionless system is obtained as follows. The translationtorsion equations of motion for the low speed stage's sun gear are expressed as, The translation-torsion equations of motion for the low speed stage's planet carrier are expressed as, The translation-torsion equations of motion for the low speed stage's ring gear are expressed as, The translation-torsion equations of motion for the low speed stage's nth planet gear are expressed as, The translation-torsion equations of motion for the intermediate stage's sun gear are expressed as, c s 2 p 2,nδ s 2 p 2,n cos(ϕ s 2 p 2,n ) The translation-torsion equations of motion for the intermediate stage's carrier are expressed as, c r 2 p 2,nδ r 2 p 2,n The translation-torsion equations of motion for the intermediate stage's ring gear are expressed as, c r 2 p 2,nδ r 2 p 2,n sin(ϕ r 2 p 2,n ) c r 2 p 2,nδ r 2 p 2,n cos(ϕ r 2 p 2,n ) The translation-torsion equations of motion for the intermediate stage's nth planet gear are expressed as, 1 ω d c s 2 p 2,nδ s 2 p 2,n sin(ϕ s 2 p 2,n ) + 1 ω 2 d k r 2 p 2,n f (δ r 2 p 2,n ) sin(ϕ r 2 p 2,n ) + 1 ω d c r 2 p 2,nδ r 2 p 2,n sin(ϕ r 2 p 2,n ) 1 ω d c r 2 p 2,nδ r 2 p 2,n cos(ϕ r 2 p 2,n ) The translation-torsion equations of motion for fixed-axis gear are expressed as, where M i = I i /r 2 i is the equivalent mass, X i are the length of the time-varying friction lever arm, c pg = 2ξ k pg /(1/M p + 1/M g ), where ξ is the meshing damping ratio. In order to perform the dimensionless treatment of the system, the dimensionless time parameter ω d and the dimensionless length parameter b c are introduced, where ω d = k s p n /( Thus, the dimensionless equation of motion in the general coordinates can be given in the matrix form as, where q is the generalized coordinate vector, M is the mass matrix, C is the damping matrix, K is the stiffness matrix, and F is the force vector. Here, the matrix M, K and C are (42 × 42) square matrices.

Calculation of time-varying meshing stiffness
The meshing stiffness k is a time-varying parameter, as the number of engaged teeth varies according to the contact ratio. The meshing of gear pair is divided into internal meshing and external meshing in the model of gear system.

Calculation of external meshing stiffness
The time-varying meshing stiffness of standard involute spur gear transmission system is solved based on the potential energy method of analytical method in this paper. According to the principle of potential energy, when the meshing force F N loads on the tooth, the tooth will experience Hertz contact deformation and bending, compression and shear deformation. Hence, the four elastic potential energies, including the Hertz contact potential energy U h , bending potential energy U b , axial potential energy U a and shear potential energy U s can be expressed as, where K h , K b , K a , K s denote the Hertz stiffness, bending stiffness, radial compression stiffness and shear stiffness corresponding to the deformation of the tooth in the direction of the meshing line under the action of the meshing force F N . Based on Hertz contact theory, Hertzian contact stiffness can be written as, where E, ν and W are Young's modulus, Poisson's ratio and width of tooth, respectively. Based on the beam theory, the potential energy stored in a meshing gear tooth can be expressed as, where α 1 is angle between meshing force and gear thickness direction, d is the height of the meshing force from the tooth root circle, and G is the shear modulus of the material. I x is inertia moment of sectional parts with height h x from meshing force. A x is area of section, , where K h , K b , K a and K s can be obtained according to the above formula. These four stiffnesses together with the stiffness K f corresponding to the deformation of the gear on the meshing line constitute all components of the gear time-varying meshing stiffness.
where the parameters of δ f are the same as that in Ref. [18]. Given K h , K b , K a and K s and K f , τ (τ = 1, 2) is the number of meshing tooth pair, the time-varying meshing stiffness of external meshing can be calculated as,

Calculation of internal meshing stiffness
The calculation method of time-varying meshing stiffness of internal gear tooth is the same as that of external gear tooth. Unlike the calculation of external gear, the cantilever beam cannot be equivalent to the base circle or the tooth root circle for integration. According to the geometric relations, x can be obtained, where z is the number of meshing internal gear teeth, α r is the pressure angle of meshing internal gear. The geometric relationship of θ f , α x , α f is shown in Ref. [18]. The moment of inertia of the internal gear cross section x point and its cross section area can be calculated by reference to the external meshing gear. Distribution force F m of gear teeth perpendicular to the center line, parallel to the center line of gear tooth and bending moment of gear tooth, can be expressed as, Since the internal gear does not consider the deformation of elementary body, the bending stiffness, radial compression stiffness and shear stiffness of the internal gear meshing tooth can be calculated as, The time-varying meshing stiffness of internal meshing can be obtained as,  Meshing damping ratio, ξ 0.07 Rated input speed, n/min 18 Coupling stiffness, N/m k s1c2 = 5 × 10 9 k s2 g1 = 5 × 10 9 Dimensionless nominal length, b c (µm) 1 0

Verification for the calculation method of TVMS
To test the validity of the above provided method(PM), the parameters of gear pairs are the same as that in Ref. [4]. The calculation results with two methods are displayed in Fig. 2a. For the value of single-tooth meshing and double-tooth meshing of TVMS, the relative errors are less than 2%, which are similar to that in Ref. [4] for traditional method(TM). Therefore, the proposed method is accurate to calculate TVMS of the gear system.
Without considering the influence of meshing phase on time-varying meshing stiffness of gear, the timevarying meshing stiffness of the wind turbine's gear system in two complete cycles is obtained. Results are shown in Fig. 2b.
The above time-varying meshing stiffness is expanded by Fourier series and the time-varying meshing stiffness in the form of Fourier series is obtained by deleting the high-order term, wherek j ,k j , ω, ϕ j are the equivalent average meshing stiffness, the time-varying meshing stiffness amplitude, the meshing frequency and the meshing phase angle, respectively.

An improved method for time-varying meshing friction calculation
In the process of gear transmission, the friction coefficient changes significantly with the relative sliding velocity, load distribution, temperature and tooth surface profile. In fact, the lubrication form of the wind turbine's gear meshing pairs is mixed EHL, which includes full oil film EHL and boundary lubrication [19]. Figure 3 depicts the moment of friction force and the internal reaction forces acting on the driving gear ( p) and driven gear(g). The geometric relationship of meshing friction force in gear meshing process is shown in Fig. 3a. Based on the Coulomb friction law, the magnitude of friction force F τ f,v (v = p, g) for the τ th (τ = 1, 2) gear tooth of pinion and gear is expressed as, where μ is the friction coefficient of mixed EHL, F τ N ,v is the dynamic meshing force of the gear pair. λ is the direction of the friction force. The normal forces acting on the driving gear is given by, The parameters have the same meaning as the wind turbine's gear system. The changes of friction force during meshing process can be divided into three period, which is capable to be found in Fig. 3. Figure 3a-b are approaching contact, no friction and leaving contact process, respectively. λ is expressed as, where ω v (v = p, g) is the meshing frequency, X τ p (t), X τ g (t) stand for the length of the time-varying friction lever arm of the τ th (τ = 1, 2) gear tooth of pinion and gear, respectively and are given by, where r v (v = p, g) are reference radius, r vb (v = p, g) are base radius, and r va (v = p, g) are addendum radius. The meshing period T m is 2π/(ω p z p ).
The friction coefficient of mixed EHL can be introduced as, Fig. 9 The dynamic characteristic curve of the system at = 0.98 The boundary lubrication friction coefficient μ B measured through experiments is generally 0.1-0.2. The full oil film EHL coefficient μ E H L can be expressed as [9], where η 0 is the absolute viscosity at oil inlet temperature, ph is the normal pressure, S R is the slide-roll ratio of the gear pair, v r and v s are the relative sliding velocity and rolling velocity at contact point, and α k1 is flare angle. b i (i = 1, 2, ..., 9) are constant coefficients that are dependent on the lubricant type. For ISO-VG-220 gear oil, The friction coefficient of mixed EHL ζ(0 < ζ < 1) is related to the load distribution percentage function. Based on a large number of experiments and numerical calculations, the load distribution percentage coefficient be expressed as, where h min is the minimum oil film thickness, σ is the root mean square (RMS) of roughness, α p is the Fig. 10 The dynamic characteristic curve of the system at = 1.15 viscosity coefficient, v m is the average speed of tooth meshing, andw is the load coefficient.

Verification for the calculation method of TVMF
To test the validity of the above method, the parameters of gear pairs used are the same as those in Ref. [18]. Comparison of calculation results of two Methods from Fig. 4a, the calculation results ITVMF of the above method are consistent with those in Ref. [18] for TVMF. Therefore, the proposed method can accurately calculate the friction coefficient of wind turbine's spur gears.
In order to further study the influence of other parameters on the friction coefficient, Fig. 4b, c gives the effect of surface roughness and input torque. As Fig. 4b shows, as absolute values of slide-to-roll ratio and sur-face roughness increasing, the coefficient of friction increases, surface roughness has great influence on friction coefficient. As could be deduced from Eq. (40), the coefficient is directly related to meshing force F N , which is determined by input torque (T c 1 ). Figure 4c indicates the TVMF under different (T c 1 ). Thus, with the increase of T c 1 , the friction coefficient increases.

Nonlinear dynamic response of the CGTS
The parameters of the CGTS are listed in Table 1, and other calculation parameters are shown in Table 2. The coupling dynamic equations are solved by the Runge-Kutta method. For better analysis the dynamic characteristics of the CGTS, numerical computations are carried out and the dynamic responses of gear sys- Fig. 11 The dynamic characteristic curve of the system at = 2 tem are obtained. In order to study the effect of friction on the dynamic response of the system, the time domain curves of dimensionless relative vibration displacement of internal and external meshing gear pairs are obtained without and with the effect of friction, as shown in Fig. 5. It can be seen from Fig. 5 that the dimensionless relative vibration displacement amplitude decreases when the effect of friction is considered. The friction force decreases the peak amplitude of the dimensionless relative vibration displacement amplitude and makes the periodicity of the system complex, which is consistent with the conclusion in the studies [17]. The correctness of the model and simulation results can be verified. Because the vibration characteristics of the internal and external meshing pairs are consistent, the relative dis-placement of the external meshing pair δ sp is analyzed for the dynamic characteristics as an example.

Effect of excitation frequency
The speed of the wind turbine's gear transmission system is the most important external excitation of the gear system, which corresponds to the excitation frequency ( ) in the gear system. When the dimensionless excitation frequency is considered as the bifurcation parameter, the bifurcation diagram of the wind turbine's gear transmission system is shown in Fig. 6. It can be seen from Fig. 6 that the wind turbine's gear transmission system has rich bifurcation and chaos behaviors, including single-period, multi- period, quasi-period and chaotic motion. The system undergoes periodic motion in the low-frequency region. Then, it enters chaotic motion. There exist two periodic windows (3-T and 4-T) in chaotic region as illustrated. Finally, system evolves from chaotic motion to single-period motion through inverse doubling route. To display the dynamic characteristics more clearly, time region diagrams, Poincare maps, phase diagrams and amplitude-frequency spectrums are obtained to analyze the nonlinear characteristics of the wind turbine's gear transmission system in detail.
When is less than 0.4, it can be seen from Fig. 6 that the bifurcation diagram presents a single point. Time domain diagram exhibits regular singleamplitude sine wave, and phase diagram shows one close loop trajectory. Poincare map has only one point, and amplitude-frequency spectrum has one-peak amplitude, as shown in Fig. 7. Obviously, the gear system is periodic-1 motion under low-frequency excitation. Corresponding to the wind turbine at low rotational speed, the operation of the gear system is relatively stable.
When is increased to 0.42, it can be seen from Fig. 8 that time domain diagram appears two peaks in one period, there are two closed circles in phase diagram, and Poincare map shows two points. Moreover, the amplitude-frequency spectrum has two peaks amplitude; the system state is changed from periodic-1 to periodic-2 motion. Corresponding to the wind turbine's gear system, there is a multi-period window in the single-periodic region.  When is changed from 0.43 to 0.97, the system's motion is periodic-1 motion. With increases to 0.98, the onset of chaotic motion is observed. As shown in Fig. 9, time domain diagram appears irregular vibration and phase diagram shows irregular countless closed motion, and the chaotic motion can be verified by irregular discrete points in projection of Poincare map and amplitude-frequency spectrum displays multiple amplitudes.
When increases to 1.12, the system appears periodic window in chaotic region. For example, when = 1.15, the time domain diagram has three amplitude changes in a period. There are three discrete points on the Poincare map. The phase diagram is a closed pattern formed by three periods of phase tra-jectories. Poincare map is three fixed points. There are three peaks in the amplitude spectrum, which appear at one-third and two-thirds of the fundamental frequency, respectively. Obviously, the system is periodic-3 motion, as shown in Fig. 10.
When increases to 1.18, it is can be seen from Fig. 6 that the system turns to chaotic motion. When continues to increase to 1.52, the systems appear periodic window again, including periodic-4 motion. When reaches to 1.66, the time domain diagram, phase diagram, Poincare map and amplitude-frequency spectrum corresponding to the chaotic response are presented in Fig. 11. Obviously, the system is in chaotic state under this excitation frequency. In the working process of the wind turbine, the operation under the Fig. 16 The dynamic characteristic curve of the system at = 2 excitation frequency of this region should be avoided as far as possible.
With the increase of , the system turns to periodic-2 motion which remains the state from 3.05 to 3.76. When is equal to 3.15, time domain diagram appears two peaks in one period, there are two closed circles in phase diagram, and Poincare map shows two points. The amplitude-frequency spectrum has two peaks amplitude, which are distributed at discrete points of /2. as illustrated in Fig. 12. When continues to increases to 3.77, the system finally returns to periodic-1 motion. the time domain diagram, phase diagram, Poincare map and amplitudefrequency spectrum corresponding to periodic-1 response are presented in Fig. 13. In engineering practice, the excitation frequency is determined by the speed of the wind turbine's gear transmission system. In order to improve stability and reliability and prolong the working cycle life, wind turbine should avoid the speed range of chaotic motion and the critical speed of state changing.

Effect of friction on bifurcation characteristics
The bifurcation characteristics is also one of the main nonlinear dynamic characteristics of gear system. Thus, it is necessary to study the influence of friction on the bifurcation characteristics of gear system. Figure 14 shows the bifurcation diagram of the dimensionless displacement δ sp using the excitation frequency as It can be seen from Fig. 14a, b that the bifurcation behavior of the system is basically unchanged with the increase of RMS of roughness in the low-frequency region, indicating that friction has little effect on the low-frequency region.
When is equal to 1.15, phase diagrams infinitely loop in the enclosed area, but never repeat. The return points in Poincare maps form a geometrically fractal structure, and amplitude-frequency spectrum is continuous. The system's motion is in state of chaos, as shown in Fig. 15. Thus, when is changed from 1.12 to 1.18, the periodic window of the chaotic region disappears.
When is equal to 2, the fluctuation of time domain curve tends to be smooth; discrete points of Poincare maps are more concentrated, as shown in Fig. 16.
Through comparing Figs. 11 and 16, the peak number of amplitude-frequency spectrum decreases. The chaotic characteristics of system attenuate. It can be seen from Fig. 14a, b that with the increase of RMS of roughness, the periodic window of inclusions in the chaotic region of the system disappears in the highfrequency region.
When is equal to 3.15, as shown in Fig. 17, the quasi-periodic motion of the system changes into chaotic motion, resulting in an increase in the chaotic region. The motion state of the system changes; the bifurcation behavior becomes complicated. When is larger than 3.27, the system enters periodic-3 motion, as shown in Fig. 18, which indicates that friction has a great influence on the highfrequency region. From the above analysis, the multiple periodic motion of the system decreases, the chaotic motion region increases, and the critical frequency of entering the chaotic motion is advanced in the highfrequency region when RMS of roughness increases.

Conclusion
A 42-DOF translation-torsion coupling dynamic model of the wind turbine's gear transmission system considering TVMS, TVMF, meshing damping, meshing error and backlash is established. The dynamic characteris-tics of the wind turbine's gear transmission system are studied, which is solved by the Runge-Kutta numerical method, and the conclusions are as follows, • Considering the TVMS and TVMF, the system frequency is an important parameter. For some values of the frequency ratio (i.e., 2 < < 3.15), the system gets into multi-period and chaos motion state. In these cases, the vibration state of the system is disordered, and such vibration frequency should be avoided as far as possible. • Mixed EHL is more in line with the engineering practice, which includes EHL and boundary lubrication. When friction is considered, the friction has little effect on the bifurcation behavior of the system in the low-frequency region. In the highfrequency region, the bifurcation behavior of the system becomes complicated and the chaotic region increases. • With the increase of friction coefficient, the periodic window of the system disappears in the lowfrequency region, the multiple periodic motion regions decrease in the high-frequency region, and the chaotic motion interval of the system increases.