Bifurcation Analysis for a Simple Dual-Rotor System with Nonlinear Intershaft Bearing Based on the Singularity Method

*is paper aims to classify bifurcation modes for two interrelated primary resonances of a simple dual-rotor system under double frequency excitations. *e four degree-of-freedom (4DOF) dynamic equations of the system considering the nonlinearity of the intershaft bearing can be obtained by using the assumed mode method (AMM) and Lagrange’s equation. A simplified method for dynamic equations is developed due to the symmetry of rotors, based on which the amplitude frequency equations for two interrelated primary resonances are obtained by using the multiple scales method. Furthermore, the validity of the simplified method for dynamic equations and the amplitude frequency equations solved by the multiple scales method are confirmed by numerical verification. Afterwards, the bifurcation analysis for two interrelated primary resonances is carried out according to the two-state-variable singularity method.*ere exist a total of three different types of bifurcation modes because of double frequency excitations of the dual-rotor system and the nonlinearity of the intershaft bearing.*e second primary resonance is more prone to have nonlinear dynamic characteristics than the first primary resonance. *is discovery indicates that two interrelated primary resonances of the dual-rotor system may have different bifurcation modes under the same dynamic parameters.


Introduction
With the development of aero-engines towards higher thrust-weight ratio and better aerodynamic stability [1], the dual-rotor structure has been widely used in the field of aero-engines. As a vital part between the lower pressure (LP) rotor and the higher pressure (HP) rotor in the dualrotor system especially in aero-engines of fighters, the nonlinearity of the intershaft bearing [2] has achieved considerable attention [3].
ere exist double frequency excitations, i.e., unbalanced excitations of LP and HP rotors, in the dual-rotor system due to the complexity of its structure. Moreover, under the nonlinear influence of the intershaft bearing, the dynamic behaviors of dual-rotor system will become much more complex and difficult to predict. erefore, to clarify bifurcation modes of the dualrotor system under different dynamic parameters is very contributive to avoid some unexpected nonlinear dynamic behaviors [4,5], so as to improve the operation safety of aero-engines. ere indeed exists some literature about the dynamic behaviors of the dual-rotor system in recent decades. Hibner [6] proposed a unique method, namely, the transfer-matrix method, to estimate vibratory responses of a dual-rotor system with the nonlinear damping coefficient in an aeroengine. Gupta et al. [7] developed a dual-rotor test rig to simulate the two spool aero-engine dynamically and compared experimental results and theoretical results for unbalance responses of the dual-rotor system. Ferraris et al. [8] presented a simple but realistic nonsymmetric coaxial dualrotor system and focused on the dynamic behaviors of the dual-rotor system under corotation and counterrotation. Different from Gupta et al. [7], a dual-rotor test rig composed of two coaxial shafts was developed by Guskov et al. [9], where two shafts are linked by an intershaft bearing and rotate individually. Zhang et al. [10,11] successively put forward the whole-beat correlation method and the nonwhole beat correlation method to identify unbalance responses of a dual-rotor system with a slight rotation speed difference. Yang et al. [12] established a dual-rotor system with imbalance-rubbing coupled fault to study dynamic responses of aero-engines with fixed-point rubbing. Wang et al. [13,14] also concentrated on the dynamic responses a dual-rotor structure of a whole aero-engine with bladecasing rubbing. Sun et al. [15,16] investigated the nonlinear dynamical behaviors of a dual-rotor aero-engine with rub impact by means of an analytic method, namely, MHB-AFT (multiharmonic balance-alternating frequency/time domain). Nevertheless, the nonlinearity of the intershaft bearing is not considered in these dual-rotor models, and none of the above literature have classified bifurcation modes of the dual-rotor system. e singularity theory is a mathematical method to study the morphology and classification of differentiable mapping near singular points [17]. Until 1985, Golubitsky and Schaeffer [18] elaborated the mathematical deductions and applications about the singularity theory in their books, based on which the singularity theory became a powerful mathematical tool to investigate bifurcation modes of the nonlinear dynamic system. Chen et al. [19,20] applied the steady-state bifurcation theory of the dynamic system to investigate topological structures of bifurcation modes for a nonlinear dynamic system, and they proposed a new method, namely, C-L method, to study bifurcation modes for the nonlinear dynamic system of a single bifurcation parameter. Furthermore, Qin et al. [21,22] developed the single state-variable singularity method into the two-state-variable singularity method through rigorous mathematical proof and gave mathematical formulas of the transition sets, i.e., the bifurcation set, the hysteresis set, and the double limit set. is research applied the two-state-variable singularity method to investigate bifurcation modes for the nonlinear dual-rotor system. e purpose of this paper is to investigate bifurcation modes for two interrelated primary resonances of a dual-rotor system based on the two-state-variable singularity method, where the double frequency excitations of the dual-rotor system and the nonlinearity of the intershaft bearing are taken into consideration. A simplified method for dynamic equations is developed due to the symmetry of rotors, and the bifurcation parameters of the dual-rotor system are also reduced because of the amplitude frequency equations solved by the multiple scales method. Afterwards, the numerical integral is performed to confirm the validity of the simplified method and the amplitude frequency equations. Furthermore, the bifurcation analysis for two interrelated primary resonances of the dual-rotor system could be carried out based on the twostate-variable singularity method. It shows the second primary resonance is more prone to have nonlinear dynamic characteristics than the first primary resonance.

Dynamic Equations of the Dual-Rotor System
2.1. Dual-Rotor System. A simple two-disk dual-rotor system supported by four points with an intershaft bearing [15,23] is depicted in Figure 1. e system is consisted of one LP rotor, one HP rotor, and an intershaft bearing. e LP rotor comprises a flexible shaft and a disk, while the HP rotor comprises a rigid rotor and a disk. e intershaft bearing, located between the LP rotor and HP rotor, is considered as a nonlinear elastic spring, which is the nonlinear source of the dual-rotor system. Herein, the subscript 1 denotes the LP rotor and the subscript 2 denotes the HP rotor without special instructions; k, c, ω, and O are the stiffness coefficient, the damping coefficient of the shaft, the rotation speed of the rotor, and the geometric center of the disk; l i (i � 1∼4) and l are the lengths of shafts. Figure 2 illustrates the lateral vibration of the dual-rotor system. e LP rotor is assumed as a flexible rotor while the HP rotor is assumed as a rigid rotor [24], since the LP shaft is much more slender than the HP shaft. Herein, x 1 and y 1 represent vertical and horizontal displacements of the LP disk; x 2 and y 2 represent vertical and horizontal displacements of the HP disk. e LP shaft and HP shaft are considered as the Euler-Bernoulli beam, in which the deformations of shear and torsion are ignored. According to the assumed mode method [25], the LP shaft is a flexible shaft; thus, the modal shape function is assumed as a sine function, as At any portion z 1 , the vertical displacement u 1 and horizontal displacement v 1 of the LP shaft can be represented by vertical displacement x 1 and horizontal displacement y 1 of the LP disk due to the modal shape function equation (1), as e rotational angles of the LP shaft around x-axis θ x and y-axis θ y are both very small, and thus, they can also be represented by x 1 and y 1 , as Higher pressure rotor Intershaft bearing Lower pressure rotor Figure 1: A simple two-disk dual-rotor system supported by four points with an intershaft bearing.
where m 1 is the mass of the LP rotor; J p1 and J d1 are the LP disk's polar and diameter moment of inertia; r 1 and ρ 1 are the radius and the density of the LP shaft (solid shaft); and the LP shaft's area moment of inertia is I 1 � (πr 4 1 /4). e HP shaft is a stubby shaft, which can be considered as a rigid rotor; thus, the modal shape function is At any portion z 2 , the vertical displacement u 2 and the horizontal displacement v 2 of the HP shaft can be represented by vertical displacement x 2 and horizontal displacement y 2 of the HP disk due to the modal shape function equation (5), as e rotational angles of HP shaft around x-axis φ x and yaxis φ y are both very small; thus, they can also be represented by x 2 and y 2 , as where h 2 (z 2 ) � (df 2 (z 2 )/f 2 (l 3 )dz 2 ) � (1/(l 3 − l 2 )). e HP disk is fixed with the HP shaft at z 2 � l 3 ; thus, displacements and rotational angles of HP disk are just equal to displacements and rotational angles of HP shaft at z 2 � l 3 . e kinetic energy of the HP rotor contains two parts, one is the kinetic energy of the HP disk T d2 and the other is the kinetic energy of the HP shaft T s2 , as x 2 x 1 y 1 y 2 Figure 2: Lateral vibration of the dual-rotor system.

Shock and Vibration 3
where m 2 is the mass of the HP rotor; J p2 and J d2 are the HP disk's polar and diameter moment of inertia; and r 2 , r 3 , and ρ 2 are the inner and the outer radiuses, and the density of the HP shaft (hollow shaft); the HP shaft's area moment of inertia is I 2 � (π(r 4 3 − r 4 2 )/4). us, the kinetic energy of the system can be obtained as e potential energy of the system is where k 1 and k 2 are the stiffness coefficients of LP shaft and HP shaft. Rayleigh's dissipation energy of the system is where c 1 and c 2 are the damping coefficients of LP shaft and HP shaft. e generalized force virtual work of the system is where ; F x and F y represent vertical and horizontal nonlinear restoring forces of the intershaft bearing; e 1 and e 2 represent unbalances of LP and HP rotors; and δ• (x 1 , y 1 , x 2 , y 2 ) represents the corresponding virtual displacement in the corresponding generalized coordinate.

Nonlinear Restoring Forces of the Intershaft Bearing.
e inner race and the LP rotor are interference fit at z 1 � l 4 ; thus, vertical displacement x i and horizontal displacement y i of the inner race are just equal to vertical displacement and horizontal displacement of the LP shaft at z 1 � l 4 ; substituting z 1 � l 4 into equation (2), we obtain e outer race and the HP rotor are interference fit at z 2 � l 4 ; thus, vertical displacement x o and horizontal displacement y o of the outer race are just equal to vertical displacement and horizontal displacement of the HP shaft at z 2 � l 4 ; substituting z 2 � l 4 into equation (6), we obtain e intershaft bearing can be seen as a nonlinear elastic spring [26,27], which is consisted of a linear stiffness k s and a cubic nonlinear stiffness k n . erefore, vertical and horizontal nonlinear restoring forces of the intershaft bearing are where k s denotes the linear stiffness of the intershaft bearing while k n denotes the cubic nonlinear stiffness of the intershaft bearing.
Comparing equations (15a) and (15b) with conventional cubic nonlinear elastic spring, the symmetry of the intershaft bearing in vertical and horizontal directions is taken into consideration for the nonlinear restoring forces of the intershaft bearing.

Dynamic Equations of the Dual-Rotor System.
e 4DOF dynamic equations of the dual-rotor system can be obtained with the aid of Lagrange's equation of the second kind, as where e dimensionless dynamic equations of the dual-rotor system are are the dimensionless parameters. e dimensionless nonlinear restoring force of the intershaft bearing is e parameters of the dual-rotor system [23]

Dynamic Equations Transferring.
Considering the symmetry of LP and HP rotors in vertical and horizontal directions, the 4DOF equations (17a)-(17d) can be simplified into a 2DOF equation with the help of a complexnumber coordinate [28]. Let e nonlinear force of the intershaft bearing in the complex-number coordinate (deduction process in Appendix) is Substitute the nonlinear force equation (21) into equations (20a) and (20b) and transpose terms as Shock and Vibration (22b) Equations (22a) and (22b) in a matrix form is where Eigenvectors and eigenvalues of the matrix K are where eig(·) is a function to compute eigenvectors and eigenvalues of the matrix. Let eigenvalues be Ω � (1/ω 2 1 ) where Ω 1 and Ω 2 are the first-and second-order critical speedof the system. And let eigenvectors be Φ � (1/(p 1 − p 2 )) and q 2 � (p 1 /(p 1 − p 2 )) . Let both sides of equation (23) be premultiplied by Φ − 1 and then expanded into a 2DoF equation as where R 1 � q 1 r 1 + (1 − q 2 )r 2 and R 2 � − q 1 r 1 + q 2 r 2 . ere are some assumptions as follows: (1) e damping of the system is weak damping, and the nonlinearity of the system is weak nonlinearity (2) Ω 2 ≫ Ω 1 (Ω 2 � ηΩ 1 , η ≈ 5.10), i.e., second-order critical speed is far greater than first-order critical speed. is paper focuses on primary resonances, i.e., first-order critical speed e dual-rotor system is excited by double frequency excitations, i.e., the unbalances of LP rotor and HP rotor. erefore, first-order critical speed contains two interrelated primary resonances, one of which is excited by the unbalance of the HP rotor and the other is excited by the unbalance of the LP rotor [23]. When the dual-rotor system rotates at a constant rotation speed ratio λ � (ω 2 /ω 1 ), in general, λ > 1 because the HP rotor rotates faster than the LP rotor. us, with the increase of the rotation speed, the first primary resonance is excited by the unbalance of the HP rotor at ω 1 ≈ (Ω 1 /λ)(ω 2 ≈ Ω 1 ) while the second primary resonance is excited by the unbalance of the LP rotor at ω 1 ≈ Ω 1 .

Stability Analysis.
e Floquet theory [29] is applied to confirm the stability of results solved by the multiple scales method, and the Floquet index can be figured out according to Hsu's method [30,31]. F(X(τ), τ), and X * denotes equilibriums of equations (25a) and (25b). Take a perturbation ΔX to X * , as (42) e stability of X * can be confirmed according to the Jacobi matrix A(τ, X * (τ)).
e approximate monodromy matrix can be obtained according to Hsu's method, as where I represents the identity matrix, T denotes the time period of the solution, which is divided into N n subintervals, ΔT is the length of one subinterval, A n is a constant matrix which represents the time-varying matrix A(τ, X * (τ)), and N j represents the number of terms in the approximation of the constant matrix A n exponential. e Floquet index of matrix M is available by calculating eigenvalues, and the stability of X * can be confirmed.

Numerical Verification.
In order to verify the validity of the simplified method for dynamic equations and the amplitude frequency equations solved by the multiple scales method, the multiple scales solutions of the simplified 2DOF equations (22a) and (22b) are compared with Runge-Kutta solutions of original 4DOF equations (17a) and (17b) around the first-order critical speed as shown in Figure 3. e ode45 function in MATLAB (Runge-Kutta method) is applied to solve original 4DOF equations (17a) and (17b). e initial states of motion for the first rotation speed ω 1 � 600 rad/s are set as where j is the number of points in one period and X j , Y j are the average values of corresponding response.
In Figure 3, it can be seen that the amplitude frequency curve of the LP rotor is basically the same as that of the HP rotor in appearance, but the amplitude of the LP rotor is much greater than that of the HP rotor, which is because the LP rotor is much more slender than the HP rotor.
ere are two primary resonances A and B caused by double frequency excitations, which are controlled by amplitude frequency equations equations (33) and (41). Furthermore, the multiple scales solutions almost coincide with the Runge-Kutta solutions, which verify the validity of the simplified method for dynamic equations proposed in this paper and provide a theoretical fundament for the multiple scales method or other analytical methods to solve the nonlinear dynamic equations of the dual-rotor system. e Runge-Kutta solutions are slightly larger than the multiple scales solutions, because dynamic responses solved by the multiple scales method only contain one frequency excitation which plays the most important role while others are ignored. Moreover, the unstable solutions can be attained by this method, while these cannot be attained by the Runge-Kutta method. It can save a lot of calculation time and cost. e results solved by this method can well reveal the inherent mechanism of some unexpected nonlinear phenomena, such as the vibration jump and the resonance hysteresis. All in all, the validity of the simplified method for dynamic equations and the amplitude frequency equations solved by the multiple scales method has been confirmed, based on which the bifurcation analysis can be carried out afterwards.

Bifurcation Analysis
e bifurcation characteristics for two interrelated primary resonances caused by double frequency excitations of the dual-rotor system are investigated in this section. Nevertheless, the local bifurcation characteristics for two primary resonances mostly depend on dynamic parameters of the dual-rotor and intershaft bearing systems, in which the linear stiffness and nonlinear stiffness of the intershaft bearing play the most important roles. erefore, this section focuses on the effect of the intershaft bearing's linear stiffness and nonlinear stiffness on the bifurcation characteristics for two interrelated primary resonances. Considering the complexity of the amplitude frequency equations, a two-state-variable singularity method [21,22,33] is applied to describe bifurcation modes for two primary resonances of the dual-rotor system in an engineering way. Amplitude frequency equations equations (33) and (41) are directly taken as the engineering unfolding, instead of the universal unfolding. e engineering unfolding may not be logical in mathematics, but it is of very important significance for engineering.
Substitute σ 1 � λ 2 − σ into equation (33) and σ 2 � 1 − σ into equation (41), and take σ � (Ω 2 1 /ω 2 1 ) as the bifurcation parameter and α 1 and α 4 as the unfolding parameters; the left side of equations (33) and (41) could be viewed as the engineering unfolding for two primary resonances of the dualrotor system, as 8 Shock and Vibration (44b) e derivatives of F 1 and F 2 can be expressed as Transition sets are consisted of three parts, which are bifurcation set, hysteresis set, and double limit set. e bifurcation set indicates that unfolding equations (44a) and (44b) contain a singular point. e hysteresis set indicates that the unfolding makes at least quadratic contact with the vertical plane of a constant bifurcation parameter. e double limit set indicates that the unfolding exists for two different limit points for the same bifurcation parameter [21]. According to the two-variable singularity method, transition sets for the bifurcation equations can be expressed as follows: Bifurcation set is as follows: Hysteresis set is as follows: Double limit set is as follows: Shock and Vibration 9 According to above calculation procedure equations (44a)-(46c), the transition sets for two primary resonances of the dual-rotor system can be attained. e transition sets contain two hysteresis set curves H 1 and H 2 , as shown in the k n − k s parameters plane in Figure 4 with k s from 0 to 15 (10 8 N/m) and k n from 0 to 150 (10 23 N/m 3 ). ere are three regions divided by H 1 and H 2 in the k n − k s parameters plane, which indicates that two primary resonances of the dual-rotor system have three different bifurcation modes.
ree different bifurcation modes for three regions I, II, and III can be obtained, as shown in Figures 5-7. Moreover, two critical bifurcation modes for two hysteresis set curves H 1 and H 2 can also be obtained, as shown in Figures 8 and 9. Figure 5 shows the bifurcation mode of region I in Figure 4, in which the linear stiffness k s is rather large and the nonlinear stiffness k n is rather small. It can found that two primary resonances A and B have a linear characteristics. e dynamic characteristics of the dual-rotor system are almost the same with the linear system in this case. Figure 6 shows the bifurcation mode of region II in Figure 4, in which the linear stiffness is slightly smaller or the nonlinear stiffness is slightly larger. It can be found that the second primary resonance B shows vibration jump and resonance hysteresis while the first primary resonance A does not, which means B shows the hardening characteristics while A does not. is phenomenon indicates that B is more prone to show nonlinear phenomena than A. Figure 7 shows the bifurcation mode of region III in Figure 4, in which the linear stiffness is rather small and the nonlinear stiffness is rather large. It can be found that both primary resonances A and B show vibration jump and resonance hysteresis, which indicates that both A and B show the hardening characteristics. e dynamic characteristics of the dual-rotor system have essential differences with the linear system in this case.
Furthermore, two critical bifurcation modes for two hysteresis set curves H 1 and H 2 are investigated. Figure 8 shows the bifurcation mode of hysteresis set curve H 1 , in   which the second primary resonance B just begins to show vibration jump in this case, while the first primary resonance A has linear characteristics. Figure 9 shows the bifurcation mode of hysteresis set curve H 2 , in which the first primary resonance A just begins to show vibration jump in this case, while the second primary resonance B has nonlinear characteristics.

Conclusions
In this paper, a simplified method for dynamic equations under double frequency excitations of a simple dual-rotor system considering the nonlinearity of the intershaft bearing has been proposed, based on which bifurcation modes for two interrelated primary resonances caused by double frequency excitations have been investigated theoretically. e multiple scales method has been applied to calculate amplitude frequency equations for two interrelated primary resonances, which have been taken as the engineering unfolding directly to investigate bifurcation modes of the dual-rotor system based on the two-state-variable singularity method. e results show that the linear stiffness and the nonlinear stiffness of the intershaft bearing have a significant influence on bifurcation modes of the system. e transition set, including two hysteresis set curves, divides the parameters plane k n − k s into three regions, indicating that there exists three different types of bifurcation modes for two interrelated primary resonances, which are, respectively, for both primary resonances having linear characteristics and for second primary resonance having hardening characteristics while first primary resonance having linear characteristics and for both primary resonances having hardening characteristics.
e future work will concentrate on the experimental study of the dual-rotor system. Appendix e nonlinear force of the intershaft bearing in the complexnumber coordinate is } is formulated as us, the nonlinear force of the intershaft bearing in the complex-number coordinate is Data Availability e dynamic simulation data used to support the findings of this study are included within the article.

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