Meshing Stiffness Parametric Vibration of Coaxial Contrarotating Encased Differential Gear Train

Planetary gears are widely used in mechanical transmission systems, but the vibration and noise affect their reliability and life. In this paper, the torsional dynamic model of an encased differential planetary gear with coaxial contrarotating outputs considering the time-varying meshing stiffness, damping, and phase difference of all gear pairs is established. By solving the equations of the derived system, three types of natural frequencies with different multiplicities of the system are obtained.)emultiscale method is used to study the parametric vibration stability caused by the time-varying meshing stiffness, and the results are verified by numerical simulation.)e dynamic characteristics of elastic meshing force are analyzed from time domain and frequency domain. )e variation of the dynamic load factor of each gear pair with input speed and the relationship between its peak position and the natural frequency of the derived system are discussed.)e results show that there is an unequal coupling phenomenon of meshing frequency between the meshing forces of different planetary sets. In the absence of external excitation, the meshing stiffness parameters not only excite the main resonance response of the system but also cause superharmonic resonance, subharmonic resonance, and combined resonance.


Introduction
Due to the existence of an encased structure, the coaxial contrarotating encased differential gear train can transmit torque through multiple paths, so it has stronger load-carrying capacity, more compact structure, and wider transmission ratio range than the ordinary multistage planetary gear train. It is widely used in the fields of helicopter transmission, marine power, and aeroengine. However, it is difficult to design the parameters of the system due to the requirement of constant speed reverse output of the internal and external output shafts, and the vibration and noise caused by internal excitations such as time-varying stiffness and errors also challenge its application. Howe and Mckibbin [1,2] studied various gear drive system configurations that can realize the coaxial counterrotating of the fan rotors, and, on this basis, the lubrication and efficiency of differential gear train configuration are analyzed.
Parametric stability of gear systems caused by time-varying meshing stiffness has been extensively studied. Lin and Parker [3,4] used multiscale method to study the effects of meshing stiffness parameters on the instability of the two-stage gear system and planetary gear system and derived a design expression for controlling the instability region by adjusting the stiffness parameters. Liu et al. [5] used the perturbation method to derive the expression of the system parameters instability with both meshing stiffness and bearing stiffness fluctuations, which is verified by Floquet theory. Qiu et al. [6,7] established a lumped parameter model of a planetary gear considering timevarying meshing stiffness, damping, and input speed fluctuations and studied the effects of damping and input speed fluctuations on instability by the perturbation analysis and numerical method. Canchi and Parker [8] studied the parameter excitation of the rotating ring gear under time-varying stiffness and obtained the conditions of the ring-planet mesh phasing and contact ratio, which can suppress the parameter instability of the flexible ring gear in the planetary gear system. Based on the hybrid continuous discrete model, Parker [9,10] obtained the parameter expression of the instability boundary of the planetary gear system with elastic continuous ring gear by using the structural modal characteristics and multiscale method.
Scholars have carried out many studies on the dynamic characteristics of meshing forces and multipath load distribution caused by time-varying meshing stiffness and error of planetary gear system. Inalpolat and Kahraman [11] established a nonlinear dynamic model including periodic time-varying mesh stiffnesses and gear tooth separation that can predict the sidebands of planetary gear set having manufacturing errors and rotating carrier and compared the predicted meshing force spectrum with the measured meshing force spectrum. Sondkar and Kahraman [12] used a three-dimensional dynamic model of the herringbone planetary gear train to analyze the variation of the maximum meshing force amplitude with the meshing frequency under different staggering phase angles of the left and right gears. Sun et al. [13] and Dong et al. [14] established the nonlinear dynamic model of planetary gear train and studied the influence of different load-sharing structures on the load-sharing performance.
At present, the research mainly focuses on single-stage or multistage planetary gear train, and the research on the dynamic coupling characteristics of the compound planetary gear train is not enough. In the coaxial contrarotating encased differential gear train, the inner and outer shafts are output in the opposite direction, and there are many output components that affect each other, so the dynamic coupling characteristics of the system are more complicated.
In this paper, a multi-degree-of-freedom torsional dynamics model considering the time-varying meshing stiffness and mesh phasing of the encased differential gear train, which is widely used in coaxial twin-rotor helicopters, is established. e multiscale method is used to analyze the vibration characteristics of the system caused by the fluctuation of meshing stiffness. e elastic meshing force is obtained by solving the differential equations using the Runge-Kutta method. e coupling characteristics of meshing frequency of different planetary set and the fluctuation characteristics of the dynamic load factor with the rotating speed are analyzed.

Dynamic Model of the Encased Differential Gear Train
Coaxial contrarotating encased differential gear train is shown in Figure 1. e encased stage consists of the sun gear s1, stepped planets ai and bi, and ring gear r1, and the differential stage consists of the sun gear s2, planets pj, and ring gear r2. e input power is transmitted to the ring gear and the carrier, respectively, through the encased stage and the differential stage to realize the reverse double output. By reasonable design of system parameters, the output speeds of ring gear r2 and carrier c2 can be equal [15]. Define the torsional deformation of each component of the gear system as generalized coordinates: x � u s1 , u s2 , u r1 , u r2 , u c2 , u a1 , . . . , u aM , u b1 , . . . , u bM , u p1 , . . . , u pN T , where M and N are the numbers of stepped planet gears and planet gears, respectively. For simplicity, the following assumptions are made: the planets of each planetary set are evenly arranged and have equal mass and moment of inertia. e torsional stiffness of each stepped planet is equal. e friction and error of the system are not considered. e torsional dynamic model of the coaxial contrarotating gear system is shown in Figure 2. In the figure, k s and k x represent the torsional stiffness of the sun gear shaft and the stepped planet gear shaft, respectively, and k r is the torsional stiffness between the rings r1 and r2. k s1ai , k r1bi , k s2pj , and k r2pj are the meshing stiffnesses between the sun s1 and the stepped planets ai, the ring r1 and the stepped planets bi, the sun s2 and the planets pj, and the ring r2 and the planets pj, respectively. Damping symbols are similar.

Stiffness.
e Fourier series expression of the timevarying meshing stiffness considering the mesh phasing is [4] In the formula, k a m and 2k v m are the mean component and fluctuation value of the mesh stiffness of the gear pair m, ω m is the meshing frequency, ε m is the contact ratio, and c m is the mesh phasing.

Damping.
e meshing damping of gear pair m is [16] c m � 2ξ m In the formula, ξ m is the meshing damping coefficient, I m1 and I m2 are the moments of inertia of the gears, and r m1 and r m2 are the base radii of the gears. e torsional damping of the shaft z is where ξ z is the relative torsional damping coefficient of the shaft z, k z is the torsional stiffness, and I z1 and I z2 are the moments of inertia of the components connected to the shaft.

The Equations of Motion for the Gear Train
According to Newton's law, the equations of motion for the encased differential gear train are obtained as follows: In the above equations, the relative displacement of gear mesh is where I h (h � s1, s2, r1, r2, c2, m, n, p) are the moments of inertia of the component h, r bh is the base circle radius of the gear h, r c2 is the distance from the center of planet to the center of the carrier, m p is the mass of planet, α s2p ′ and α r2p ′ are the working pressure angles of the sun-planet and ringplanet mesh, T in is the input torque, and T ou1 and T ou2 are the load torques applied to the carrier and ring gear, respectively. Substituting equation (14) into (6)- (13), the matrix form of the system equations can be obtained as follows: where M, C,K, and Q are generalized mass matrix, damping matrix, stiffness matrix, and external load array, respectively, as shown in Appendix A. Disregarding damping and external forces, the equation for the free vibration of the system is In the equation, K a is the mean stiffness matrix, including torsional stiffness and comprehensive meshing stiffness, and K v is the variable stiffness matrix.

Dynamic Equations Decoupled and Stability Analysis
Let K v in equation (16) be zero; then the derived system equation is e eigensolutions can be obtained by solving formula (17).
where φ (i) are the vibration modes of the derived system and ω i are natural frequencies. Let where Φ � [φ (1) , φ (2) , . . . , φ (n) ] is the main modal matrix; q � q 1 , q 2 , . . . , q n T . Substituting equation (19) into equation (16) and multiplying Φ T by left, the regular modal equation can be obtained as Expansions of equation (20) are Using multiscale method, time variables representing different scales are introduced: e solution of equation (21) is expressed as Substituting equation (23) into equation (21), collecting terms of the same power in ε on both sides of the equation yields where A n (T 1 ) is the amplitude in complex form, ε � max ε s1a , ε r1b , ε s2p , ε r2p , and cc is the complex conjugate of preceding items. Substitution of equation (24) into equation (25) yields It can be seen from the above that, in addition to the fact that the meshing frequency is close to the natural frequency of the system to excite the main resonance, when the frequency lω m is close to the natural frequency combination ω n ± ω r , there will be a secular term in the equation and an unstable region of parameter vibration in the gear system.

Meshing Stiffness Parametric Vibration
e parameters of the coaxial contrarotating encased differential gear train are shown in Table 1.
e meshing frequency of each gear pair of the system is [15] f s1ai � n in z s1 60 , where n in is the input speed (r/min).      f/Hz     (a)   Table 2. It can be divided into three types according to its multiplicity: eight natural frequencies always have multiplicity 1, two groups of M−1 multiples (f 4 /f 5 , f 12 /f 13 ), and one group of N-1 multiples (f 8 /f 9 ), which, respectively, correspond to the torsional vibration mode, stepped planet mode, and planet mode of the system. Frequency of 0 represents the rigid body motion of the system, and the magnitude of multiple frequencies is independent of the number of planets [17].

Dynamic Characteristics of Elastic Meshing Force.
e elastic meshing force of gear pair m in the encased differential gear train is expressed as F m (t) � k m (t) δ m . It is defined that input torque T in � 1283.5 N·m is constant, and the (c)  Figure 5: Spectrum of meshing forces (n s1 � 3000 r/min). (a)  (c)   Mathematical Problems in Engineering output torque T ou1 � T ou2 � iT in /2, where i is the gear system transmission ratio. e Runge-Kutta method is used to solve the numerical solution of the differential equations.

Time-Domain Analysis of Meshing Force.
When the torque is constant and the input speed is ns1 � 4000 r/min, the meshing force of gear pairs is shown in Figure 3. It can be seen from the figure that the dynamic meshing force fluctuates around the static meshing force value shown by the dotted line. e meshing force F s1a1 fluctuates with its own meshing period. e meshing forces F r1b1 , F s2p1 , and F r2p1 mainly fluctuate in their own periods, and they also have high-frequency fluctuations of small amplitude. is is caused by the excitation of the time-varying stiffness parameters such as mesh phasing and contact ratio.

Frequency-Domain Analysis of Meshing Force.
e frequency spectrum of system meshing force at different input speeds and constant torque is shown in Figures 4-6. It can be seen from the figures that, in addition to the meshing frequency of its own planetary set, the meshing frequency of other planetary sets and natural frequencies f 2 and f 3 also appear in the meshing force, which indicates that the system has multifrequency coupling vibration.
As shown in Figure 4, when the input speed is n in � 4000 r/min, the meshing frequency f r1b1 appears in the frequency spectrum of the meshing force F s1a1 of the encased stage, and the encased stage meshing frequency f s1a1 appears in the spectrum of differential stage meshing forces F s2p1 and F r2p1 , but there is no differential stage meshing frequency in the encased stage meshing force spectrum, which indicates that the meshing frequency   coupling relationship between the meshing forces of different planetary sets is not equal. As shown in Figures 5  and 6, when the rotating speed is n in � 3000 r/min and n in � 2000 r/min, in addition to the above coupling effect, there is also coupling of differential stage meshing frequency in the encased stage meshing force spectrum. is is because the torsional internal force of the connecting shaft couples the vibration displacements of different frequencies together and then acts on the two connected parts in the form of acting force and reaction force, respectively. e connected coupling parts will return back the vibration displacement with new frequency to the meshing force of the planetary set in the form of equivalent meshing line deformation, resulting in the intercoupling between the meshing forces of different planetary sets [18].

Analysis of the Dynamic Load Characteristics.
In the compound planetary gear transmission, different meshing states will cause difference between the multiple external and internal meshing forces of the same planetary set and will also form the different dynamic characteristics of the meshing forces of each planetary set. In this paper, the dynamic load factor is used to evaluate the dynamic characteristics of meshing forces. e dynamic load factor G m is the ratio of the maximum dynamic load to the static load on the gear teeth during a meshing period; that is [19], where F m ′ is the static meshing force. It can be seen from Figure 7 that, in the range of the speed variation, the dynamic load factor of the gear pair s1ai is significantly larger than that of other meshing pairs, which is because the gear pair s1ai is a high-speed gear train and its vibration is relatively severe. e load distribution among the planets in each planetary set of the system is very uniform.
According to Figures 7(a) and 7(b), the system has a strong main resonance response when the meshing frequency f s1ai is close to the natural frequency f 2 , and the dynamic load factors G s1ai and G r1bi of the encased stage gear pair are very large. Due to the time-varying meshing stiffness, the meshing frequency f s1ai causes the superharmonic responses and 1/2 subharmonic response of the encased stage, and the meshing frequencies f r1bi and f s2pj also cause the superharmonic responses.
From Figures 7(c) and 7(d), in addition to the main resonance response, superharmonic response, and subharmonic response, the dynamic load factor of the differential stage also peaks at f s2pj � f 2 + f 3 , corresponding to the combined resonance of the additive type.

Conclusions
In this paper, the torsional dynamic characteristics of the coaxial contrarotating encased differential gear system are studied, and the natural frequency of the system is obtained. e time-domain and frequency-domain characteristics of the elastic meshing force are analyzed, and the variation of dynamic load factor with the input rotating speed is studied.
(1) Natural frequencies of the encased differential gear train can be divided into 3 types according to its multiplicity, which correspond to the torsional vibration mode, stepped planet mode, and planet mode. (2) e dynamic meshing force fluctuates around the static meshing force value. e meshing force of the gear pair s1a1 fluctuates with its own meshing period. e meshing forces of gear pairs s1b1, s2p1, and r2p1 fluctuate mainly in their own periods, and they also have high-frequency fluctuations of small amplitude. (3) In addition to the meshing frequency of its own planetary set, the meshing force spectrum also includes the meshing frequency of other planetary sets and the natural frequency of the system. (4) e dynamic load factor of the gear pair s1ai is significantly larger than those of the other meshing pairs. e dynamic load factor of the gear pair of the system peaks not only at the main resonance response but also at the superharmonic response, subharmonic response, and combined resonance response.

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 there are no conflicts of interest regarding the publication of this paper.