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

Authors: Peng Gao; Yushu Chen; Lei Hou
Journal: Shock and Vibration
 (2020)

---

## Abstract

This paper aims to classify bifurcation modes for two interrelated primary resonances of a simple dual-rotor system under double frequency excitations. The 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. There 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. The second primary resonance is more prone to have nonlinear dynamic characteristics than the first primary resonance. This discovery indicates that two interrelated primary resonances of the dual-rotor system may have different bifurcation modes under the same dynamic parameters.


---

## Body

## 1. 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 dual-rotor system especially in aero-engines of fighters, the nonlinearity of the intershaft bearing [2] has achieved considerable attention [3]. There 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. Therefore, to clarify bifurcation modes of the dual-rotor 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.There 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 aero-engine. 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 dual-rotor 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 non-whole 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 blade-casing 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.The 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. This research applied the two-state-variable singularity method to investigate bifurcation modes for the nonlinear dual-rotor system.The 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 two-state-variable singularity method. It shows the second primary resonance is more prone to have nonlinear dynamic characteristics than the first primary resonance.


## 2. 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. The system is consisted of one LP rotor, one HP rotor, and an intershaft bearing. The LP rotor comprises a flexible shaft and a disk, while the HP rotor comprises a rigid rotor and a disk. The 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; li (i = 1∼4) and l are the lengths of shafts.Figure 1
A simple two-disk dual-rotor system supported by four points with an intershaft bearing.Figure2 illustrates the lateral vibration of the dual-rotor system. The 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, x1 and y1 represent vertical and horizontal displacements of the LP disk; x2 and y2 represent vertical and horizontal displacements of the HP disk.Figure 2
Lateral vibration of the dual-rotor system.The 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(1)f1z1=sinπz1l,z1∈0,l.At any portionz1, the vertical displacement u1 and horizontal displacement v1 of the LP shaft can be represented by vertical displacement x1 and horizontal displacement y1 of the LP disk due to the modal shape function equation (1), as(2)u1z1,t=f1z1f1l1x1t=f1z1f1l1x1,v1z1,t=f1z1f1l1y1t=f1z1f1l1y1.The rotational angles of the LP shaft aroundx-axis θx and y-axis θy are both very small, and thus, they can also be represented by x1 and y1, as(3)θxz1,t≈−∂v1z1,t∂z1=−df1z1f1l1dz1y1t=−h1z1y1,θyz1,t≈∂u1z1,t∂z1=df1z1f1l1dz1x1t=h1z1x1,where h1z1=df1z1/f1l1dz1.The LP disk is fixed with the LP shaft atz1=l1; thus, displacements and rotational angles of LP disk are just equal to displacements and rotational angles of LP shaft at z1=l1. The kinetic energy of the LP rotor contains two parts, one is the kinetic energy of the LP disk Td1 and the other is the kinetic energy of the LP shaft Ts1, as(4a)Td1=12m1u˙12l1+v˙12l1+12Jd1θ˙x2l1+θ˙y2l1+12Jp1ω12−2ω1θxl1θ˙yl1=12m1x˙12+y˙12+12Jd1h12l1x˙12+y˙12+12Jp1ω12+2ω1h12l1x˙1y1,(4b)Ts1=ρ1πr122∫0lu˙12+v˙12dz1+ρ1I12∫0lθ˙x2+θ˙y2dz1+ρ1I1∫0lω12−2ω1θxθ˙ydz1=ρ1πr122f12l1x˙12+y˙12∫0lf12dz1+ρ1I12x˙12+y˙12∫0lh12dz1+ρ1I1ω12l+2ω1x˙1y1∫0lh12dz1,where m1 is the mass of the LP rotor; Jp1 and Jd1 are the LP disk’s polar and diameter moment of inertia; r1 and ρ1 are the radius and the density of the LP shaft (solid shaft); and the LP shaft’s area moment of inertia is I1=πr14/4.The HP shaft is a stubby shaft, which can be considered as a rigid rotor; thus, the modal shape function is(5)f2z2=z2−l2,z2∈l2,l4.At any portionz2, the vertical displacement u2 and the horizontal displacement v2 of the HP shaft can be represented by vertical displacement x2 and horizontal displacement y2 of the HP disk due to the modal shape function equation (5), as(6)u2z2,t=f2z2f2l3x2t=f2z2f2l3x2,v2z2,t=f2z2f2l3y2t=f2z2f2l3y2.The rotational angles of HP shaft aroundx-axis φx and y-axis φy are both very small; thus, they can also be represented by x2 and y2, as(7)φxz2,t≈−∂v2z2,t∂z2=−df2z2f2l3dz2y2t=−h2z2y2,φyz2,t≈∂u2z2,t∂z2=df2z2f2l3dz2x2t=h2z2x2,where h2z2=df2z2/f2l3dz2=1/l3−l2.The HP disk is fixed with the HP shaft atz2=l3; thus, displacements and rotational angles of HP disk are just equal to displacements and rotational angles of HP shaft at z2=l3. The kinetic energy of the HP rotor contains two parts, one is the kinetic energy of the HP disk Td2 and the other is the kinetic energy of the HP shaft Ts2, as(8a)Td2=12m2u˙22l3+v˙22l3+12Jd2φ˙x2l3+φ˙y2l3+12Jp2ω22−2ω2φxl3φ˙yl3=12m2x˙22+y˙22+12Jd2h22l3x˙22+y˙22+12Jp2ω22+2ω2h22l3x˙2y2,(8b)Ts2=ρ2πr32−r222∫l2l4u˙22+v˙22dz2+ρ2I22∫l2l4φ˙x2+φ˙y2dz2+ρ2I2∫l2l4ω22−2ω2φxφ˙ydz2=ρ2πr32−r222f22l3x˙22+y˙22∫l2l4f22dz2+ρ2I22x˙22+y˙22∫l2l4h22dz2+ρ2I2ω22l4−l2+2ω2x˙2y2∫l2l4h22dz2,where m2 is the mass of the HP rotor; Jp2 and Jd2 are the HP disk’s polar and diameter moment of inertia; and r2, r3, 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 I2=πr34−r24/4.Thus, the kinetic energy of the system can be obtained as(9)T=Td1+Ts1+Td2+Ts2.The potential energy of the system is(10)V=12k1x12+y12+12k2x22+y22,where k1 and k2 are the stiffness coefficients of LP shaft and HP shaft.Rayleigh’s dissipation energy of the system is(11)D=12c1x˙12+y˙12+12c2x˙22+y˙22,where c1 and c2 are the damping coefficients of LP shaft and HP shaft.The generalized force virtual work of the system is(12)δW=m1ω12e1cosω1tδx1+m1ω12e1sinω1tδy1+m2ω22e2cosω2tδx2+m2ω22e2sinω2tδy2−FxF1δx1−F2δx2−FyF1δy1−F2δy2,where F1=f1l4/f1l1 and F2=f2l4/f2l3; Fx and Fy represent vertical and horizontal nonlinear restoring forces of the intershaft bearing; e1 and e2 represent unbalances of LP and HP rotors; and δ• (x1, y1, x2, y2) represents the corresponding virtual displacement in the corresponding generalized coordinate.


### 2.2. Nonlinear Restoring Forces of the Intershaft Bearing

The inner race and the LP rotor are interference fit atz1=l4; thus, vertical displacement xi and horizontal displacement yi of the inner race are just equal to vertical displacement and horizontal displacement of the LP shaft at z1=l4; substituting z1=l4 into equation (2), we obtain(13)xi=u1l4,t=f1l4f1l1x1=F1x1,yi=v1l4,t=f1l4f1l1y1=F1y1.The outer race and the HP rotor are interference fit atz2=l4; thus, vertical displacement xo and horizontal displacement yo of the outer race are just equal to vertical displacement and horizontal displacement of the HP shaft at z2=l4; substituting z2=l4 into equation (6), we obtain(14)xo=u2l4,t=f2l4f2l1x2=F2x2,yo=v2l4,t=f2l4f2l1y2=F2y2.The intershaft bearing can be seen as a nonlinear elastic spring [26, 27], which is consisted of a linear stiffness ks and a cubic nonlinear stiffness kn. Therefore, vertical and horizontal nonlinear restoring forces of the intershaft bearing are(15a)Fx=ksxi−xo+knxi−xoxi−xo2−3yi−yo2=ksF1x1−F2x2+knF1x1−F2x2F1x1−F2x22−3F1y1−F2y22,(15b)Fy=ksyi−yo+knyi−yoyi−yo2−3xi−xo2=ksF1y1−F2y2+knF1y1−F2y2F1y1−F2y22−3F1x1−F2x22,where ks denotes the linear stiffness of the intershaft bearing while kn 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.


### 2.3. Dynamic Equations of the Dual-Rotor System

The 4DOF dynamic equations of the dual-rotor system can be obtained with the aid of Lagrange’s equation of the second kind, as(16a)M1x¨1+c1x˙1+ω1g1y˙1+k1x1=m1ω12e1cosω1t−F1Fx,(16b)M1y¨1+c1y˙1−ω1g1x˙1+k1y1=m1ω12e1sinω1t−F1Fy,(16c)M2x¨2+c2x˙2+ω2g2y˙2+k2x2=m2ω22e2cosω2t+F2Fx,(16d)M2y¨2+c2y˙2−ω2g2x˙2+k2y2=m2ω22e2sinω2t+F2Fy,where M1=m1+H12Jd1+ρ1πr12/f12l1∫0lf12dz1+ρ1I1∫0lh12dz1, M2=m2+H22Jd2+ρ2πr32−r22/f22l3∫l2l4f22dz2+ρ2I2∫l2l4h22dz2, g1=H12Jp1+2ρ1I1∫0lh12dz1, and g2=H22Jp2+2ρ2I2∫l2l4h22dz2. Moreover, the gyroscopic effect is reflected in the parameters M1,M2,g1,andg2.The dimensionless dynamic equations of the dual-rotor system are(17a)X1″+C1X1′+G1Y1′+k1M1ω12X1=E1cosτ−F1M1ω12FX,(17b)Y1″+C1Y1′−G1X1′+k1M1ω12Y1=E1sinτ−F1M1ω12FY,(17c)X2″+C2X2′+G2Y2′+k2M2ω12X2=λ2E2cosλτ+F2M2ω12FX,(17d)Y2″+C2Y2′−G2X2′+k2M2ω12Y2=λ2E2sinλτ+F2M2ω12FY,where τ=ω1t, λ=ω2/ω1 (rotation speed ratio), E1=m1/M1, E2=m2/M2⋅e2/e1, X1=x1/e1, Y1=y1/e1, X2=x2/e1, Y2=y2/e1, C1=c1/M1ω1, C2=c2/M2ω1, G1=g1/M1, and G2=λg2/M2 are the dimensionless parameters.The dimensionless nonlinear restoring force of the intershaft bearing is(18)FXFY=Fxe1Fye1T=ksXi−Xo+e12knXi−XoXi−Xo2−3Yi−Yo2ksYi−Yo+e12knYi−YoYi−Yo2−3Xi−Xo2=ksF1X1−F2X2+e12knF1X1−F2X2F1X1−F2X22−3F1Y1−F2Y22ksF1Y1−F2Y2+e12knF1Y1−F2Y2F1Y1−F2Y22−3F1X1−F2X22.The parameters of the dual-rotor system [23] are shown as follows:(19)m1=97.37kg,Jp1=3.6907kg⋅m2,Jd1=1.8454kg⋅m2,m2=108.3kg,Jp2=4.0114kg⋅m2,Jd2=2.0060kg⋅m2,k1=k2=8×107N/m,c1=c2=600N⋅s/m,ρ1=ρ2=7800kg/m3,e1=1μm,e2=3μm,l=1.5913m,l1=0.2250m,l2=0.3555m,l3=0.8675m,l4=1.4918m,r1=20mm,r2=30mm,r3=35mm,λ=1.2,ks=3×108N/m,kn=1×1023N/m3.


## 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. The system is consisted of one LP rotor, one HP rotor, and an intershaft bearing. The LP rotor comprises a flexible shaft and a disk, while the HP rotor comprises a rigid rotor and a disk. The 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; li (i = 1∼4) and l are the lengths of shafts.Figure 1
A simple two-disk dual-rotor system supported by four points with an intershaft bearing.Figure2 illustrates the lateral vibration of the dual-rotor system. The 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, x1 and y1 represent vertical and horizontal displacements of the LP disk; x2 and y2 represent vertical and horizontal displacements of the HP disk.Figure 2
Lateral vibration of the dual-rotor system.The 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(1)f1z1=sinπz1l,z1∈0,l.At any portionz1, the vertical displacement u1 and horizontal displacement v1 of the LP shaft can be represented by vertical displacement x1 and horizontal displacement y1 of the LP disk due to the modal shape function equation (1), as(2)u1z1,t=f1z1f1l1x1t=f1z1f1l1x1,v1z1,t=f1z1f1l1y1t=f1z1f1l1y1.The rotational angles of the LP shaft aroundx-axis θx and y-axis θy are both very small, and thus, they can also be represented by x1 and y1, as(3)θxz1,t≈−∂v1z1,t∂z1=−df1z1f1l1dz1y1t=−h1z1y1,θyz1,t≈∂u1z1,t∂z1=df1z1f1l1dz1x1t=h1z1x1,where h1z1=df1z1/f1l1dz1.The LP disk is fixed with the LP shaft atz1=l1; thus, displacements and rotational angles of LP disk are just equal to displacements and rotational angles of LP shaft at z1=l1. The kinetic energy of the LP rotor contains two parts, one is the kinetic energy of the LP disk Td1 and the other is the kinetic energy of the LP shaft Ts1, as(4a)Td1=12m1u˙12l1+v˙12l1+12Jd1θ˙x2l1+θ˙y2l1+12Jp1ω12−2ω1θxl1θ˙yl1=12m1x˙12+y˙12+12Jd1h12l1x˙12+y˙12+12Jp1ω12+2ω1h12l1x˙1y1,(4b)Ts1=ρ1πr122∫0lu˙12+v˙12dz1+ρ1I12∫0lθ˙x2+θ˙y2dz1+ρ1I1∫0lω12−2ω1θxθ˙ydz1=ρ1πr122f12l1x˙12+y˙12∫0lf12dz1+ρ1I12x˙12+y˙12∫0lh12dz1+ρ1I1ω12l+2ω1x˙1y1∫0lh12dz1,where m1 is the mass of the LP rotor; Jp1 and Jd1 are the LP disk’s polar and diameter moment of inertia; r1 and ρ1 are the radius and the density of the LP shaft (solid shaft); and the LP shaft’s area moment of inertia is I1=πr14/4.The HP shaft is a stubby shaft, which can be considered as a rigid rotor; thus, the modal shape function is(5)f2z2=z2−l2,z2∈l2,l4.At any portionz2, the vertical displacement u2 and the horizontal displacement v2 of the HP shaft can be represented by vertical displacement x2 and horizontal displacement y2 of the HP disk due to the modal shape function equation (5), as(6)u2z2,t=f2z2f2l3x2t=f2z2f2l3x2,v2z2,t=f2z2f2l3y2t=f2z2f2l3y2.The rotational angles of HP shaft aroundx-axis φx and y-axis φy are both very small; thus, they can also be represented by x2 and y2, as(7)φxz2,t≈−∂v2z2,t∂z2=−df2z2f2l3dz2y2t=−h2z2y2,φyz2,t≈∂u2z2,t∂z2=df2z2f2l3dz2x2t=h2z2x2,where h2z2=df2z2/f2l3dz2=1/l3−l2.The HP disk is fixed with the HP shaft atz2=l3; thus, displacements and rotational angles of HP disk are just equal to displacements and rotational angles of HP shaft at z2=l3. The kinetic energy of the HP rotor contains two parts, one is the kinetic energy of the HP disk Td2 and the other is the kinetic energy of the HP shaft Ts2, as(8a)Td2=12m2u˙22l3+v˙22l3+12Jd2φ˙x2l3+φ˙y2l3+12Jp2ω22−2ω2φxl3φ˙yl3=12m2x˙22+y˙22+12Jd2h22l3x˙22+y˙22+12Jp2ω22+2ω2h22l3x˙2y2,(8b)Ts2=ρ2πr32−r222∫l2l4u˙22+v˙22dz2+ρ2I22∫l2l4φ˙x2+φ˙y2dz2+ρ2I2∫l2l4ω22−2ω2φxφ˙ydz2=ρ2πr32−r222f22l3x˙22+y˙22∫l2l4f22dz2+ρ2I22x˙22+y˙22∫l2l4h22dz2+ρ2I2ω22l4−l2+2ω2x˙2y2∫l2l4h22dz2,where m2 is the mass of the HP rotor; Jp2 and Jd2 are the HP disk’s polar and diameter moment of inertia; and r2, r3, 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 I2=πr34−r24/4.Thus, the kinetic energy of the system can be obtained as(9)T=Td1+Ts1+Td2+Ts2.The potential energy of the system is(10)V=12k1x12+y12+12k2x22+y22,where k1 and k2 are the stiffness coefficients of LP shaft and HP shaft.Rayleigh’s dissipation energy of the system is(11)D=12c1x˙12+y˙12+12c2x˙22+y˙22,where c1 and c2 are the damping coefficients of LP shaft and HP shaft.The generalized force virtual work of the system is(12)δW=m1ω12e1cosω1tδx1+m1ω12e1sinω1tδy1+m2ω22e2cosω2tδx2+m2ω22e2sinω2tδy2−FxF1δx1−F2δx2−FyF1δy1−F2δy2,where F1=f1l4/f1l1 and F2=f2l4/f2l3; Fx and Fy represent vertical and horizontal nonlinear restoring forces of the intershaft bearing; e1 and e2 represent unbalances of LP and HP rotors; and δ• (x1, y1, x2, y2) represents the corresponding virtual displacement in the corresponding generalized coordinate.


## 2.2. Nonlinear Restoring Forces of the Intershaft Bearing

The inner race and the LP rotor are interference fit atz1=l4; thus, vertical displacement xi and horizontal displacement yi of the inner race are just equal to vertical displacement and horizontal displacement of the LP shaft at z1=l4; substituting z1=l4 into equation (2), we obtain(13)xi=u1l4,t=f1l4f1l1x1=F1x1,yi=v1l4,t=f1l4f1l1y1=F1y1.The outer race and the HP rotor are interference fit atz2=l4; thus, vertical displacement xo and horizontal displacement yo of the outer race are just equal to vertical displacement and horizontal displacement of the HP shaft at z2=l4; substituting z2=l4 into equation (6), we obtain(14)xo=u2l4,t=f2l4f2l1x2=F2x2,yo=v2l4,t=f2l4f2l1y2=F2y2.The intershaft bearing can be seen as a nonlinear elastic spring [26, 27], which is consisted of a linear stiffness ks and a cubic nonlinear stiffness kn. Therefore, vertical and horizontal nonlinear restoring forces of the intershaft bearing are(15a)Fx=ksxi−xo+knxi−xoxi−xo2−3yi−yo2=ksF1x1−F2x2+knF1x1−F2x2F1x1−F2x22−3F1y1−F2y22,(15b)Fy=ksyi−yo+knyi−yoyi−yo2−3xi−xo2=ksF1y1−F2y2+knF1y1−F2y2F1y1−F2y22−3F1x1−F2x22,where ks denotes the linear stiffness of the intershaft bearing while kn 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.


## 2.3. Dynamic Equations of the Dual-Rotor System

The 4DOF dynamic equations of the dual-rotor system can be obtained with the aid of Lagrange’s equation of the second kind, as(16a)M1x¨1+c1x˙1+ω1g1y˙1+k1x1=m1ω12e1cosω1t−F1Fx,(16b)M1y¨1+c1y˙1−ω1g1x˙1+k1y1=m1ω12e1sinω1t−F1Fy,(16c)M2x¨2+c2x˙2+ω2g2y˙2+k2x2=m2ω22e2cosω2t+F2Fx,(16d)M2y¨2+c2y˙2−ω2g2x˙2+k2y2=m2ω22e2sinω2t+F2Fy,where M1=m1+H12Jd1+ρ1πr12/f12l1∫0lf12dz1+ρ1I1∫0lh12dz1, M2=m2+H22Jd2+ρ2πr32−r22/f22l3∫l2l4f22dz2+ρ2I2∫l2l4h22dz2, g1=H12Jp1+2ρ1I1∫0lh12dz1, and g2=H22Jp2+2ρ2I2∫l2l4h22dz2. Moreover, the gyroscopic effect is reflected in the parameters M1,M2,g1,andg2.The dimensionless dynamic equations of the dual-rotor system are(17a)X1″+C1X1′+G1Y1′+k1M1ω12X1=E1cosτ−F1M1ω12FX,(17b)Y1″+C1Y1′−G1X1′+k1M1ω12Y1=E1sinτ−F1M1ω12FY,(17c)X2″+C2X2′+G2Y2′+k2M2ω12X2=λ2E2cosλτ+F2M2ω12FX,(17d)Y2″+C2Y2′−G2X2′+k2M2ω12Y2=λ2E2sinλτ+F2M2ω12FY,where τ=ω1t, λ=ω2/ω1 (rotation speed ratio), E1=m1/M1, E2=m2/M2⋅e2/e1, X1=x1/e1, Y1=y1/e1, X2=x2/e1, Y2=y2/e1, C1=c1/M1ω1, C2=c2/M2ω1, G1=g1/M1, and G2=λg2/M2 are the dimensionless parameters.The dimensionless nonlinear restoring force of the intershaft bearing is(18)FXFY=Fxe1Fye1T=ksXi−Xo+e12knXi−XoXi−Xo2−3Yi−Yo2ksYi−Yo+e12knYi−YoYi−Yo2−3Xi−Xo2=ksF1X1−F2X2+e12knF1X1−F2X2F1X1−F2X22−3F1Y1−F2Y22ksF1Y1−F2Y2+e12knF1Y1−F2Y2F1Y1−F2Y22−3F1X1−F2X22.The parameters of the dual-rotor system [23] are shown as follows:(19)m1=97.37kg,Jp1=3.6907kg⋅m2,Jd1=1.8454kg⋅m2,m2=108.3kg,Jp2=4.0114kg⋅m2,Jd2=2.0060kg⋅m2,k1=k2=8×107N/m,c1=c2=600N⋅s/m,ρ1=ρ2=7800kg/m3,e1=1μm,e2=3μm,l=1.5913m,l1=0.2250m,l2=0.3555m,l3=0.8675m,l4=1.4918m,r1=20mm,r2=30mm,r3=35mm,λ=1.2,ks=3×108N/m,kn=1×1023N/m3.


## 3. Multiple Scales Solution

### 3.1. 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 complex-number coordinate [28]. Let r1=X1+iY1, r2=X2+iY2, and Fb=FX+iFY (i2=−1); equations (17a)–(17d) are transformed into(20a)r1″+C1−iG1r1′+k1M1ω12r1=E1eiτ−F1M1ω12Fb,(20b)r2″+C2−iG2r2′+k2M2ω12r2=λ2E2eiλτ+F2M2ω12Fb.The nonlinear force of the intershaft bearing in the complex-number coordinate (deduction process in Appendix) is(21)Fb=FX+iFY=ksri−ro+e12knr¯i−r¯o3=ksF1r1−F2r2+e12knF1r¯1−F2r¯23.Substitute the nonlinear force equation (21) into equations (20a) and (20b) and transpose terms as(22a)r1″+k1+F12ksM1ω12r1−F1F2ksM1ω12r2=E1eiτ−C1−iG1r1′−F1e12knM1ω12F1r¯1−F2r¯23,(22b)r2″−F1F2ksM2ω12r1+k2+F22ksM2ω12r2=λ2E2eiλτ−C2−iG2r2′+F2e12knM2ω12F1r¯1−F2r¯23.Equations (22a) and (22b) in a matrix form is(23)r″+Kr=Q,where r=r1r2, K=1/ω12k1+F12ks/M1−F1F2ks/M1−F1F2ks/M2k2+F22ks/M2, and Q=Q1Q2=E1eiτ−C1−iG1r1′−F1e12kn/M1ω12F1r¯1−F2r¯23λ2E2eiλτ−C2−iG2r2′+F2e12kn/M2ω12F1r¯1−F2r¯23.Eigenvectors and eigenvalues of the matrixK are(24)Φ,Ω=eigK,where eig⋅ is a function to compute eigenvectors and eigenvalues of the matrix.Let eigenvalues beΩ=1/ω12Ω1200Ω22, where Ω1 and Ω2 are the first- and second-order critical speedof the system. And let eigenvectors be Φ=p1p211; the inverse matrix of eigenvectors is Φ−1=q11−q2−q1q2, where q1=1/p1−p2 and q2=p1/p1−p2 .Let both sides of equation (23) be premultiplied by Φ−1 and then expanded into a 2DoF equation as(25a)R1″+Ω12ω12R1=q1Q1+1−q2Q2,(25b)R2″+Ω22ω12R2=−q1Q1+q2Q2,where R1=q1r1+1−q2r2 and R2=−q1r1+q2r2.There are some assumptions as follows:(1)
The 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. This paper focuses on primary resonances, i.e., first-order critical speedThe dual-rotor system is excited by double frequency excitations, i.e., the unbalances of LP rotor and HP rotor. Therefore, 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. Thus, 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.


### 3.2. First Primary Resonance atω1≈Ω1/λω2≈Ω1

Time scalesT1=τ and T2=ετ, and the tuning parameter σ1 satisfies(26)Ω12ω12=λ2−εσ1.Assume that solutions of equations (25a) and (25b) are(27a)R1=R11+εR12,(27b)R2=R21+εR22.Substitute equations (26)–(27b) into equations (25a) and (25b) and equate coefficients of ε0 and ε1 to zero, as ε0 is(28a)D12R11+λ2R11=q1E1eiT1,(28b)D12R21+λη2R21=−q1E1eiT1.Also, ε1 is(29a)D12R12+λ2R12=−2D1D2R11+σ1R11+1−q2λ2E2eiλT1−q1C1−iG1p1D1R11+p2D1R21−1−q2C2−iG2D1R11+D1R21+γ1p1F1−F2R¯11+p2F1−F2R¯213,(29b)D12R22+λη2R22=−2D1D2R21+σ1η2R21−q2λ2E2eiλT1+q1C1−iG1p1D1R11+p2D1R21−q2C2−iG2D1R11+D1R21+γ2p1F1−F2R¯11+p2F1−F2R¯213,where D1=∂⋅/∂T1, D2=∂⋅/∂T2, γ1=1−q2F2/M2ω12−q1F1/M1ω12e12kn, γ2=q2F2/M2ω12+q1F1/M1ω12e12kn, and η=Ω2/Ω1.Assume that general solutions of equations (28a) and (28b) are(30a)R11=A1eiλT1+φ1+cc+a1eiT1,(30b)R21=A2eiληT1+φ2+cc+a2eiT1,where a1=q1E1/λ2−1 and a2=−q1E1/λη2−1.Substitute equations (30a) and (30b) into equations (29a) and (29b), and then, eliminate secular terms and separate real and imaginary parts, as(31a)D2A1=−12λλα4A1+λ21−q2E2sinφ1,(31b)A1D2φ1=−12λα1A13+α2A1A22+λα3+σ1A1+λ21−q2E2cosφ1,(31c)D2A2=−12β4A2,(31d)A2D2φ2=−12ληA2β1A22+β2A12+λβ3+ησ1η,where α1=3γ1p1F1−F23, α2=6γ1p1F1−F2p2F1−F22, α3=−q1p1G1−1−q2G2, α4=q1p1C1+1−q2C2, β1=3γ2p2F1−F23, β2=6γ2p1F1−F22p2F1−F2, β3=q1p2G1−q2G2, and β4=−q1p2C1+q2C2.Equate right sides of equations (31a)–(31d) to zero; we obtain A2=0 due to equation (31c). Since A2 is the amplitude of the harmonic term with frequency Ω2 at ω1≈Ω1/λ, A2=0 is correct when Ω2≫Ω1 (assumption 2). Therefore, equations (31a)–(31d) are simplified as(32a)D2A1=−12λλα4A1+λ21−q2E2sinφ1,(32b)A1D2φ1=−12λα1A13+λα3+σ1A1+λ21−q2E2cosφ1.Finally, eliminateφ1 in equations (32a) and (32b); the amplitude frequency equation for the first primary resonance at ω1≈Ω1/λ is(33)α12A3+2α1λα3+σ1A2+λα3+σ12+λ2α42A−λ41−q22E22=0,where A=A12.


### 3.3. Second Primary Resonance atω1≈Ω1

Time scalesT1=τ and T2=ετ, and the tuning parameter σ2 satisfies(34)Ω12ω12=1−εσ2.Assume that solutions of equations (25a) and (25b) are(35a)R1=R11+εR12,(35b)R2=R21+εR22.Substitute equations (34)–(35b) into equations (25a) and (25b) and equate coefficients of ε0 and ε1 to zero, as ε0 is(36a)D12R11+R11=1−q2λ2E2eiλT1,(36b)D12R21+η2R21=q2λ2E2eiλT1.Also ε1 is(37a)D12R12+R12=−2D1D2R11+σ2R11+q1E1eiT1−q1C1−iG1p1D1R11+p2D1R21−1−q2C2−iG2D1R11+D1R21+γ1p1F1−F2R¯11+p2F1−F2R¯213,(37b)D12R22+η2R22=−2D1D2R21+σ2η2R21−q1E1eiT1+q1C1−iG1p1D1R11+p2D1R21−q2C2−iG2D1R11+D1R21+γ2p1F1−F2R¯11+p2F1−F2R¯213.Assume that general solutions of equations (36a) and (36b) are(38a)R11=B1eiT1+φ1+cc+b1eiλT1,(38b)R21=B2eiηT1+φ2+cc+b2eiλT1,where b1=1−q2λ2E2/1−λ2 and b2=q2λ2E2/η2−λ2.Substitute equations (38a) and (38b) into equations (37a) and (37b), and then, eliminate secular terms and separate real and imaginary parts, as(39a)D2B1=−12α4B1+q1E1sinφ2,(39b)B1D2φ1=−12α1B13+α2B1B22+α3+σ2B1+q1E1cosφ2,(39c)D2B2=−12β4B2,(39d)B2D2φ2=−12ηB2β1B22+β2B12+β3+ησ2η.Equate right sides of equations (39a)–(39d) to zero; we obtain B2=0 due to equation (39c). Since B2 is the amplitude of the harmonic term with frequency Ω2 at ω1≈Ω1, B2=0 is correct when Ω2≫Ω1 (assumption 2). Therefore, equations (39a)–(39d) are simplified as(40a)D2B1=−12α4B1+q1E1sinφ2,(40b)B1D2φ2=−12α1B13+α3+σ2B1+q1E1cosφ2.Finally eliminatingφ2 in equations (40a) and (40b), the amplitude frequency equation for the second primary resonance at ω1≈Ω1 is(41)α12B3+2α1α3+σ2B2+α3+σ22+α42B−E12q12=0,where B=B12.


### 3.4. Stability Analysis

The 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].LetX=R1,R2,R1′,R2′T=X1,X2,X3,X4T; transform equations (25a) and (25b) into X′τ=FXτ,τ, and X∗ denotes equilibriums of equations (25a) and (25b). Take a perturbation ΔX to X∗, as(42)ΔX′≈∂Fτ,X∗∂X∗ΔX=Aτ,X∗τΔX.The stability ofX∗ can be confirmed according to the Jacobi matrix Aτ,X∗τ.The approximate monodromy matrix can be obtained according to Hsu’s method, as(43)M=ΨT=∏n=1NneAnΔT=∏n=1NnI+∑j=1NjAnΔTjj!,where I represents the identity matrix, T denotes the time period of the solution, which is divided into Nn subintervals, ΔT is the length of one subinterval, An is a constant matrix which represents the time-varying matrix Aτ,X∗τ, and Nj represents the number of terms in the approximation of the constant matrix An exponential. The Floquet index of matrix M is available by calculating eigenvalues, and the stability of X∗ can be confirmed.


### 3.5. 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. The ode45 function in MATLAB (Runge–Kutta method) is applied to solve original 4DOF equations (17a) and (17b). The initial states of motion for the first rotation speed ω1=600rad/s are set as X1X1′Y1Y1′X2X2′Y2Y2′=00000000, but for other rotation speeds, the initial states of motion are set as the final states of motion for the previous rotation speed. Integration duration should be long enough to ensure steady-state responses; for equations (17a) and (17b), integration duration are set as τ=0,2000π. The root mean square (RMS) [32] is applied to represent the amplitude of the vibration solved by the Runge–Kutta method. The amplitudes of LP and HP rotor are rj=∑i=1N^Xji−X¯j2+Yji−Y¯j2/N,j=1,2, where j is the number of points in one period and X¯j, Y¯j are the average values of corresponding response.Figure 3
Amplitude frequency curves solved by the multiple scales method and the Runge–Kutta method. (a) For the LP rotor. (b) For the HP rotor. Solid lines represent stable solutions, and dotted lines represent unstable solutions.
(a)(b)In Figure3, 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. There 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. The 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. The 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.


## 3.1. 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 complex-number coordinate [28]. Let r1=X1+iY1, r2=X2+iY2, and Fb=FX+iFY (i2=−1); equations (17a)–(17d) are transformed into(20a)r1″+C1−iG1r1′+k1M1ω12r1=E1eiτ−F1M1ω12Fb,(20b)r2″+C2−iG2r2′+k2M2ω12r2=λ2E2eiλτ+F2M2ω12Fb.The nonlinear force of the intershaft bearing in the complex-number coordinate (deduction process in Appendix) is(21)Fb=FX+iFY=ksri−ro+e12knr¯i−r¯o3=ksF1r1−F2r2+e12knF1r¯1−F2r¯23.Substitute the nonlinear force equation (21) into equations (20a) and (20b) and transpose terms as(22a)r1″+k1+F12ksM1ω12r1−F1F2ksM1ω12r2=E1eiτ−C1−iG1r1′−F1e12knM1ω12F1r¯1−F2r¯23,(22b)r2″−F1F2ksM2ω12r1+k2+F22ksM2ω12r2=λ2E2eiλτ−C2−iG2r2′+F2e12knM2ω12F1r¯1−F2r¯23.Equations (22a) and (22b) in a matrix form is(23)r″+Kr=Q,where r=r1r2, K=1/ω12k1+F12ks/M1−F1F2ks/M1−F1F2ks/M2k2+F22ks/M2, and Q=Q1Q2=E1eiτ−C1−iG1r1′−F1e12kn/M1ω12F1r¯1−F2r¯23λ2E2eiλτ−C2−iG2r2′+F2e12kn/M2ω12F1r¯1−F2r¯23.Eigenvectors and eigenvalues of the matrixK are(24)Φ,Ω=eigK,where eig⋅ is a function to compute eigenvectors and eigenvalues of the matrix.Let eigenvalues beΩ=1/ω12Ω1200Ω22, where Ω1 and Ω2 are the first- and second-order critical speedof the system. And let eigenvectors be Φ=p1p211; the inverse matrix of eigenvectors is Φ−1=q11−q2−q1q2, where q1=1/p1−p2 and q2=p1/p1−p2 .Let both sides of equation (23) be premultiplied by Φ−1 and then expanded into a 2DoF equation as(25a)R1″+Ω12ω12R1=q1Q1+1−q2Q2,(25b)R2″+Ω22ω12R2=−q1Q1+q2Q2,where R1=q1r1+1−q2r2 and R2=−q1r1+q2r2.There are some assumptions as follows:(1)
The 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. This paper focuses on primary resonances, i.e., first-order critical speedThe dual-rotor system is excited by double frequency excitations, i.e., the unbalances of LP rotor and HP rotor. Therefore, 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. Thus, 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.


## 3.2. First Primary Resonance atω1≈Ω1/λω2≈Ω1

Time scalesT1=τ and T2=ετ, and the tuning parameter σ1 satisfies(26)Ω12ω12=λ2−εσ1.Assume that solutions of equations (25a) and (25b) are(27a)R1=R11+εR12,(27b)R2=R21+εR22.Substitute equations (26)–(27b) into equations (25a) and (25b) and equate coefficients of ε0 and ε1 to zero, as ε0 is(28a)D12R11+λ2R11=q1E1eiT1,(28b)D12R21+λη2R21=−q1E1eiT1.Also, ε1 is(29a)D12R12+λ2R12=−2D1D2R11+σ1R11+1−q2λ2E2eiλT1−q1C1−iG1p1D1R11+p2D1R21−1−q2C2−iG2D1R11+D1R21+γ1p1F1−F2R¯11+p2F1−F2R¯213,(29b)D12R22+λη2R22=−2D1D2R21+σ1η2R21−q2λ2E2eiλT1+q1C1−iG1p1D1R11+p2D1R21−q2C2−iG2D1R11+D1R21+γ2p1F1−F2R¯11+p2F1−F2R¯213,where D1=∂⋅/∂T1, D2=∂⋅/∂T2, γ1=1−q2F2/M2ω12−q1F1/M1ω12e12kn, γ2=q2F2/M2ω12+q1F1/M1ω12e12kn, and η=Ω2/Ω1.Assume that general solutions of equations (28a) and (28b) are(30a)R11=A1eiλT1+φ1+cc+a1eiT1,(30b)R21=A2eiληT1+φ2+cc+a2eiT1,where a1=q1E1/λ2−1 and a2=−q1E1/λη2−1.Substitute equations (30a) and (30b) into equations (29a) and (29b), and then, eliminate secular terms and separate real and imaginary parts, as(31a)D2A1=−12λλα4A1+λ21−q2E2sinφ1,(31b)A1D2φ1=−12λα1A13+α2A1A22+λα3+σ1A1+λ21−q2E2cosφ1,(31c)D2A2=−12β4A2,(31d)A2D2φ2=−12ληA2β1A22+β2A12+λβ3+ησ1η,where α1=3γ1p1F1−F23, α2=6γ1p1F1−F2p2F1−F22, α3=−q1p1G1−1−q2G2, α4=q1p1C1+1−q2C2, β1=3γ2p2F1−F23, β2=6γ2p1F1−F22p2F1−F2, β3=q1p2G1−q2G2, and β4=−q1p2C1+q2C2.Equate right sides of equations (31a)–(31d) to zero; we obtain A2=0 due to equation (31c). Since A2 is the amplitude of the harmonic term with frequency Ω2 at ω1≈Ω1/λ, A2=0 is correct when Ω2≫Ω1 (assumption 2). Therefore, equations (31a)–(31d) are simplified as(32a)D2A1=−12λλα4A1+λ21−q2E2sinφ1,(32b)A1D2φ1=−12λα1A13+λα3+σ1A1+λ21−q2E2cosφ1.Finally, eliminateφ1 in equations (32a) and (32b); the amplitude frequency equation for the first primary resonance at ω1≈Ω1/λ is(33)α12A3+2α1λα3+σ1A2+λα3+σ12+λ2α42A−λ41−q22E22=0,where A=A12.


## 3.3. Second Primary Resonance atω1≈Ω1

Time scalesT1=τ and T2=ετ, and the tuning parameter σ2 satisfies(34)Ω12ω12=1−εσ2.Assume that solutions of equations (25a) and (25b) are(35a)R1=R11+εR12,(35b)R2=R21+εR22.Substitute equations (34)–(35b) into equations (25a) and (25b) and equate coefficients of ε0 and ε1 to zero, as ε0 is(36a)D12R11+R11=1−q2λ2E2eiλT1,(36b)D12R21+η2R21=q2λ2E2eiλT1.Also ε1 is(37a)D12R12+R12=−2D1D2R11+σ2R11+q1E1eiT1−q1C1−iG1p1D1R11+p2D1R21−1−q2C2−iG2D1R11+D1R21+γ1p1F1−F2R¯11+p2F1−F2R¯213,(37b)D12R22+η2R22=−2D1D2R21+σ2η2R21−q1E1eiT1+q1C1−iG1p1D1R11+p2D1R21−q2C2−iG2D1R11+D1R21+γ2p1F1−F2R¯11+p2F1−F2R¯213.Assume that general solutions of equations (36a) and (36b) are(38a)R11=B1eiT1+φ1+cc+b1eiλT1,(38b)R21=B2eiηT1+φ2+cc+b2eiλT1,where b1=1−q2λ2E2/1−λ2 and b2=q2λ2E2/η2−λ2.Substitute equations (38a) and (38b) into equations (37a) and (37b), and then, eliminate secular terms and separate real and imaginary parts, as(39a)D2B1=−12α4B1+q1E1sinφ2,(39b)B1D2φ1=−12α1B13+α2B1B22+α3+σ2B1+q1E1cosφ2,(39c)D2B2=−12β4B2,(39d)B2D2φ2=−12ηB2β1B22+β2B12+β3+ησ2η.Equate right sides of equations (39a)–(39d) to zero; we obtain B2=0 due to equation (39c). Since B2 is the amplitude of the harmonic term with frequency Ω2 at ω1≈Ω1, B2=0 is correct when Ω2≫Ω1 (assumption 2). Therefore, equations (39a)–(39d) are simplified as(40a)D2B1=−12α4B1+q1E1sinφ2,(40b)B1D2φ2=−12α1B13+α3+σ2B1+q1E1cosφ2.Finally eliminatingφ2 in equations (40a) and (40b), the amplitude frequency equation for the second primary resonance at ω1≈Ω1 is(41)α12B3+2α1α3+σ2B2+α3+σ22+α42B−E12q12=0,where B=B12.


## 3.4. Stability Analysis

The 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].LetX=R1,R2,R1′,R2′T=X1,X2,X3,X4T; transform equations (25a) and (25b) into X′τ=FXτ,τ, and X∗ denotes equilibriums of equations (25a) and (25b). Take a perturbation ΔX to X∗, as(42)ΔX′≈∂Fτ,X∗∂X∗ΔX=Aτ,X∗τΔX.The stability ofX∗ can be confirmed according to the Jacobi matrix Aτ,X∗τ.The approximate monodromy matrix can be obtained according to Hsu’s method, as(43)M=ΨT=∏n=1NneAnΔT=∏n=1NnI+∑j=1NjAnΔTjj!,where I represents the identity matrix, T denotes the time period of the solution, which is divided into Nn subintervals, ΔT is the length of one subinterval, An is a constant matrix which represents the time-varying matrix Aτ,X∗τ, and Nj represents the number of terms in the approximation of the constant matrix An exponential. The Floquet index of matrix M is available by calculating eigenvalues, and the stability of X∗ can be confirmed.


## 3.5. 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. The ode45 function in MATLAB (Runge–Kutta method) is applied to solve original 4DOF equations (17a) and (17b). The initial states of motion for the first rotation speed ω1=600rad/s are set as X1X1′Y1Y1′X2X2′Y2Y2′=00000000, but for other rotation speeds, the initial states of motion are set as the final states of motion for the previous rotation speed. Integration duration should be long enough to ensure steady-state responses; for equations (17a) and (17b), integration duration are set as τ=0,2000π. The root mean square (RMS) [32] is applied to represent the amplitude of the vibration solved by the Runge–Kutta method. The amplitudes of LP and HP rotor are rj=∑i=1N^Xji−X¯j2+Yji−Y¯j2/N,j=1,2, where j is the number of points in one period and X¯j, Y¯j are the average values of corresponding response.Figure 3
Amplitude frequency curves solved by the multiple scales method and the Runge–Kutta method. (a) For the LP rotor. (b) For the HP rotor. Solid lines represent stable solutions, and dotted lines represent unstable solutions.
(a)(b)In Figure3, 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. There 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. The 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. The 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.


## 4. Bifurcation Analysis

The 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. Therefore, 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. The 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 σ=Ω12/ω12 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 dual-rotor system, as(44a)F1=α12A3+2α1λα3+λ2−σA2+λα3+λ2−σ2+λ2α42A−λ41−q22E22,(44b)F2=α12B3+2α1α3+1−σB2+α3+1−σ2+α42B−E12q12.The derivatives ofF1 and F2 can be expressed as(45a)F1A=3α12A2+4α1λα3+λ2−σA+λα3+λ2−σ2+λ2α42,(45b)F1AA=6α12A+4α1λα3+λ2−σ,(45c)F1σ=−2α1A2−2λα3+λ2−σA,(45d)F2B=3α12B2+4α1α3+1−σB+α3+1−σ2+α42,(45e)F2BB=6α12B+4α1α3+1−σ,(45f)F2σ=−2α1B2−2α3+1−σB,(45g)F1B=F2A=F1AB=F2AB=0.Transition sets are consisted of three parts, which are bifurcation set, hysteresis set, and double limit set. The bifurcation set indicates that unfolding equations (44a) and (44b) contain a singular point. The hysteresis set indicates that the unfolding makes at least quadratic contact with the vertical plane of a constant bifurcation parameter. The 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:(46a)Bi=α1,α4∈R2∃A,B,σs.t.F1=0,F2=0,F1AF2B−F1BF2A=0,F1AF2σ−F1σF2A=0.Hysteresis set is as follows:(46b)Hy=α1,α4∈R2∃A,B,σs.t.F1=0,F2=0,F1AF2B−F1BF2A=0,F1AA′−F1BB′=0,F1Af2−F2Af1=0,W′≠0,W′=A′,B′,f1=F1AAA′2+2F1ABA′B′+F1BBB′2,f2=F2AAA′2+2F2ABA′B′+F2BBB′2.Double limit set is as follows:(46c)DL=α1,α4∈R2∃W1,W2,σs.t.F1=0,F2=0,F1AF2B−F1BF2A=0,W1≠W2,W=A,B.According to above calculation procedure equations (44a)–(46c), the transition sets for two primary resonances of the dual-rotor system can be attained. The transition sets contain two hysteresis set curves H1 and H2, as shown in the kn−ks parameters plane in Figure 4 with ks from 0 to 15 (108 N/m) and kn from 0 to 150 (1023 N/m3). There are three regions divided by H1 and H2 in the kn−ks parameters plane, which indicates that two primary resonances of the dual-rotor system have three different bifurcation modes. Three 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 H1 and H2 can also be obtained, as shown in Figures 8 and 9.Figure 4
Transition sets inkn−ks parameters plane.Figure 5
Amplitude frequency curves for the bifurcation mode of the region I in Figure4 (kn=1×1023N/m3,ks=12×108N/m). (a) For the LP rotor. (b) For the HP rotor. Solid lines represent stable solutions, and dotted lines represent unstable solutions.
(a)(b)Figure 6
Amplitude frequency curves for the bifurcation mode of the region II in Figure4 (kn=30×1023N/m3,ks=9×108N/m). (a) For the LP rotor. (b) For the HP rotor. Solid lines represent stable solutions, and dotted lines represent unstable solutions.
(a)(b)Figure 7
Amplitude frequency curves for the bifurcation mode of the region III in Figure4 (kn=80×1023N/m3,ks=8×108N/m). (a) For the LP rotor. (b) For the HP rotor. Solid lines represent stable solutions, and dotted lines represent unstable solutions.
(a)(b)Figure 8
Amplitude frequency curves for the bifurcation mode of the hysteresis set curveH1 in Figure 4 (kn=20×1023N/m3,ks=10×108N/m). (a) For the LP rotor. (b) For the HP rotor. Solid lines represent stable solutions, and dotted lines represent unstable solutions.
(a)(b)Figure 9
Amplitude frequency curves for the bifurcation mode of the hysteresis set curveH2 in Figure 4 (kn=45×1023N/m3,ks=10×108N/m). (a) For the LP rotor. (b) For the HP rotor. Solid lines represent stable solutions, and dotted lines represent unstable solutions.
(a)(b)Figure5 shows the bifurcation mode of region I in Figure 4, in which the linear stiffness ks is rather large and the nonlinear stiffness kn is rather small. It can found that two primary resonances A and B have a linear characteristics. The 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. This 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. The 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 curvesH1 and H2 are investigated. Figure 8 shows the bifurcation mode of hysteresis set curve H1, 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 H2, 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.


## 5. 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. The 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. The results show that the linear stiffness and the nonlinear stiffness of the intershaft bearing have a significant influence on bifurcation modes of the system. The transition set, including two hysteresis set curves, divides the parameters planekn−ks 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. The future work will concentrate on the experimental study of the dual-rotor system.



---

*Source: 7820635-2020-05-20.xml*