An Improved Rigid Multibody Model for the Dynamic Analysis of the Planetary Gearbox in a Wind Turbine

This paper proposes an improved rigid multibody model for the dynamic analysis of the planetary gearbox in a wind turbine. The improvements mainly include choosing the inertia frame as the reference frame of the carrier, the ring, and the sun and adding a new degree of freedom for each planet. An element assembly method is introduced to build the model, and a time-varying mesh stiffness model is presented. A planetary gear study case is employed to verify the validity of the improved model. Comparisons between the improvementmodel and the traditional model show that the natural characteristics are very close; the improvedmodel can obtain the right equivalentmoment of inertia of the planetary gear in the transient simulation, and all the rotation speeds satisfy the transmission relationships well; harmonic resonance and resonance modulation phenomena can be found in their vibration signals. The improved model is applied in a multistage gearbox dynamics analysis to reveal the prospects of the model. Modal analysis and transient analysis with and without time-varying mesh stiffness considered are conducted. The rotation speeds from the transient analysis are consistent with the theory, and resonance modulation can be found in the vibration signals.


Introduction
The planetary gearboxes are widely used in automotive and aeronautical applications because of their advantages, such as high power density and compactness.The structure of the planetary gear is quite complex: several flexible gear pairs mesh dynamically and the axes of the planet gears rotate at the same time, which leads to dynamics that are nonlinear and nonstationary.Accurate dynamic models for the planetary gear have great benefits for gearbox design, operation, and maintenance.Much study has been applied to this research area.
Dynamic models can be categorized into rigid multibody models (also called lumped parameter models), flexible multibody models, and rigid-flexible coupled models according to their handling methods of body flexibility.Rigid multibody models consist of rigid bodies that are linked by joints and stiffness elements.The degrees of flexibility of the components, including teeth and bearings, are lumped together and modeled by stiffness elements.Saada and Velex [1] proposed a rigid multibody model for determining the critical frequencies for tooth loading on spur and helical gear planetary trains.Lin and Parker [2] built a similar analytical model and revealed that the vibration modes of the spur planetary gear can be classified into rotational, translational, and planet modes.Contact and force coupling between mating flanks have a great effect on dynamics and they are usually modeled by a spring in the rigid multibody model, such as [1,2].For the mesh stiffness, Umezawa et al. [3] used the contact length to calculate the time-varying meshing stiffness externally.Velex and Maatar [4] developed an original lumped parameter model in which there is a three-dimensional mesh model based on total flank deviation that accounts for nonlinear time-varying mesh stiffness, gear errors, and tooth shape modifications.The torsional model can be seen as a special case of a rigid multibody model, in which all degrees of freedom (DOFs) except for the torsional DOFs are constrained.This type of model is usually used in the power train dynamic analysis, such as the work in [5].In a more complicated multiphysical and multicomponent dynamic system model, the planetary gearbox may even be reduced to a lumped mass-spring model, such as the wind turbine gearbox model in FAST [6].
Information inside the body, such as local deformations and corresponding stresses, cannot be calculated from the rigid multibody model; in these scenarios, a flexible multibody model should be adopted, which includes additional generalized DOFs to represent the deformations of individual bodies.The finite element technique [7] is dominant in flexible component modeling.Valco [8] built a pure finite element model of a planetary system, in which both the gears and gear mesh were modeled by finite elements.In this model, a large number of DOFs are needed inside the contact zones to make the contact equations well-conditioned, which increases the computational complexity greatly.Generalpurpose finite element programs are not customized for the special needs of gear contact analysis, so the practical use of pure finite element model is limited.A combination finite element and semianalytical method for a three-dimensional contact problem was proposed by Vijayakar [9].In this model, for points within the contact zone, semianalytical techniques were used to compute the relative deformations and stresses.The "near-field" semianalytical solution and the "far-field" finite element solutions are matched at a "matching surface."Kahraman et al. [10,11] used this model to study the effect of ring flexibility on the quasistatic behavior and presented guidelines for the design of a planetary internal gear.Ambarisha and Parker [12] studied the nonlinear dynamics of planetary gears using the rigid multibody model and this finite element model, and the comparisons validated the effectiveness of the lumped parameter model in analyzing the dynamics of planetary gears.
Flexible bodies and rigid bodies can also be assembled in a rigid-flexible coupled model.Abousleiman and Velex [13] developed a hybrid model for quasistatic and dynamic analyses of planetary gear sets.Ring gears were modeled by 3D finite elements and a modal condensation technique was used to decrease the number of DOFs.Chapron et al. [14] used a dynamic model, consisting of 3D two-node gear elements connected to deformable ring-gears discretized into beam elements, to analyze optimum profile modifications (PMs) in planetary gears (PGTs) with regard to dynamic mesh forces.Peeters [15,16] used a rigid-flexible multibody model for wind turbine drive train dynamics, in which all gears are modeled as deformable parts, the component mode synthesis (CMS) technique is employed for the reduction of the FE models, and the models for the bearings and the tooth contact forces are the same as in the rigid multibody model.
In addition to the flexibility, many other factors also have an effect on the mesh force and make the dynamics more complicated, including the layout of the planet gears [17], contact ratio [18], manufacturing and assembly errors of bearing and gears such as clearances, misalignments, and eccentricities [13,19], and defects on the gear flanks and bearings [20,21].
Compared with the flexible multibody model, the rigid multibody model can obtain an acceptable result with a much lower time cost, so it is widely used in gearbox dynamics analysis.Typical applications include load-sharing and fatigue [22,23], vibration and noise control [22,23], and fault diagnosis [17,24].The gearbox is one of the key components of a wind turbine.It is useful to apply the rigid multibody model in the wind turbine gearbox dynamic analysis.The gearbox usually consists of several planetary stages and parallel stages, and it is coupled with a rotor and generator to form the drive train [5,25].To avoid modulations between the meshing pulsations and the carrier angular velocity, the kinematic equations of traditional rigid multibody model are written relative to rotating frames fixed to the carrier center and planet centers.This characteristic of the traditional model that uses a rotating frame as the reference frame leads to the following problems: (1) The rotation speed of the housing needs to be constrained to the current rotation speed of the carrier.
Most studies simply use a fixed constraint for an approximate treatment, resulting in an incorrect rotation speed for the planetary gears, as will be further illustrated in Section 5.
(2) The models for the planetary gears cannot be directly assembled into a gearbox model because each planetary gear stage has its own reference rotating frame.
These drawbacks make the traditional rigid multibody model hard to use in wind turbine gearbox dynamics research.The objective of this research is to make some modifications to the traditional rigid multibody model to enable its direct assembly with the dynamic models of other components to build a dynamic model for a wind turbine gearbox and drive train.The improvements mainly include the following: (1) the DOFs of the carrier, ring, and sun gear are defined with respect to the inertial frame and (2) each planet has a revolution angular displacement DOF that equals the rotational angular displacement of carrier.A simplified time-varying stiffness model is also proposed.Comparisons between the improved rigid multibody model and the traditional model will be conducted by a planetary gear study case to verify the validity of the improved model, with a modal analysis, sensitivity analysis of the natural frequencies, and transient analysis with the time-varying mesh stiffness considered.To reveal the prospects of the improved model in wind turbine gearbox dynamics analysis, it will be applied for the dynamic analysis of a gearbox that consists of two planetary gear stages and one helical parallel gear stage, including modal analysis and transient analysis both with and without the time-varying mesh stiffness considered.

Structure of Planetary Gear Stage.
The structure diagram of a planetary gear stage with three planets is shown in Figure 1.The ring is fixed and the connection between the ring and housing is modeled as a bearing.Both of the outer rings of the bearings of the carrier and the ring are fixed to the housing.Denote the rotation angular displacements as   ,  = ℎ, , , , ,  = 1, . . ., , where ℎ, , , ,  represent the housing, the carrier, the ring, the sun, and the th planet, respectively.The corresponding rotation speeds are θ  .The radius is denoted as   , and the tooth number of each gear is   .The transmission ratio of the planetary gear stage is The rotation speed of the sun and the planets is For the rigid multibody model, the tooth numbers in (1) and (2) should be replaced by the corresponding radii.The tooth mesh frequency is

Kinematics Analysis of Planetary Gear.
Figure 2 shows the system coordinates, the DOFs of each gear, and the joints.
The reference frame for the carrier, the ring, and the sun is the inertial frame -, where , , and  represent the horizontal, vertical, and axial directions, respectively.The reference frame   -      for th planet is a moving frame rotating around its own axis, while it revolves around the sun gear at the same time; both revolution and rotation speeds are equal to carrier rotation speed;   ,   , and   represent the radial, tangential, and axial directions.
In Figure 2, , ,  represent the displacements in the , , and  direction,  represents the bend angle and its subscripts ,  represent the bend direction,   represents the angular displacement and its subscript  represents the rotation about the  direction, and the subscript  = ℎ, , , ,  ( = 1, . . ., ).  represents the revolution angular displacement of the th planet.  represents the angle from the -axis of th planet gear to the centerline between  it and the sun gear.  and   represent the sun-planet mesh stiffness and the ring-planet stiffness, respectively.  ,   ,   represent the stiffness of the bearings.
The generalized DOFs for the planets are shown as in (5).There are seven DOFs; the first six are the same as in (4), and the seventh is the revolution angular displacement around the sun gear: q p = [  ,   ,   ,   ,   ,   ,   ] ,  = 1, . . ., . (5) For the planets, the kinematics of the axial and bending movements are the same as in the traditional model.The transient position of the planets in     plane r p is where i p , j p are the base vectors of   and   axes, respectively.Then, the acceleration can be derived as Because the planet has a revolution angular displacement DOF, we define a nominal base vector in the revolution direction as k p , according to the relationship between the revolution angular displacement and its equivalent tangential displacement    =     , so k p = j p /  .Then, the acceleration of the planet can be transformed as

Dynamic Model of Planetary Gear
The dynamic model in matrix form is shown as (9), which can be derived from either Newton's second law or Hamilton's principle.Newton's second law is employed in this paper: 3.1.Force Analysis of Planetary Gear.Referring to Figures 1  and 2, the forces sustained by the carrier include the input torque, the carrier bearing contact force, and the force from the planet bearings.The forces acting on the ring consist of its bearing contact force and the tooth contact forces from the planets.The forces sustained by the sun consist of its bearing force and the tooth contact forces from the planets.The forces sustained by each planet consist of its bearing force, the tooth contact forces from the ring, and the sun and the nonconservative inertial forces.
The bearing in our work is modeled as several independent linear springs, the stiffness of each of which can be calculated by Hertz theory or finite element analysis [26,27].
The deflection is independent of the reference frame, so the deformation coordinate relationship for the planet bearings and mesh springs can be obtained in the rotating frame referring to the traditional model in [1,2,16].For the planet bearings, we define the generalized bearing contact force as where  and  represent the force and moment, respectively; subscripts 1 and 2 represent the carrier and the planet gear, respectively; subscripts , , and  represent directions.For example,  1 represents the horizontal force on the carrier caused by the bearing contact force.Only the deflections in the  plane,   and   , are different from those in the traditional model.The deflections can be calculated as The planet torque in the revolution direction is where   and   are the instantaneous stiffness of the bearing in the horizontal and vertical directions.
The deformation of the sun-planet mesh spring is where subscript 1 represents the sun and 2 means the planet.  =   −  , where   ,  = 1, 2, represents the angle from the  axis of gear  to the centerline between the sun and the planet and   is the pressure angle in the transverse plane.  is the helix angle if the gear is helical; otherwise it is zero.Define the generalized sun-planet mesh force as where the meanings of the notations are similar to those in (10).The first twelve elements can be obtained by referring to [1,2,16].The revolution torque on the planet gear caused by the sun-planet mesh force is Because the reference frame of the planet gear is a rotating frame, and  2 is not varying with time, let  2 = 0.The planet moves around the sun, so  1 =   +  0 , where  0 is the initial circumferential position of the planet gear.It can be found that  1 in the improved model varies with the planet circumferential position, while it is a constant value in the rotating frame.
Similar to the ring-planet mesh deformation in (13), the ring-planet mesh deformation is where subscripts 1 and 2 represent the ring gear and the planet gear, respectively, and the other notations are similar to those in (13).Define the generalized mesh force between the ring gear and the planet gear as The first twelve elements can be obtained referring to [1,2,16].The last element is (18)

Dynamic Model Based on the Element Assembly Method.
The element assembly method is a widely used method in finite element analysis to build the total matrices M, C, and K [7], and it can also be used to build the rigid multibody model [1].The element assembly method can be divided into three substeps: (1) numbering the DOFs, (2) assembling the elements, and (3) processing the constraints.There are 4 +  nodes in this model, which are the housing, the carrier, the ring, the sun, and  planets.Each node except for the planet nodes has six DOFs, as shown in (4), and each planet node has seven DOFs, as shown in (5).The total number of DOFs is, therefore, 24 + 7.The elements involved in this dynamic model include the mass elements for the carrier, the sun, the ring, and the planets, the centrifugal stiffness matrix and Coriolis damping matrix for the planets, the stiffness element for the bearings of the carrier, the sun, and the ring, the planetary bearing stiffness element, the sun-planet mesh stiffness element, and the ring-planet mesh stiffness element.All these elements can be derived from the force model described in Section 3.1.
The mass element for the carrier, the sun, and the ring is a 6 × 6 matrix, expressed as where diag(⋅) represents a diagonal matrix,  is the mass of the gear,   and   are the bending inertias of the gear about the horizontal and vertical axes, and   is the rotary inertia about its central axis.The corresponding generalized DOFs are given in (4).The mass matrix for the planets can be derived according to (8) as The first six diagonal elements have the same meaning as the elements in the carrier mass matrix, while the seventh diagonal element is the planetary rotary inertia about the carrier central axis.The corresponding generalized DOFs are given in (5).
The stiffness element for the bearings of the carrier, ring, and sun is a 12 × 12 matrix.The corresponding generalized DOFs are the combinations of the DOFs of the gear and housing.In this paper, the bearings are modeled as six independent springs, so the stiffness matrix is where  radial ,  axial ,  tilt ,  torsion represent the radial, axial, tilt, and torsional stiffness of the bearings, respectively.The centrifugal stiffness matrix K cp and Coriolis damping matrix C cp for the planets are shown as in (22), which can be derived based on (8).Their corresponding generalized DOFs are The moment of inertia of the wind turbine drive train is usually very large, and the rotation speeds are less than 2000 rpm, so it is reasonable to ignore these two elements in the wind turbine gearbox dynamic model.
The planetary bearing stiffness element, the sun-planet mesh stiffness element, and the ring-planet mesh stiffness element can be calculated by where q be is the combination of the DOFs of the two bodies connected by the bearing and q mie consists of the DOFs of the two gears.The shape of all these stiffness matrices is 13 × 13.Embedding all these matrices according to the planetary gear structure, the total matrix without constraints can be obtained.
There are two constraints in the improved model: (1) the housing is fixed and (2) the planet revolution speeds are equal to the carrier rotation speed.For the first constraint, the rows and columns corresponding to the housing DOFs are deleted from the total matrix.For the second constraint, the rows and columns corresponding to the revolution speed DOF of the planets are first added to the rows and columns corresponding to the carrier rotation speed DOF and then deleted from the total matrix.Thus, the total matrices M, C, and K in (9) are obtained.
The total matrix K is a function of   ,   , and   ( = 1, . . ., ).From the physical perspective,   in ( 13) and ( 16) reflects the effect of the planet gear transient position on the gearbox vibration.This is one of the advantages of the improved model compared with the traditional model.In the following study case, we will omit the effect of   on the mesh force.Details about the time-varying mesh stiffness (  and   ) will be given in Section 4.
The damping matrix C can be built by the element assembly method if the damping coefficients of all the joint elements are known.In this way, the Coriolis damping matrix can also be added into the total damping matrix.Considering that the damping coefficient is hard to obtain and has a slight effect on the research objective of this paper, matrix C is built by the Rayleigh damping model for simplification.A similar usage can be found in [28].The Rayleigh damping model is widely used in elastic bodies defined as where coefficients  and  can be calculated based on the modal damping ratios [29].In this paper, coefficients  and  are chosen by experience.Matrix K in ( 24) is replaced by matrix K mean in practice; see (25).Because the stiffness matrix K is nonlinear, to improve the stability of the transient simulation, the dynamic model is transformed into the form of (25), with all the time-varying elements in K moved to the right side: M q + C q + K mean q = F eq () , where F eq () is called the equivalent load and K mean is constant and obtained by taking the time-varying mesh stiffness with the synthesizing stiffness.
Benefiting from the DOFs of the carrier and the sun gear being defined in the inertial frame, the improved planetary gear stage model can be regarded as a super element and directly assembled with other dynamic models to build the gearbox model or coupled with the rotor and generator to build the wind turbine drive train model.This will be illustrated in the following sections.

A Model for the Time-Varying Mesh Stiffness
The time-varying mesh stiffness has a great effect on the planetary gear dynamics, such as causing parametric excitation of the transmission error [29], affecting the vibration instability, and controlling the dynamic tooth loads [28,30], so it should be considered in the dynamic model.In this paper, the time-varying mesh stiffness model is categorized into predefined (static) and nonlinear (dynamic).A predefined model means that the stiffness is a function of time and other parameters independent of the current state, so it acts just like a predefined external load.Gu et al. [31] proposed an approximate formula for the time-varying mesh stiffness function for ideal solid spur and helical gears.The results compare very well with those obtained by using twodimensional finite element models.The predefined stiffness mothed has been widely used in the dynamics of gears with defects [32,33].For a nonlinear model, the transient mesh stiffness is a result of history movements and the current state.This method dynamically considers the effect of factors, such as the occasionally contact loss caused by tooth separation.Velex and Maatar [4] developed a three-dimensional mesh model based on the total flank deviation that accounts for the nonlinear time-varying mesh stiffness, gear errors, and tooth shape modifications.In this paper, a simplified dynamic timevarying mesh stiffness model is built, in which transmission error is omitted and mesh phase relationships among the sunplanet and ring-planet meshes are considered.For a meshing spur tooth pair, the contact length is constant and the contact position is changing with time.For each tooth, define , , and  as the position of first contact, the position of last contact, and the position with maximum mesh stiffness, respectively.Define the relative contact positions of  and  as   = 0 and   = , where  is the transverse contact ratio.Suppose the stiffness of  and  as   =   =  min .The stiffness of point  can be calculated by the synthesizing stiffness   =  mean /(0.75  + 0.25), where   is the transverse contact ratio.Suppose  is the midpoint between  and , so   = /2.Suppose the stiffness  min =   /, where  = 1.5 for an external gear mesh and  = 1.67 for an internal gear mesh.A parabola is used to fit the functional relationship between the mesh stiffness and the relative contact position  based on these three points (0,  min ), (/2,   ), and (,  min ), shown as follows: For a gear pair, several teeth will mesh at the same time.When there is no transmission error in the gear pairs, the following method can be used to calculate the tooth pairs in meshing and the corresponding contact positions.
For one of the gears in a fixed axis meshing gear pair, suppose one of its teeth in meshing is tooth  0 with the relative contact position equaling 0 when the gear rotation angular displacement is  0 ,  = 1, 2. For each gear, the teeth are numbered starting from one along its rotational direction.The sun-planet mesh and ring-planet mesh gear pairs can be transformed to fixed axis gear pairs in the rotating frame that is fixed to the carrier.Then, when the gear rotation angular displacements are   , one of the meshing teeth on the first gear can be calculated by the following: where   is the tooth number,   is the relative contact position, and   is the relative rotation angular displacement: This model cannot be used directly to calculate the corresponding meshing tooth of the other gear, because the real rotation angular displacement does not satisfy the no transmission error assumption.Substituting  2 =  1  1 / 2 into (27) can obtain the corresponding meshing tooth  2 and the contact position  2 .It can be found that  1 =  2 = .This meshing tooth pair is denoted as ( 1 ,  2 , ).All the meshing tooth pairs can be calculated by {( 1 − ,  2 − ,  + ) |  +  < ,  = 0, 1, 2, . ..} .(29) Combining ( 26) and ( 29), the total meshing stiffness can be calculated by where   1 , 2 (⋅) is the mesh stiffness function between tooth  1 of the first gear and tooth  2 of the second gear.Tooth defects that affect the meshing stiffness can be simulated by this model.Phase relationships among the sun-planet and ringplanet meshes are very important for the planetary gear.Parker and Lin gave a complete analysis in [34].Phase relationships determine the relationships between the initial meshing tooth   0 ,  = 1, 2;  = , ;  = 1, . . ., , and the relative contact position   0 (in ( 27)) of these meshing gear pairs.Here, subscript 2 represents the planet gear for all the gear pairs.Superscript  represents the gear pair of the sun gear and the th planet gear, and  represents gear pair of the ring gear and the th planet gear.For all the gears, the initial rotation angular displacement is zero,   0 = 0.For simplification, suppose   0 = 0. Define  1 10 = 0,  1 10 = 0,  This time-varying mesh stiffness model can also be applied for helical gears by replacing the contact position with the contact length and modifying the mesh stiffness function () of one tooth pair.If a high-precision mesh stiffness model is needed, the three-dimensional mesh model in [4] is recommended.
Figure 3 shows the time-varying stiffness of the planetary gear used in Section 5 simulated by this model.The rotation speed of the carrier is constant, θ  = 120 rpm.It can be found that there is a phase shift between the sun-planet mesh stiffness and the ring-planet mesh stiffness.

Verification Based on a Planetary Gear Study Case
In this section, comparisons between the improved rigid multibody model and the traditional model will be conducted to verify the validity of the improved model.Lin's model is chosen as the traditional model and the rotation angular displacement of the housing in their model is treated with fixed constraint.Because the gears in the model for comparison are spur gears, all the axial and bending movement DOFs in the improved model are fixed.

Parameters of the Planetary Gear Model.
Table 1 shows the parameters of the planetary gear we analyzed [2,16].All the gears are spur gears.The number of planets is  = 3.   show that this finding also applies for the improved model.Table 2 shows a comparison between the natural frequencies calculated by the two models.It can be found that, while the translational natural frequencies are the same, the rotational natural frequencies are slightly different.Except for the first rotational natural frequency, the maximum relative error is 9.7%.These errors may be caused by the improper treatment of the housing rotation angular displacement DOF in the traditional model.

Modal
Figures 4 and 5 show the comparison of the modal energy distributions between the improved model and the model of Lin and Parker [2].Modal energy is calculated by the strain energy, which can be found in [35].There are four bars from 1 to 4 in each subplot, which represent the rotational energy of all the components except for the planets, the translational energy of all the components except for the planets, the rotational energy of all the planets, and the translational energy of all the planets.It can be found that  all the translational mode distributions are the same.For the second rotational mode (the fourth subplot), there is more energy located in the planet translational movements in our results than in the results of Lin and Parker [2].
Figure 6 shows the natural mode shapes of the second rotational mode.It can be found that the planet movements in the two results are very similar, while the planet rotational movements in the results of Lin and Parker [2] are bigger  than our results.This is consistent with the modal energy distribution analysis.

Modal Frequency Sensitivity Analysis.
Figure 7 shows the changes in the planetary gear natural frequencies along with the sun-planet stiffness.The stiffness varies from 17 to 19 Nm.It can be found that all the sensitivity lines between the two results are very similar, with the major difference being that the fifth and sixth rotational natural frequencies of the improved model are larger than those of Lin.The natural frequency sensitivity to the sun-planet mesh stiffness can be calculated by the modal strain energy distribution [35].
Figure 8 shows the normalized strain energy distributions calculated by the improved model corresponding to the vertical lines in Figure 7. Numbers 1 to 6 on the  axis represent the six rotational natural frequencies from low to high, and numbers 7 to 18 represent the twelve translational natural frequencies in ascending order.It can be found that as the stiffness increases all the sensitivities concenter on the highest rotational and the highest translational (including horizontal and vertical) natural frequencies.[2].For each column, the subplots from top to bottom are the carrier rotation speed, the sun rotation speed, and first planet rotation speed.

Transient Dynamic Analysis.
There are several methods that can be used to solve (9); the implicit Newmark algorithm is adopted in this paper.A description of the Newmark method can be found in [29].Considering a nonlinear mesh stiffness, an implicit algorithm [36] is used to enhance the stability of the solving process.
The initial state of all the DOFs is zero.The driving torque on the carrier is constant,  1 () = 10 Nm, and the load torque  2 () on the sun gear is where ratio represents the planetary gear transmission ratio, calculated by (1); here, ratio = 4.56. is a constant rotation speed; here, let  = 120 rpm. 0 is the time when the rotation speed of the carrier reaches the speed  for the first time; it is obtained by the solver dynamically.Qualitative analysis shows that the rotation speeds of the carrier, the sun, and the planets will increase linearly until a time larger than  0 , at which point the load starts to switch in and the rotation speed of the carrier θ  will ultimately fluctuate at approximately 120 rpm at last.The total inertia cast on the sun gear side is  = 0.124, and the time  0 can be estimated by  0 ≈ ( ⋅ )/( ⋅ ratio) = 0.156 s.The coefficients of Rayleigh damping are chosen by experience as  = 0.01256,  = 2.12 − 7. The total simulation time is 5 seconds, the time step is 1 − 6 s, and the data sampling rate is 40 kHz.The time-varying mesh stiffness model in Section 4 depends on the rotation speed of the planets.For the model of Lin and Parker [2], the rotation speed of the carrier is used to estimate the instantaneous rotation speed of the planets.
The rotation speeds in the first 2 s calculated by the two models are given in Figure 9.The results show that each of these speeds is first increased linearly and then converges to a certain value, which is consistent with the theory.Taking the average value of the last 1 s as the equilibrium speed, both of the simulated carrier rotation speeds are approximately 119 rpm, while the theoretical value without damping is 120 rpm.This slight error is acceptable.
Theoretical rotation speeds of the sun and the planets can be calculated based on the carrier rotation speed and the transmission ratio.The time  0 and the rotation speeds are listed in Table 3.It can be found that the time  0 calculated by the improved model is consistent with the theoretical value, which means that the equivalent moment of inertia of the planetary gear stage is right; the planet rotation speed calculated by Lin and Parker [2] is less than the theoretical value, and the planet tangential velocity is larger than the theoretical value, while all the values calculated  by the improved model coincide well with the theoretical values.This comparison verifies the validity of the improved model.
The vertical vibration acceleration of the ring and the rotational acceleration of the sun after 1 s are used for signal analysis.Figure 10 shows the waveform and spectrum of the ring vertical acceleration vibration by the two models.Figure 11 shows the spectrum details by local upscaling.It can be found that the highest amplitude appears to be approximately 1102 Hz, and the next ones are at 743 Hz and 1891 Hz.These frequencies are translational natural frequencies; refer to Table 2. Complicated modulation exists in the high frequency region.These findings validate that the resonant modulation phenomenon exists in the vibration signal of the planetary gear.Comparisons reveal that the spectra from these two methods are very similar, while the modulation in the high frequency region is more obvious by the improved model.
The 5000 to 10000 Hz frequency band is chosen for the demodulated resonance analysis (DRA). Figure 12 shows the results of the DRA and its partial enlarged views.The meshing frequency 192 Hz and its harmonics can be found.
No significant information can be found in its low frequency region.
Figures 13 and 14 show the spectrum analysis of sun rotation acceleration and the results of the envelop spectrum analysis.Comparisons reveal that their high frequency ranges are very similar, while the component at approximately 1330 Hz of the improved model is very outstanding.This component is the result of harmonic resonance, since 1334 Hz is the second rotational natural frequency of the improved model, and it is close to seven times of the current meshing frequency.Compared with the analysis results of the ring vertical vibration acceleration, it can be concluded that the harmonic resonance is more obvious in the rotational vibration than in the translational vibration.

An Application in Planetary Gearbox Dynamics
To show the prospects of the improved planetary gear model in the wind turbine gearbox dynamic analysis, a dynamic analysis of a gearbox consisting of two planetary gear stages and one helical parallel gear stage is conducted.The parameters of the two planetary gear stages are listed in Table 1.The first planetary stage has 4 planets while the second stage has 3 planets.The parameters of the parallel gear stage are listed in Table 4 [16].The total transmission ratio is ratio = 20.73 and the total rotation inertia cast on the input shaft side is  = 5.21.The coupling between the planetary stages is implemented as a rigid link between the rotational DOF of the first stage's sun gear and the rotational DOF of the second stage's carrier.The coupling between the second planetary stage and the third stage is implemented as a torsional spring with stiffness 15 Nm/rad connecting the rotational DOF of the sun and the rotational DOF of the parallel stage's input gear.

Modal Analysis.
The natural frequencies calculated by the improved model are categorized based on the mode  energy distribution and listed in Table 5, where ( * 2) means that the frequency appears two times.Figure 16 shows two of the modal energy distributions.
The second stage of the gearbox is the same as that described in Section 5, and Table 6 shows the comparison between their natural frequencies.It can be found that the rotational natural frequencies are different when only the planetary gear stage is analyzed and embedded in the gearbox.This is caused by the rotation coupling among the stages.Because there are no translational and bending DOFs involved in the couplings, the translational natural frequencies are the same.

Transient Dynamic Analysis.
The driving torque put on the input shaft is  1 = 100 Nm, and the load torque put on the output shaft follows (31).Here, ratio = 20.73, = 60 rpm, and  0 can be estimated as  0 ≈ 0.3273 s.The coefficients of Rayleigh damping are chosen by experience as  = 0.01256 and  = 2.12 − 7. The simulation time is 4 seconds, the time step is 1 − 6 s, and the data sampling rate is 40 kHz.Both models are solved with and without the time-varying mesh stiffness considered.
Figure 17 shows the simulation rotation speeds in the first two seconds with time-varying mesh stiffness considered.It can be found that the changing processes of the rotational speeds agree well with expectations.The simulation equilibrium rotation speeds take the average value of the last one second of data, and the theoretical rotation speeds of the gears are calculated based on the simulation of the input shaft rotation speed.All these speeds and the time  0 are listed in Tables 7 and 8. Comparisons show that the time  0 calculated by the simulation is consistent with the theoretical value, which means that the equivalent moment of the gearbox is right, and the simulation rotation speeds are well consistent with the theoretical results in both cases.So, this model can be directly coupled with the rotor model and the generator model for the wind turbine drive train dynamic analysis.
Figure 18 shows the spectrum of the vertical vibration acceleration of the ring of the first planetary stage.Figure 19 shows the spectrum of the rotation acceleration of the output gear.No matter whether the mesh stiffness is constant or time-varying, the resonance phenomenon can be found in both spectra.The spectrum with the time-varying mesh stiffness is more complicated and resonance demodulation can be found in its high frequency region.
In this study, the effect of the planet circumferential position on the vibration is omitted, and the time-varying mesh stiffness model is simplified.In future work, more factors should be taken into account to evaluate their effects and obtain a more accurate result.

Conclusions
An improved rigid multibody model and a time-varying mesh stiffness model for a planetary gear are proposed.Compared with the traditional model, the improved model for the planetary gear stage can be regarded as a super element, directly assembled with other dynamic models to build the gearbox model or coupled with the rotor and generator to build a wind turbine drive train model.
The improved model is effective for planetary gear dynamics.Comparisons between the results are calculated by the improved model and the traditional model, showing that (1) the translational natural frequencies and the sensitivity of the natural frequencies are the same; (2) the rotational natural frequencies have a small difference, and the relative error of the second rotational natural frequency (reaching 9.71%) is the largest, which may be caused by the improper treatment of the housing rotation speed constraint in the traditional model; (3) in transient analysis, the rotation speed of the planets and the equivalent moment of inertia of the planetary

Figure 1 :
Figure 1: Structure diagram of a planetary gear stage, where number 1 represents the carrier, 2 represents the ring, 3 represents the sun, 4 represents one planet, 5 represents the carrier bearing, 6 represents the sun bearing, 7 represents the planet bearing, 8 represents the inner mesh between the ring and planet, 9 represents the outer mesh between the sun and the planet, 10 represents the ring bearing, and 11 represents the housing.

Figure 2 :
Figure 2: The system coordinates and the DOFs of each gear.

Figure 3 :
Figure 3: The time-varying stiffness in a planetary gear stage: (a) sun-first planet, (b) sun-second planet, and (c) ring-first planet.

Figure 4 :Figure 5 :
Figure 4: Modal energy distributions by the improved model.

Figure 6 :
Figure 6: Natural mode shapes of the second rotational mode: (a) by the improved model and (b) by the model of Lin and Parker [2].Dashed lines are the equilibrium positions, solid lines are the deflected positions, and dots represent the component centers.

Figure 7 :
Figure 7: Natural frequency versus the sun-planet stiffness: (a) by the improved model and (b) by the model of Lin and Parker [2].

Figure 9 :
Figure 9: Instantaneous rotation speeds: (a), (c), and (e) by the improved model and (b), (d), and (f) by the model of Lin and Parker[2].For each column, the subplots from top to bottom are the carrier rotation speed, the sun rotation speed, and first planet rotation speed.

Figure 10 :
Figure 10: Waveform and spectrum of the vertical vibration acceleration of the ring: (a) and (c) from the improved model and (b) and (d) from Lin's model [2].

Figure 11 :
Figure 11: Partial enlarged views of the spectrum of the ring vertical vibration acceleration: (a) and (c) from the improved model and (b) and (d) from Lin's model [2].

Figure 12 :
Figure 12: DRA of the ring vertical vibration acceleration with frequency band from 5000 Hz to 10000 Hz: (a) and (c) are the spectrum of DRA and its partial enlarged view from the improved model and (b) and (d) are from Lin's model [2].

Figure 13 :Figure 14 :
Figure 13: Waveform and spectrum of the sun rotational acceleration: (a) and (c) from the improved model and (b) and (d) from Lin's model [2].

Figure 15 :
Figure 15: Structure of the gearbox: 1 represents the first planetary gear stage, 2 represents the second planetary gear stage, 3 represents the parallel gear stage, 4 represents the input shaft, and 5 represents the output shaft.

Figure 16 :
Figure 16: Modal energy distributions: (a), (c), and (e) are distributions of the second rotational natural frequency of the second stage; (b), (d), and (f) correspond to the third rotational natural frequency of the second stage.The 3 bars in the subplots in the first row represent the energy of the three stages.The 3 bars in the subplot in the second row represent the rotational energy of all the components except for the planets, the translational energy of all the components except for the planets, and the energy of all the planets, respectively; in the subplots in the third row, the first four bars and the second four bars represent the energies in the first stage and the second stage (refer to Figure 4), bar 9 represents the rotational energy of the third stage, and bar 10 represents the other energy of the third stage.

Table 1 :
Parameters of the planetary gear model.

Table 2 :
Comparison between the natural frequencies calculated by the two models.

Table 3 :
Equilibrium speeds from simulations and theory.

Table 5 :
Natural frequencies of the gearbox.

Table 6 :
Comparison between the natural frequencies calculated by the two models.

Table 7 :
Simulation rotation speeds with constant mesh stiffness.