Nonlinear dynamic characteristics induced by inter-shaft rub-impact in a counter-rotational dual-rotor system of aeroengine

Abstract

The influence and mechanism of inter-shaft rub-impact have not been given a wide attention in exciting literature. This paper aims to insight into the coupled bending-torsional analysis of a counter-rotational dual-rotor system nonlinear vibration characteristics under inter-shaft rubbing fault. Firstly, a dual-rotor finite element model is established, which is based on Timoshenko beam element considering torsional deformations. In particular, the impact force is described by the nonlinear cylindrical collision model. Secondly, the modal characteristics are clarified in detail by the Campbell diagram and modal shapes diagram. Finally, the Newmark-β method is applied to solve the dual-rotor system governing equations with inter-shaft rub-impact numerically. Bifurcation diagram, three-dimensional spectrum, whirl orbit, spectrum plot, and Poincare map are used to preliminarily demonstrate the rotor system dynamic behaviors with inter-shaft rub-impact and the influences of friction coefficient and elastic modulus of the coating. Results reveal the vibration coupling effect and complicated vibration characteristics caused by inter-shaft rub-impact.

Keywords

Dual-rotor system inter-shaft rub-impact nonlinear vibration analysis bifurcation

Introduction

As the core component of multiple advanced aeroengines, dual-rotor structure has been considered as a common configuration due to its greater advantages than single rotor in achieving better aerodynamic stability. In order to shorten the loads span, the inter-shaft bearing is used to connect the LPR (low pressure rotor) and HPR (high pressure rotor), which causes the coupling effect between rotors.¹ Another form of dual-rotor structure, that is, shared support-rotor system whose rotors are supported independently without inter-shaft bearing, is considered as one of the ways to lighten the engine.² The introduction of inter-shaft bearing or shared support-rotor system constitutes a more complex engine structure, which leads to an indirect mutual influence between the rotors. Under the effect combining multiple excitation sources and coupling vibrations, complicated dynamic behaviors of dual-rotor system could be observed. In addition, various faults may occur in the aeroengine, such as the blade-casing rub-impact,³ rotor to stator rub-impact,⁴ misalignment,⁵ and bearing defects.⁶ These faults will affect rotor dynamic properties and lead to aeroengine unstable operation. In severe cases, these faults will cause damage and failure of the equipment. Therefore, it is of practical significance for investigating the dual-rotor systems complicated vibration characteristics and reducing the vibrations amplitude of the whole aeroengine subjected to various faults.

The ever-changing analytical techniques and numerical models provide effective approaches to analyze the dual-rotor system dynamic response besides experiments. Thanks to decades efforts of many scholars, the lumped-mass method,⁷ transfer matrix method,⁸ and finite element method (FEM)⁹ are verified fit for describing the dynamics of dual-rotor. The lumped-mass method is applicable to qualitative analysis of rotor, and the transfer matrix method is suitable for simple analysis, such as calculating the critical speed. The first two methods used to be applied; furthermore, the FEM has been widely used recently due to its better applicability and higher accuracy. Chen¹⁰ introduced a new rotor-ball-bearing-casing coupling dynamic model by using FEM, and the modeling method was verified by experiments and fault simulations. Ma¹¹ proposed a FE model of a shaft–disk–blade system, and the model vibration subjected to rubbing between the blade-tip and casing was studied. Kang¹² established a finite element model of a dual-rotor system considering the intermediate bearing nonlinearity and studied the parameter influence on the back whirling stability through the experimental and simulation method. Additionally, Wang,¹³ Lu,¹⁴ and Ma¹⁵ et al. adopted FEM to build respective dual-rotor systems and studied the dynamic properties. These studies demonstrated that the FEM can be considered as a powerful tool applied to the dual-rotor system nonlinear dynamic analysis.

Reducing rotor-stator clearance provides a common way to reduce specific fuel consumption and improve engine efficiency dramatically, resulting in probability of rubbing between rotating components and stationary parts. Over the past few decades, many researchers have focused primarily on the nonlinear dynamical phenomenon of rub-impact fault which is caused by other rotor system faults and has an important effect on the aeroengine reliability. Chu et al.¹⁶ established the Jeffcott rub-impact model and found there were three paths leading to chaos with the increase of rotation speed. Patel et al.¹⁷ proposed a rotor-stator rub-impact dynamic model by using beam element containing torsional degree of freedom, and the rubbing force was characterized by contact stiffness, damping, and Coulomb friction. Yang et al.¹⁸ simulated the rotor system dynamic model with imbalance-fixed point rubbing coupling faults, and the influence of the material of the casing convex point on rub-impact was studied. Herisanu, N et al.¹⁹ explored the nonlinear vibration of Jeffcott rotor with rub-impact by deriving the governing equations of the rotor and establishing an approximate analytical solution. Zhang et al.²⁰ developed a novel indicator for rotor system with rub-impact based on the vibration power flow theory, and the accuracy of the theoretical was validated by a corresponding test. In addition, some key phenomena and conclusions were summarized by Muszynska,²¹ Jacquet-Richardet,²² Jiang,²³ and Ma.²⁴

The researches mentioned above are related to the single rotor system rub-impact. Due to more complicated components, when rubbing is aroused in the dual-rotor system, it will lead to various nonlinear characteristics. In recent years, the dual-rotor rub-impact faults were revealed by many scholars through various means including theoretical analysis, numerical simulation, and experiment. Considering the nonlinear hysteretic force model, Yang et al.^.25,26 developed the fixed point rubbing motion equations of the dual-rotor system, and several combination frequencies were found as the particular fault frequencies in the case of fix point rubbing. Zhang et al.²⁷ studied the dual-rotor system dynamic vibrations with rub-impact considering bending-torsion coupling by using the Lagrange equation, the simulation results showed that the torsional vibration exhibits richer characteristic frequency. Jin et al.^28,29 established a dual-rotor model considering misalignment and blade-casing rubbing, the nonlinear vibration characteristics verified by numerical simulation and experiment were revealed, and a strong self-excited vibration occurred for a small blade-tip clearance. Liu et al.³⁰ established an engine dual-rotor nonlinear support structure model, in which the inter-shaft bearing, squeeze film damper, casing, and nonlinear Hertz contact force from rolling bearing were considered, and nonlinear characteristics of the rubbing fault were analyzed in detail. Ma et al.³¹ established a typical aeroengine dual-rotor system analytical model, and the influence of various factors, including the unbalanced position, the rotor speed, rubbing stiffness, clearance, and friction coefficient on the dual-rotor vibration response was studied. Prabith et al.³² proposed a dual-rotor aeroengine numerical procedure for capturing the nonlinear dynamical responses and obtained substantial progress on the dual-rotor system response characteristics undergoing multi-disk rub-impact. The simulation results showed that the multi-disk rub-impact nonlinear behaviors were more complicated, and a very large amplitude was caused by the dry friction backward whirl.

Sealing structure is used in dual-rotor engines to prevent the fluid flow in adjacent paces, which may lead to interaction of rotors at the sealing installation position.³³ As a new vibration transmission way between rotors, the inter-shaft rub-impact will result in new rotors coupling modes and produce complex nonlinear characteristics. In the previous researches, a few scholars paid attention to the rubbing fault vibration characteristics. The research of Chen³⁴ on inter-shaft rub-impact indicated that periodic, quasi-periodic, and multi-periodic appeared in the dual-rotor system with increase of rotating speed. Based on the dual-rotor FE model considering inter-shaft rub-impact, Yu et al.^35,36 compared the differences of modal characteristics and vibration when the rotor system rotated in the same and opposite directions and studied the effects of rotational speed, contact stiffness, and friction coefficient on dynamic characteristics. Generally speaking, the actual rotor–rotor contact is cylindrical contact. However, in the previous researches on inter-shaft rub-impact, impact force was all expressed by the point contact model, such as the piecewise linear model, which is difficult to describe the inter-shaft rub-impact force. Moreover, the inter-shaft rub-impact will produce friction torque on both shafts, which will affect the rotor system vibration response. Therefore, it is informative to carry out this kind of dynamical analysis of dual-rotor system for rotor system fault diagnoses and aeroengine running stability.

In this paper, a counter-rotation aeroengine dual-rotor system simulation model is established for studying nonlinear vibration with inter-shaft rub-impact. The paper structure is organized as follows. In section 2, the finite element model of a simplified dual-rotor system is proposed and the nonlinear impact force is considered. In section 3, the modal shapes, bifurcation phenomenon, and vibration responses under inter-shaft rub-impact are investigated, and the effects of different friction coefficients and elastic modulus are studied. Some conclusions are offered in section 4.

Dynamic model with inter-shaft rub-impact

Dual-rotor system finite element model

An ordinary aeroengine dual-rotor system composition is shown in Figure 1, which includes a LPR supported by bearings No. 1 and 4 in 1-0-1 manner, and a HPR supported by bearings No. 2 and 3 in 0-2-0 manner. Since the inter-shaft bearing is not adopted, so the coupling between rotors will not occur under the ideal working condition of the engine. The finite element model which can reflect the dual-rotor system nonlinear dynamic response is established. Meanwhile, the following simplification are made

(1) The Timoshenko beam element model considering gyroscopic moments, rotary inertias, bending, torsional, and shear deformations is applied to simulate rotors.

(2) The lumped-mass model is used to represent blades and disks, and the linear spring and elastic damping are adopted to simulate supports.

(3) The LPR and HPR consist of two rigid turntables and two supports, respectively. The relative position of the rotor and the disk is shown in Figure 1.

(4) The seal is fixed on the LPR, where the inter-shaft rub-impact occurs and the rubbing nodes are marked in Figure 1.

(5) An abrasion-resistant coating is painted on the rubbing node surface of HPR, and the thermal effect in rub-impact is not considered.

Figure 1.

The dual-rotor system finite element model.

According to the simplified structure, the system is divided into 22 beam elements and 24 nodes, and the rubbing force is considered to be applied to the rubbing nodes shown in the Figure 1. Geometric parameters of simplified rotor system are listed in Table 1, and material parameters of rotor system are listed in Table 2.

Table 1.

The geometric parameters for rotors.

Parameter	The value between nodes
Length of LPR (mm)	L₁₂	L₂₃	L₃₄	L₄₅	L₅₆	L₆₇
	40.3	39.5	33.7	96.6	96.6	90
	L₇₈	L₈₉	L₉₁₀	L₁₀₁₁	L₁₁₁₂	L₁₂₁₃
	90	22.7	86.5	86.5	43	299
Length of HPR (mm)	L₁₄₁₅	L₁₅₁₆	L₁₆₁₇	L₁₇₁₈	L₁₈₁₉	L₁₉₂₀
	35.4	53.3	49.9	49.8	69	69
	L₂₀₂₁	L₂₁₂₂	L₂₂₂₃	L₂₃₂₄
	68	65.4	42.4	23.7
Diameter of LPR (mm)	D₁₂	D₂₃	D₃₄	D₄₅	D₅₆	D₆₇
	30	35	35	54	54	47.6
	D₇₈	D₈₉	D₉₁₀	D₁₀₁₁	D₁₁₁₂	D₁₂₁₃
	47.6	37.6	27.6	27.6	27.6	22
Inner diameter of HPR (mm)	ID₁₄₁₅	ID₁₅₁₆	ID₁₆₁₇	ID₁₇₁₈	ID₁₈₁₉	ID₁₉₂₀
	61.2	61.2	61.2	61.2	53.6	53.6
	ID₂₀₂₁	ID₂₁₂₂	ID₂₂₂₃	ID₂₃₂₄
	53.6	31.6	31.6	31.6
Outer diameter of HPR (mm)	OD₁₄₁₅	OD₁₅₁₆	OD₁₆₁₇	OD₁₇₁₈	OD₁₈₁₉	OD₁₉₂₀
	70	70	72.6	72.6	65	65
	OD₂₀₂₁	OD₂₁₂₂	OD₂₂₂₃	OD₂₃₂₄
	65	42	39	39

Table 2.

Material parameters of rotor system.

Parameters	Values
Elastic modulus of shaft	$E_{b}$
(GPa)	210
Density of shaft	$ρ$
(kg/m³)	7850
Poisson’s ratio	0.3
Mass of disk	$m_{d 1}$	$m_{d 2}$	$m_{d 3}$	$m_{d 4}$
(kg)	5.378	8.8852	3.9933	3.3437
Diametric inertial moment	$J_{d 1}$	$J_{d 2}$	$J_{d 3}$	$J_{d 4}$
(kg·m²)	0.00164	0.038	0.0146	0.012
Polar inertial moment	$J_{p 1}$	$J_{p 2}$	$J_{p 3}$	$J_{p 4}$
(kg·m²)	0.00328	0.076	0.0292	0.024
Support stiffness	$k_{s 1}$	$k_{s 2}$	$k_{s 3}$	$k_{s 4}$
(N/m)	5 × 10⁷	1.1× 10⁷	8×10⁷	5 × 10⁶
Support damping	$c_{s 1}$	$c_{s 2}$	$c_{s 3}$	$c_{s 4}$
(N·s/m)	200	100	100	200
Eccentricity	$e_{1}$	$e_{2}$	$e_{3}$	$e_{4}$
(10^-5m)	2	1	1	2
Elastic modulus of coating	$E_{c}$
(MPa)	350

Beam element matrix

Generally, the homogeneous assumptions are made for the materials used in the rotors. As can be seen in Figure 2, the Timoshenko beam element adopted in this paper includes 2 nodes, each node contains 5-DOFs including lateral vibration and rotation in z direction. The generalized displacement vector of shaft element is defined as $q_{e} = {[\begin{array}{l} w_{x 1} & w_{y 1} & θ_{x 1} & θ_{y 1} & θ_{z 1} & w_{x 2} & w_{y 2} & θ_{x 2} & θ_{y 2} & θ_{z 2} \end{array}]}^{T}$ , where $w_{x i}$ ( $i$ = 1, 2) denotes the two translations in the x direction, $w_{y i}$ ( $i$ = 1, 2) denotes the two translations in the y direction, $θ_{x i}$ ( $i$ = 1, 2) denotes the two rotations around the x axis, $θ_{y i}$ ( $i$ = 1, 2) denotes the two rotations around the y axis, and $θ_{z i}$ ( $i$ = 1, 2) denotes the two rotations around the z axis.

Figure 2.

Discrete elements of the finite element model: (a) Timoshenko beam element and (b) disc element.

Several physical parameters of beam, including second moment of area $I$ , radius of gyration $r_{g}$ , and shear deformation coefficient $ϕ$ , are obtained to simplify the expressions, and these physical parameters are given as

I = \frac{π}{64} (D^{4} - d^{4})

(1)

r_{g} = \sqrt{\frac{I}{A}}

(2)

ϕ = \frac{12 E I}{κ G A l^{2}} κ = {\begin{cases} 6 (1 + ν) / (7 + 6 ν) \\ 2 (1 + ν) / (4 + 3 ν) \end{cases} \begin{array}{l} for solid beam \\ for hollow beam \end{array}

(3)

where geometry parameters

D

d

, and

l

denote the outer diameter, inner diameter, and length of the shaft-beam element; material property parameters of the shaft-beam element

E

G

ν

, and

ρ

are the Young’s modulus, shear modulus, Poisson ratio, and density, respectively. Here, the mass matrix of shaft element is obtained as³⁷

M_{e} = ρ A l [\begin{array}{c} m_{1} \\ 0 & m_{1} \\ 0 & - m_{2} & m_{5} \\ m_{2} & 0 & 0 & m_{5} & s y m m e t r i c \\ 0 & 0 & 0 & 0 & 2 m_{7} \\ m_{3} & 0 & 0 & m_{4} & 0 & m_{1} \\ 0 & m_{3} & - m_{4} & 0 & 0 & 0 & m_{1} \\ 0 & m_{4} & m_{6} & 0 & 0 & 0 & m_{2} & m_{5} \\ m_{4} & 0 & 0 & m_{6} & 0 & - m_{2} & 0 & 0 & m_{5} \\ 0 & 0 & 0 & 0 & m_{7} & 0 & 0 & 0 & 0 & 2 m_{7} \end{array}]

(4)

where

m_{1} = \frac{\frac{13}{35} + \frac{7}{10} ϕ + \frac{1}{3} ϕ^{2} + \frac{6}{5} {(\frac{r_{g}}{l})}^{2}}{{(1 + ϕ)}^{2}}

m_{2} = \frac{\frac{9}{70} + \frac{3}{10} ϕ + \frac{1}{6} ϕ^{2} - \frac{6}{5} {(\frac{r_{g}}{l})}^{2}}{{(1 + ϕ)}^{2}}

m_{3} = \frac{\frac{11}{210} + \frac{11}{120} ϕ + \frac{1}{24} ϕ^{2} + (\frac{1}{10} - \frac{1}{2} ϕ) {(\frac{r_{g}}{l})}^{2}}{{(1 + ϕ)}^{2}} l

m_{4} = \frac{\frac{13}{420} + \frac{3}{40} ϕ + \frac{1}{24} ϕ^{2} - (\frac{1}{10} - \frac{1}{2} ϕ) {(\frac{r_{g}}{l})}^{2}}{{(1 + ϕ)}^{2}} l

m_{5} = \frac{\frac{1}{140} + \frac{1}{60} ϕ + \frac{1}{120} ϕ^{2} + (\frac{1}{30} + \frac{1}{6} ϕ - \frac{1}{6} ϕ^{2}) {(\frac{r_{g}}{l})}^{2}}{{(1 + ϕ)}^{2}} l^{2}

m_{6} = \frac{\frac{1}{105} + \frac{1}{60} ϕ + \frac{1}{120} ϕ^{2} + (\frac{2}{15} + \frac{1}{6} ϕ + \frac{1}{3} ϕ^{2}) {(\frac{r_{g}}{l})}^{2}}{{(1 + ϕ)}^{2}} l^{2}

m_{7} = \frac{r_{g}^{2}}{6 ρ}

The stiffness matrix is obtained as³⁷

K_{e} = K [\begin{array}{c} k_{1} \\ 0 & k_{1} \\ 0 & - k_{2} & k_{3} \\ k_{2} & 0 & 0 & k_{3} & s y m m e t r i c \\ 0 & 0 & 0 & 0 & k_{5} \\ - k_{1} & 0 & 0 & - k_{2} & 0 & k_{1} \\ 0 & - k_{1} & k_{2} & 0 & 0 & 0 & k_{1} \\ 0 & - k_{2} & k_{4} & 0 & 0 & 0 & k_{2} & k_{3} \\ k_{2} & 0 & 0 & k_{4} & 0 & - k_{2} & 0 & 0 & k_{3} \\ 0 & 0 & 0 & 0 & - k_{5} & 0 & 0 & 0 & 0 & k_{5} \end{array}]

(5)

where

K = \frac{E I}{(1 + ϕ) l^{3}}

k_{1} = 12, k_{2} = 6 l, k_{3} = (4 + ϕ) l^{2}

k_{4} = (2 - ϕ) l^{2}, k_{5} = \frac{2 G I}{l K}

The gyroscopic matrix is obtained as³⁷

G_{e} = 2 Ω ρ A l [\begin{array}{c} 0 \\ - p & 0 \\ - q & 0 & 0 \\ 0 & - q & - s & 0 & s y m m e t r i c \\ 0 & 0 & 0 & 0 & 0 \\ 0 & - p & - q & 0 & 0 & 0 \\ p & 0 & 0 & - q & 0 & - p & 0 \\ - q & 0 & 0 & o & 0 & q & 0 & 0 \\ 0 & - q & - o & 0 & 0 & 0 & q & - s & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{array}]

(6)

where

Ω

is the rotating angular speed,

p = \frac{\frac{6}{5} {(\frac{r_{g}}{l})}^{2}}{{(1 + ϕ)}^{2}}, q = \frac{- (\frac{1}{10} - \frac{1}{2} ϕ) {(\frac{r_{g}}{l})}^{2}}{{(1 + ϕ)}^{2}} l

s = \frac{(\frac{2}{15} + \frac{1}{6} ϕ + \frac{1}{3} ϕ^{2}) {(\frac{r_{g}}{l})}^{2}}{{(1 + ϕ)}^{2}} l^{2}, o = \frac{- (\frac{1}{30} + \frac{1}{6} ϕ - \frac{1}{6} ϕ^{2}) {(\frac{r_{g}}{l})}^{2}}{{(1 + ϕ)}^{2}} l^{2}

The specific expressions of mass matrix $M_{d}$ and the gyroscopic matrix $G_{d}$ of disks unit simulated by lumped-mass model are obtained as³⁷

M_{d} = d i a g [\begin{array}{l} m_{d} & m_{d} & J_{d} & J_{d} & J_{p} \end{array}]

(7)

G_{d} = Ω [\begin{array}{l} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & J_{p} & 0 \\ 0 & 0 & - J_{p} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{array}]

(8)

Assuming that the support is symmetrical in the x and y directions, the bearing is simulated with the linear spring-damping model to obtain satisfactory results and simplify the calculation. The element stiffness matrix and damping matrix are given by

K_{s} = d i a g [\begin{array}{l} k_{s} & k_{s} & 0 & 0 & 0 \end{array}]

(9)

C_{s} = [\begin{array}{l} c_{s} & 0 & 0 & 0 & 0 \\ 0 & c_{s} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{array}]

(10)

The element matrices of disk and rotating shaft are assembled according to the corresponding node numbers to form the overall mass matrix $M$ , overall stiffness matrix $K$ , and overall gyroscopic matrix $G$ . Rayleigh damping matrix $C$ is given by

C = ζ M + η K

(11)

ζ = \frac{4 π f_{n 1} f_{n 2} (ξ_{1} f_{n 2} - ξ_{2} f_{n 1})}{(f_{n 2}^{2} - f_{n 1}^{2})}, η = \frac{ξ_{2} f_{n 2} - ξ_{2} f_{n 1}}{π (f_{n 2}^{2} - f_{n 1}^{2})}

(12)

where

f_{n 1}

and

f_{n 2}

denote the first and second natural frequencies, respectively, and

ξ_{1}

and

ξ_{2}

are corresponding modal damping ratios. In this paper,

ξ_{l 1}

= 0.02 and

ξ_{l 2}

= 0.015 for LPR, and

ξ_{h 1}

= 0.014 and

ξ_{h 2}

= 0.005 for HPR.

Rubbing model

The accuracy of the contact force model is a crucial factor affecting the analysis results in the simulation of dynamical systems.³⁸ Obviously, unlike blade-casing rub-impact in the point contact framework, inter-shaft rub-impact is mechanically closer to the cylindrical contact. Meanwhile, the contact model that has been proved to be effective and widely used in point contact, such as the Lankarani-Nikravesh model, is simply a coarse approximation to estimate the contact forces between two cylindrical bodies.³⁹ Based on the shortcomings of point contact model and existing cylindrical contact model, a simple cylindrical contact force model is summarized by Liu et al.,⁴⁰ where contact force can be directly calculated by given penetration

F_{n} = \frac{1}{2} π δ^{*} E^{*} \sqrt{\frac{δ^{*}}{2 (Δ R + δ^{*})}}

(13)

where

\frac{1}{E^{*}} = \frac{1 - ν_{1}^{2}}{E_{1}} + \frac{1 - v_{2}^{2}}{E_{2}}

E_{i}

and

ν_{i}

(

i

= 1, 2) are the Young’s modulus and the Poisson’s ratio of two cylinders,

Δ R = R_{2} - R_{1}

represents the radial clearance of cylinders, and

δ^{*}

is the penetration of cylinders.

For revealing the inter-shaft rub-impact dynamic characteristic, the rubbing model between shafts is established as shown in Figure 3. The rubbing between LPR and HPR is simplified as the contact between two circles. $δ = \sqrt{{(x_{l} - x_{h})}^{2} + {(y_{l} - y_{h})}^{2}}$ is defined as the center distance of the two circles and $δ_{0}$ is defined as the initial clearance between shafts. When $δ$ is greater than the $δ_{0}$ with the two shafts rotating, it is considered that rub-impact occurs. Moreover, the inter-shaft rubbing force expression is derived from the force analysis combining the cylindrical contact model and Coulomb model. Therefore, the rubbing forces in x and y directions and the torque in the direction of rotation around the shaft are as follows

{\begin{cases} f_{x l} = H F_{n} (1 - \frac{δ_{0}}{δ}) (- (x_{l} - x_{h}) + s i g n (v_{r e l}) μ (y_{l} - y_{h})) \\ f_{y l} = H F_{n} (1 - \frac{δ_{0}}{δ}) (- s i g n (v_{r e l}) μ (x_{l} - x_{h}) - (y_{l} - y_{h})) \\ f_{t l} = μ H F_{n} R_{l} \end{cases}

(14)

{\begin{cases} f_{x h} = - f_{x l} \\ f_{y h} = - f_{y l} \\ f_{t h} = - f_{t l} \end{cases}

(15)

where

f_{x i}

and

f_{y i}

(

i

l

h

) denote the three components of the rubbing forces in the lateral, vertical, and rotation directions, respectively.

H

denotes

H e a v i s i d e

function,

μ

is the frictional coefficient, and

s i g n

denotes symbolic function. Relative speed of rotors at the rubbing point is represented by

v_{r e l}

, and the calculation formula is shown as

v_{r e l} = (Ω_{l} R_{l} + ω_{l} r_{l}) - (Ω_{h} R_{h} + ω_{h} r_{h})

ω_{i (i = l, h)} = \frac{1}{x^{2} + y^{2}} (x \dot{y} - \dot{x} y)

(16)

where

ω_{l}

and

ω_{h}

denote the rotors whirl speeds;

r_{l}

and

r_{h}

denote the rotors radii, respectively.

Ω_{l}

and

Ω_{h}

are the rotors rotation speeds;

R_{l}

and

R_{h}

denote the displacements of the rotors.

Figure 3.

Inter-shaft rub-impact force model: (a) Undeflected rotors and (b) rotors with inter-shaft rub-impact.

Global equation of motion

The external forces on the rotor are mainly composed of rub-impact forces and unbalance forces of the disk. The general dual-rotor system motion equation can be given by

M \ddot{q} + (C + G) \dot{q} + K q = F_{u} + F_{r}

(17)

where unbalance force vector is presented as

F_{u} = {[\begin{array}{l} 0 & \dots & \underset{u n b a l a n c e n o d e s}{\underset{⏟}{\begin{array}{l} f_{u x} & f_{u y} & 0 & 0 & f_{u t} \end{array}}} & \dots & 0 \end{array}]}^{T}

(18)

{\begin{cases} f_{u x} = m e ({\dot{φ}}^{2} \cos φ + \ddot{φ} \sin φ) \\ f_{u y} = m e ({\dot{φ}}^{2} \sin φ - \ddot{φ} \cos φ) \\ f_{u t} = m e (\ddot{x} \sin φ - \ddot{y} \cos φ) \end{cases}

(19)

where

m

and

e

denote the unbalanced mass and eccentricity of disk,

φ = Ω t + θ_{z} + θ_{0}

is the rotation angle and the initial phase

θ_{0}

= 0 for LPR and HPR, and

\ddot{x}

and

\ddot{y}

denote the lateral and vertical acceleration of disk, respectively. Rubbing force is presented as

F_{r} = [\begin{array}{l} 0 & \dots & \underset{r u b b i n g n o d e s f o r L P}{\underset{⏟}{\begin{array}{l} f_{x l} & f_{y l} & 0 & 0 & f_{t l} \end{array}}} & 0 & \dots & \underset{r u b b i n g n o d e s f o r H P}{\underset{⏟}{\begin{array}{l} f_{x h} & f_{y h} & 0 & 0 & f_{t h} \end{array}}} & 0 & \dots \end{array}]

(20)

In this article, the bending-torsional coupling is reflected in two aspects: firstly, the torsion angle coordinate and the lateral displacement coordinate are coupled in the unbalanced force according to equation (19); secondly, the friction torque is a function of the lateral displacement and torsion angle according to equation (13) and equation (14).

Numerical results and discussions

Influence on inherent characteristics

In this subsection, the dual-rotor system modes and critical speeds without or with inter-shaft rub-impact are compared to further study the modal changes caused by rub-impact. The speed ratio is defined as $η = Ω_{h} / Ω_{l}$ , and $η$ is equal to −1.2 in this paper, the negative sign represents the counter-rotation. In order to calculate the natural frequency of the system, homogeneous form of equation (17) is exhibited:

M \ddot{q} + (C + G) \dot{q} + K q = 0

(21)

By introducing column vectors, the order of equation (21) is reduced, and the second-order linear differential equation is expressed as the first-order linear differential equation:

\dot{v} = A v

(22)

where

A = [\begin{array}{c} - M^{- 1} (C + G) & - M^{- 1} K \\ I & 0 \end{array}]

(23)

v = [\begin{array}{l} \dot{q} \\ q \end{array}]

(24)

The eigenvalues of matrix $A$ are the rotor system natural frequencies, which are in the form of $λ_{k} = α_{k} \pm j ω_{k}$ . Due to the influence of gyroscopic moments, $ω_{k}$ will change with the rotor speed $Ω$ , and the intersections of the change curve and equation $ω_{k} = Ω$ are the critical speeds of each order.

When HPR and LPR move independently without inter-shaft rub-impact, the Campbell diagram and the first few mode shapes for each shaft are calculated separately at $Ω_{l}$ = 14325 rpm, and results are presented in Figures 4 and 5. For the mode diagram, “+” represents the rotors forward whirl (FW) vibration mode.

Figure 4.

Campbell diagram without inter-shaft rub-impact of different shafts: (a) LPR and (b) HPR.

Figure 5.

Modal shapes for $Ω_{l}$ = 14325 rpm without inter-shaft rub-impact: (a) n = 1+, (b) n = 2+, (e) n = 3+ for the LPR, (c) n = 1+, and (d) n = 2+ for the HP rotor.

The dynamic stiffness of rotors is changed by the introduction of rubbing force, which means adding a support at the rub-impact node. In general, the support stiffness $k_{r}$ reflects the severity of rubbing and is measured in terms of Hertz contact stiffness. In this subsection, $k_{r}$ is set as $5 \times 10^{5} N / m$ , and $k_{x} = k_{y} = k_{r}$ . Ignoring the nonlinear effect of friction, the linear effect of inter-shaft rub-impact can be regarded as a spring connecting LPR and HPR, and the stiffness matrix of the spring is

K_{c} = [\begin{array}{r} \underset{(5 o - 4, 5 o - 4)}{\underset{⏟}{k_{x}}} & 0 & \underset{(5 o - 4, 5 n - 4)}{\underset{⏟}{- k_{x}}} & 0 \\ 0 & \underset{(5 o - 3, 5 o - 3)}{\underset{⏟}{k_{y}}} & 0 & \underset{(5 o - 3, 5 n - 3)}{\underset{⏟}{- k_{y}}} \\ \underset{(5 n - 4, 5 o - 4)}{\underset{⏟}{- k_{x}}} & 0 & \underset{(5 n - 4, 5 n - 4)}{\underset{⏟}{k_{x}}} & 0 \\ 0 & \underset{(5 n - 3, 5 o - 3)}{\underset{⏟}{- k_{y}}} & 0 & \underset{(5 n - 3, 5 n - 3)}{\underset{⏟}{k_{y}}} \end{array}]

(25)

where

o

and

n

denote the rubbing nodes on the LPR and HPR, respectively.

The Campbell diagrams and vibration modes corresponding to the first five modal characteristics at $Ω_{l}$ = 14325 rpm are exhibited in Figures 6 and 7, respectively. The line of blue dots represents the position of the shaft without vibration. (a), (b), and (e) in Figure 7 are the mode shapes of the LPR, and (c) and (d) are the HPR mode shapes. The critical speeds before and after rubbing are listed in Table 3.

Figure 6.

Campbell diagram considering rubbing as elastic spring, $k_{r}$ = 5 × 10⁵ N/m.

Figure 7.

Modal shapes considering rubbing as elastic spring (a) n = 1+, (b) n = 2+, (c) n = 3+, (d) n = 4+, and (e) n = 5+ for $Ω_{l}$ = 14325 rpm.

Table 3.

Table of critical speeds with rubbing and without rubbing.

Order	Critical speeds		Excitation source
Order	With rubbing	Without rubbing	Excitation source
1	6601	6303	LPR
2	12770	12224	LPR
3	13007	12606	HPR
4	22130	22012	HPR
5	32423	32422	LPR

It can be found from Figure 7 that the rotor mode shapes correspond to the n = 1+-5+modes shown in Figure 5. Moreover, the vibration of the LPR can be excited by HPR as exhibited in Figure 7(c) and (d). The LPR vibration is transmitted to the HPR and induces it to generate weak vibration as shown in Figure 7(a) and (b). These phenomena show that the energy of one shaft can be transferred to the other shaft by inter-shaft rub-impact. Compared the critical speeds with and without rubbing in Table 2, it is found that the rotor critical speeds increase due to inter-shaft rub-impact. This result originates from the fact that the introduction of the support stiffness at the rubbing point improves the system stiffness. It is observed that the critical speed increase of the LPR third-order FW critical speed is not obvious. Combining with Figures 5 and 7, this phenomenon appears because rubbing node of LPR is located at the vibration node as shown in Figure 7(e), which causes the poor effect of adding stiffness.

Bifurcation and nonlinear response analysis

In this section, the dual-rotor system nonlinear characteristics with inter-shaft rub-impact are conducted. The governing equation based on the equation (17) derived in section 2 is implemented in MATLAB and solved by the Newmark-β method. The optimal calculation results and the convergence of equation iteration are ensured by setting $Δ t$ as 10⁻⁴ s. The total integration time is set as 4 s, which is long enough to ensure the stability of the vibration response solution. The case in which the calculation reaches the steady state is discussed. The eccentricities of disks and the elastic modulus of coating are list in Table 2. The clearance between the two shafts of rubbing position is $δ_{0}$ = 2 × 10⁻⁵ m and friction coefficient $μ$ = 0.1. The corresponding phases are all 0^o. The phase difference between disk 1 and disk 4 on the LPR is taken as the torsion angle.

The radial displacement at rubbing section can be used to distinguish the rubbing types according to the rubbing force mechanical expression given by equation (14). The variation curves of radial displacement maximum values $δ_{\max}$ and radial displacement minimum values $δ_{\min}$ with increasing speed are shown in Figure 8. When $δ_{\max}$ < $δ_{0}$ , there is no rub-impact between rotors. When $δ_{\min}$ < $δ_{0}$ < $δ_{\max}$ , the rubbing states are intermittent. When $δ_{\min}$ > $δ_{0}$ , continuous rub-impact occurs.

Figure 8.

The relative displacement of rubbing nodes for $μ$ = 0.1 and $E_{c}$ = 3.5 × 10⁸ Pa.

As shown in Figure 8, continuous rub-impact occurs when 5539 rpm< $Ω_{l}$ <8308 rpm and 11937 rpm< $Ω_{l}$ <13460 rpm. Additionally, the rubbing types can be summarized as follows: intermittent rub-impact - continuous rub-impact - intermittent rub-impact - continuous rub-impact - intermittent rub-impact.

The bifurcation diagram and three-dimensional spectrum of bending and torsional deformations for disk 4, are obtained as shown in Figure 9. It can be seen from Figure 9(a) that the nonlinear characteristic of disk 4 varies with the rotating speed. When $Ω_{l}$ <5252 rpm, rub-impact has not occurred, which can be verified from Figure 8, and LPR stabilizes in the 1T-periodic motion. When 5252 rpm< $Ω_{l}$ <6016 rpm, the motion of disk 4 enters into 5T-period. When 6016 rpm< $Ω_{l}$ <7353 rpm, the quasi-periodic motion appears. When 7353 rpm< $Ω_{l}$ <12319 rpm, 5T-period appears again. When 12319 rpm< $Ω_{l}$ <13752 rpm, 5T-periodic motion is found. When $Ω_{l}$ >13752 rpm, disk 4 returns to 1T-periodic motion. Combined with Figure 8, it is found that quasi-periodic can only emerge in the continuous rub-impact state.

Figure 9.

Dual-rotor system dynamic response for $μ = 0.1$ and $E_{c}$ = 3.5 × 10⁸ Pa: (a) bifurcation diagram, (b) lateral vibration three-dimensional spectrum, and (c) torsional vibration three-dimensional spectrum.

According to the phenomenon shown in Figure 9(b), there are obvious frequencies $f_{l}$ and $f_{h}$ which represent the rotation speed of rotors and multiple frequency. Moreover, the rotation speed region of multiple frequency appearing is approximately 6016 rpm–7353 rpm, which corresponds to the interval where the quasi-periodic motion appears for the first time in Figure 9(a). Figure 9(c) shows that different combined frequencies occur, including $2 f_{l}$ , $2 f_{h}$ , $f_{h} - f_{l}$ , $f_{l} + f_{h}$ , and $4 f_{l}$ . It can be observed that torsional vibration response is more complicated.

To further reveal the nonlinear responses of inter-shaft rub-impact, the axis orbit, time-domain waveform, lateral frequency spectrum, Poincare map of disk 4 lateral displacement, Poincare map of LPR torsion angle, and rubbing force at different speeds are obtained as shown in Figure 10. The rotor whirl status could be judged by comparing the order of the red point, corresponding to vibration at time $t$ , and cyan point, corresponding to vibration at time $t + 3 Δ t$ , marked in the orbits.

Figure 10.

Response of disk 4 at different rotation speeds for $μ$ = 0.1 and $E_{c}$ = 3.5 × 10⁸ Pa: (a) $Ω_{l}$ = 5252 rpm, (b) $Ω_{l}$ = 7162 rpm, (c) $Ω_{l}$ = 9550 rpm, and (d) $Ω_{l}$ = 12892 rpm.

It can be seen from Figure 10(a) that the whirl orbit is regular and waveform approximates harmonic curves, and rotation frequencies $f_{l}$ and $f_{h}$ appear in the frequency spectrum. This phenomenon corresponds to the three-dimensional spectrum and demonstrates the interaction of rotors because of rub-impact. Moreover, five points come forth in the lateral displacement Poincare map of Figure 10(a), which indicates that the LPR undergoes a 5T-periodic motion at $Ω_{l}$ = 5252 rpm. When $Ω_{l}$ = 7162 rpm, 3 $f_{l}$ can be observed in spectrum plot and its amplitude is relatively small, which indicates that the time-varying rubbing forces is much smaller than the unbalanced force. Additionally, a closed circle appears in the lateral displacement Poincare map, which indicates that the system enters quasi-periodic motion. At $Ω_{l}$ = 9550 rpm, Figure 10(c) shows that some repetitive cross curves appear in orbits and five separate points are observed in the lateral displacement Poincare map. Therefore, it can be concluded that the rotor system is 5T-periodic motion. Figure 10(d) shows that the dual-rotor response remains 5T-periodic motion at $Ω_{l}$ = 12892 rpm. And the system rotates counterclockwise, that is, it is at FW status. Besides, it can be seen that the system is 5T-periodic at $Ω_{l}$ = 5252 rpm, $Ω_{l}$ = 9550 rpm, and $Ω_{l}$ = 12892 rpm; when $Ω_{l}$ = 7162 rpm, the system is quasi-periodic based on the torsional vibration Poincare maps. Compared with the lateral vibration results, it is proved that the torsional vibration shows the same periodic motion characteristics with lateral vibration.

According to the simulation results above, it can be seen that the system response experiences period-five motion and quasi-periodic motion in the speed region 5252 rpm< $Ω_{l}$ <13752 rpm. Moreover, it can be concluded that the appearance of 3 $f_{l}$ will lead to quasi-periodic motion while the amplitude of frequency $2 f_{l} + f_{h}$ is relatively more obvious.

Effect of friction coefficient

In this section, the research conducted focused on the role of the friction coefficient on the dual-rotor system dynamic properties, the displacement curve, bifurcation diagram, and three-dimensional spectrum under friction coefficient $μ$ = 0.3 are obtained and shown in Figure 11 and Figure 12, respectively.

Figure 11.

The relative displacement of rubbing nodes for $μ$ = 0.3 and $E_{c}$ = 3.5 × 10⁸ Pa.

Figure 12.

Dual-rotor system dynamic response for $μ$ = 0.3 and $E_{c}$ = 3.5 × 10⁸ Pa: (a) bifurcation diagram, (b) lateral vibration three-dimensional spectrum, and (c) torsional vibration three-dimensional spectrum.

As shown in Figures 8 and 11, there is no fundamental difference in rubbing types and discipline of change exhibited by two diagrams. A noticeable discrepancy is that the maximum $δ_{\max}$ = 0.108 mm occurs at $Ω_{l}$ = 6971 rpm in Figure 8, while the maximum $δ_{\max}$ = 0.209 mm occurs at $Ω_{l}$ = 13465 rpm in Figure 11. The maximum vibration amplitude ascends prominently. The interval widths of continuous rubbing close to the second-order critical speed of LPR are also widened.

Compared with Figure 9(a), it can be found from Figure 12 that the system response evolution process is approximately similar to that of $μ$ = 0.1 when $Ω_{l}$ <7353 rpm. However, it should be noted that the rotational speed ending quasi-periodic motion, $Ω_{l}$ = 7831 rpm, is different. When 10791 rpm< $Ω_{l}$ <13370 rpm, quasi-periodic motion occurs again. The interval of quasi-periodic motion corresponds to the speed region where continuous rub-impact occurs in Figure 11. These results show that higher friction coefficient will widen the system rotating speed region for appearing the quasi-periodic motion obviously.

It can be observed from Figure 12(b) that the rotors rotation frequencies $f_{l}$ and $f_{h}$ and combined frequency 3 $f_{l}$ appear in three-dimensional spectrum. Meanwhile, the FW natural frequency $f_{n 1 +}$ and $f_{h} - 2 f_{n 1 +}$ can be found in speed region 10790–13273 rpm. The peaks of $f_{n 1 +}$ are higher than that of $f_{l}$ , indicating that the natural frequency plays a leading role in dual-rotor system and inter-shaft rub-impact induces the self-excited rubbing motion. By comparing with the response under lower friction coefficient, it can be found that the amplitude of disk 4 significantly increases, suggesting that the increase of friction coefficient could aggravate the vibration level. In addition, it can note that the speed region $f_{n 1 +}$ occurring corresponds to the continuous rub-impact region in Figure 11 and quasi-periodic motion region in Figure 12(a), indicating that $f_{n 1 +}$ may be the characteristic frequency that the above-mentioned rotor motion emerges. The reason for the difference may lie in which the speed interval of bifurcation diagram is 95.5 rpm, and that of three-dimensional spectrum is 190 rpm. Figure 12(c) shows richer characteristic frequencies than Figure 12(b) and Figure 9(c), including $2 f_{l}$ , $f_{l} + f_{h}$ , $f_{h} - f_{n 1 +}$ , 2 $f_{n 1 +}$ , $f_{l} + f_{n 1 +}$ , and $f_{h} + f_{n 1 +}$ . This phenomenon suggests that with the increase of friction coefficient, the more abundant frequency contents can be observed, which means the rub-impact is more severe.

To evaluate the effect of the friction coefficient in-depth, the axis orbit, time-domain waveform, lateral frequency spectrum, Poincare map of disk 4 lateral displacement Poincare map of LPR torsion angle, and rubbing force are obtained as shown in Figure 13. At $Ω_{l}$ = 12892 rpm, it can note that the orbits show an ellipse composed of complex lines and the frequency spectrum consists of frequency $f_{l}$ , $f_{h}$ , $f_{n 1 +}$ , and $f_{h} - 2 f_{n 1 +}$ . Meanwhile, one can observe that the rotation direction of disk 4 is clockwise, revealing that the LPR is at the backward whirl (BW) status. Figure 11 shows that the system is continuous rub-impact at $Ω_{l}$ = 12892 rpm. In addition, it can be seen that the amplitude of other frequency components is far less than that of $f_{n +}$ , indicating that the natural mode dominates this BW motion.

Figure 13.

Response of disk 4 at different rotation speeds for $μ$ = 0.3 and $E_{c}$ = 3.5 × 10⁸ Pa: (a) $Ω_{l}$ = 5252 rpm, (b) $Ω_{l}$ = 7162 rpm, (c) $Ω_{l}$ = 9550 rpm, and (d) $Ω_{l}$ = 12892 rpm.

Effect of elastic modulus of the coating painted on seal

According to equation (13), the materials of coatings painted on HPR will significantly affect its stiffness and further affect the rub-impact response. For $E_{c}$ = 1.1 × 10⁹ Pa and $μ$ = 0.3, the displacement curve, bifurcation diagram, and three-dimensional spectrum of disk 4 with rub-impact are shown in Figures 14 and 15.

Figure 14.

The relative displacement of rubbing nodes for $μ$ = 0.3 and $E_{c}$ = 1.1 × 10⁹ Pa.

Figure 15.

Dual-rotor system dynamic response for $μ$ = 0.3 and $E_{c}$ = 1.1 × 10⁹ Pa: (a) bifurcation diagram, (b) lateral vibration three-dimensional spectrum, and (c) torsional vibration three-dimensional spectrum.

Compared with Figure 11, when the elastic modulus increases under the same friction coefficient, the rubbing types remain similar; thus, intermittent rub-impact and continuous rub-impact appear alternately. However, there are obvious differences between two diagrams in speed region of rubbing motion. The dual-rotor system first enters continuous rubbing at $Ω_{l}$ = 5921 rpm, which is 382 rpm larger than the speed exhibited in Figure 11. Meanwhile, a relatively narrow continuous rubbing speed region 10982 rpm< $Ω_{l}$ <11555.5 rpm can be observed in Figure 14.

As shown in Figure 15, the dual-rotor system nonlinear characteristics become more apparent. At $Ω_{l}$ = 7162 rpm, a jump phenomenon can be observed in the response of disk 4. When 10982 rpm< $Ω_{l}$ <11555 rpm, the system motion is quasi-periodic; when $Ω_{l}$ = 10218 rpm and $Ω_{l}$ = 11842 rpm, the system is in chaotic motion as shown in Poincare diagrams in Figure 14(a), which indicates that when 8404 rpm< $Ω_{l}$ <10982 rpm or 11555 rpm< $Ω_{l}$ <14038 rpm, chaotic periodic motion and quasi-periodic motion appear alternately. Compared with Figure 14, the interval of quasi-periodic corresponds to continuous rub-impact region and the interval of Chaos and Q-P corresponds to intermittent rub-impact region.

From the simulation results, we can see that abundant combined frequencies occur in Figure 15(b), including $f_{h} - 2 f_{n 1 +}$ , $2 f_{n 1 +} - f_{h}$ , $f_{n 1 +} - f_{h} + f_{l}$ , $f_{h} + f_{l} - f_{n 1 +}$ , and $2 f_{h} - f_{n 1 +}$ , $3 f_{l}$ . Moreover, one can also find two speed regions [8498–10982] rpm and [11555–14037] rpm where the combined frequencies of natural frequency and rotation frequencies manifest, which can correspond well to the rotating speed region identified in Figures 14 and 15(a). The appearance of these frequencies may indicate that the system enters intermittent rub-impact motion and Chaos and Q-P motion. It can note from Figure 15(c) that some other frequency contents appear, such as $f_{h} - f_{n 1 +}$ , $2 f_{n 1 +}$ , $f_{l} + f_{n 1 +}$ , $f_{h} + f_{n 1 +}$ , and $f_{l} + f_{h}$ . Compared with simulation results in the case of $E_{c}$ = 3.5 × 10⁸ Pa, one can see that the forms of combined frequencies of natural frequency and rotating frequencies are more diverse. These behaviors show that the inter-shaft rubbing responses are sensitive to the elastic modulus of the seal.

To further explore the dual-rotor system nonlinear characteristics, the axis orbit, time-domain waveform, lateral frequency spectrum, Poincare map of disk 4 lateral displacement Poincare map of LPR torsion angle, and rubbing force are extracted. It can be seen from Figure 16(a) that a circle consisting of five intersecting curves appears in the axis orbit, and five points appear in Poincare diagram, indicating that the rotor of LPR is 5T-periodic motion. At $Ω_{l}$ = 6685 rpm, a noisy one-circle orbit can be observed in orbit, and the system response is quasi-periodic. This conclusion is drawn because a closed loop formed by the points presented in the Poincare map, as shown in Figure 16(b). At $Ω_{l}$ = 9550 rpm and $Ω_{l}$ = 12892 rpm, the orbit and time-domain waveform are irregular and the responses of disk 4 exhibit very complicated features. The displacements exhibited in Figure 14 indicate that the system is intermittent rub-impact motion at these two speeds. The whirl direction is clockwise, indicating the system is at BW status. Meanwhile, many frequency components shown in frequency spectrum validate the analysis correctness about three-dimensional spectrum. It can be found from Figure 16(c) that some smooth circles constructed by discrete points come forth, meaning that the system motion is quasi-periodic. In addition, Figure 16(d) shows the system enters quasi-periodic, which can be determined by one closed petal curve.

Figure 16.

Response of disk 4 at different rotation speed for $μ$ = 0.3 and $E_{c}$ = 1.1 × 10⁹ Pa: (a) $Ω_{l}$ = 5252 rpm, (b) $Ω_{l}$ = 7162 rpm, (c) $Ω_{l}$ = 9550 rpm, and (d) $Ω_{l}$ = 12892 rpm.

Conclusion

In this study, based on the Timoshenko beam element model including torsional degrees of freedom, the dual-rotor system dynamics model under nonlinear inter-shaft rubbing force has been established by FEM after a series of simplifications. Considering the limitation of the point impact model in analyzing cylindrical contact, the cylindrical contact model and the Coulomb model are used to form the expression of rub-impact force. The dual-rotor system dynamic equation has been numerically solved using the Newmark-β method. The Campbell diagram, modal shape, three-dimensional spectrum, whirl orbit, etc., have been obtained according to the simulation results under different rubbing parameters. Then, the rotor system nonlinear dynamic behaviors have been investigated. Some conclusions are summarized and listed as follows:

(1) The inter-shaft rub-impact has a positive effect on the increase of critical speed and will lead to the coupled vibration of the rotors according to the comparison results of the Campbell diagrams and the mode shapes. Meanwhile, the LPR/HPR vibration can be excited by the vibration of HPR/LPR, which reveals the behavior of energy transfer between the two rotors by inter-shaft rub-impact.

(2) It can be concluded from the displacement diagrams that when the dual-rotor system suffers from inter-shaft rubbing force, the intermittent rub-impact motion, and continuous rub-impact can be observed.

(3) A more unstable system motion form will be emerged under higher friction coefficient and higher coating elastic modulus. As the rotation speed gradually increases in a certain range, the typical nonlinear dynamic behaviors including multi-periodic, quasi-periodic, and chaotic motions are revealed. Moreover, the periodicity of the system always corresponds to specific rubbing state and characteristic frequencies.

(4) The characteristic frequency of torsional vibration is not exactly likely that of bending vibration, and they can be combined with each other for inter-shaft rub-impact fault diagnosis. And the torsional vibration shows the same periodic motion characteristics with the lateral vibration.

(5) The friction coefficient and rubbing stiffness can significantly affect the dual-rotor system response, including stimulating the reverse whirling state and making the rub-impact responses complicated. Many frequencies can be observed in three-dimensional spectrums including rotation frequencies, natural frequency, and combination frequencies.

List of symbols.

$w_{x 1}$ , $w_{y 1}$ , $w_{x 2}$ , and $w_{y 2}$	Translations of shaft element node
$θ_{x 1}$ , $θ_{y 1}$ , $θ_{z 1}$ , $θ_{x 2}$ , and $θ_{y 2}$ , $θ_{z 2}$	Rotations of shaft element node
$q_{e}$	Generalized displacement vector of shaft element
$D$ and $d$	Outer diameter and inner diameter of shaft element
$l$ , $ν$ , and $ρ$	Length, Poisson ratio, and density of shaft element
$E$ and $G$	Elastic modulus and shear modulus of shaft element
$I$	Second moment of area
$r_{g}$	Radius of gyration
$ϕ$	Shear deformation coefficient
$M_{e}$ , $K_{e}$ , and $G_{e}$	Mass, stiffness, and gyroscopic matrices of shaft element
$M_{d}$ and $G_{d}$	Mass and gyroscopic matrices of disk
$m_{d}$ , $J_{d}$ , and $J_{p}$	Mass, diametric inertial moment, and polar inertial moment of disk
$K_{s}$ and $C_{s}$	Stiffness and damping matrices of support
$k_{s}$ and $c_{s}$	Support stiffness and support damping
$M$ , $K$ , $C$ , and $G$	Mass, stiffness, damping, and gyroscopic matrices of rotor system
$F_{n}$	Radial impact force
$δ_{0}$ , $δ$ , and $δ^{*}$	Initial clearance, axial distance, and penetration between seal and HPR
$E_{b}$	Elastic modulus of LPR and HPR
$R_{l}$ and $R_{h}$	Curvature radius of LPR and HPR
$ν_{1}$ and $ν_{2}$	Poisson ratio of coatings painted on HPR and seal
$E$	Composite modulus of two colliding cylinders
$H$ and $s i g n$	Heaviside function and symbolic function
$x_{l}$ and $y_{l}$	Lateral displacement and vertical displacement of LPR
$x_{h}$ and $y_{h}$	Lateral displacement and vertical displacement of HPR
$μ$	Frictional coefficient
$v_{r e l}$	Relative speed
$\dot{x}$ and $\dot{y}$	Lateral velocity and vertical velocity
$Ω_{l}$ and $Ω_{h}$	Whirl speed of LPR and whirl speed of HPR
$ω_{l}$ and $ω_{h}$	Rotation speed of LPR and rotation speed of HPR
$r_{l}$ and $r_{h}$	Complex displacement of LPR and complex displacement of HPR
$F_{u}$ and $F_{r}$	Unbalance force vector and rubbing force vector
$m$ and $e$	Mass and eccentricity of disk
$Δ t$	Time step
$E_{c}$	Elastic modulus of coatings of seal

Highlights

• Inter-shaft rub-impact will occur due to the use of sealing structures.

• The FEM model of dual-rotor system is established based on the Timoshenko beam considering torsional deformations.

• The nonlinear cylindrical collision model is applied to describe the impact forces.

• The coupled bending-torsional vibration responses of the dual-rotor system with inter-shaft rub-impact are investigated.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

This work has been carried out with the supports of the National Science and Technology Major Project (J2019-II-0011-0031) and the National Science and Technology Major Project (J2019-II-0016-0037).

References

Xi-guan

Gui-huo

Tang

, et al. Modeling method and dynamic characteristics of high-dimensional counter-rotating dual rotor system. J Aero Power 2014; 29: 585–595.

Bing-long

Chao

Kang

. Coupling vibration characteristics analysis and experiment of shared support-rotors system. J Aero Power 2020; 35: 2293–2305.

Zeng

Zhao

, et al. Rubbing dynamic characteristics of the blisk-casing system with elastic supports. Aero Sci Technol 2019.

Chang-Jian

Chen

. Non-linear dynamic analysis of rub-impact rotor supported by turbulent journal bearings with non-linear suspension. Int J Mech Sci 2008; 50: 1090–1113. DOI: 10.1016/j.ijmecsci.2008.02.003.

Jin

Huang

, et al. The applications of POD method in dual rotor-bearing systems with coupling misalignment. Mech Syst Signal Process 2021; 150: 107236. DOI: 10.1016/j.ymssp.2020.107236.

Cheng

Zhang

, et al. Mechanical characteristics and nonlinear dynamic response analysis of rotor-bearing-coupling system. Appl Math Model 2021; 93: 708–727. DOI: 10.1016/j.apm.2020.12.041.

Wang

Hong

Jiang

. Characteristics of stick-slip oscillations in dry friction backward whirl of piecewise smooth rotor/stator rubbing systems. Mech Syst Signal Process 2020; 135: 106387. DOI: 10.1016/j.ymssp.2019.106387.

Maharathi

Dash

Beheraf

. Dynamic behavior analysis of a dual-rotor system using the transfer matrix method. Int J Acoust Vib 2004; 9: 115–128.

Zhong

Chen

, et al. Nonlinear response analysis for a dual-rotor system supported by ball bearing. Int J Non Lin Mech 2021; 128: 103627. DOI: 10.1016/j.ijnonlinmec.2020.103627.

10.

Chen

. A coupling dynamic model for whole aero-engine vibration and its verification. J Aero Power 2012; 27: 241–254.

11.

, et al. Vibration response analysis of a rotational shaft-disk-blade system with blade-tip rubbing. Int J Mech Sci 2016; 107: 110–125. DOI: 10.1016/j.ijmecsci.2015.12.026.

12.

Kang

Cao

Hou

, et al. Analysis of backward whirling characteristics of a dual-rotor system caused by unbalance. Measurement 2022; 203: 111982. DOI: 10.1016/j.measurement.2022.111982.

13.

Wang

Liu

Jiang

, et al. Casing vibration response prediction of dual-rotor-blade-casing system with blade-casing rubbing. Mech Syst Signal Process 2019; 118: 61–77. DOI: 10.1016/j.ymssp.2018.08.029.

14.

Zhong

Chen

, et al. Modeling and dynamic characteristics analysis of blade-disk dual-rotor system. Complexity 2020; 2020: 1–13. DOI: 10.1155/2020/2493169.

15.

Zhai

Wang

, et al. Unbalance vibration characteristics and sensitivity analysis of the dual-rotor system in aeroengines. J Aero Eng 2021; 34. DOI: 10.1061/(asce)as.1943-5525.0001197.

16.

. Bifurcation and chaos in a rub-impact Jeffcott rotor system. J Sound Vib 1998; 201: 1–18.

17.

Patel

Zuo

Zhao

. Nonlinear lateral-torsional coupled motion of a rotor contacting a viscoelastically suspended stator. Nonlinear Dynam 2012; 69: 325–339.

18.

Yang

Cao

. Study on imbalance-fixed point rubbing coupling faults of dual-disc cantilever rotor. J Aero Power 2016; 31: 307–316.

19.

Herisanu

Marinca

. Analytical study of nonlinear vibration in a rub-impact jeffcott rotor. Energies 2021; 14: 8298. DOI: 10.3390/en14248298.

20.

Zhang

Yang

, et al. A novel diagnosis indicator for rub-impact of rotor system via energy method. Mech Syst Signal Process 2023; 185: 109825. DOI: 10.1016/j.ymssp.2022.109825.

21.

Muszynska

. Rotor-to-stationary element rub-related vibration phenomena in rotating machinery -- literature suryey. Shock Vib Digest 1989; 21: 3–11.

22.

Jacquet-Richardet

Torkhani

Cartraud

, et al. Rotor to stator contacts in turbomachines. Review and application. Mech Syst Signal Process 2013; 40: 401–420. DOI: 10.1016/j.ymssp.2013.05.010.

23.

Jiang

Chen

. Advances in the research on nonlinear phenomena in rotor/stator rubbing systems. Advances in mechanics 2013; 43: 132–148.

24.

Yin

Guo

, et al. A review on dynamic characteristics of blade-casing rubbing. Nonlinear Dynam 2016; 84: 437–472. DOI: 10.1007/s11071-015-2535-x.

25.

Yang

Cao

, et al. Prediction of dynamic characteristics of a dual-rotor system with fixed point rubbing-Theoretical analysis and experimental study. Int J Mech Sci 2016; 115-116: 253–261. DOI: 10.1016/j.ijmecsci.2016.07.002.

26.

Yang

Cao

Wang

, et al. Fixed-point rubbing characteristic analysis of a dual-rotor system based on the Lankarani-Nikravesh model. Mech Mach Theor 2016; 103: 202–221. DOI: 10.1016/j.mechmachtheory.2016.05.007.

27.

Tiancheng

Shuqian

. Analysis and experiment of coupled bending and torsional vibration of a rub-impact dual-rotor system. J Aero Power 2019; 34: 643–655.

28.

Jin

Liu

Yang

, et al. Nonlinear vibrations of a dual-rotor-bearing-coupling misalignment system with blade-casing rubbing. J Sound Vib 2021; 497: 115948. DOI: 10.1016/j.jsv.2021.115948.

29.

Jin

Liu

. Fault simulation of blade-casing rubbing for dual-rotor system of aero-engines. Acta Aeronautica Astronautica Sinica 2022; 43: 498–511.

30.

Liu

Fei

, et al. Nonlinear vibration response of a complex aeroengine under the rubbing fault. Nonlinear Dynam 2021; 106: 1869–1890. DOI: 10.1007/s11071-021-06717-4.

31.

Shan

Zhai

, et al. Vibration characteristics and simulation verification of the dual-rotor system for aeroengines with rub-impact coupling faults. Shock Vib 2021; 2021: 1–24. DOI: 10.1155/2021/6622065.

32.

Prabith

Krishna

IRP

. Response and stability analysis of a two-spool aero-engine rotor system undergoing multi-disk rub-impact. Int J Mech Sci 2022; 213. DOI: 10.1016/j.ijmecsci.2021.106861.

33.

Chen

Hong

. Structure of gas turbine engine. Beijing: Beihang University Press, 2010.

34.

Song-ting

Zhi-qiang

. Rubbing vibration analysis for a counter-rotating dual-rotor systemv. J Vib Shock 2012; 43: 498–511.

35.

Wang

Hou

, et al. Dynamic characteristics of an aeroengine dual-rotor system with inter-shaft rub-impact. Mech Syst Signal Process 2022; 166: 108475. DOI: 10.1016/j.ymssp.2021.108475.

36.

Hou

Wang

, et al. Insights into the nonlinear behaviors of dual-rotor systems with inter-shaft rub-impact under co-rotation and counter-rotation conditions. Int J Non Lin Mech 2022; 140: 103901. DOI: 10.1016/j.ijnonlinmec.2021.103901.

37.

Chen

Gao

, et al. The coupled bending-torsional dynamic behavior in the rotating machinery: modeling, simulation and experiment validation. Mech Syst Signal Process 2022; 178: 109306. DOI: 10.1016/j.ymssp.2022.109306.

38.

da Silva

Marques

da Silva

, et al. A compendium of contact force models inspired by Hunt and Crossley's cornerstone work. Mech Mach Theor 2022; 167: 20. DOI: 10.1016/j.mechmachtheory.2021.104501.

39.

Pereira

Ramalho

Ambrosio

. Applicability domain of internal cylindrical contact force models. Mech Mach Theor 2014; 78: 141–157. DOI: 10.1016/j.mechmachtheory.2014.03.010.

40.

Liu

Zhang

Yang

. The FEM analysis and approximate model for cylindrical joints with clearances. Mech Mach Theor 2007; 42: 183–197. DOI: 10.1016/j.mechmachtheory.2006.02.006.