Abstract
Considering time-varying meshing stiffness, comprehensive errors, and piecewise backlash nonlinearities of gear and spline, a torsional nonlinear dynamic model of star gear-rotor coupling transmission system of (Geared Turbofan Engine) GTF aeroengine is established. By using the Runge–Kutta numerical integration method, the dynamic responses are solved, analyzed, and illustrated with the bifurcation parameters including input rotational speed, gear backlash, damping ratio, and comprehensive meshing errors. The motions of the star gearing system and diverse nonlinear dynamic characteristics are identified through global bifurcation, FFT spectra, Poincaré map, and the phase diagram. The results reveal that the star gear-rotor system exhibits abundant torsional nonlinear behaviors, including multiperiodic, quasi-periodic, and chaotic motions. Furthermore, the roads to chaos via quasi-periodicity, period-doubling scenario, and mutation are demonstrated. These results provide an understanding of undesirable torsional dynamic motion for the GTF transmission system and provide a reference for the design and control of gear system.
1. Introduction
The most advanced GTF aeroengine currently uses a star gearing system as the main decelerator for fan, allowing the low-pressure rotor to operate at high efficiency and high speed to match the optimum speed of the fan. The star gearing drive is shared by multiple fixed-axis star gears. Fixed-axis gear train is used in the entire transmission system; therefore, the strength, rigidity, working life, and reliability of the system are improved; however, harsh operating condition may bring diverse nonlinear dynamic behaviors, which are investigated in this work.
The gear pair with multiple nonlinear parameters were studied in many kinds of research; the effects caused by nonlinearity were detected. Nadolski and Pielorz [1] investigated the distribution nonlinear loads on gear teeth theoretically. Chang-Jian and Chang [2] conducted dynamic analysis on the spur gear pair of single-degree-of-freedom and compared the vibration response with and without nonlinear suspension. He et al. [3, 4] established a new analytical spur gear model incorporating sliding friction and realistic time-varying stiffness, and the effects were detected. Yang et al. [5] analyzed nonlinear time-varying dynamic model for right-angle gear pairs, and the effects of backlash and asymmetric mesh were considered with an enhanced multiterm harmonic balance method. Wang et al. [6] proposed a new method for double helical gears based on modified optimization to minimize vibration amplitude while raising transmission efficiency; the optimal modification parameters were obtained through the minimum amplitude of transmission error. Wang et al. [7] studied the gear pairs related to the friction, backlash, and time-varying gear meshing stiffness. The bifurcation, chaos, and corresponding maximum Lyapunov exponent of the gear system are studied, and the key parameters are determined. Yang et al. [8] conducted vibration analysis for spur gear system with the clearance nonlinearity to detect tooth crack detection. Wang et al. analyzed nonlinear dynamics of spur gear system for railway locomotive [9], established multi-degree-of-freedom spur gear system [10], and investigated the nonlinear characteristics of the system with bending-torsional coupling vibration.
Then the nonlinear dynamic analysis on gear-rotor-bearing system was investigated, in which the nonlinear behavior became more complex due to the nonlinearity brought with bearings. Gurkan and Nevzat Özgüven [11] studied the effects of backlash and bearing clearance in a geared flexible rotor and the interactions between these two nonlinearities. Baguet and Jacquenot [12] analyzed the interactions between gears, shafts, and hydrodynamic journal bearings of gear drives, and the time-varying properties and nonlinearity of gears and bearings were introduced. Siyu et al. [13] established a nonlinear gear rotor transmission model and investigated the effect of gear friction and backlash on the system; the results indicated that the vibration magnitude enlarged with the dynamic backlash, which may bring the system into chaotic motion. Chang-Jian [14] investigated the gear-rotor-bearing system with nonlinear effects including nonlinear rub-impact force, nonlinear meshing force, nonlinear couple-stress fluid film force, and nonlinear suspension effect. Han et al. [15] studied a geared rotor system with slant cracked shaft under multiple internal and external excitations and found that torsional excitations have significant influence on the forced response spectra. Zhou et al. [16] established MDOF lumped parameter dynamic model of the spur gear rotor bearing system with gear backlash, transmission error, eccentricity, and gravity and also considered the coupled lateral-torsional vibration and the coupled multibody dynamics.
As gearbox has various applications in wind turbine and aero industry, the nonlinear and dynamic analysis for the gearbox, especially for those composed of planetary gears, were studied. Guo and Parker [17] studied the nonlinear behavior of tooth wedging related to planet bearing forces by analyzing the dynamic response of the planetary gears in wind turbine. Kim et al. [18] proposed a new dynamic model for the gear set considering torsional motion with the time-varying pressure angle and contact ratio, and the dynamic responses between the modified and previous models were demonstrated. Bahk and Parker [19] conducted dynamic analysis on the planetary gears, where the nonlinear behaviors caused by the time-varying mesh stiffness at resonance exhibited jump phenomena and subharmonic vibration. Rao et al. [20] performed nonlinear torsional analysis of multistage gear system with flexible shafts. Xiang et al. [21, 22] investigated the nonlinear dynamic of multistage gears, in which the nonlinear parameter multiclearances were considered. Hou and Cao [23] built up pure rotational model of Planetary Gears-Rotor System in Geared Turbofan Engines and carried out nonlinear dynamic analysis to obtain the torsional vibration responses between gear and rotor. Shuai et al. [24] built precise model for complex tooth surface microtopography and carried out nonlinear friction dynamics analysis for high performance face gear. Mo et al. [25, 26] also conducted analytical investigation of load sharing characteristics for herringbone planetary gear train with flexible support and multipower face gear split flow system, respectively. Wang et al. [27] established the model of gear-rotor-bearing system of GTF gearbox with flexible support of planet carrier and studied the load sharing behavior of the system.
According to our extensive search throughout the existing literature, there are a few studies specialized in nonlinear dynamics analysis of the transmission system for GTF gearbox; furthermore, the effect of interaction between multiple nonlinear parameters (gear backlash, spline flank clearance) on dynamic response of the transmission system is rarely seen in those studies. Therefore, it is essential to take all the nonlinear factors into the model for the more in-depth research.
In this paper, a nonlinear torsional dynamic model of transmission system of GTF gearbox is proposed, in which the influence of crucial nonlinear parameters, including input rotational speed, gear backlash, comprehensive meshing errors, and meshing damping on the vibration response of the transmission system, is explored and quantified through time zone and frequency zone analysis, and the nonlinear behavior is studied by analyzing the states of motion, mechanisms of bifurcation, and roads to chaos.
2. Dynamic Model of Transmission System of GTF Gearbox
2.1. GTF Gearbox Transmission System
Figure 1(a) shows the physical model of GTF gearbox (highlighted in red color), and Figure 1(b) shows the schematic diagram of GTF gearbox transmission system with elastic support base of planet carrier. Based on the elastic support model of the planet carrier in the GTF gearbox, the fan-driven gearbox adopts a five-way shunting herringbone gear transmission structure, which mainly includes the sun gear, star gear, ring gear, planet carrier, and input and output shaft. The sun gear is a floating part, which is splined with the input shaft and meshes with five circumferential evenly distributed star gears. Star gears are internally supported by bearing and meshed with the ring gear, which is a semi-floating component and taken as the output end of the gear train, connected to the output shaft by the bolts and the gravity of the entire gearbox, and the torque generated by the gear meshing is carried through planet carrier supported by elastic base.

(a)

(b)
2.2. Dynamic Model of Spline Joint
In the GTF gearbox transmission system, the external spline on input shaft is connected with the internal spline inside the sun gear, and the input torque is transmitted to the star gear train through the external spline of the input shaft and the spline inside the sun gear, the schematic diagram of spline connection is shown in Figure 2(a), and the curve of piecewise linear function of flank clearance for spline joint is shown in Figure 2(b).

(a)

(b)
When the spline joint is subjected to the torque, due to the existence of the spline flank clearance, as shown by and in the Figure 2(a), when the spline pair is stressed, the spline teeth with small backlash are contacted first, and then the teeth with bigger clearance; therefore, the actual number of teeth engaged is less than the number of teeth designed [28]; the piecewise function of the spline flank clearance for calculation is shown as equation (1).
In equation (1), the is the original minimum clearance at each side of spline joint pair, is the number of teeth of spline, and is the equivalent rotational displacement of spline joint pair.
2.3. Dynamic Model of Star Gearing System
Figure 3(a) shows the physical model of star gearing system; Figure 3(b) shows the schematic diagram of the dynamic model of the star gearing system. Since the gears of the star gearing system are characterized by fixed-axis rotation, an absolute coordinate system is established for the system. Consider the torsional displacement of each component. The subscripts , , represent the sun gear, star gears, and ring gear, respectively.

(a)

(b)
In the dynamic model of star gearing system, , , , represent the meshing stiffness and damping of sun-star meshing pair and star-ring meshing pair. The value of stiffness can be obtained by adopting the calculation formula proposed by Matter and Velex [29]. In addition, the value of damping can be obtained from the following formulas.
For the meshing damping of gear pair, we can use formula (2) to calculate. and are the meshing damping ratio of gear pair; in general, the value ranged from 0.03 to 0.1. , , and are the radius of base circle for sun gear, star gear, and ring gear. , , and are rotational inertia for sun gear, star gear, and ring gear.
For the torsional damping of gear pair, we can use formula (3) to calculate. is the torsional damping ratio of gear pair, in general, the value ranged from 0.005 to 0.075, is the torsional stiffness of driving gear, and is the torsional stiffness of driven gear, is the rotational inertia of driving gear, and is the rotational inertia of driven gear.
The elastic meshing force and the viscous meshing force of the internal and external meshing pairs can be expressed as
In the previous formula, and are the nonlinear function of gear backlash, which can be described by the following uniform form.
In the previous formula, indicates the relative displacement of external or internal meshing pairs and represents the half gear backlash of external or internal meshing pairs.
Due to the nonlinear function of gear backlash, there is transition between single-sided contact, double-sided contact, and no contact during operation of the system. Figure 4 shows the piecewise nonlinear function ; it can be drawn from the figure that the impact is not observed in a gear system when the displacement lies in the region and , which is shown in Figure 4 as zone I. Double-sided impact exists when the displacement lies in the region and , shown as zone II, while the displacement lies in the region and or and ; the system presents single-sided impact, shown as zone III.

2.4. Dynamic Model of Transmission System of GTF Gearbox
Figure 5 shows the dynamic model of the GTF gearbox transmission system; a pure torsional nonlinear dynamic model of the star gearing system and input and output shaft with spline joints connection is established by the lumped mass method.

By integrating the input and output rotor, spline joints, and star gearing system, the overall coupled dynamics model of the GTF gearbox can be obtained. The model has 9 + N degrees of freedom. The generalized displacement array of the system can be expressed aswhere represents the torsional displacement of each component in the transmission system and is the number of star gears.
3. Mathematical Model of Transmission System of GTF Gearbox
The kinetic and potential energy of the GTF transmission system are expressed regarding the generalized coordinates; then the equations of motion can be derived using Lagrange’s equation; all the components of the system are assumed to have torsional motion; therefore, the kinetic energy is rotational kinetic energy. The total kinetic energy of the transmission system is composed of the kinetic energy of sun gear, star gears, ring gear, spline connections, and shafts, which can be expressed by
The potential energy of the transmission system can be expressed aswhere each term is the contributions of each spline meshing pair and gear meshing pair, respectively.
When the input shaft is subjected to only an external moment in the counterclockwise direction, all the generalized forces are zero, except the generalized force corresponding to the coordinate .
The equations of motion for the GTF transmission system are derived by using Lagrange’s equation, as follows:
Substitution of equations (7), (8), and (9) into equation (10) leads to the equations of motions (dynamic differential equation) of the transmission system.where is the moment of inertia, is the angle of rotation, is the input or output torque, is the pitch radius of the outer or inner spline, is the base circle radius of members in star gearing system, is helix angle of base circle, and represent torsional damping and contact damping, respectively, and are torsional stiffness and contact stiffness, and represents the piecewise function of spline flank clearance.
In order to eliminate the rigid body displacement, the following variables were introduced.
Substituting the relative displacements into equation (14) yielded the governing equation.
The dimensionless time parameter and the dimensionless displacement parameter were used for nondimensionalized equations with and . The dimensionless parameters, including the displacement, velocity, acceleration, gear backlash, and meshing error, can be defined as
4. Analysis of the Nonlinear Dynamic Responses
The nonlinear behaviors of the model can be obtained by focusing on the relative torsional vibration responses of external meshing pair and spline joint pair. The influence of key parameters, including the rotational speed of input shaft, gear backlash, comprehensive meshing error, and damping ratio of meshing pair, on the nonlinear dynamic response was studied. The basic parameters of the star gearing system used in this paper are shown in Table 1. The dimensionless ordinary differential equations of the system were solved by the fourth-order Runge–Kutta numerical integration method. The influence rules of excitation parameters on the system were characterized by means of bifurcation diagrams, Poincaré maps, phase diagrams, and FFT spectrum. There are two spline joints in the fan drive gearbox. One is used to transmit torque between the input shaft and sun gear, and the other is used to transmit torque between the output shaft and fan. The basic parameters of the spline are as shown in Table 2.
Figure 6 shows the mesh stiffness of external and internal meshing pairs; it can be drawn from the stiffnesses in a period that internal mesh stiffness is larger than that of external mesh, which is caused by the increase of overall contact ratio for internal meshing pairs.

(a)

(b)
4.1. Effect of Input Rotational Speed on the Dynamic Response of the System
In this section, the vibration response of the GTF transmission system will be conducted as the variation of input rotational speed. Keeping other parameters the same, the bifurcation diagram of relative displacement of external meshing pair is shown in Figure 7. Here, the input rotational speed (r/min) is treated as the control parameter varying in the range of . With the increase of input rotational speed, it can be observed that the system undergoes periodic motion, quasi-periodic motion, and chaotic motion subsequently. The gear system exhibits T-periodic motion when for relative external meshing displacement; with the increase of input rotational speed from to , the periodic motion is replaced by quasi-period motion; when the speed increases, the gear system exhibits chaotic motion through quasi-periodicity until the speed arrives at 8700 rpm. As the control parameter is further increased, the system turns back to T-periodic motion by inverse bifurcation while the input rotational speed varies from to , then the system shows quasi-periodic motion as the control parameter is in the range from to , and finally the system performs T-periodic motion as . In order to accurately determine the motion form, the phase plane diagrams, Poincaré maps, time domain response diagrams, and amplitude-frequency spectrums are plotted to illustrate dynamic response of the system in Figure 8.


In this section, the vibration response of the gear system in the relative displacement of external meshing pair will be analyzed at three input rotational speeds of ,, and . As shown in Figure 9, the gear system shows different vibration features in different speed. When the control parameter , the first-order meshing frequency is the dominant response and the amplitude of first-order meshing frequency is much bigger than the response of other orders of meshing frequency in the frequency spectrum under logarithmic representation, besides first five orders of meshing frequency, the side frequency is not obvious, the Poincaré map corresponding to has several points, and the phase diagram reveals single line. Thus, it indicates that the system performs nT-periodic motion. As the control parameter , although first-order meshing frequency is still the dominant response, the de-multiplication and multiplication frequency components can be observed. The Poincaré map shows a disordered point set and the phase diagram is a chaotic attractor with a random crew, which means that the system exhibits a chaotic behavior. As the control parameter arrives at , the side frequencies are not that complex but also evident, the Poincaré map becomes a phase-locked loop, and the phase diagram has a trajectory torus, which demonstrates that the gear system changes into quasi-periodic motion.

4.2. Effect of Gear Backlash on the Dynamic Response of the System
As a primarily nonlinear excitation in gear system, the effect of gear backlash on the response of the dynamic system is analyzed. Figures 10(a) and 10(b) illustrate the bifurcation diagram of relative displacement of spline joint pair of input shaft and external gear meshing pair, respectively, with the half gear backlash as control parameter. It can be drawn that the system exhibits T-periodic motion when ; as increases while in the range from 5 to 9 , the T-periodic motion is replaced by the quasi-periodic motion; as the backlash further increases beyond 9 , the system turns into chaotic motion. The road to chaos is via period-doubling scenario, and the vibration amplitude becomes larger as chaotic motion appears.

(a)

(b)
Figure 8 presents the vibration response for , , and , respectively. It can be drawn from the frequency spectrum that as the control parameter increases, although the first-order meshing frequency is the dominant response, the amplitude of higher order meshing frequency becomes larger. The system undergoes T-periodic motion when with only two points at the Poincaré map and a single trajectory line at the phase diagram. As the backlash increases to , the Poincaré map shows phased-locked loops and the phase trajectory repeat themselves every period with different adjacent trajectory, which illustrate the quasi-periodic motion. As the gear backlash , the chaotic motion is demonstrated by the appearance of a cloud of points at the Poincaré map and a typical strange attractor of chaos with the random phase trajectories at the phase diagram.
4.3. Effect of Meshing Damping Ratio on the Dynamic Response of the System
The bifurcation diagrams of relative displacement of spline meshing pair of input shaft and external gear meshing pair of the system are shown in Figure 11, where the damping ratio of gear meshing was taken as control parameter. From the bifurcation diagram, chaotic motion can be observed within damping ratio ; then it turned into quasi-periodic motion through inverse bifurcation while the control parameter ; as the damping ratio further increased beyond , the system turned into T-periodic motion, which indicates that increase of damping ratio of gear meshing could effectively suppress chaos motion.

(a)

(b)
Figure 12 shows vibration response for , , and ; from the frequency spectrums, it can be drawn that, with the increase of damping ratio , the amplitude of higher orders of meshing frequencies shows descending trend as the first order remains the dominant response. When , a chaotic attractor with random phase trajectory exists at the phase diagram and randomly positioned points in the Poincaré map, which demonstrates that the system undergoes chaotic motion. As the control parameter approaches , the system turns into quasi-period motion through T-periodic motion, the Poincaré map reveals a phase-locked loop, and the phase diagram presents a trajectory torus. As the damping ratio arrives at , the system experiences a process of inverse bifurcation and finally enters T-periodic motion, which can be proved by the several points in Poincaré map and phase-locked loop in phase diagram.

4.4. Effect of Comprehensive Meshing Error on the Dynamic Response of the System
The comprehensive meshing error in the star gearing system is more complicated on account of the correlation of external and internal meshing. Meanwhile, this error brings nonlinear excitation because of the rotation and revolution of star gears. Figure 13 shows the bifurcation diagram of relative displacement of external meshing pairs as the comprehensive meshing error is treated as the control parameter. It can be concluded from the diagram that as the meshing error reaches (), the motion of the system transfers from T-periodic motion to quasi-periodic motion and to chaotic motion when . With the further increase of comprehensive meshing error, the chaotic motion is replaced by quasi-periodic motion through T-periodic motion when ; then the motion turns back to chaotic motion when ; finally the motion becomes T-periodic motion when the control parameter is beyond .

Figure 14 shows vibration response for , , and , respectively. For the spectrum, the higher orders of meshing frequency become more obvious and more frequency components appear when the T-periodic motion is replaced by quasi-periodic motion and chaotic motion. The Poincaré map has only several points, and the phase diagram has single trajectory at , implying the T-periodic motion. When the control parameter reaches , the return points become a fractal structure of strange attractor at the Poincaré map and there are random chaotic trajectories at the phase diagram, which indicate that the system is in chaos. As for , the Poincaré maps exhibit phase-locked loop and the phase diagrams reveal trajectory torus; thus, the star gearing system presents quasi-periodic motion.

5. Conclusion
In this paper, the coupling star gear-rotor dynamic model with spline connection for GTF aeroengine was established, based on which the pure rotational nonlinear dynamic analysis was conducted. The influences of pivotal parameters, including the input rotational speed, gear backlash, comprehensive meshing error, and damping ratio of gear meshing, on the torsional nonlinear response were studied. The dynamic characteristics were identified and illustrated by bifurcation diagrams, Poincaré maps, phase diagrams, and spectrum. The conclusions drawn from the study are summarized as follows:(1)The increase of input rotational speed could bring the system into chaotic motion; however, as the speed increases, the system undergoes T-periodic motion, quasi-periodic motion, and chaotic motion and finally turns into T-periodic motion through inverse bifurcation. The nonlinear dynamic characteristics such as continuous frequency components can be observed due to the increase of input rotational speed; therefore, higher input speed should be avoided, especially for certain intervals of speed range, which could cause chaotic motion for the gear system.(2)With the increase of gear backlash, the amplitude of relative displacement of spline joint pair and external meshing pair becomes higher, and the system undergoes T-periodic motion, quasi-periodic motion, and chaotic motion successively; therefore, the gear backlash should be constrained as small as possible.(3)When the damping ratio of gear meshing is relatively small, the system is prone to entering into chaos motion; while the damping ratio increases, the chaotic motion is replaced by the quasi-periodic motion; with further increase of the damping ratio, the system turns into T-periodic motion, which means higher damping ratio could constrain the vibration amplitude of relative displacement of spline joint pair and meshing pair; thus it inhibits the system form getting into chaotic motion.(4)The effect of comprehensive meshing error on the response of the system is complex; there are multiple transitions from quasi-periodic motion and chaotic motion; nevertheless, when the comprehensive error is too small or too large, the response of system remains in T-periodic motion. Therefore, the factors leading to the meshing error should be controlled in order not to bring the system into chaotic motion.
Data Availability
All the data studied in this article are derived from a specific type of GTF gearbox transmission system, ensuring that the research content is close to the actual work situation.
Conflicts of Interest
The authors declare that they have no potential conflicts of interest.
Acknowledgments
This research was supported by the National Natural Science Foundation of China (Grant no. 51775265).