Modeling and dynamic characteristic analysis of dual rotor-casing coupling system with rubbing fault

Abstract

With the rapid development of aero-engine manufacturing technology, the dual-rotor system has been employed in part of turbofan engine in order to improve the working performance of aircraft more efficiently. In this study, taking the counter-rotation dual-rotor as the research object, the dynamic model of dual rotor-casing coupling system is established by the aid of MATLAB. The dynamic frequency curves are in good agreement with the results in references and calculated by FEM method, that shows the validity and feasibility of the model. The local rub-impact dynamic model of dual rotor-casing coupling system is established, and rubbing analysis is carried out using Newmark-β method. The effects of rotating speed and speed ratio on local rub-impact response are deeply discussed. The results show that with the increase of rotating speed, combined frequencies and frequency multiplication components are more significant. In addition, speed ratio has a great influence on the periodic motion of the system. With the increase of the absolute value of the speed ratio, the whirl radius of the outer rotor and the normal rubbing force increase dramatically.

Introduction

The dual rotor system is widely used in aeroengine. In order to reduce vibration amplitudes and improve the performance of the whole machine, such as reliability, stability, and so on, the research of dual rotor structure is deeply studied. In the dual rotor system, a special intershaft-bearing structure connects the inner rotor and outer rotor, that causes the coupling of the vibration behavior of the inner and outer rotors. ¹ The dynamic characteristics of the dual rotor system become more complex compared with single rotor system.

With less variables and equations, the transfer matrix method is more efficient to solve rotor dynamic problems. Considering the gyroscopic moment, shear deformation, and bearing characteristics, Bansal et al. ² established the analytical structure matrix of the dual rotor system by the transfer matrix method. Using the Muller quadratic interpolation method, they obtained the complex eigenvalues and corresponding modes and analyzed the stability of the dual rotor bearing system. Adopting the transfer matrix method to establish the motion equation of the multiaxial rotating machinery, Hibner et al.^3–5 calculated the critical speed and nonlinear damping response of the system and discussed the stability of the squeeze film damper bearing through theoretical and experimental analyses.

The modal synthesis method is widely used in calculating the eigenvalue of the system by computing the eigenvalues of multiple substructures. By the aid of the modal synthesis method, Glasgow et al. ⁶ established the motion differential equation of the dual rotor system and calculated the critical speed and mode shape of the system. They found that the modal truncation number had a great influence on the accuracy of the positive and negative precession frequency of the system. Without considering the damping of the system, Li and Gunter^7,8 analyzed the application of the substructure mode synthesis method in the calculation of the natural characteristics of the dual rotor system. Through discussing the influence of the modal truncation error on the dynamic characteristics of the system, they depicted the reasonable range of the modal truncation number.

Considering the gyroscopic effect and bearing damping, Childs ⁹ established a typical dynamic model of the dual rotor system and analyzed the inherent characteristics on the basis of research of the traditional Jeffcott rotor. Considering the gyroscopic moment, bending, and shear deformation of rotating shafts, Chiang et al.^10,11 adopted finite element method to study the dynamic characteristics of single axis- and dual axis-bearing turbomachinery systems. They analyzed the inherent characteristics of the system and discussed the influence of speed ratio of high- to low-pressure rotor on the critical speed.

Some scholars had carried out the experimental research on the dual rotor system. Gupta et al.^12,13 calculated the critical speed, mode shape, and unbalanced response of the dual rotor system through the transfer matrix method. The accuracy of the results was validated through comparing with the test results obtained by a set of dual rotor test equipment. The influence of the coupling support between the inner and outer rotors was studied by the numerical simulation and the experimental method. In addition, Shafei, ¹⁴ Rajan, ¹⁵ Hu, ¹⁶ and Shiau ¹⁷ had also conducted relevant research on dual rotor systems.

In aeroengine field, rotor–stator rubbing is one of the most frequent faults. ¹⁸ Many scholars had studied the single rotor-casing rub-impact and achieved abundant results.^19–26 Considering the influence of coating deposited on the outer surface of disc and the inner surface of casing on rub-impact, Yang et al.^27,28 established a dynamic model of aeroengine in the case of the coupling fault of mass unbalance and local rub-impact. Based on Lankarani–Nikravesh model and Runge–Kutta numerical algorithm, they deduced a new rubbing force model and solved the dynamic rubbing response of system. Finally, the validity of the dynamic model was verified by comparing with the experimental results.

Considering the influence of gyroscopic effect, Sun et al.^29,30 adopted multiharmonic balance combined with the alternating frequency/time-domain technique (MHB-AFT) to obtain the nonlinear dynamic equation of an 8-DOF dual rotor with rub-impact fault. The complex nonlinear phenomena caused by rub impact, such as combined harmonic vibration, hysteresis, and resonance peak shift, were solved by the aid of arc-length continuation method. The stability of numerical solution was discussed using Floquet theory. In order to further understand the steady-state response characteristics of the dual rotor system under blade-casing rubbing, they discussed the effects of different control parameters such as mass eccentricity, intershaft bearing stiffness, and rotating speed on the dynamic characteristics of the system.

Regarding the rub-impact plate as an elastic diaphragm, Xu et al. ³¹ proposed a new finite element rubbing model. Considering the elastic deformation, contact penetration, and nonlinear spring damping support, the transient response of the system under rubbing condition was solved based on the wheel rim-elastic diaphragm rubbing force model. The rub-impact experiments were carried out on the dual rotor test-bed, and the experimental results were very consistent with those of the theoretical simulation.

Considering the coupling effect of the high- and low-pressure rotor through the intershaft bearing and the disc deflection caused by gyroscopic effect, Jia et al. ³² established a new axial–radial coupling dynamic rubbing model for the counter-rotation dual rotor system. The detailed expressions of the bearing coupling forces were given according to the two different forms of intershaft bearing.

Based on the rigid-flexible coupling multibody model, Han et al. ³³ established the local rub-impact model of the dual rotor system using MSC. By using ADAMS software, they analyzed the nonlinear dynamic characteristics of the system. Considering the flexibility of the casing, Wang et al. ³⁴ established a collision force model between a disk and a fixed stopper using the Lankarani–Nikravesh model. Using numerical integration method, they discussed the influence of speed ratio, initial clearance, the mass eccentricity, and the bearing stiffness between the shafts on the nonlinear characteristics of the system under rubbing condition.

Lu et al. ³⁵ presented a review of proper orthogonal decomposition (POD) techniques applied to order reduction in a variety of research areas. The applications of POD method in high-dimensional nonlinear dynamic systems were prospected in detail in order to provide direct guidance for researchers in various fields of engineering. Taking into account the nonlinearities of all of the supporting rolling element bearings, Jin et al. ³⁶ established a dynamic model of a complex dual rotor-bearing system of an aero-engine based on the finite element method with rigid disc, cylindrical beam element, and conical beam elements. The nonlinear vibration responses of the dual rotor-bearing system were studied in the case of different clearances of the intershaft bearing. Lu et al. ³⁷ applied the POD method to the dimension reduction of dual rotor-bearing experiment rig for the first time. The dynamical characteristics of the reduced-order model were compared with the full order system and experimental results, which verified that the POD method presented higher computational efficiency and accuracy.

On the basis of above research studies, this study establishes the coupling dynamic model of dual rotor-bearing-casing system, and the inherent characteristics of the system are analyzed, the accuracy of which is validated by comparing with the results in references. Using this model, the effect of rotating speed, speed ratio, and mass unbalance on the rubbing responses is studied. The evolution rules of the axis trajectory, the periodic motion, and bifurcation behavior of the system along with the variation of speed ratio are studied. The difference of dynamic responses between considering the rub impact or not is analyzed, and the influences of rub impact on axis trajectory, time-domain response, and amplitude spectrum are pointed out.

Dual rotor-casing system

Modeling

Figure 1 is a schematic diagram of the dual rotor-casing coupling system of the aeroengine, which consists of a high-pressure rotor (outer rotor), a low-pressure rotor (inner rotor), a casing coupled with the high- and low-pressure rotor, and some supports. Only the transverse and bending vibrations of the rotors are considered, while the translation along and torsional vibration around z-axis are ignored.

Figure 1.

Schematic of a dual rotor-casing coupling system.

Timoshenko beams with constant cross-section are adopted to establish the model of the two rotors and the casing. The linear spring-damping model is used to simulate the casing-foundation, the inner rotor-casing, and outer rotor-casing supports. k_bx1, k_by1, k_bx4, and k_by4 are the support stiffness of the inner rotor-casing supports at the left and right ends in x- and y-directions, respectively. c_bx1, c_by1, c_bx4, and c_by4 are the damping of the inner rotor-casing supports at the left and right ends in x- and y-directions, respectively. k_bx2, k_by2, c_bx2, and c_by2 are the stiffness and damping of outer rotor-casing supports at the left end in x- and y-directions, respectively. k_bx3, k_by3, c_bx3, and c_by3 are the stiffness and damping of inner–outer rotor intershaft bearing in x- and y-directions. k_bx5, k_by5, k_bx6, k_by6, c_bx5, c_by5, c_bx6, and c_by6 are the stiffness and damping of casing-foundation supports at the left and right ends in x- and y-directions, respectively.

Parameters of shafts, casing, disks, and bearings are excerpted from Ref., ³⁸ shown in Tables 1 and 2, respectively. The parameters of the shaft segment elements of the inner and outer rotors and the casing are shown in Table 1, while those of the disks and supports are shown in Table 2. Four discs are arranged on the inner and the outer rotor shaft, located at nodes 3, 11, 16, and 19, respectively. ³⁸

Table 1.

The element data of dual-rotor and casing system.

	Element number	Segment length (mm)	Inner diameter (mm)	Outer diameter (mm)	Material parameters
Inner axis	1	50	0	22	ρ=7800 kg/m3E=210 GPaυ=0.3
	2	50	0	22
	3	50	0	22
	4	50	0	22
	5	50	0	22
	6	50	0	22
	7	50	0	22
	8	50	0	22
	9	50	0	22
	10	50	0	22
	11	50	0	22
	12	50	0	22
Outer axis	13	35	40	52
	14	40	40	52
	15	50	40	52
	16	50	40	52
	17	50	40	52
	18	40	40	52
	19	35	40	52
Casing	20	30	200	204
	21	50	200	204
	22	50	200	204
	23	50	200	204
	24	35	200	204
	25	40	200	204
	26	50	200	204
	27	50	200	204
	28	50	200	204
	29	40	200	204
	30	35	200	204
	31	50	200	204
	32	50	200	204
	33	50	200	204
	34	30	200	204

Table 2.

Parameters of disks and bearings.

	Geomotry parameters	Material parameters
Disk 1 and 2	m1 = 9.683 kg, Rd = 100 mm, hd = 40 mm, Jp = 0.04898 kg·m2, Jd = 0.0258 kg·m2	ρ = 7800 kg/m3
Disk 3 and 4	m2 = 9.139 kg, Rd = 100 mm, hd = 40 mm, Jp = 0.04878 kg·m2, Jd = 0.02560 kg·m2	E = 210 GPa
Bearings	kbx1 = kby1 = 2.6 × 107 N/m, cbx1 = cby1 = 1.1 × 103 N·s/m	υ = 0.3
	kbx3 = kby3 = 8.75 × 10 N/m, cbx3 = cby3 = 1.1 × 103 N·s/m
	kbx2 = kby2 = kbx4 = kby4 = 1.75 × 106 N/m, cbx2 = cby2 = cbx4 = cby4 = 1.1 × 103 N·s/mkbx5 = kby5 = kbx6 = kby6 = 2.5 × 109 N/m, cbx5 = cby5 = cbx6 = cby6 = 2.1 × 103 N·s/m

The expression of element stiffness matrix K ^e is written as ³⁹ :

The equations of motions

Dual rotor-casing matrix

Based on the idea of finite element method, the dual rotor-casing coupling system is modeled using MATLAB software. The two nodes spatial beam element is adopted to establish the model of the rotating shafts. For ordinary spatial beam element, each node has six degrees of freedom, which are translation along and rotation around the x-, y-, and z-directions, respectively. The element is shown in Figure 2, where x _A , y _A , z _A , x _B , y _B , and z _B are translations along x-, y-, and z-directions while θ _Ax , θ _Ay , θ _Az , θ _Bx , θ _By , and θ _Bz are rotations around x-, y- and z-directions for nodes A and B, respectively.

Figure 2.

Finite element model of shaft unit.

Without considering the translation along and rotation around the z- direction, the degree of freedom of the beam element is

$${\mathbf{q}}^e={\left(x_A{y}_A{\theta }_{Ax}{\theta }_{Ay}{x}_B{y}_B{\theta }_{Bx}{\theta }_{By}\right)}^{\text{T}}$$

(1)where superscript e represents parameter for the element. The expression of shear influence factor is given by ³⁹

$$\phi =\frac{12EI}{\kappa A_{s}G{l}^2},\kappa =\{\begin{array}{l}6\left(1+\upsilon \right)/\left(7+6\upsilon \right)\hspace{1em}\text{for solid}\mathsf{ }\text{beam}\mathsf{ }\text{element}\\ 2\left(1+\upsilon \right)/\left(4+3\upsilon \right)\hspace{1em}\text{for hollow}\mathsf{ }\text{beam}\mathsf{ }\text{element}\end{array}$$

(2)where, E, G, ρ, υ, and A _s are the elastic modulus, shear modulus, density, Poisson’s ratio of the material, and the cross-sectional area of the element, respectively.

The formula of element mass matrix M ^e is given by ³⁹

$${\mathbf{M}}^e=\rho A_{s}l\left[\begin{array}{cccccccc}a& & & & & & & \\ 0& a& & & & & & \\ 0& -c& g& & & \text{symmetric}& & \\ c& 0& 0& g& & & & \\ b& 0& 0& d& a& & & \\ 0& b& -d& 0& 0& a& & \\ 0& d& f& 0& 0& c& g& \\ -d& 0& 0& f& -c& 0& 0& g\end{array}\right]$$

(3)where

$$a=\frac{\frac{13}{35}+\frac{7}{10}\phi +\frac{1}{3}{\phi }^2+\frac{6}{5}{\left(r_g/l\right)}^2}{{\left(1+\phi \right)}^2},b=\frac{\frac{9}{70}+\frac{3}{10}\phi +\frac{1}{6}{\phi }^2-\frac{6}{5}{\left(r_g/l\right)}^2}{{\left(1+\phi \right)}^2}$$

$$c=\frac{\left[\frac{11}{210}+\frac{11}{120}\phi +\frac{1}{24}{\phi }^2+\left(\frac{1}{10}-\frac{1}{2}\phi \right){\left(r_g/l\right)}^2\right]l}{{\left(1+\phi \right)}^2}$$

$$d=\frac{\left[\frac{13}{420}+\frac{3}{40}\phi +\frac{1}{24}{\phi }^2-\left(\frac{1}{10}-\frac{1}{2}\phi \right){\left(r_g/l\right)}^2\right]l}{{\left(1+\phi \right)}^2}$$

$$f=\frac{\left[\frac{1}{140}+\frac{1}{60}\phi +\frac{1}{120}{\phi }^2+\left(\frac{1}{30}+\frac{1}{6}\phi -\frac{1}{6}{\phi }^2\right){\left(r_g/l\right)}^2\right]l^2}{{\left(1+\phi \right)}^2}$$

$$g=\frac{\left[\frac{1}{105}+\frac{1}{60}\phi +\frac{1}{120}{\phi }^2+\left(\frac{2}{15}+\frac{1}{6}\phi +\frac{1}{3}{\phi }^2\right){\left(r_g/l\right)}^2\right]l^2}{{\left(1+\phi \right)}^2}$$

$${\mathbf{K}}^e=\left[\begin{array}{cccccccc}h& & & & & & & \\ 0& h& & & & & & \\ 0& -i& j& & & \text{symmetric}& & \\ i& 0& 0& j& & & & \\ -h& 0& 0& -i& h& & & \\ 0& -h& i& 0& 0& h& & \\ 0& -i& k& 0& 0& i& j& \\ i& 0& 0& k& -i& 0& 0& j\end{array}\right]$$

(4)where $h=\frac{12EI}{l^3(1+\phi )}$ , $i=\frac{6EI}{l^2(1+\phi )}$ , $j=\frac{(4+\phi )EI}{l(1+\phi )}$ , and $k=\frac{(2-\phi )EI}{l(1+\phi )}$

The formula of element gyroscopic matrix G ^e is depicted as ³⁹

$${\mathbf{G}}^e=2\Omega \rho A_{s}l\left[\begin{array}{cccccccc}0& & & & & & & \\ -p& 0& & & & & & \\ -q& 0& 0& & & \text{symmetric}& & \\ 0& -q& -s& 0& & & & \\ 0& -p& -q& 0& 0& & & \\ p& 0& 0& -q& -p& 0& & \\ -q& 0& 0& w& q& 0& 0& \\ 0& -q& -w& 0& 0& q& -s& 0\end{array}\right]$$

(5)where

$$p=\frac{6r_g^2/5}{l^2{(1+\phi )}^2}$$

$$q=\frac{-\left(\frac{1}{10}-\frac{1}{2}\phi \right)r_g^2}{l{\left(1+\phi \right)}^2}$$

$$s=\frac{\left(\frac{2}{15}+\frac{1}{6}\phi +\frac{1}{3}{\phi }^2\right)r_g^2}{{\left(1+\phi \right)}^2}$$

$$w=\frac{-\left(\frac{1}{30}+\frac{1}{6}\phi -\frac{1}{6}{\phi }^2\right)r_g^2}{{\left(1+\phi \right)}^2}$$

where l is the unit length, J is the polar moment of inertia of the section, I is the moment of inertia of the section, $\Omega =\dot{\theta }$ is the rotating angular speed, and $r_g=\sqrt{I/A_s}$ is the radius of gyration.

The linear spring-damping model is chosen for describing the dynamic behavior of the bearing. The element stiffness matrix and damping matrix of rotor-casing and casing-foundation supports are denoted by K _b1, K _b2, K _s1, K _s2, C _b1, C _b2, C _s1, and C _s2, respectively, while the specific expressions of the elements are shown as follows ³⁸

$$\begin{array}{l}\{\begin{array}{l}{\mathbf{K}}_{\text{b}1}=\text{diag}\left(k_{\text{b}x1}{k}_{\text{b}y1}00\right)\\ {\mathbf{K}}_{\text{b}2}=\text{diag}\left(k_{\text{b}x2}{k}_{\text{b}y2}00\right)\end{array}\\ \{\begin{array}{l}{\mathbf{K}}_{\text{s}1}=\text{diag}\left(k_{\text{s}x1}{k}_{\text{s}y1}00\right)\\ {\mathbf{K}}_{\text{s}2}=\text{diag}\left(k_{\text{s}x2}{k}_{\text{s}y2}00\right)\end{array}\end{array}$$

(6)

$$\begin{array}{l}\{\begin{array}{l}{\mathbf{C}}_{\text{b}1}=\text{diag}\left(c_{\text{b}x1}{c}_{\text{b}y1}00\right)\\ {\mathbf{C}}_{\text{b}2}=\text{diag}\left(c_{\text{b}x2}{c}_{\text{b}y2}00\right)\end{array}\\ \{\begin{array}{l}{\mathbf{C}}_{\text{s}1}=\text{diag}\left(c_{\text{s}x1}{c}_{\text{s}y1}00\right)\\ {\mathbf{C}}_{\text{s}2}=\text{diag}\left(c_{\text{s}x2}{c}_{\text{s}y2}00\right)\end{array}\end{array}$$

(7)

The lumped mass model is used for the disk. The specific expressions of the element mass matrix M _d and the gyroscopic matrix G _d of the disk are given by ³⁸

$${\mathbf{M}}_{\text{d}}=\text{diag}\left[m_{\text{d}}{m}_{\text{d}}{J}_d{J}_d\right]$$

(8)

$${\mathbf{G}}_{\text{d}}=\mathrm{\Omega }\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& J_p\\ 0& 0& -J_p& 0\end{array}\right]$$

(9)

Matrix assembling of dual rotor-casing system

The general differential equations of motion for dual rotor-casing system are given by ³⁹

$$M\ddot{q}+\left(C+G\right)\dot{q}+Kq=F$$

(10)where M , G , C , K , q , and F are mass matrix, damping matrix (including supporting damping matrix and gyroscopic matrix of disk and shaft), Rayleigh damping matrix, stiffness matrix, generalized displacement vector, and external force vector, respectively. Rayleigh damping matrix C is given by ³⁹

$$\mathbf{C}=\zeta \mathbf{M}+\eta \mathbf{K}$$

(11)

$$\zeta =\frac{4\text{π}{f}_{\text{n}1}{f}_{\text{n}2}\left({\xi }_1{f}_{\text{n}2}-{\xi }_2{f}_{\text{n}1}\right)}{\left(f_{\text{n}2}^2-f_{\text{n}1}^2\right)},\eta =\frac{{\xi }_2{f}_{\text{n}2}-{\xi }_1{f}_{\text{n}1}}{\text{π}\left(f_{\text{n}2}^2-f_{\text{n}1}^2\right)}$$

(12)where f_n1 and f_n2 are the first and second natural frequencies of the system (Hz), respectively, and ξ₁ and ξ₂ are corresponding modal damping ratios.

The assembling of stiffness matrix of the dual rotor-casing coupling system K is shown in Figure 3. The assembling method of mass matrix M and damping matrix G is similar to that of stiffness matrix, so the unnecessary details are ignored here.

Figure 3.

Schematic diagram of assembled stiffness matrix for dual rotor-casing coupling system.

Assuming that the intershaft bearing is located at the i-th node of the inner rotor and the j-th node of the outer rotor. The expressions of the stiffness matrix K ^cb and the damping matrix C ^cb for the intershaft bearing are shown as follows(13)

Model verification

Table 3 shows the comparison of natural frequencies using MATLAB and ANSYS under static state. It can be found that the results of ANSYS are slightly higher than those of MATLAB; however, the maximum error is 1.587%.

Table 3.

Natural frequencies using MATLAB and ANSYS (Ω = 0 r/min).

Order	ANSYS (Hz)	MATLAB (Hz)	Error (%)
1	35.102	35.002	−0.028
2	35.102	35.002	−0.028
3	94.265	94.122	−0.015
4	94.265	94.122	−0.015
5	197.899	198.001	0.005
6	197.899	198.001	0.005
7	267.036	262.797	−1.587
8	267.036	262.797	−1.587
9	298.801	294.105	−1.572
10	298.801	294.105	−1.572
11	315.103	310.182	−1.562
12	315.103	310.182	−1.562
13	683.184	680.297	−0.004
14	683.184	680.297	−0.004

Figures 4(a),(b) are Campbell diagrams of the dual rotor-casing system under the co-rotation and counter-rotation conditions. It can be found that:

(1) The dynamic frequency lines calculated by MATLAB are basically consistent with those computed by ANSYS, and the first three order forward and backward whirl frequencies are in good agreement, which shows the validity and feasibility of the model established in this article.

(2) With the increase of the rotating speed, the counter-rotation weakens the influence of the gyroscope torque more significantly. In Figure 4(b), the distinctions of the forward and backward precession frequencies for each order under counter-rotation condition are obviously smaller than those under co-rotation, that is the reason why the counter-rotation dual rotor system is more widely used in the aeroengine.

Figure 4.

Campbell diagram of dual rotor-casing system: (a) co-rotation and (b) counter-rotation.

The critical speed diagram of the dual rotor under co-rotation is shown in Figure 5. The solid line and dotted line in the diagram represent the critical speeds under the excitation applied on the low-pressure rotor and on the high-pressure rotor, respectively. The intersection points of the curves with the dash-dotted line, implying the condition that η = 1.5, numbered as 1, 2, 3, 4, 5, and 6, are the critical speeds of the dual rotor-casing system during the acceleration process. The results agree well with those in Ref., ¹ and the maximum error is only 0.75% (see Table 4).

Figure 5.

Critical speed diagram using the model of Ref. ¹ (a) and present model (b).

Table 4.

Critical speed distribution and error.

Order	Present model	The model of Ref. 1	Error (%)
1	2310	2310	0
2	2430	2426	0.16
3	5895	5914	0.32
4	6010	6025	0.25
5	11,940	12,030	0.75
6	11,980	12,060	0.66

In order to further validate the accuracy of the present model, the transient displacements of the center of Disk 4 are studied and compared with those in Ref. ¹ under the exponential law Ω = 10,000 [1−exp (−1.05 t)] rpm and the linear law Ω = $\{\begin{array}{c}2500t,0\le t<4\\ \mathrm{10000,4}\le t\le 6\end{array}$ rpm. The specific simulation parameters are as follows: total time: 6 s, time step size, 0.0004 s, speed ratio η = ω₂/ω₁ = 1.5, and model damping ratio ξ₁ = ξ₂ = 0.001. The simulation results using the model of Ref. ¹ and present model under the two different speed-up laws are shown in Figure 6. It can be seen from the figure that the time-domain waveforms of the two models under different laws are quite consistent, which further proves the accuracy of the model in this article.

Figure 6.

Transient displacements of the center of Disk 4 using the model of Ref. ¹ under the linear law (a), the exponential law (b), using present model under the linear law, (c) and the exponential law (d).

Rubbing analysis of dual rotor-casing system

Newmark-β method is adopted to solve the response of dual rotor-casing rubbing. With the characteristics of high efficiency and numerical stability, the differential form of Newmark-β integration method is ¹

$$\{\begin{array}{l}{\dot{u}}_{t+\mathrm{\Delta }t}={\dot{u}}_t+\left[\left(1-\delta \right){\ddot{u}}_t+\delta {\ddot{u}}_{t+\Delta t}\right]\mathrm{\Delta }t\\ u_{t+\mathrm{\Delta }t}=u_t+{\dot{u}}_t\mathrm{\Delta }t+\left[\left(0.5-\alpha \right){\ddot{u}}_t+\alpha {\ddot{u}}_{t+\mathrm{\Delta }t}\right]\mathrm{\Delta }{t}^2\end{array}$$

(14)where δ and α are parameters that can be adjusted according to the requirements of accuracy and numerical integration stability, and Δt is the integration step length. According to equation (14), the following equation can be derived

$$\{\begin{array}{l}{\ddot{u}}_{t+\mathrm{\Delta }t}=\frac{1}{\alpha \mathrm{\Delta }{t}^2}\left(u_{t+\mathrm{\Delta }t}-u_t\right)-\frac{1}{\alpha \mathrm{\Delta }t}{\dot{u}}_t-\left(\frac{1}{2\alpha }-1\right){\ddot{u}}_t\\ {\dot{u}}_{t+\mathrm{\Delta }t}=\frac{\delta }{\alpha \mathrm{\Delta }t}\left(u_{t+\mathrm{\Delta }t}-u_t\right)+\left(1-\frac{\delta }{\alpha }\right){\dot{u}}_t+\mathrm{\Delta }t\left(1-\frac{1}{2\alpha }\right){\ddot{u}}_t\end{array}$$

(15)

By substituting equation (15) into equation (10), we can get

$$M^{\prime }{\ddot{u}}_{t+\mathrm{\Delta }t}+C^{\prime }{\dot{u}}_{t+\mathrm{\Delta }t}+K^{\prime }{u}_{t+\mathrm{\Delta }t}={F^{\prime }}_{t+\mathrm{\Delta }t}$$

(16)

Equation (16) can be simplified as

$$K^{\ast }{u}_{t+\mathrm{\Delta }t}=F_{t+\mathrm{\Delta }t}^{\ast }$$

(17)where

$$\{\begin{array}{l}{K}^{\ast }=\frac{1}{\alpha \mathrm{\Delta }t}{M}^{\prime }+\frac{\delta }{\alpha \mathrm{\Delta }t}{C}^{\prime }+K^{\prime }\\ F_{t+\mathrm{\Delta }t}^{\ast }={F^{\prime }}_{t+\mathrm{\Delta }t}+\left[\frac{1}{\alpha \mathrm{\Delta }{t}^2}{u}_t+\frac{1}{\alpha \mathrm{\Delta }t}{\dot{u}}_t+\left(\frac{1}{2\alpha }-1\right){\ddot{u}}_t\right]M^{\prime }+\left[\frac{\delta }{\alpha \mathrm{\Delta }t}{u}_t+\left(\frac{\delta }{\alpha }-1\right){\dot{u}}_t+\left(\frac{\delta }{2\alpha }-1\right)\mathrm{\Delta }t{\ddot{u}}_t\right]C^{\prime }\end{array}$$

(18)

When the Newmark integral parameters meet the following conditions, the result of the equation will be unconditionally stable

$$\alpha =\frac{1}{4}{\left(1+{\gamma }_1\right)}^2,\hspace{1em}\delta =\frac{1}{2}+{\gamma }_1,\hspace{1em}{\gamma }_1\ge 0,\hspace{1em}\left({\gamma }_1=\text{0}\hspace{0.17em}\text{in}\hspace{0.17em}\text{this}\hspace{0.17em}\mathrm{paper}\right)$$

(19)

Local rubbing force model

Figure 7 is a simplified local-rubbing schematic of dual rotor-casing system. When the system is in a static state, the rotor center is located at O_d, and the casing center O_c coincides with the coordinate origin O. R_d and R_c are radius of the disk and inner-wall of the casing, respectively. The static parallel misalignment e between the two centers is preset before rubbing. For calculation convenience, the condition that x_m = y_m is used. The minimum clearance c between the rotor and the casing is expressed as

$$c=R_{\text{c}}-R_{\text{d}}-\sqrt{x_{\text{m}}^2+y_{\text{m}}^2},e=\sqrt{x_{\text{m}}^2+y_{\text{m}}^2}$$

(20)

Figure 7.

Local-rubbing schematic of dual rotor-casing system.

In the actual working process, due to the existence of installation error, rub impact, and other fault factors, the rotor center and casing center deviate from the original position, located at O_d’ and O_c’, respectively. In this case, the radial relative displacement ξ is given by

$$\xi =\sqrt{{\left(x_{\text{d}}-x_{\text{c}}+x_{\text{m}}\right)}^2+{\left(y_{\text{d}}-y_{\text{c}}+y_{\text{m}}\right)}^2}$$

(21)

When ξ < e + c, there is no rubbing; when ξ ≥ e + c, rubbing phenomenon occurs. Setting δ_c as the penetration depth between rotor and casing, then

$${\delta }_{\text{c}}=\{\begin{array}{l}\xi -e-c,\xi \ge e+c\\ 0,\xi \text{<}e+c\end{array}$$

(22)

Normal and tangential rubbing forces can be written as

$$\{\begin{array}{l}{F}_{\text{N}}\left(x,y\right)=k_{\text{r}}{\delta }_{\text{c}},\xi \ge e+c\\ F_{\text{T}}\left(x,y\right)=\mu F_{\text{N}}\end{array}$$

(23)where μ and k_r are friction coefficient and rubbing stiffness, respectively.

Effect of rotating speed on rubbing response

Node 3 on the inner rotor (at the location of left disk) and the node 25 on the casing are set as the rub-impact points. The specific simulation parameters are as follows: the unbalance magnitude of four disks f_u = 156 × 10⁻⁶ kg m, speed ratio η = ω₂/ω₁ = −1.2, initial clearance c = 1×10⁻⁴ m, rubbing stiffness k_r = 8 × 10⁶ N/m, the static parallel misalignment e = 8 × 10⁻⁵ m, friction coefficient μ = 0.3, and model damping ratio ξ₁ = ξ₂ = 0.08.

Figure 8 is a three-dimensional spectrum of rub-impact response of the center of disk 1 at different speeds. Here, the operating speed range is 800∼5800 r/min. While the rotating speed range is 800∼2000 r/min, the number of times that rubbing occurs is relatively fewer. The spectrum mainly shows the corresponding rotating frequency of inner and outer rotors, denoted by f₁ and f₂, respectively. When the rotating speed Ω = 1700 r/min, amplitude amplification phenomenon occurs on the spectrum line corresponding to the outer rotor, this is because the rotating frequency of the outer rotor arouses the first natural frequency of the system. When the rotating speed reaches the range of 2000∼4000 r/min, besides rotating frequencies, frequency multiplication, frequency division, and the combination frequency of the rotating frequencies appear, such as f₂ − f₁, 1/2f₂, 2f₁ − f₂, 2f₁, f₁ + f₂, and 2f₂. When the rotating speed of inner rotor Ω=2300 r/min, the first natural frequency of the system is excited and obvious amplitude amplification phenomenon occurs (point B in Figure 8). When the rotating speed reaches the range of 4000∼5800 r/min, the rubbing between the rotor and the casing becomes more vigorous, and some high frequency components and combined frequency appear. The amplitudes are increasing, and the frequency multiplication components are much more significant. The second-order natural frequencies of the system are excited by the inner and outer rotors (points C and D in Figure 8), and the amplitude amplification phenomenon is more highlighted under the corresponding speed.

Figure 8.

Three-dimensional spectrum of inner rotor rubbing under different rotating speeds.

Effect of speed ratio on rubbing response

Figure 9 is a three-dimensional spectrum of rub-impact response under different speed ratios. It can be seen from the figure that the spectrum lines present a capital “N” shape, and the high-frequency component in the spectrum is not prominent. There are mainly five spectral lines, corresponding to the frequencies of f₂ − f₁, f₁, f₂, 2f₁, and f₁ + f₂. Under the condition that η = −1, the rotating speeds of inner and outer rotors are in the same value, resulting in the amplitude superposition at the corresponding spectral lines, so that the amplitude at point B in the figure is significantly higher than those of other places. Furthermore, amplitude amplifications occur at points A and C. This is because the combined frequency f₂ − f₁ excites the first natural frequency of the system and resonance occurs at point A. What is more, the rotating frequency of outer rotor f₂ excites the second natural frequency of the system at point C. With the increase of the absolute value of the speed ratio η, amplitude variations corresponding to the rotating frequency of inner rotor are not notable, while those corresponding to outer rotor appear an increasing trend. It can be found that the factors which cause the resonance of the system are rather complicated. On the one hand, both of the inner and outer rotor rotating frequencies can excite the natural frequencies of the system. Due to the effect of speed ratio, the excitations of natural frequencies by outer rotor are prior to those by inner rotor (point C in the figure); on the other hand, under the fixed speed condition, different values of speed ratio arouse the resonance of the system at different combination frequencies. Here, the resonance peaks are not only limited to the rotating frequencies of the inner and outer rotors but also located at combined frequency points (point A in the figure).

Figure 9.

Three-dimensional spectrum of inner rotor rubbing under different speed ratios.

The influence of unbalanced force on the response of dual rotor system can be intuitively reflected by axis orbit diagram under co-rotation and counter-rotation conditions. The rotating speed is set as Ω = 2300 r/min, and Figures 10, 11 show the axis orbit diagram of node 3 on the inner rotor under counter-rotation and co-rotation conditions, respectively. It can be seen that the axis orbit under counter-rotation condition is simpler and more regular than under co-rotating condition. Most of the axis orbits under counter-rotation condition appear “petal” shape, and the number of “petals” is different under different speed ratios; while axis orbits under co-rotating condition appear the shape of big circle covering small circle. By comparing with the former, the axis orbits are much more complex. At the condition that η = 1 or η = −1, the axis orbits under counter-rotation condition present the shape of inclined ellipses, while those under co-rotation condition present the shape of circles. At the condition that η = 1.5 or η = −1.5, there are “sudden changes” in the axis orbits of the inner rotor, at this time, period-2 motion occurs in both cases.

Figure 10.

Rotor orbit of counter-rotation dual-rotor under different speed ratios.

Figure 11.

Rotor orbit of co-rotation dual-rotor under different speed ratios.

In addition, the axis orbits at the condition that 10η is odd are obviously more complex than those at the condition that 10η is even. The inner and outer rotors under counter-rotation condition rotate at reverse directions, that weakens the gyroscopic effect and reduces the amplitude of the unbalanced response, so that the counter-rotation dual rotor system is more stable than the co-rotation dual rotor system. It can also be seen from the comparison of the two figures. Therefore, most scholars focus on the dynamic characteristics of the counter-rotating dual rotor system.

Figure 12 is a bifurcation diagram of the displacement of the inner rotor under different speed ratios with rubbing fault. The region (−2, −1) is selected as the speed ratio range, and the interval is set as −0.02. The rotating speed is Ω = 3000 r/min. Different from the traditional single rotor system, the response mainly carries out period-2, period-5, and period-10 motions when only one digit after the decimal point of the speed ratio is considered, such as −1.5, −1.7, and −1.9. When the speed ratios are integers, such as η = −1 and η = −2 in the diagram, the response carries out single period motion. When two decimals are considered, the response appears multi-period motions. It shows that the speed ratio has a great influence on the periodic motion of the response.

Figure 12.

Bifurcation diagram of inner rotor under different speed ratios.

The two characteristic speed ratios η = −1.5 and η = −1.7 corresponding to period 2 and period 10 motions are selected to study the rubbing responses. Figure 13 shows the rubbing responses of node 3 on inner rotor and node 25 on the casing (η = −1.5). It can be seen that the rubbing level is relatively lower. The wave clipping phenomenon of the inner rotor is not significant, also there is almost no change at response amplitude of the casing before and after rub impact, that can be proved by the axis orbit diagram. During the rubbing process, the rotating frequency of inner rotor f₁ and outer rotor f₂ and the frequency component of f₂ − f₁ are mainly excited. There are many high frequency components in the response of the casing, but the amplitudes are very small. It can be seen from the spectrum in Figure 13(b) that the influence of the rotating frequency of outer rotor f₂ on the response of the casing is much larger than that of inner rotor f₁.

Figure 13.

System rubbing responses with η = −1.5: (a) node 3 on inner rotor, (b) node 25 on casing. (Note: the figures from left to right are time-domain waveform, rotor orbit and amplitude spectrum).

Figure 14 shows the rub-impact response of the inner rotor and the casing (η = −1.7). Wave clipping phenomenon occurs after rubbing, that is also reflected in the axis orbit diagrams. The amplitude amplification phenomenon appears at combined frequency f₂ − f₁=35 Hz, which is close to the first natural frequency of the system. A lot of combined frequencies are excited on the frequency spectrum of the inner rotor, such as 2f₁ − f₂, f₂ − f₁ and f₁ + f₂.

Figure 14.

System rubbing responses with η = −1.7: (a) node 3 on inner rotor and (b) node 25 on casing. (Note: the figures from left to right are time-domain waveform, rotor orbit, and amplitude spectrum).

Figure 15 shows the normal rubbing forces under different speed ratios. The response presents period-2 motion, period-10 motion, and period-10 motion at the condition that η = −1.5, η = −1.7, and η = −1.9, respectively. With the increase of the absolute value of the speed ratio, the whirl radius of the outer rotor and the normal rubbing forces increase dramatically.

Figure 15.

Normal rubbing force of casing node 25 under different speed ratios: (a) η = −1.5, (b) η = −1.7, and (c) η = −1.9.

Effect of mass unbalance on rubbing response

The impact of mass unbalance magnitudes on the rubbing response is discussed in this section. The specific simulation parameters are as follows: the rotating speed Ω 2300 rpm, speed ratio η = ω₂/ω₁ = −1.2, initial clearance c = 1 × 10⁻⁴ m, rubbing stiffness k_r = 8 × 10⁶ N/m, the static parallel misalignment e = 8 × 10⁻⁵ m, friction coefficient μ = 0.3, and model damping ratio ξ₁ = ξ₂ = 0.08.

Figures 16, 17 show the rubbing response of node 3 on the inner rotor under the condition that mass unbalance magnitude f_u = 200 × 10⁻⁶ kg m and f_u = 550 × 10⁻⁶ kg m, respectively. It can be seen from the figure that the whirl trajectory values and response amplitudes of the inner rotor increase with the growth of mass unbalance magnitude. When the rubbing occurs, the wave clipping phenomenon of the inner rotor becomes more significant under higher value of mass unbalance magnitude. Due to the coupling effect of the inner rotor, outer rotor, and the casing, the axis orbit of the inner rotor becomes more complicated, and the spectrum amplitude is larger with the increase of the mass unbalance magnitude. The response frequency components mainly concentrate on the rotating frequency of the inner rotor f₁, the outer rotor f₂, combined frequency f₂ − f₁, and multiple rotating frequencies.

Figure 16.

System rubbing responses with f_u = 200 × 10⁻⁶ kg m: (a) time-domain waveform, (b) rotor orbit, and (c) amplitude spectrum.

Figure 17.

System rubbing responses with f_u = 550 × 10⁻⁶ kg m: (a) time-domain waveform, (b) rotor orbit, and (c) amplitude spectrum.

Figure 18 shows the response of node 3 on the inner rotor under the condition that mass unbalance magnitude f_u = 550 × 10⁻⁶ kg m without considering the impact of rubbing. The response of the inner rotor presents cyclostationary state during the time-domain history under the condition that no rubbing occurs [see Figure 18(a)]. In the rotor orbit as shown in Figure 18(b), the order of magnitude is the same as that in Figure 17(b), but in a more regular form. The frequency components are concentrated on the rotating frequency of the inner rotor f₁ and the outer rotor f₂, without any combined frequencies or multiple frequency components [see Figure 18(c)].

Figure 18.

System responses with f_u = 200 × 10⁻⁶ kg m neglecting the rubbing impact: (a) time-domain waveform, (b) rotor orbit, and (c) amplitude spectrum.

Conclusion

The model of dual rotor-casing coupling system is established, and then the local rubbing response is analyzed in this article. The effects of rotating speed and speed ratio on rub-impact response are studied. Some detailed conclusions are listed as follows:

(1) The dynamic frequency curves computed by the model in this article are in good agreement with the results in references and calculated by FEM method, that shows the validity and feasibility of the model established in this article.

(2) The counter-rotation weakens the influence of the gyroscope torque, that is the reason why the counter-rotation dual rotor system is more widely used in the aeroengine.

(3) As the rotating speed increases, the rubbing between the rotor and the casing becomes more vigorous, and some high frequency components and combined frequency appear. Meanwhile, the frequency multiplication components are more significant.

(4) Speed ratio has a great influence on the periodic motion of the system. In addition, with the increase of the absolute value of the speed ratio, the whirl radius of the outer rotor and the normal rubbing force increase dramatically.