Time-Domain Dynamic Modeling and Analysis of Complex Heavy-Duty Gearbox Considering Floating Effect

Featured Application: The proposed time-domain dynamic modeling method can be applied to analysis and design of heavy-duty gearboxes in wind turbines, cranes, drilling rigs and so on. Abstract: Expert insights into the time-domain dynamic behavior of heavy-duty gearboxes form the foundations of design evaluation and improvement. However, in the existing lateral–torsional coupling (LTC) modeling method for gearboxes that is normally used for frequency-domain dynamic behavior, the meshing forces are modeled as spring dampers with ﬁxed acting points on the meshing gears to simulate only the transient LTC effect, and thus the steady state characteristic in the time domain cannot be obtained due to the unrealistic distortion of positions and orientations as the gear angles increase. In this paper, a novel and generally applicable LTC modeling method for heavy-duty gearboxes, mainly planetary gear sets with ﬂoating components, is proposed by using space-ﬁxed spring dampers with ﬂoating acting points on the meshing gears to study the time-domain dynamic response and to support the dynamic design of heavy-duty gearboxes. Based on the proposed method, a LTC model of a 2 megawatt (MW) wind turbine gearbox with ﬂoating components considering the time-varying meshing stiffness, bearing stiffness, torsional stiffness, and ﬂoating effect was established. The simulated results of representative components were in accordance with experimental results on a test rig, and dynamic behavior was calculated.


Introduction
Gear systems are the most significant and extensively adopted mechanical transmission devices in modern applications such as automobile, locomotive, mining, metallurgy, electric power, petroleum and chemical industries. As mechanical equipment tends to develop in the direction of larger scale and higher performance, gear systems with high power, high speed and heavy load are in great demand. The working environment for gear transmission systems is extremely complex, especially for the heavy-duty gearboxes adopted in large ships, offshore platforms, wind power machines, etc. There exist not only external incentives introduced by the variable load condition and power plant, but also internal incentives from the time-varying meshing stiffness, the gear transmission error, the meshing impact, etc. Taking the wind turbine gearbox (WTG) as an example, as the drivetrain that connects the rotor hub with the generator converts the wind energy into electrical energy, the wind field turbulence through the hub and the electricity grid disturbance through the generator become the external excitations [1]. Meanwhile, since most WTGs contain planetary gear sets with floating components because of the advantages of compact structure, high power density, and coaxial arrangement of transmission shafts [2], the varying mesh stiffness, meshing incentive, impact of floating components, and composition error generate the internal excitations [3,4]. These make the service condition of WTGs quite unique and rigorous. The service life of a wind turbine is usually designed to last more than 20 years [5]. However, internal WTGs have the highest failure rate of all wind turbine components [6,7], and the main cause of wind turbine standstill is the vibration issue [8]. Hence, the dynamic behavior and working performance of gearboxes have a significant impact on the whole machine, and research on the vibration characteristics is critical to the design and redesign of gearboxes.
Considerable efforts have been made to give comprehensive insights into the dynamic behavior and vibration issues of gearboxes, especially of planetary gear sets [9,10]. Kahraman et al. [11] derived a time-varying dynamic model of a planetary transmission including the manufacturing error, assembly variation and time-varying stiffness. Lin and Parker et al. [12] developed an analytical lumped-parameter model of a planetary gear set and used it to investigate the natural frequencies and vibration modes. Abousleiman et al. [13,14] investigated the vibration behavior of a planetary gear set with an elastic ring gear and an elastic planet carrier. Chen and Shao [15] derived a mesh stiffness model of an internal gear pair with a tooth root crack in the ring gear based on the potential energy principle. Christopher et al. [16] proposed a finite element formulation for the dynamic response of gear pairs. Maláková et al. [17] analyzed the main influencing factors on the meshing stiffness, as changes in the stiffness are a significant source of noise and vibration. She also designed a non-circular gear transmission system with continuously changing gear ratio and studied its kinematical characteristics [18]. Girsang et al. [19] studied the dynamic response of a WTG by enhancing the capability of FAST (a wind turbine CAE tool) [20] through the integration of a gearbox model built using Simscape, which took into account excitations from both the wind field and the generator. Zhu et al. [21,22] used the lumped parameter method to study the dynamic behavior of a WTG with flexible pins and proved that the flexible pins can improve the load sharing ability. Park et al. [23] investigated the influence of non-torque loads on the dynamic behavior of planetary gears, including the loading distribution and loading sharing, and concluded that considering the non-torque load during the design stage is quite important to accurately determine the design load that will guarantee the service life of a WTG.
In large-scale mechanical systems, heavy-duty gearboxes with planetary gear sets are adopted extensively due to their unique advantages of high power density and high transmission ratio; however, this yields complex configurations and characteristic serving circumstances. Thus, more stringent design criteria are desired than for common industrial gearboxes, and the dynamic forces become the major limiting factor to achieving better performance [24]. In previous studies concerning gearboxes, the assessment of dynamic behavior is usually centered on the frequency domain. Hence, a study of dynamic behavior in the time domain becomes significant, as it can confirm the dynamic loading factors, which are usually chosen by experience from a manual according to the traditional design code for common gearboxes. In this paper, an appropriate and generally applicable modeling approach of heavy-duty gearboxes, mainly the planetary gear set, is proposed to investigate dynamic behavior in the time domain, and a LTC model of a 2 MW WTG was built in which the time-varying meshing stiffness, bearing stiffness, torsional stiffness, and floating factor were considered. The model was then verified by experiments on a test rig, and dynamic characteristics of the WTG were finally revealed on the verified model.

LTC Modeling Approach of Heavy-Duty Gearbox
To establish a dynamic model of a heavy-duty gearbox, modeling of a planetary gear set with floating components serves as the most basic as well as the most significant portion. This is mainly because common gear transmissions can be simplified or derived from planetary gear sets. Therefore, a planetary gear set model will be highlighted and introduced as a representative basis in this section.

LTC Models of Planetary Gear Set
Interactive incentives between the gear engagements occur mainly through the meshing forces. The existing LTC model of a planetary gear set is shown in Figure 1a. The translation flexibility is simplified as linear spring dampers in the XY plane, while the flexibility between gear engagement, which is mainly caused by bending deflections and contact deformations of gear teeth, is simplified as a spring damper with equivalent average stiffness along the line of engagement with fixed acting points [12,15,[25][26][27][28][29][30]. Thus, each component is modeled with three degrees of freedom (DOF), and the bearing stiffness and time-invariant meshing stiffness are considered. This method has an inherent defect in that only small angular displacements of gears are permitted, and thus only the transient vibration characteristic of gear pairs can be solved and obtained.

LTC Models of Planetary Gear Set
Interactive incentives between the gear engagements occur mainly through the meshing forces. The existing LTC model of a planetary gear set is shown in Figure 1a. The translation flexibility is simplified as linear spring dampers in the XY plane, while the flexibility between gear engagement, which is mainly caused by bending deflections and contact deformations of gear teeth, is simplified as a spring damper with equivalent average stiffness along the line of engagement with fixed acting points [12,15,[25][26][27][28][29][30]. Thus, each component is modeled with three degrees of freedom (DOF), and the bearing stiffness and time-invariant meshing stiffness are considered. This method has an inherent defect in that only small angular displacements of gears are permitted, and thus only the transient vibration characteristic of gear pairs can be solved and obtained. This paper proposes a novel method for LTC modeling. The elasticity and viscidity of gears, carriers, shafts, bearings, etc., are represented by time-varying stiffness-damping forces, while the bodies of gears, carriers and shafts are all modeled as rigid elements. Thus, a final multibody system of the gear set with consideration of time-varying stiffness, friction, damping, etc. is established. Based on Newton's third law, the interaction of two bodies can be represented by an action force and an opposite reaction force that are functions related to the corresponding characteristic parameters of the two bodies. According to the basic law on parallel translation of force, forces on the contact point are equal to the combined effects of component force in XY plane and component torque around the Z axis in the center. With these proper transformations, the meshing effect can be modeled as a normal pressure force, represented by a time-varying stiffness spring damper force, and a tangential force, represented by a synthesis of a viscous damping force and a coulomb friction force (Figure 1b). These two kinds of force can eventually be replaced by time-varying component forces and torques with respect to the generalized coordinates on gear centers of the input and output gears respectively. Hence, large angular displacements of gears are allowed, and the entire coupling model can be utilized for the calculation of system dynamic behavior in the time domain.
The proposed method operates as space-fixed spring dampers with floating acting points, which guarantees that the positions and orientations of spring dampers coincide with the actual meshing forces regardless of changes in the centroids or the angles of the meshing gears. The proposed method is compared with the existing method in Table 1. This paper proposes a novel method for LTC modeling. The elasticity and viscidity of gears, carriers, shafts, bearings, etc., are represented by time-varying stiffness-damping forces, while the bodies of gears, carriers and shafts are all modeled as rigid elements. Thus, a final multibody system of the gear set with consideration of time-varying stiffness, friction, damping, etc. is established. Based on Newton's third law, the interaction of two bodies can be represented by an action force and an opposite reaction force that are functions related to the corresponding characteristic parameters of the two bodies. According to the basic law on parallel translation of force, forces on the contact point are equal to the combined effects of component force in XY plane and component torque around the Z axis in the center. With these proper transformations, the meshing effect can be modeled as a normal pressure force, represented by a time-varying stiffness spring damper force, and a tangential force, represented by a synthesis of a viscous damping force and a coulomb friction force (Figure 1b). These two kinds of force can eventually be replaced by time-varying component forces and torques with respect to the generalized coordinates on gear centers of the input and output gears respectively. Hence, large angular displacements of gears are allowed, and the entire coupling model can be utilized for the calculation of system dynamic behavior in the time domain.
The proposed method operates as space-fixed spring dampers with floating acting points, which guarantees that the positions and orientations of spring dampers coincide with the actual meshing forces regardless of changes in the centroids or the angles of the meshing gears. The proposed method is compared with the existing method in Table 1.
in which R is the ring gear, P is the planet gear, S is the sun gear, and C is the planet carrier.
A schematic diagram of a force analysis of the whole gear set is shown in Figure 2. The reaction forces on the ith planet gear include N i RP and f i RP , the normal pressure force and tangential force that the ring gear exerts; N i SP and f i SP , the normal pressure force and tangential force that the sun gear exerts; and F i CP and T i CP , the reaction force and reaction torque from the bearing. Force analyses of the ring gear and the sun gear can be presented in the same way, which consist of the meshing forces and the reaction bearing forces if they exist.

Detailed Formulations
Due to the unique status of the planetary gear set in heavy-duty gearboxes, detailed formulations of LTC dynamic equations of a planetary gear set with floating components are emphasized. Each object in the planetary gear set is modeled with three DOFs. The translational displacements and rotational angle are chosen as the generalized coordinates: ( , , ) , , in which R is the ring gear, P is the planet gear, S is the sun gear, and C is the planet carrier.
A schematic diagram of a force analysis of the whole gear set is shown in Figure 2. The reaction forces on the ith planet gear include

Force Derivation of the Engagement of the Sun and ith Planet Gears
The composite deformation in normal and tangent directions of a meshing can be obtained from the generalized displacements; thus, the normal force N and tangent force f are derived. As shown in Figure 3, two main influencing factors exist: the difference between the actual and the ideal angle, and the time-varying center distance between the sun and planet gear that is mainly introduced by the floating components.
Under the ideal condition where the deformations of bearings and gear teeth are not considered, the movement of the ith planet and the sun gear can be regarded as fixed-axis rotation. The angles fit a transmission ratio equation as follows:

Force Derivation of the Engagement of the Sun and ith Planet Gears
The composite deformation in normal and tangent directions of a meshing can be obtained from the generalized displacements; thus, the normal force N and tangent force f are derived. As shown in Figure 3, two main influencing factors exist: the difference between the actual and the ideal angle, and the time-varying center distance between the sun and planet gear that is mainly introduced by the floating components. , and Z denotes the tooth number.
As a result, the angle difference The linear deformation caused by this angle difference is where ( ) is the pitch radius of the ith planet gear, which is timevarying during the meshing process and can be derived from the center distance , and i PS α in the equation is the engagement angle: As shown in Figure 3b, , which denotes the center distance variation of i PS d , is another crucial factor in calculating the composite deformation. The linear deformation caused by the variation of center distance is The normal composite deformation of the equivalent spring damper between the planet and sun gears is therefore given by The tangent composite deformation that is mainly caused by variation in center distance can be represented as Furthermore, during the gear meshing process, the stiffness of an equivalent spring damper varies with changes in the gear meshing pair number, which toggles between 1 and 2 for a spur gear and is larger for a helical gear. According to the parameter identification of gear meshing stiffness, the relationship between the single pair and double pair Under the ideal condition where the deformations of bearings and gear teeth are not considered, the movement of the ith planet and the sun gear can be regarded as fixed-axis rotation. The angles fit a transmission ratio equation as follows: where γ i PS = arctan y i P − y S / x i P − x S is the time-varying position angle between the center connection line O i P O S and X axis, However, when considering the flexibility of the gear body, the deflection of gear teeth and the deformation of gear contact, the preceding equation does not fit anymore. The ideal angle θ i P−PS of the planet gear under the condition of a rotation angle θ S and a position angle ∆γ i PS is As a result, the angle difference ∆θ i P−PS with respect to the ith planet gear is The linear deformation caused by this angle difference is where is the pitch radius of the ith planet gear, which is time-varying during the meshing process and can be derived from the center distance , and α i PS in the equation is the engagement angle: , which denotes the center distance variation of d i PS , is another crucial factor in calculating the composite deformation. The linear deformation caused by the variation of center distance is The normal composite deformation of the equivalent spring damper between the planet and sun gears is therefore given by Appl. Sci. 2021, 11, 6876 6 of 19 The tangent composite deformation that is mainly caused by variation in center distance can be represented as Furthermore, during the gear meshing process, the stiffness of an equivalent spring damper varies with changes in the gear meshing pair number, which toggles between 1 and 2 for a spur gear and is larger for a helical gear. According to the parameter identification of gear meshing stiffness, the relationship between the single pair and double pair gear meshing stiffness can be represented by a coefficient ξ (generally 1 < ξ < 2). Starting from the moment that a single meshing begins, the time-varying meshing stiffness in one cycle is where θ is the meshing phase angle, and θ T is the cycle phase angle.
In addition, exchanges between the single and double pair of different planet-sun gear pairs are asynchronous during the process of loading sharing. Thus, the initial phase angles ϕ i representing the meshing state are introduced. Supposing that the total number of planets is n, and under the premise that design, manufacturing and assembly errors are ignored, the initial phase angle of an ith planet gear that is different from the adjacent one is Thus, the time-varying meshing stiffness of the ith planet-sun gear meshing pair is Based on the above derivations of the time-varying stiffness and composite deformations in the normal and tangent directions, the normal force exerted on the ith planet gear by the sun gear in Figure 2 is calculated as The direction of this force is tangent to the base circle of the planet gear. The angle between it and the X axis is The tangent friction force between the planet and sun gear is composed of the viscous damping and coulomb friction forces: In this equation, µ PS is the coulomb friction coefficient, c PS is the viscous damping coefficient, and sgn (x) is the sign function.
The tangent friction force is perpendicular to the normal force, direction of which is opposite to the direction of the relative sliding velocity. The angle between it and the X axis is

Force Derivation of the Engagement of the Ring and ith Planet Gear
Analysis procedures to derive the force between the ring and ith planet gears are similar to the analysis procedures for that between the sun and ith planet gears. The difference is mainly concentrated in a small fluctuation in meshing stiffness between the ring and planet gear due to the high contact ratio, such that meshing stiffness can be regarded as a constant in some situations, which will be proved via parameter identification results in a later section.
The center distance of the ith planet and ring gear is The time-varying position angle between the center connection line and the X axis is Under the ideal condition, the rotational angles fit a transmission ratio equation as follows: The ideal angle θ i P−PR of a planet gear under the condition of rotation angle of θ R and As a result, the angle difference ∆θ i P−PR with respect to the ith planet gear is The linear deformation caused by the angle difference is where α i PR is the engagement angle and r i P−PR is the time-varying pitch radius of the ith planet gear, which can be derived from the center distance d i PR , r i P−PR = d i PR · Z P /(Z R − Z P ). The linear deformation caused by the variation in center distance is Thus, the equivalent normal composite deformation of the spring damper between the ring and ith planet gears is As in the previous section, the tangent deformation mainly caused by variation in center distance between the ring and planet gear is k PR denotes the equivalent meshing stiffness of the ring and planet gears, and the fluctuation its value can be ignored. Based on the previous derivation of the deformation, the normal force on the ith planet gear can be computed as Its direction is tangent to the base circle of the ith planet gear, and the angle between it and the X axis is The tangent friction force f i RP between the ith planet gear and ring gear is composed of the viscous damping and coulomb friction forces: in which the definitions of µ PR , c PR and the function sgn(x) are almost the same as in the previous subsection. The tangent friction force is perpendicular to the normal force, direction of which is opposite to the direction of relative sliding velocity. The angle between it and the X axis is

Bearing Force Analysis between Carrier and ith Planet Gear
The interaction between the carrier and the planet gear can be represented using an equivalent spring damper force and a friction torque. R P denotes the distribution radius of the planets, and γ i CP0 is the initial distribution angle, namely the angle between the center line of the ith planet and the carrier along the X axis. As such, the deformation of the ith planet relative to the carrier is The rotational angle of the ith planet relative to the carrier is δ i CP−θ = θ i P − θ C . Therefore, the reaction force and torque on the center of the ith planet are computed as follows: where k CPB is the bearing stiffness, d CPB is the nominal diameter of bearings, µ CPB is the coulomb friction coefficient, and c CPB is the viscous damping coefficient.

System Dynamic Equations
To construct the sun gear dynamic equations, supposing that k SB and c SB are the bearing stiffness and the rotational viscous damping coefficient of the sun gear bearing if it exists, respectively, and that there exist n planet gears: .. For the ring gear, the dynamic equations derivation process is similar. k RB and c RB are assumed to be the ring bearing stiffness and the rotational viscous damping coefficient if they exist, respectively.
Plugging in the corresponding terms, the equations can be rewritten as M R ..
For the ith planet gear, as the applied forces on it include the engaging forces from the ring and sun gear and the bearing force, its dynamic equations are given by M P ..
For the planet carrier, similarly, dynamic equations can be expressed as

Dynamic Modeling of a WTG Test Rig
A full test rig which was used to test and evaluate the dynamic performance of a 2 MW WTG is introduced in this section. A LTC model considering specific structures of this WTG was constructed using the proposed modeling approach of a planetary gear set with floating components.

Structure of the Full Test Rig
The 2 MW WTG modeled is depicted in Figure 4. It primarily consists of three gear sets: a low-speed planetary gear set (LSP), a high-speed planetary gear set (HSP) and a parallel-shaft helical gear set (PSH). To reach the design objective of power dividing, the LSP had four transverse fixed planet gears and a floating sun gear, while the HSP had a fixed ring gear and a floating carrier with three planet gears affixed with bearings. The wind field excitation acted on the ring gear of the LSP, and the grid load from the electric generator acted on the output axis of the PSH. The intricate dynamics of the entire floating component, including the sun gear from the LSP, the carrier and three planet gears from the HSP, made this WTG quite distinct from common designed ones.
A full test rig which was used to test and evaluate the dynamic performance of a 2 MW WTG is introduced in this section. A LTC model considering specific structures of this WTG was constructed using the proposed modeling approach of a planetary gear set with floating components.

Structure of the Full Test Rig
The 2 MW WTG modeled is depicted in Figure 4. It primarily consists of three gear sets: a low-speed planetary gear set (LSP), a high-speed planetary gear set (HSP) and a parallel-shaft helical gear set (PSH). To reach the design objective of power dividing, the LSP had four transverse fixed planet gears and a floating sun gear, while the HSP had a fixed ring gear and a floating carrier with three planet gears affixed with bearings. The wind field excitation acted on the ring gear of the LSP, and the grid load from the electric generator acted on the output axis of the PSH. The intricate dynamics of the entire floating component, including the sun gear from the LSP, the carrier and three planet gears from the HSP, made this WTG quite distinct from common designed ones.  Figure 5 shows the 3D geometry model of the entire test rig. Part 1 is a loading motor, the speed of which was controlled, and part 2 is an auxiliary reduction gearbox. A combination of the two was utilized to imitate the wind input. Part 5 is a loading generator, the torque of which was controlled. Part 4 is the tested 2 MW WTG, which is connected with Part 2 through a universal coupling and with Part 5 through a flexible flange.

Dynamic Model of the Full Test Rig
A LTC model of the full test rig was constructed based on the proposed method in order to analyze the dynamic loading factors and dynamic behavior.  Figure 5 shows the 3D geometry model of the entire test rig. Part 1 is a loading motor, the speed of which was controlled, and part 2 is an auxiliary reduction gearbox. A combination of the two was utilized to imitate the wind input. Part 5 is a loading generator, the torque of which was controlled. Part 4 is the tested 2 MW WTG, which is connected with Part 2 through a universal coupling and with Part 5 through a flexible flange.
A full test rig which was used to test and evaluate the dynamic performance of a 2 MW WTG is introduced in this section. A LTC model considering specific structures of this WTG was constructed using the proposed modeling approach of a planetary gear set with floating components.

Structure of the Full Test Rig
The 2 MW WTG modeled is depicted in Figure 4. It primarily consists of three gear sets: a low-speed planetary gear set (LSP), a high-speed planetary gear set (HSP) and a parallel-shaft helical gear set (PSH). To reach the design objective of power dividing, the LSP had four transverse fixed planet gears and a floating sun gear, while the HSP had a fixed ring gear and a floating carrier with three planet gears affixed with bearings. The wind field excitation acted on the ring gear of the LSP, and the grid load from the electric generator acted on the output axis of the PSH. The intricate dynamics of the entire floating component, including the sun gear from the LSP, the carrier and three planet gears from the HSP, made this WTG quite distinct from common designed ones.  Figure 5 shows the 3D geometry model of the entire test rig. Part 1 is a loading motor, the speed of which was controlled, and part 2 is an auxiliary reduction gearbox. A combination of the two was utilized to imitate the wind input. Part 5 is a loading generator, the torque of which was controlled. Part 4 is the tested 2 MW WTG, which is connected with Part 2 through a universal coupling and with Part 5 through a flexible flange.

Dynamic Model of the Full Test Rig
A LTC model of the full test rig was constructed based on the proposed method in order to analyze the dynamic loading factors and dynamic behavior.

Dynamic Model of the Full Test Rig
A LTC model of the full test rig was constructed based on the proposed method in order to analyze the dynamic loading factors and dynamic behavior. Figure 6 shows the schematic of the simplified LTC dynamic model of the full test rig. Not only the meshing effect but also the bearing flexibility and torsional elasticity of shafts and joints were considered. In the model, the gear meshing effects of the LSG and HSG gear sets were considered using the proposed method. The PSH was regarded as a spur gear set with high contact ratio under the condition that axial flexibility was not considered. The imitated wind input that combined a loading motor and a reduction gearbox was considered an excitation with rotational inertia, and the output generator was considered a torque loading with rotational inertia. In addition, torsional models of the spline connection, the universal coupling and the flange joint were built, and the radial stiffness of different bearings was included.
HSG gear sets were considered using the proposed method. The PSH was regarded as a spur gear set with high contact ratio under the condition that axial flexibility was not considered. The imitated wind input that combined a loading motor and a reduction gearbox was considered an excitation with rotational inertia, and the output generator was considered a torque loading with rotational inertia. In addition, torsional models of the spline connection, the universal coupling and the flange joint were built, and the radial stiffness of different bearings was included. This proposed LTC model in Figure 6 was implemented with the aid of the MS. AD-AMS software due to its superiority in calculation of multibody systems. In the simulated model, translational displacements and rotational angles of all components were acquired in real time. The time-varying deformations and relative velocities of the elastic elements in Figure 6 were calculated using the derived equations. The meshing forces of gears, torques of shafts, bearing forces, etc. were then also computed using the equations in Section 2, and finally acted on the corresponding components using the force function of ADAMS. It should also be noted that parameters like meshing stiffness were also time-varying according to the related switching conditions derived in Section 2. Tables 2 and 3 show the main design parameters of the LSP, HSP and PSH gear sets, based on which basic parameters of gears and gear pairs could be derived and calculated.  As shown in Figure 7, parameter identifications of the gear meshings, the connections and the shafts were conducted by means of Finite Element Analysis (FEA). In the process, the loading was increased gradually until reaching the service load, and the stiffness value This proposed LTC model in Figure 6 was implemented with the aid of the MS. ADAMS software due to its superiority in calculation of multibody systems. In the simulated model, translational displacements and rotational angles of all components were acquired in real time. The time-varying deformations and relative velocities of the elastic elements in Figure 6 were calculated using the derived equations. The meshing forces of gears, torques of shafts, bearing forces, etc. were then also computed using the equations in Section 2, and finally acted on the corresponding components using the force function of ADAMS. It should also be noted that parameters like meshing stiffness were also time-varying according to the related switching conditions derived in Section 2. Tables 2 and 3 show the main design parameters of the LSP, HSP and PSH gear sets, based on which basic parameters of gears and gear pairs could be derived and calculated.  As shown in Figure 7, parameter identifications of the gear meshings, the connections and the shafts were conducted by means of Finite Element Analysis (FEA). In the process, the loading was increased gradually until reaching the service load, and the stiffness value was calculated as the slope of the nonlinear force-displacement curve on the operating point. For gear meshing, both the single and double pair contact situations were calculated. Detailed identification results are shown in Tables 3-5.

Parameter Identification of the Dynamic Model
Parameter identification of the radial bearing stiffness was conducted in a similar way. The inner or outer ring was fixed while the other one was loaded until reaching the service load, and the bearing stiffness was calculated through the nonlinear relationship between loading force and relative radial displacement. Table 6 shows the identification results. With these identified parameters, dynamic behavior of the 2 MW WTG could be simulated.   Parameter identification of the radial bearing stiffness was conducted in a similar way. The inner or outer ring was fixed while the other one was loaded until reaching the service load, and the bearing stiffness was calculated through the nonlinear relationship between loading force and relative radial displacement. Table 6 shows the identification results. With these identified parameters, dynamic behavior of the 2 MW WTG could be simulated.

Experimental Validation and Dynamic Analysis of the WTG
To demonstrate the validity of the novel proposed modeling method of heavyduty gearboxes, experiments were conducted on the 2 MW WTG test rig specified in Figures 4 and 5, and the experimental results were compared with the simulation results from the LTC model of the entire test rig under loading conditions identical to the rated working environment. Furthermore, dynamic loading factors of different gear sets, which are crucial to the design of WTG system, were calculated based on the LTC model that was verified by experiments.

Configuration of the Full Test Rig
The experiments were conducted in a workshop of SANY (a wind power equipment company). The test rig was loaded according to the rated working condition of the 2 MW WTG: the rated rotating speed of the imitated wind input that combined the loading motor and the reduction gearbox was 15 rpm and the rated loading torque from the generator, acting on the PSH output axis, was 1.06 × 10 4 N·m. Single-axis vibration acceleration sensors (PCB 352C34) were utilized to measure the vibration behavior of the WTG system, and the detailed position arrangement of sensors is shown in Figure 8. The data acquisition bandwidth was set to be 8192 Hz. Sensor signals were imported into a LMS SCADAS III data acquisition system, and the accompanying software LMS Test was employed for data processing and analysis. The simulation LTC dynamic model was loaded on the input equivalent moment of inertia progressively to a constant rotating speed of 15 rpm and on the output equivalent moment of inertia with an identical torque, to make sure that the loading conditions of the simulation and experiment WTG were identical. The variable stiffness coefficients of the three gear sets in the simulation model were set to be 1.1, 1.4 and 1.05, respectively, based on the identification results.

Comparison of the Simulated and Experimental Results
Based on the identified parameters and loading conditions, theoretical meshing frequencies of the tested WTG were computed through multiplying the rotating frequency by the tooth number. The results are listed in Table 7, and were 29.5 Hz, 158.2 Hz and 628.3 Hz, respectively.

Comparison of the Simulated and Experimental Results
Based on the identified parameters and loading conditions, theoretical meshing frequencies of the tested WTG were computed through multiplying the rotating frequency by the tooth number. The results are listed in Table 7, and were 29.5 Hz, 158.2 Hz and 628.3 Hz, respectively. Acceleration signals from the X and Y directions were synthesized to reveal the lateral vibration behavior of the full test rig. Figures 9 and 10 are the comparison diagrams for the front chassis and the output shaft respectively, and Table 8 shows the comparison of the main frequency components of the experimental and simulated results from different locations on the entire WTG system.        Due to the time-varying stiffness, the gear meshing effect can be regarded as inner excitations to some extent, and thus the vibration characteristic of the entire system became complex, which can be identified from the abundant frequency components of both vibration spectra. Experimental frequencies from the front chassis mainly concentrated on 62.          On this basis, the dynamic loading factors between gears, which are significant for the design stage, were calculated to be 1.12, 1.37 and 1.27 and are listed in Table 9. The results were different from the chosen factors of 1.02, 1.08, 1.07 obtained from an experiential manual. Under larger dynamic loadings that exceed the design limits, gear cracks may occur. Further dynamic examinations of the whole system in both the frequency and On this basis, the dynamic loading factors between gears, which are significant for the design stage, were calculated to be 1.12, 1.37 and 1.27 and are listed in Table 9. The results were different from the chosen factors of 1.02, 1.08, 1.07 obtained from an experiential manual. Under larger dynamic loadings that exceed the design limits, gear cracks may occur. Further dynamic examinations of the whole system in both the frequency and time domains can be conducted using the proposed method, while these examinations are beyond the scope of the existing method.

Conclusions
In this paper, a novel and generally applicable LTC modeling approach for calculating the dynamic behavior of heavy-duty gearboxes in the time domain with consideration of the floating effect, time-varying meshing stiffness, bearing stiffness, torsional stiffness, etc., is proposed. This approach is more comprehensive and apparently gives clearer insights into the time-domain dynamics of heavy-duty gearboxes, such as WTGs, than the existing method. Based on the proposed method, a LTC model of a full test rig was built, which was used to test a 2 MW WTG consisting of two planetary gear sets and a parallel-shaft helical gear set. Dynamic experiments were conducted on the test rig. Validity of the proposed method and the built model were demonstrated, as the simulated results were in accordance with the experimental results. According to the time-domain results of meshing forces in the verified model, larger dynamic loading factors than the designed ones chosen by a traditional empirical method were found. This means that gear damage may occur in actual working conditions and the built 2 MW WTG system must be redesigned and improved. In conclusion, the proposed method can give a new perspective on heavy-duty gearboxes that is preferable to the existing method, and can provide sufficient time-domain information to contribute to the design, evaluation and redesign of large-scale mechanical systems as represented by WTG systems. Rotational angle of the ith planet relative to the carrier k CPB Translational stiffness of the bearings between carrier and planet gear d CPB Nominal diameter of bearings between carrier and planet gear µ CPB Coulomb friction coefficient of the bearings between carrier and planet gear c CPB Viscous damping coefficient of the bearings between carrier and planet gear k SB , c SB Bearing stiffness and rotational viscous damping coefficient of the sun gear bearing k RB , c RB Bearing stiffness and rotational viscous damping coefficient of the ring gear bearing M a , I a , G a (a=S,P,R,C) Mass, inertia and gravity of the sun gear, planet gear, ring gear and carrier respectively T a (a=S,R,C) Torques on the sun gear, ring gear and carrier respectively