Analysis of Dynamic Characteristics of the Multistage Planetary Gear Transmission System with Friction Force

Multistage planetary gear transmission system has been widely utilized in engineering practice due to the salient characteristics, such as high bearing load and large speed ratio. (is paper addresses a two-stage planetary gearbox and establishes a system coupling torsional dynamical model which considers the time-varying mesh stiffness, friction forces, and interstage coupling factors. Meanwhile, the friction and lubrication states are classified to comprehensively analyze the calculation of friction coefficients under different conditions. Considering the time-varying influence of friction on the tooth surface under the condition of fluid lubrication, the vibration response under parametric excitation is solved by a numerical method. A multistage planetary transmission test bench is built in the back-to-back form so as to test the vibration of the two-stage planetary gearbox. It shows that the simulation results of the dynamical model are consistent with the test data. Consideration of the calculation of friction on the tooth surface and the friction coefficients is helpful for the establishment of the more accurate dynamical model and lays the foundation for the structural design, fault diagnosis, and dynamic optimization of the multistage planetary gear transmission system.


Introduction
e multistage planetary gear transmission system is mainly used in complex machines such as large wind power generator, mining machinery, and shield machines. It is popularized in the industry due to advantages including compact structure, large transmission ratio, high load, and smooth transmission. However, because of the numerous structures, complex dynamic coupling, and abundant nonlinear factors of multistage planetary gear devices, the accuracy of the kinetic model can hardly be guaranteed. Indepth research on dynamic characteristics is still challenging, which restricts the development and application of multistage planetary gears. Traditionally, the analysis of vibration and noise of multistage planetary gears is based on dynamics research. Vibration and noise not only affect the working environment but also reduce the normal operating life of equipment. erefore, it is of great importance to study the dynamic characteristics of the multistage planetary gear transmission system in both engineering application and theoretical research. Planetary gear dynamics has been an important and active research topic in the field of mechanical transmission in recent years. Kahraman [1] established the systematic dynamical model and studied the inherent characteristics of planetary gear transmission systems. Lin and Parker [2] constructed a pure torsional dynamical model to study the instability of planetary gear parameters caused by the variation of mesh stiffness. Sun and Hu [3] established a dynamical model of bendingtorsional coupling and analyzed the nonlinear frequency response characteristics of planetary gear transmission. Song et al. [4] presented the modified torsional dynamics model of the 2K-H planetary gear transmission device which considered the support stiffness of the bearing. Parker et al. [2,5,6] proposed a pure torsional dynamical model and analyzed the dynamic response of the meshing phase. Zhu et al. [7] investigated a translation-torsional dynamical model of the single-stage planetary gear to compare and analyze the influence of the tooth surface friction on dynamic characteristics. Bank and Parker [8] studied the analytic solution of nonlinear dynamics of planetary gear transmission. Li et al. [9] analyzed the periodic orbits of the planetary gear nonlinear transmission system by the improved PNF method. Kim et al. [10] suggested a bendingtorsional coupling dynamical model under the consideration of the time-varying coincidence degree and the time-varying meshing angle. Qin et al. [11] addressed the transmission system of the shield machine and established a dynamical model of the multistage planetary gear coupling to compare and analyze the simulation data and the experimental data. Liu et al. [12] analyzed the time-domain and frequencydomain response of the system and verified the effectiveness of the modular modeling method by modularizing the multistage planetary transmission system. Zhang et al. [13] presented a translational-rotational dynamic model of a two-stage closed-form planetary gear set under the consideration of the rotational and translational displacements to investigate the dynamic response and to avoid resonance. Liu et al. [14] addressed a new elastic-discrete model of the planetary gear, in which the sliding friction and elastic ring gears are considered simultaneously. Wei et al. [15] proposed a comprehensive, fully coupled, dynamic method by applying a virtual equivalent shaft element and used the proposed method to guide the design reliability and lowvibration planetary gear systems. Chen et al. [16] established a novel dynamic model of the planetary gear transmission by using Newton's theory, in which some key factors such as time-variant meshing stiffness, phase relationships, and tooth contact characteristics are considered. Li et al. [17] established a coupled dynamic model to predict the modulation sidebands of a two-stage planetary gear set and analyzed the influence of crack propagation on the mesh stiffness. Mo et al. [18,19] established the dynamic models for the multipower face gear split flow system and herringbone planetary gears, respectively, by the lumped-parameter method and studied the load-sharing characteristics of the transmission system. Ambarisha and Parker [20] proposed two independent models: one a lumped-parameter mathematical model and the other a finite element model, with different formulations, and different mesh modeling assumptions are used to analyze the nonlinear dynamics of planetary gears. e above research mainly focuses on the dynamics of single-stage planetary gears, while the dynamics of multistage planetary gears has rarely been investigated. Moreover, the excitation factors in these studies are simplified, and the calculation method of the friction on the tooth surface and friction coefficient is not considered in most papers. However, Howard et al. [21] explored the effect of friction on the resultant gear case vibration by using a simplified gear dynamic model and found that the diagnostic techniques worked clearly when friction was included in the model. Inspired by these observations, this paper will focus on the two-stage planetary gear transmission system and establish the system coupling torsional dynamical model by considering the calculation methods of the friction force on the tooth surface and friction coefficient, the time-varying meshing stiffness, and the interstage coupling. is model is more precise, which can be effectively used to analyze the vibration and predict the dynamic force of multistage planetary gears. It provides some basis for the dynamic design and fault diagnosis of the multistage planetary gear transmission system. Accurate analytical modeling is needed to estimate the relative gear vibration and predict dynamic forces in industrial applications.
is paper is organized as follows: first, a new coupling dynamical model is presented by using the lumpedparameter method on the basis of the Lagrange equation and geometric relation analysis. en, the vibration response under parametric excitation is solved by a numerical method, and the dynamic characteristics and parameter influence rules of the key components are obtained. Finally, the experimental results and simulation data are compared and analyzed by carrying out the bench test of the planetary gearbox, and the validity of the analytical model is verified.

Coupling Dynamical Model of the Two-Stage Planetary Gear Transmission System
In this paper, two planetary gear drive systems are considered, which are used to establish the dynamical model of the two-stage model, and there are a total of three planet wheels at each stage. e coupling between the two gears is considered as the torsional spring. e transmission coupling between stage I and stage II is expressed by the torsional springs and dampers. e internal transmission structure and the schematic diagram of the studied system are as shown in Figures 1(a) and 1(b), respectively.

Dynamical Model.
In the torsional dynamical mode as shown in Figure 2, n � 1, 2, . . ., N is the number of planet gears. For the dynamical model of the single-stage drive system with N planet gears, the torsional dynamics model includes N + 3 degree of freedom when only the torsional degree of freedom is considered. Let j � c, r, s, n, where k jt is the torsional support stiffness of components, u j is the torsional line displacement of components, k sn (t), c sn , and e sn (t) are the time-varying meshing stiffness, damping, and transmission error of the sun gear and nth planet, k rn (t), c rn , and e rn (t) are the time-varying meshing stiffness, damping, and transmission error of the ring gear and nth planet gear, φ 2 is the planet gear phase angle, and k e 12 and c e 12 are the interstage coupling stiffness and coupling damping.
Let gear meshing compression be positive, and the torsional model of the two-stage planetary gear train includes 2N + 6 degrees of freedom. e generalized column coordinate x is Equation (2) shows the relative displacement of the sun gear, ring gear, and planetary gear: 2 Shock and Vibration

Shock and Vibration
where z s and z r are the working pressure angle of the sun and ring gear. e coupled torsional dynamic equation of the two-stage planetary gear transmission system with frictional force is formulated as the following equation: where superscripts I and II denote the 1st and 2cd stage of the two-stage gear system; J j are the moments of inertia of each component; r j are the equivalent radii of the rotating parts; r e is the radii of the shaft coupling; m p is the mass of a planet gear; F sn , F rn are the nth planet-ring and planet-sun meshing forces; f sn , f rn are the friction forces; L sn , L rn are the arm of friction forces; i 1 , i 2 are the transmission ratios of the 1st and 2cd stage; and T in , T out are the input torque and load torque, respectively. e dynamic equations of the coupled torsional system are arranged as the following equation: where x is the generalized coordinate vector; M is the mass matrix; C is the damping coefficient matrix; K jt is the torsional support stiffness matrix; K (t) is the time-vary meshing stiffness matrix; and F 1 and F 2 are the external and internal loads, respectively.

Transmission Error and Damping
Coefficient. e meshing error is simplified according to the sine curve as the following equation: where E is the amplitude of the meshing error and ϕ is the initial phase angle.
According to the definition of critical damping and vibration theory, the meshing damping coefficient is given as the following equation: where ζ v is the meshing damping ratio in the range of 0.02∼0.17; I 1 and I 2 are the moments of inertia of the two meshing gears; k 0 is the mean meshing stiffness; and r b1 and r b2 are the base radius of gear wheels.

Time-Varying Meshing
Stiffness. e meshing stiffness can usually be obtained by the energy method or finite element method. e literature and previous work show that the time-varying meshing stiffness usually presents the periodic characteristics of trapezoidal waves [20]. So, here, the time-varying mesh stiffness is expressed by the Fourier series as the following equation: where k 1 � (2(k m1 − k m2 )/nπ)sin[(2 − ε)nπ], k 0 � k m1 (ε − 1) + k m2 (2 − ε), ϕ � (2 − ε)nπ, and k m1 , k m2 are the meshing stiffness of the double teeth and single tooth; ε is the contact ratio; k 0 is the mean meshing stiffness; k 1 is the harmonic component of the meshing stiffness; ω m is the meshing frequency; and ϕ is the harmonic phase of the meshing stiffness. 4 Shock and Vibration Figure 3 shows an external meshing process diagram (the internal engagement process diagram is not given in this paper). N 1 N 2 is the theoretical length of the meshing line, B 1 B 2 is the actual length of the meshing line, r s and r n are the base radius of the sun wheel and planet wheel, r as and r an are the addendum circle radius of the sun wheel and planet wheel, and zs is the meshing angle of the sun wheel and planet wheel, where i � 1, 2 is meshing tooth pair 1 and meshing tooth pair 2.

Frictional Force Arm of the Tooth Surface.
Ignoring the meshing-in and meshing-out deformation of the gear tooth contact, according to Figure 3, the frictional force arm between the sun wheel and the planet wheel is expressed as the following equation: where L sn,i is the friction force arm of the friction force arm of the meshing point relative to the sun gear in the ith sunplanet gear pair; L ns,i is the friction force arm of the meshing point relative to the planet gear in the the ith sun-planet gear pair; and p sn is the normal pitch between the sun gear and the planetary gear. Similarly, the frictional force arm between the inner gear and the planetary gear is expressed as the following equation: where L nr,i the friction force arm of the meshing point relative to the planet gear in the the ith ring-planet gear pair; L rn,i is the friction force arm of the meshing point relative to the ring gear in the the ith ring-planet gear pair; and p rn is the normal pitch between the inner gear and the planetary gear.

Friction on the Tooth Surface.
Friction on the tooth surface can be expressed as the following equation: where δ is the coefficient of friction direction, varying with the relative speed of the driving gear and driven gear, valued as 1, 0, and − 1, μ is the friction coefficient, and F is the normal load on the tooth surface. According to tribology and lubrication theory, the Stribeck curve is used to describe the lubrication state of the friction pair surface. e Stribeck curve uses the friction coefficient as the ordinate and the expression of viscosity, speed, and load as the abscissa to illustrate several forms of lubrication. Equation (11) expresses the lubrication state of the gear is very complicated and depends on the film thickness ratio: where h min is the minimum oil film thickness and R p and R g are the roughness of gear surface 1 and gear surface 2.
According to the ratio (λ) of the thickness of the oil film to the overall roughness of the surface of the friction pair, the lubrication states are usually divided into three types: boundary lubrication, mixed lubrication, and hydrodynamic lubrication. Figure 4 shows the gear meshing under different lubrication states. When λ <1, the gear teeth are in the boundary lubrication state, and the Buckingham formula can be used. When 1 ≤ λ ≤ 3, the gear teeth are in a mixed lubrication state, that is, the gear teeth are in the partial elastic hydrodynamic lubrication state, and the Kelley and Lemanski formulas can be used in the calculation. When λ >3, the gear teeth are in the full elastohydrodynamic lubrication state, and the full elastohydrodynamic lubrication theory is used to calculate the friction coefficients.
Buckingham semiempirical formula [7] is given as the following equation: where ] is the relative sliding velocity.
Kelley-Lemanski formula is given as the following equation: where R is the surface roughness; ω is the load on the unit width; V s is the sliding speed; V 1 denotes the total tooth speed; V r defines the rolling speed; η 0 is the absolute viscosity; and r p , r g are the base radii. e formula for full elastohydrodynamic lubrication [7] is described as the following equation: where are the tangential velocity of the main wheel and the driven wheel; α is the viscosity coefficient; h 0 is the oil film thickness at the rupture of the dynamic oil film; η is the dynamic viscosity of oil; P is the normal load; and ϕ, β, and k are the coefficients of friction coefficient calculation.

Dynamic Response Analysis
According to the aforementioned structure, the input component is the sun gear, and the output unit is the planet carrier in the planetary gearbox. e planetary gears are evenly spaced, and the carrying capacity is the same. e parameters of each component of the transmission system Shock and Vibration 5 are as shown in Table 1. e meshing stiffness is calculated by the elastic theory and Ishikawa method. According to the experimental equipment configuration, the load torque is as follows: no load, 504 Nm, and 1016 Nm, and the three operating conditions at the rotating speed of 500 r/min are tested as shown in Table 2.
In the simulation, the load is 504 Nm, and the input speed is 500 r/min. e friction coefficients are calculated according to the formula under the boundary lubrication condition. Under this operating condition, the factors such as the tooth surface friction, the time-varying meshing stiffness, and the transmission error are taken into account, and the dynamic response of the transmission system is obtained by using the Runge-Kutta method to solve the differential equations. Figure 5 shows the time-domain response of vibration displacement of planetary gears at different levels, and the peak-to-peak values of vibration displacement of first and second stages are 19.593 μm and 12.315 μm, respectively. Figure 6 shows the time-domain profile of vibration displacement of sun gears at different levels, and the peakto-peak values of vibration displacement of first and second stages are 31.019 μm and 14.316 μm, respectively. Figure 7 shows the vibration displacement-time domain performance of the ring gear at different levels, and the peakto-peak values of vibration displacement of the first and second stages are, respectively, 13.949 μm and 9.626 μm.
According to the peak value of each component and the effective RMS value, the vibration of the sun gear is the most obvious one with the same vibration, and the vibration of the ring gear is small under the vibration at the same level. However, the vibration of the first stage is greater than that of the second stage in the vibration of different levels. e frequency spectrum is shown in Figure 8. It can be seen that the peak frequency of the vibration of the ring gear occurs mainly in the first stage of gear frequency (167.3 Hz), 3f 1 frequency doubling, and the second stage of gear frequency f 2 (39.3 Hz) in the range of 1000 Hz. e peak frequency of the vibration of the second ring gear occurs mainly in the first stage of gear frequency, 4f 1 frequency doubling, and the second stage of gear frequency f 2 . e result shows that the vibration of the first-stage planetary gears is much larger than the second-stage vibration. e main reason that the second-stage frequency doubling is not obvious is the low meshing frequency of the second stage, and the vibration of the structure is not easily excited. e vibration caused by the first-stage planetary gears' meshing frequency has a larger vibration contribution on the first-stage and second-stage gear rings. e vibration caused by the second-stage meshing frequency has a larger contribution in the second stage, but the contribution in the first stage is relatively small.

Vibration Test of the Two-Stage
Planetary Gearbox

Experimental Setup.
Considering the large transmission of multistage planetary gears, a back-to-back planetary gear system is adopted for better control of operating conditions. e test setup used in the experiments is shown in Figure 9, including a driving motor, torque speed sensors, a test gearbox, a slave gearbox, and a load motor. e test setup and measurement points of the planetary gearbox test stand are shown in Figure 10. e measurement points are arranged at the position with obvious vibration. No. 1～2 are vertical and horizontal directions of the box

Test Data Analysis.
e sampling frequency is 10,240 Hz, and the resolution is 1 Hz. e vibration acceleration data of each measuring point are collected by the        acceleration sensors, and the velocity and displacement response curves are calculated by numerical integration. e vibration acceleration curve of the ring gear (measuring points 3 and 5) at each stage is provided in Figure 11, and the vibration speed curve is illustrated in Figure 12. e vibration displacement is also provided in Figure 13. e vibration peak value and the effective RMS value are summarized in Tables 3 and 4. e test acceleration spectrum is given in Figure 14. In the range of 1000 Hz, the vibration peak of the first-stage ring gear mainly shows the first-stage frequency doubling 3f 1 (167.3 Hz) and the second-stage frequency doubling 3f 2 (39.3 Hz). e vibration peak of the second-stage ring gear mainly shows the first-stage frequency doubling 4f 1 . However, the gear meshing frequency and the frequency doubling of the second stage of gear teeth f 2 and the frequency doubling of the second-stage gear are not obvious. As shown in the figures, under this test condition, the system vibration frequency is mainly concentrated at 0∼3500 Hz, and the sideband frequency characteristics are not obvious. is shows the high manufacturing accuracy of the gearbox.

Conclusion
(1) Aiming at the complex dynamic problems of multistage planetary gears, the influence of dynamic parameters is further analyzed, and a detailed dynamic model is established. Considering the factors such as the friction on the tooth surface, the calculation of friction coefficients, the time-varying meshing stiffness, and the interstage coupling, the coupled torsional dynamical model of a multistage planetary gear transmission system is established by the lumped-parameter method. (2) According to the peak value of each component and the effective value RMS, the vibration of the sun gear is the most obvious in the same vibration level, and the vibration of the ring gear is small in the vibration at the same level. However, the vibration of the first stage is greater than that of the second stage in the vibration of different levels. Moreover, the frequency peaks in the frequency domain are mainly expressed as tooth frequency and its frequency doubling of each stage. e vibration of the first stage is greater than that of the second stage, and the vibration characteristics of the second-stage tooth frequency are included. (3) e vibration data of the two-stage planetary gearbox were collected. It is found that the vibration displacement peak value and the effective value RMS of the test data are close to the simulation value. e frequency-domain components are concise with the simulation data. is shows that the calculation of the friction force and friction coefficient contributes to the refinement of the dynamical model. e research results lay the foundation for the accurate dynamic modeling and structural optimization of the multistage planetary gear transmission system.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.