Abstract
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 improvement model and the traditional model show that the natural characteristics are very close; the improved model can obtain the right equivalent moment 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.
1. 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. General-purpose 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.
2. Modeling and Equations of Motion
2.1. 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 , , 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
2.2. 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 . 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 carrier, sun, and ring are shown as in (4). There are six DOFs, which are translational motion in , , and , bending along and , and rotational motion
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:
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 is where , 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 , according to the relationship between the revolution angular displacement and its equivalent tangential displacement , so . Then, the acceleration of the planet can be transformed as
3. 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.
For the convenience of equation expression, we define displacements , , , and , where , . It can be found that they have the same dimension.
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 aswhere 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, 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 iswhere and are the instantaneous stiffness of the bearing in the horizontal and vertical directions.
The deformation of the sun-planet mesh spring iswhere subscript 1 represents the sun and 2 means the planet. , where , , 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 aswhere 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 is not varying with time, let . The planet moves around the sun, so , where is the initial circumferential position of the planet gear. It can be found that 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 iswhere 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
3.2. 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 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, . 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 matrix, expressed as where 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 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 , , , represent the radial, axial, tilt, and torsional stiffness of the bearings, respectively.
The centrifugal stiffness matrix and Coriolis damping matrix for the planets are shown as in (22), which can be derived based on (8). Their corresponding generalized DOFs are Consider
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 bywhere is the combination of the DOFs of the two bodies connected by the bearing and consists of the DOFs of the two gears. The shape of all these stiffness matrices is . 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 (). 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 aswhere coefficients and can be calculated based on the modal damping ratios [29]. In this paper, coefficients and are chosen by experience. Matrix in (24) is replaced by matrix 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:where is called the equivalent load and 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.
4. 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 two-dimensional 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 time-varying mesh stiffness model is built, in which transmission error is omitted and mesh phase relationships among the sun-planet 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 and , where is the transverse contact ratio. Suppose the stiffness of and as . The stiffness of point can be calculated by the synthesizing stiffness , where is the transverse contact ratio. Suppose is the midpoint between and , so Suppose the stiffness , where for an external gear mesh and 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 , , and , 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 with the relative contact position equaling when the gear rotation angular displacement is , . 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 into (27) can obtain the corresponding meshing tooth and the contact position . It can be found that . This meshing tooth pair is denoted as (). All the meshing tooth pairs can be calculated by
Combining (26) and (29), the total meshing stiffness can be calculated bywhere is the mesh stiffness function between tooth of the first gear and tooth of the second gear. Tooth defects that affect the meshing stiffness can be simulated by this model.
Phase relationships among the sun-planet and ring-planet 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 , ; ; , and the relative contact position (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, . For simplification, suppose . Define , , , , and all the other initial meshing tooth numbers can be calculated based on these values:(i)For the planet gears, , .(ii)For the sun gear, , .(iii)For the ring gear, , .
and , ; , depend on the tooth number of the gears and the layout of the planet gears. If the tooth number of planet is even, then . If the planets are equally spaced, then .
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, rpm. It can be found that there is a phase shift between the sun-planet mesh stiffness and the ring-planet mesh stiffness.
(a)
(b)
(c)
5. 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.
5.1. 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 .
5.2. Modal Analysis
According to the research of Lin and Parker [2], the planetary gear has six rotational modes, six translational modes, and zero planet modes. The results 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.
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.
(a)
(b)
5.3. 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 to 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.
(a)
(b)
(a)
(b)
(c)
5.4. 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, Nm, and the load torque on the sun gear is where ratio represents the planetary gear transmission ratio, calculated by (1); here, . is a constant rotation speed; here, let . 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 , 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 , and the time can be estimated by s. The coefficients of Rayleigh damping are chosen by experience as , . The total simulation time is 5 seconds, the time step is 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.
(a)
(b)
(c)
(d)
(e)
(f)
Theoretical rotation speeds of the sun and the planets can be calculated based on the carrier rotation speed and the transmission ratio. The time and the rotation speeds are listed in Table 3. It can be found that the time 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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
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.
(a)
(b)
(c)
(d)
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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
6. 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.
6.1. Parameters of the Planetary Gear Model
The structure of the gearbox for study is shown as in Figure 15.
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 and the total rotation inertia cast on the input shaft side is . 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 Nm/rad connecting the rotational DOF of the sun and the rotational DOF of the parallel stage’s input gear.
6.2. Modal Analysis
The natural frequencies calculated by the improved model are categorized based on the mode energy distribution and listed in Table 5, where () means that the frequency appears two times. Figure 16 shows two of the modal energy distributions.
(a)
(b)
(c)
(d)
(e)
(f)
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.
6.3. Transient Dynamic Analysis
The driving torque put on the input shaft is Nm, and the load torque put on the output shaft follows (31). Here, , rpm, and can be estimated as s. The coefficients of Rayleigh damping are chosen by experience as and . The simulation time is 4 seconds, the time step is 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 are listed in Tables 7 and 8. Comparisons show that the time 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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
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.
7. 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 gear stage calculated by traditional model are wrong, while those calculated by the improved model agree very well with the theory; (4) the spectra of the simulation signals are very similar, and the results validate that the resonant modulation and harmonic resonance phenomena exist in the vibration signal.
The improved model can be used for the real wind turbine planetary gearbox dynamic analysis. It has been successfully applied in the dynamic analysis of a gearbox that consists of two planetary gear stages and one helical parallel gear stage. The results show that (1) all the simulation rotation speeds satisfy the theoretical transmission ratio relations; (2) the equivalent moment of inertia of the gearbox is consistent well with the theoretical value; (3) no matter whether the mesh stiffness is constant or time-varying, the resonance phenomenon can be found in the spectrum; and (4) resonance demodulation can be found in the high frequency region of the ring vertical vibration signal.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This work was supported by the National Natural Science Foundation of China (Grant 51174273).