Analysis of Nonlinear Vibration Characteristics of the Concentric Face-Gear Split-Torque Transmission System

+e concentric face-gear split-torque transmission system (CFGSTTS) is a new type of transmission that has significant applications in helicopter main gearboxes. To study the influence of various parameters on the dynamic characteristics of the CFGSTTS, a 23-degree-of-freedom translation-torsion nonlinear dynamic model was established based on the lumped parameter theory. +e model includes tooth backlash, error excitation, time-varying meshing stiffness with meshing phase difference, meshing damping, and elastic support deformation. +e excitation conditions for the time-varying meshing stiffness of face-gear pairs were calculated based on the strain energy theory. +e bifurcation characteristics of the system with different parameters were obtained by the nonlinear dynamics numerical analysis method. +e research shows that the system exhibits rich vibration response characteristics at different rotating speeds. +e amplitude of the vibration displacement in the system bifurcation diagram increases significantly with the increase of the tooth backlash and input torque, whereas the amplitude decreases constantly with the increase of the meshing damping. +e critical rotational speed at which chaotic motion occurs increases significantly with increasing input torque and damping ratio but decreases with increasing tooth backlash. +e bearing clearance has a weak influence on the vibration displacement amplitude of the system and the speed range of chaotic motion.


Introduction
Face-gear drive has a high contact ratio and reduction ratio, and the spur gear meshed with it is not affected by the axial force, which can make the structure of the gearbox more compact. erefore, it is mostly used in high-speed, heavyduty, and high-power helicopter main transmissions. In recent years, numerous fundamental research studies [1][2][3] have been conducted to investigate the meshing theory, contact characteristics, and manufacturing methods of facegear transmission. Litvin et al. [4] studied the application of face gear in helicopter transmission in 1992 and proved that face-gear drive has the advantage of weight reduction compared with bevel gear transmission and can achieve more accurate power allocation.
Since then, the split-torque transmission of face gear was systematically investigated in the A.R.T. (Advanced Rotorcraft Transmission) research project and the RDS-21(Rotorcraft Drive Systems for the 21st Century) program [5,6]. In the early 21st century, an Enhanced Rotorcraft Drive System (ERDS) program was implemented by the U.S. Army Aviation Applied Technology Directorate in cooperation with Boeing [7]. Recently, the research on CFGSTTS has gradually been becoming a popular issue, and Li et al. [8][9][10][11] conducted studies on the load sharing performance and power flow of this configuration.
However, the research on the nonlinear dynamic response characteristics of face-gear transmission mainly focused on the single face-gear pair. Jin et al. [12] established a multi-degree-of-freedom bending-torsion-axis vibration dynamic model of an orthogonal face-gear transmission by the lumped parameter method and analyzed the effects of the main parameters on the dynamic responses, bifurcation characteristics, and load coefficients. Lin and Ran [13] studied the dynamic responses of an orthogonal face-gear drive at different rotational speeds and the effects of load and friction coefficient on system vibration. In addition, Li et al. [14][15][16] formulated a bending-torsion coupling dynamic model of a nonorthogonal face-gear transmission system and analyzed the influences of modification parameters, tooth backlash, time-varying meshing stiffness, friction, and other excitation factors on the dynamic characteristics. Hu et al. [17][18][19] developed a coupling vibration lumped parameter dynamic model based on a single face-gear pair and studied the influences of the directional rotation radius, transmission error, load, clearance, and other factors on the dynamic response. Moreover, Liu et al. [20] investigated the dynamic characteristics of 8DOF helical face-gear transmission and discussed the chaotic motion and dynamic load characteristics through bifurcation diagrams, time history diagrams, phase diagrams, Poincaré diagrams, and dynamic load coefficient diagrams. Based on the lumped parameter theory, Dong et al. [21] gave a bending-torsion coupling dynamic model of face gear and studied the mathematical relationship between various factors and the dynamic response characteristics. Feng et al. [22] studied the modeling and analysis method of the torsional vibration of a face-gear pair, established a single-input single-output torsional vibration model, and investigated the influences of various factors on the dynamic characteristics. In addition, some scholars [23,24] carried out the dynamic load and load sharing performance study of the transmission system including face gear, but there are few articles on the nonlinear response of CFGSTTS. erefore, based on existing research, the dynamic characteristics of the CGSTTS are investigated by the lumped mass method. e research object is shown in Figure 1. At the beginning of the research, the meshing stiffness is calculated by the strain energy method, which can more accurately consider the influences of gear web and other structures, and is expressed as high-order harmonic. Furthermore, the meshing phase differences between different meshing pairs of the system are analyzed, and the expression of time-varying meshing stiffness including phase difference is formed. Finally, the nonlinear dynamic model of bending-torsion-axis is established, which includes multiple nonlinear factors such as time-varying meshing stiffness, meshing phase difference, meshing damping, static transmission error, tooth backlash, bearing clearance, and so on.
e dynamic differential equations of the system are solved by the 4th and 5th order variable step Runge-Kutta method, and the nonlinear vibration characteristics of the system and the influences of various parameters on the dynamic characteristics are obtained.

Time-Varying Meshing Stiffness and Meshing
Phase Difference

Calculation of Time-Varying Meshing
Stiffness. e flexible components in the transmission system have a great impact on dynamic performance. Mo et al. [25,26] considered the flexible parts in the dynamic system. e number of teeth participating in meshing is continuously changing during the transmission process, which causes the meshing stiffness to change. An accurate calculation of time-varying meshing stiffness is the basis for studying the dynamic characteristics of the face-gear system. However, to meet the strict weight and strength requirements of aviation gear, the web structure of the gear is generally more complex. In this article, for instance, the upper face-gear web has holes and stiffeners and its thickness is thin. Also, because the web is thin, there will be a large deformation of the web during the meshing process, which will affect the meshing stiffness, so the web should be considered in the calculation. is makes the calculation of meshing stiffness more difficult. erefore, to accurately consider the flexibility of the structure, the strain energy method is used to calculate the time-varying meshing stiffness of the upper and lower face-gear pairs under the quasistatic condition. e main design parameters of the face-gear pair are shown in Table 1.
In the quasistatic loaded contact analysis of gears, considering supporting structures such as shafts and bearings will significantly increase the amount of calculation, so the multi-tooth model is adopted in this paper, as shown in Figure 2. In finite element analysis, solid elements generally have no rotational degrees of freedom. Accordingly, the reference points P R1 and P R2 are, respectively, established on the rotation axis of the face gear and the cylindrical gear. en, they are coupled with the side and inner hole surfaces of the multi-tooth model, and then the calculation is solved through three load steps, as shown below.
(1) Firstly, all DOFs of the point P R1 are constrained, and a small rotational displacement around the rotation axis of the pinion is applied to the point P R2 to eliminate the initial clearance between the tooth surfaces. (2) en, all DOFs of the point P R2 are constrained, and the rotational degree of freedom of the point P R1 around the rotation axis of the face gear is released. At last, the torque is applied on the node P R1 by multiple load substeps. (3) Finally, we maintained the boundary conditions of the point P R1 and released the rotational degree of freedom of the point P R2 around the rotation axis of the pinion. e rotational displacement is introduced by the step-by-step loading method to calculate the strain energy of the gear pair at different meshing positions.
When the face gear meshes under quasistatic conditions, the loading process is slow, and the kinetic energy and dissipated energy of the gear can be ignored. erefore, the work done by the external force is stored in the gear pair as elastic strain energy. Assume that the extrusion stiffness, shear stiffness, bending stiffness, Hertz contact stiffness of gear teeth, and elastic stiffness of the gear body are k ia , k is , k ib , k h , and k ie , respectively. e meshing strain energy of gear teeth is shown as    Here, i � f, c, and it represents face gear and pinion, respectively. Face-gear mesh stiffness can be considered as the combination of k ia , k is , k ib , k h, and k ie , and then the total stiffness K is where F n is the normal meshing force and U m is the strain energy, which can be obtained by finite element analysis. Substituting the normal meshing force F n and the strain energy U m calculated by the finite element analysis into formula (2), the meshing stiffness of the face-gear pair can be obtained as shown in Figure 3. As illustrated in Figure 3, under the same gear parameters as well as load conditions, the value of the mesh stiffness of the upper face gear is smaller than that of the lower face gear. us, the flexibility of the web has a significant impact on the meshing stiffness. erefore, the meshing stiffness considering the flexibility of the web can better reflect the actual meshing situation. e meshing stiffness changes periodically with time, so it can be expanded according to the Fourier series expansion principle, and its expression is shown in (3). e coefficients acquired by the 3rd harmonic expansion of the meshing stiffness are listed in Table 2.

Meshing Phase Differences.
If the number of pinion teeth for the CFGSTTS is even, there are no meshing phase differences between all meshing pairs. However, if the number of pinion teeth is odd, there will be phase differences between different meshing pairs. erefore, the phase difference needs to be analyzed. To simplify the description, the gears 1 to 7 denote input gear 1, input gear 2, idler gear 1, idler gear 2, tail gear, and lower and upper face gears, respectively. Five pinions and two face gears form ten meshing pairs, which are marked as Mij (i � 1, 2, 3, 4, 5, j � 6, 7). Figure 4 shows the initial meshing state of M36. e plane P 1 is the middle plane of the tooth groove of gear 3, and the plane P 2 is the middle plane of the tooth of gear 6. In the initial meshing state, the plane P 1 coincides with P 2 . Figure 5 shows the initial meshing state of M37. e planes P 3 and P 4 are the middle planes of the tooth and tooth groove of gear 3, respectively. e planes P 5 and P 6 are the middle planes of the tooth groove and tooth of gear 7, respectively. Additionally, φ is the angle between planes P 3 and P 4 . In the initial state of M37, the plane P 3 coincides with P 5 .
Two benchmarks are needed to figure out the meshing phase differences, so the gear pair M36 and its initial meshing state are defined as the reference. Taking gear pair M37 as an example, after turning gear 3 counterclockwise by an angle φ, the plane P 4 will coincide with P 6 . At this moment, the meshing state of M37 is consistent with the reference state.
is means that the gear pair M37 lags behind the gear pair M36 by an angle φ � π/z 1 , where z 1 is the number of pinion teeth. en, we define the phase difference between gear pairs M37 and M36 as c M37 M36 c M37 M36 : where T m is the meshing cycle and n 1 is the rotational speed of pinion. e phase differences of other gear pairs are related to the distribution position of pinions. e installation positions and distribution angles of pinions 1-5 relative to idler 1 in the circumferential direction of face gear are shown in Figure 6. β i represents the angle between gear i and the coordinate axis z 3 . When dec(β i z 6 /360) � 1/2, (i � 1, 2, 3, 4, 5), there are phase differences between gear pair Mi6 and the reference M36, which are c M i6 M 36 � 1/2. e symbol "dec" denotes complementation. In contrast, there are no phase differences between Mi7 and M36. When dec(β i z 6 /2π) � 0, there are no phase differences between Mi6 and M36, and the phase differences between Mi7 and M36 will be generated. e phase differences between Mij and M36 are listed in Table 3.

The Nonlinear Dynamic Differential Equations
According to the gear installation position as shown in Figure 6, the local coordinate systems of pinion and face gear are defined as follows: (1) e x-axis and y-axis of the local coordinate system of the upper face gear are the same as those of the lower face gear, and the z-axis of both face gears point from the lower face gear to the upper. (2) e x-axis direction of pinion is tangential. e zaxis direction is from the inner end of face gear to the outer end, and the y-axis is determined by the righthand rule.
Different lubrication states have a great influence on tooth surface friction [27], which makes dynamics more complicated. In order to simplify the analysis, the tooth surface friction is not considered. In addition, all gears are simulated with mass and moment of inertia, and all shafts are considered as rigid bodies. Supporting bearings are also equivalently treated as stiffness and damping acting on each gear. According to the characteristics of face-gear drive, when the pinion is a spur gear, the pinion has no axial force. erefore, the degree of freedom of the pinion along its axis is not considered. Moreover, this article does not consider the shimmy of each gear. In summary, the transmission system has 23 degrees of freedom, expressed as {x p , y p , θ p , x q , y q , z q , θ q }, where x p , y p , and θ p (p � 1, 2, . . ., 5) represent the translational vibration of the pinion and the torsional vibration around its respective rotation axis, and x q , y q , z q , and θ q (p � 6, 7) represent the translational and torsional vibration of the upper and lower face gear.
Generally, in the CFGSTTS, the power delivered from input gear 1 and input gear 2 is transferred partly to the upper face gear and partly to the lower face gear. en, the lower face gear transmits partially the power to the upper face gear through the idler gear 1 and idler gear 2, and part of the power is output through the tail gear. e schematic diagram of power flow direction is illustrated in Figure 7. erefore, according to the force characteristics of pinion and face gear, the nonlinear dynamic model is shown in Figure 8.
Due to the vibration displacement and transmission error, the relative displacements δ n6 and δ n7 between pinion n (n � 1, 2) and face gears 6 and 7 along the meshing line are given by

Mathematical Problems in Engineering
Tail gear Idler gear2 Input gear1

Input gear2
Idler gear1 x 5 x 4 Figure 6: e distribution and angles of the pinion.   Mathematical Problems in Engineering where r n is the indexing circle radius of gear n; r f is the equivalent radius of face gear; β n is the angle of gear n relative to idler gear 1; and e n6 and e n7 are the static transmission errors of the gear pair Mn6 and Mn7, respectively, which can be written as where e 0 and e r are the constant value and variable amplitude of the static transmission error, respectively, and φ r is the initial phase. e relative displacements δ k6 and δ k7 along the meshing line between pinion k (k � 3, 4, 5) and gears 6 and 7 are, respectively, e normal dynamic load can be expressed as where c pq is the meshing damping of gear p and gear q, obtained by where m p and m q represent the equivalent mass of gears p and q; ξ is the damping ratio of the meshing pair; and k av,q is the average value of the meshing stiffness.
e nonlinear function of gear backlash, f(δ pq ), between gear p and gear q can be expressed as where b h is half of the gear backlash.
Considering that there is a gap in the supporting bearing, it can also be expressed by a gap nonlinear function, f(I).
Assuming that the gap is b b , then f(I) can be given by Mathematical Problems in Engineering Here, I is the vibration displacement of the gear in its local coordinate system.
According to Newton's second law, the dynamic differential equations of pinions 1-5 in the CFGSTTS can be obtained as follows: where m p (p � 1, 2, 3, 4, 5) are the masses of gear p; J p is the moment of inertia of gear p around the rotation axis; k px and k py represent the radial stiffness of the supporting bearing of gear p, respectively; and c px , and c py represent the radial damping of supporting bearing of gear p, respectively. e dynamic differential equations of the lower and upper face gears are given as (17) and (18), respectively.

Mathematical Problems in Engineering
where m q (q � 6, 7) are the masses of gear q; J q is the moment of inertia of gear q around the rotation axis; k qx and k qy represent the radial stiffness of the supporting bearing of gear q, respectively; and c qx , and c qy represent the radial damping of the supporting bearing of gear q, respectively.

Nonlinear Dynamic Characteristics
Dynamic differential equations (12)- (18) are strongly nonlinear with variable parameters, which are difficult to solve by analytical methods. us, the 4-5 order variable step adaptive Runge-Kutta method is used to solve it in this paper. In the solution, the integration step is set to T m /600, and the integration time is taken as 1000T m . e parameters involved in the system equations are listed in Tables 1-4, respectively. Besides, some common qualitative analysis methods, including time history, phase trajectory, Poincaré section, and FFT spectrum, are used to analyze the nonlinear vibration characteristics of the system. Figure 9 shows the global bifurcation characteristics of the relative displacement δ 17 with the rotation speed Ω. It can be seen from the figure that with the increase of rotating speed, the system exhibits rich nonlinear behaviors and complex motion states, including periodic motion, multi-periodic motion, quasiperiodic motion, and chaotic motion. Quasiperiodic motion and periodic motion appear alternately. Some periodic windows are narrow, while others last for a large rotation speed range. us, the system can operate in a more stable periodic motion by adjusting the input speed under the actual operating conditions.
To further explore the system response and chaotic behavior, the motion characteristics at a specific speed are given in Figures 10-17. In the range of Ω � 1000 rpm-1575 rpm, the system behaves as a stable periodic motion, as shown in Figure 10. According to the time history, the displacement response of the system is periodic. e trajectory in the phase diagram is a noncircular closed curve, and the intersection of the phase trajectory and Poincaré section is a single discrete point. Furthermore, the spectrum diagram shows that there are amplitude components at the frequency of if 0 (i � 1, 2, 3), which indicates that there are super-harmonic components in the system response. erefore, the motion of the system is single periodic at Ω � 1500 rpm.
However, when the input speed changes in the range of 1600 rpm-2175 rpm, the system response constantly switches between chaotic, single-periodic, and multi-periodic motions. Figure 11 presents the chaotic motion when the rotational speed Ω is 2000 rpm. At this speed, the vibration displacement in the time history diagram has no obvious regularity, the phase portrait exhibits that the phase trajectory is a staggered curve filled with a certain space, and the Poincaré map is an irregular point set. e FFT spectrum has some components of frequency division and frequency doubling, and there are obvious continuous sidebands.
When the rotation speed is 2200 rpm, the 2T-periodic subharmonic response of the system appears under the action of the limit cycle attractor.
e system response diagrams are illustrated in Figure 12. e displacement curve repeats periodically, the corresponding phase trajectory is a noncircular closed curve, and two discrete points appear in the Poincaré cross section diagram. As shown in the spectrum diagram, there is an amplitude component at a frequency of f 0 /2, which indicates that the system response has a subharmonic component and its motion is 2T-periodic. After that, the system enters a stable 1T-periodic motion at 2875 rpm through reverse bifurcation. When the rotation speed is 3250 rpm, the system enters a chaotic state through the path of period-doubling bifurcation under the action of the limit toroidal attractor. Subsequently, the system in the speed range of 3250-4050 rpm has short 2T-periodic and 4T-periodic windows near Ω � 3800 rpm, although it mainly behaves as chaotic motion. As an example, Figure 13 depicts the response of the system at Ω � 3650 rpm, which is a 4T-periodic motion. As the rotating speed approaches 4075 rpm, it re-enters the 1T-periodic movement through backward bifurcation. 7000 8000 9000 10000            When the input speed is in the range of 4075 rpm-8220 rpm, the system response is mainly expressed as a stable 1T-periodic motion. As Ω is further increased to around 7625 rpm, the 1T-periodic motion converts to a 2T-periodic movement that only maintains a small range of speeds, as shown in Figure 14. Especially, there exists the jump phenomenon of response amplitude at Ω � 8300 rpm, but the motion state does not change before and after the jump. After that, the system state changes from single-periodic motion (Ω � 8300 rpm-8475 rpm) to 2Tperiodic motion at 8500 rpm and then sustains 2T-periodic motion in a wide speed range.
At higher system speed (Ω � 8950 rpm), the system response transits to single-periodic motion. However, when    the system speed increases to 9075 rpm, the system exhibits quasiperiodic response characteristics, as shown in Figure 15. e movement of the system in the time history diagram is approximately periodic. Also, the trajectory in the phase plane is a closed curve with a certain width. e Poincaré cross section diagram shows a closed curve, and the frequencies in the spectrum diagram are incommensurable. In summary, the system motion is quasiperiodic at Ω � 9100 rpm.
After that, the response amplitude of the system increases gradually and then enters the chaotic state, as shown in Figure 16 (Ω � 9425 rpm). From the time history diagram, the system movement is aperiodic.
e Poincaré map presents a set of points in a limited area. However, the system re-enters the quasiperiodic movement as the speed increases, as shown in Figure 17 (Ω � 9900 rpm). ere are periodic components in the displacement curve as observed from the time history diagram. e phase orbit is a closed curve band with a certain width, and the Poincaré map shows a closed curve. In addition, the spectral lines in the spectrum diagram are discretely distributed at the points of the combined frequency. erefore, the motion of the system at Ω � 9900 rpm is quasiperiodic.

Influence of Other Factors on Bifurcation Characteristics
Figures 18-21 describe the global bifurcation characteristics of the system with speed at different damping ratios, input torques, tooth backlashes, and bearing clearances. It is evident from Figure 18 that as the damping ratio increases, the energy dissipated by the system continues to rise, and the speed range and amplitude of the chaotic motion gradually decrease. e jumping point of the movement state at higher speed gradually moved backward, but the position of the bifurcation point at lower speed did not change significantly. Consequently, the increase of the damping ratio expands the range of rotation speed for stable periodic motion. As shown in Figure 19, with the increase of input torque, the width of the window of chaotic motion tends to decrease, but the amplitude progressively increases. At a higher speed, the critical speed of the system entering the chaotic state gradually moves backward. Notably, when the input torque is 500 Nm and 800 Nm, the system response shows a bifurcation point and jumping point around Ω � 4000 rpm.
is means that the increased torque is beneficial to system stability.  In Figure 20, when the tooth backlash is varied from 20 μm to 50 μm, the system response mainly shows a stable single-periodic motion, but the vibration displacement amplitude gradually increases. As the backlash increases to 80 μm and 120 μm, the system appears to be in chaotic motion, and the speed range of chaotic movement and the amplitude of vibration displacement increase simultaneously. As a result, the increase of tooth backlash leads to a larger speed range of unstable system motion.
As shown in Figure 21, with the increase of bearing clearance, there is no significant change in the amplitude of vibration displacement of the system as well as the location and form of the bifurcation point. In addition, a significant number of examples show that the variation of the bearing gap has a weak effect on the system response.

Conclusions
Based on the strain energy theory, a method for calculating the meshing stiffness of face gear with web structure is developed in this paper. Considering the meshing phase differences caused by the odd number of pinion teeth, a 23degree-of-freedom bending-torsion-axis nonlinear dynamic model is established. e model includes time-varying meshing stiffness, meshing damping, tooth surface error, tooth side clearance, bearing clearance, and other nonlinear factors, which can comprehensively describe the vibration characteristics of the CFGSTTS. By using the Runge-Kutta numerical analysis method, the global bifurcation characteristics of the system with speed are studied, and the effects of damping ratio, load, tooth backlash, and bearing clearance on the bifurcation characteristics are analyzed. e conclusions of this paper are as follows: (1) e dynamic response characteristics of the system are complex. With the change of speed, it presents four kinds of steady-state responses, including simple harmonic response, subharmonic response, quasiperiodic response, and chaotic response. (2) Increasing the damping ratio can gradually reduce the width and amplitude of the chaotic motion area so that the system can obtain a wider stable motion range. As the input torque increases, the width of the chaotic motion window of the system gradually decreases, while the vibration amplitude increases. By contrast, the range of unstable motion of the system is enlarged by the increase of tooth backlash. All in all, the motion state of the system can be controlled by selecting the appropriate parameters. it is proved that the influence of bearing clearance on the system motion is weak. In some cases, the bifurcation cannot be observed with the bearing clearance as the parameter.

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.