Research on modeling and dynamic characteristics of complex coaxial rotor system

Fei Wang1 , Gui-Huo Luo2 , Xi-Guan Yang3 , Hai-Tao Cui4

1, 2, 3, 4Nanjing University of Aeronautics and Astronautics, Jiangsu Province Key Laboratory of Aerospace Power System, Nanjing 210016, China

1Corresponding author

Journal of Vibroengineering, Vol. 19, Issue 3, 2017, p. 1524-1545. https://doi.org/10.21595/jve.2016.17324
Received 23 June 2016; received in revised form 30 October 2016; accepted 21 November 2016; published 15 May 2017

Copyright © 2017 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License
Table of Contents Download PDF References
Cite this article
Views 101
Reads 35
Downloads 839
Abstract.

To make up deficiency of the finite element method in predicting nonlinear dynamic characteristics of coaxial rotor systems, nonlinear dynamic model of a coaxial rotor system was established with a method combining the finite element method and the fixed interface modal synthesis method. Then an implicit time domain method was presented to solve the nonlinear equations of motion thus dynamic characteristics of the rotor system can be obtained. The computational efficiency of this method largely depends on the number of degrees of freedom with nonlinear forces acting on. With nonlinear forces of squeeze film damper and intermediate bearing considered, nonlinear dynamic response characteristics of the coaxial rotor system under multiple unbalance forces were studied in this work. The results showed that the unbalance excitation frequencies are dominant in the responses of the rotor system. Besides, due to coupling effect of the intermediate bearing some combinations of the unbalance excitation frequencies were also observed in the spectrogram. Stability and periodicity of the rotor system was investigated with bifurcation diagram, Poincare map and phase diagram. It was found that the rotor system executes multiple periods orbital motion under relatively low rotational speeds. With the increasing of rotational speed, the rotor system would execute quasi-periodic motion, chaotic motion and periodic motion again. The quasi-periodic motion and chaotic motion are closely related with the SFD. Finally, under relatively low speed, the nonlinear model was validated by comparing the simulation results with the experimental data. The proposed modeling and solving method is expected to provide theoretical and engineering basis for improving prediction of nonlinear dynamic characteristics of complex rotor systems.

Keywords: nonlinear vibration, counter-rotating coaxial rotor system, intermediate bearing, modal synthesis method, finite element method.

1. Introduction

Rotor system is the main source of vibration of aeroengines. Stability and reliability of aeroengines mostly depend on dynamic behaviors of the rotor system. In recent years, counter-rotating coaxial system has been applied to many aeroengines, such as RB211, GE120, F119 and so on. The counter-rotating technology is beneficial to improve the fuel consumption rate and thrust-weight ratio as well as reduce the gyroscopic torque of aeroengines [1]. Although advantageous for many aspects, it has made dynamic behaviors of the rotor system more complex, especially when coupled with squeeze film dampers(SFD) and intermediate bearing, both of which have been widely applied in modern aeroengines. Surveys of the dynamic issues associated with counter-rotating dual-rotor system are listed in references [2, 3].

The most commonly recurring problems in rotor dynamics are excessive steady state synchronous vibration levels and subsynchronous rotor instabilities. The first problem is caused mostly by the synchronous unbalance force and, generally, the second one is caused by nonlinear forces exerted by nonlinear components such as SFD and bearing. The steady state synchronous vibration levels may be reduced by improved balancing, or by moving critical speeds of the system out of the operating range [4], or by introducing external damping to limit peak amplitudes at critical speeds [5]. Subsynchronous rotor instabilities can be avoided by rising the natural frequency of the rotor system, or by introducing damping to increase the onset rotor speed of instability [6]. However, all these problem solving measures are based on accurate approximation of the actual rotor system.

To better approximate the rotor system, all the forces, linear and nonlinear, should be considered in the dynamic model. However, the problem can be quite complicated when nonlinear forces are considered. The problem of nonlinear response of rotor systems are discussed in the extensive investigations of references [7-15]. Glasgow D. A. [8] employed the fixed interface modal synthesis method to reduce the size of the rotor system to obtain better computational efficiency. In reference [11, 12], nonlinear response of a dual-rotor system with intermediate bearing was studied with transfer matrix method. Wang W. [14] improved the free interface modal synthesis method to analyze dynamic characteristics of a multi-rotor system. In 2004, Hsiao-Wei [15] established the mathematical model of a dual-rotor system with the finite element method to investigate the influences of speed ratio of the rotor system and stiffness of the intermediate bearing. From all the above references it can be seen that the transfer matrix method, finite element method and component modal synthesis method are commonly used in rotor dynamic analysis. In present work, nonlinear model of a counter-rotating coaxial test rig is established with a method combining the finite element method and the component modal synthesis method for accuracy and efficiency.

Notwithstanding the many investigations that there have been of the dynamic behavior of rotor systems with nonlinear forces acting on, the problem of developing efficient method for analyzing nonlinear dynamic behaviors of rotor system still needs to be solved. All these facts necessitate that more detailed models should be established and more efficient and robust numerical methods should be developed to give better approximation of dynamic behaviors of rotor systems.

In this paper, a counter-rotating coaxial system was modeled with finite element method and fixed interface modal synthesis method. Nonlinear forces exerted by the squeeze film damper and the intermediate bearing were considered. With this nonlinear model and the improved Newmark-β method, dynamic characteristics of the rotor system was obtained and then further analysis was made and compared with the experimental results.

2. Model development

2.1. Finite element formulation

Fig. 1 shows the rotor element considered in this work. Each element consists of two nodes, with four degrees of freedom at each node, namely the lateral displacement and the slope. The lateral displacements in X and Y directions are denoted with u and v respectively. The slopes are denoted with θ and φ for angular displacements in X and Y directions respectively.

The nodal displacement vector can be described as:

(1)
q e = u e 1 v e 1 θ e 1 φ e 1 u e 2 v e 2 θ e 2 φ e 2 T .

The displacement vector within the element can be interpolated as:

(2)
u e ξ ξ , t v e ξ ξ , t θ e ξ ξ , t φ e ξ ξ , t = N 1 0 0 N 2 N 3 0 0 N 4 0 N 1 - N 2 0 0 N 3 - N 4 0 0 - N ˇ 1 N ˇ 2 0 0 - N ˇ 3 N ˇ 4 0 N ˇ 1 0 0 N ˇ 2 N ˇ 3 0 0 N ˇ 4 q e .

The shape functions in Eq. (2) are:

(3a)
N 1 = 1 - 1 l e 2 + 12 g 12 g l e ξ + 3 ξ 2 - 2 l e ξ 3 , N 2 = 1 l e 2 + 12 g l e 2 + 6 g ξ - 2 l e 2 + 6 g l e ξ 2 + ξ 3 , N 3 = 1 l e 2 + 12 g 12 g l e ξ + 3 ξ 2 - 2 l e ξ 3 , N 4 = 1 l e 2 + 12 g - 6 g ξ - l e 2 - 6 g l e ξ 2 + ξ 3 ,
(3b)
N ˇ 1 = 1 l e 2 + 12 g 6 l e ξ 2 - 6 ξ , N ˇ 2 = 1 l e 2 + 12 g l e 2 + 12 g - 4 l e 2 + 12 g l e ξ + 3 ξ 2 , N ˇ 3 = 1 l e 2 + 12 g 6 ξ - 6 l e ξ 2 , N ˇ 4 = 1 l e 2 + 12 g 12 g - 2 l e 2 l e ξ + 3 ξ 2 .

Fig. 1. Diagram of beam element

 Diagram of beam element

In Eq. (3), le is the element length and g is formulated as given below:

(4)
g = E e I e κ G e A e ,

where Ee the Young’s Modulus, Ie the area moment of inertia, Ge the shear modulus, Ae the cross-sectional area. κ is the shear coefficient, and for a hollow circular shaft section [16, 17]:

(5)
κ = 6 1 + μ 2 1 + λ 2 2 7 + 12 μ + 4 μ 2 1 + λ 2 2 + 4 5 + 6 μ + 2 μ 2 λ 2 ,

where μ is the Poisson’s ratio, λ is the ratio of the inner shaft radius to the outer shaft radius. Hence for a solid shaft λ0.

With rotary inertia and shear effects considered, the kinetic energy and strain energy for a shaft element is:

(6a)
T e = 1 2 0 l e ρ e A e u ˙ e ξ 2 ξ , t + ρ e I e φ ˙ e ξ 2 ξ , t d ξ ,
(6b)
U e = 1 2 0 l e φ e ξ ξ , t ξ T E e I e φ e ξ ξ , t ξ d ξ + 1 2 κ 0 l e γ T G e A e γ d ξ ,

where ρe denotes the density of the material and γ is the transverse shear strain [18].

And thus the element mass matrix:

(7)
m e = ρ e A e 0 l e N T N d ξ + ρ e I e 0 l e g N ' ' ' + N ' T g N ' ' ' + N ' d ξ .

The element stiffness matrix:

(8)
k e = E e I e 0 l e N ' ' T N ' ' d ξ + E e I e g 0 l e N ' ' ' T N ' ' ' d ξ ,

where:

(9)
N ' = d N 1 N 2 N 3 N 4 d ξ , N ' ' = d 2 N 1 N 2 N 3 N 4 d ξ 2 , N ' ' ' = d 3 N 1 N 2 N 3 N 4 d ξ 3 .

The element gyroscopic matrix:

(10)
g e = - 2 ρ e I e 0 l e B 2 T B 1 - B 1 T B 2 d ξ ,

where:

(11)
B 1 B 2 = 0 - N ˇ 1 N ˇ 2 0 0 - N ˇ 3 N ˇ 4 0 N ˇ 1 0 0 N ˇ 2 N ˇ 3 0 0 N ˇ 4 .

For disk element, the mass and gyroscopic matrices are:

(12a)
m e d = m d i 0 0 0 0 m d i 0 0 0 0 I d i 0 0 0 0 I d i ,
(12b)
g e d = 0 0 0 0 0 0 0 0 0 0 0 I p i 0 0 - I p i 0 ,

where mdi, Idi and Ipi are the mass, diametral moment of inertia and polar moment of inertia of the disk, respectively.

2.2. Nonlinear forces of the supports

For squeeze film dampers analyzed in this work, based on the short bearing assumption and the Reynolds boundary conditions, nonlinear forces of the SFD can be expressed as [19, 20]:

(13a)
f x s f d = - μ s R L 3 c 2 x 2 + y 2 1 2 x ε ˙ I 2 + ε ψ ˙ I 1 - y ε ˙ I 1 + ε ψ ˙ I 3 , f y s f d = - μ s R L 3 c 2 x 2 + y 2 1 2 y ε ˙ I 2 + ε ψ ˙ I 1 + x ε ˙ I 1 + ε ψ ˙ I 3 ,
(13b)
ε = x 2 + y 2 / c ,
(13c)
ε ˙ = x x ˙ + y y ˙ / c x 2 + y 2 ,
(13d)
ψ ˙ = y x ˙ - x y ˙ / x 2 + y 2 ,
(13e)
tan ψ = y / x ,

where x and y are the horizontal and vertical displacements of the journal. Ij (j= 1, 2, 3) are Sommerfeld integrals. R and L are radius and length of the SFD respectively. μs and c denote the dynamic viscosity of the film and the radial clearance of the SFD.

As for the rolling ball bearing, based on pure rolling assumption and the Hertz contact theory, nonlinear force of the rolling ball bearing [21] is:

(14a)
f x b = k n j = 1 N b u θ j ξ H u θ j sin θ j , f y b = k n j = 1 N b u θ j ξ H u θ j cos θ j ,
(14b)
θ j = 2 π j - 1 / N b + ω c t ,
(14c)
u θ j = x i r - x o r sin θ j + y i r - y o r cos θ j - γ / 2 ,
(14d)
H u θ j = 0 ,                 u θ j 0 , u θ j ,         u θ j > 0 ,
(14e)
ω c = ω i n r + ω o u t R / R + r ,

where the superscript ir and or in Eq.(14c) denote inner ring and outer ring of the bearing. kn in Eq.(14a) is the contact stiffness between rollers and rings. Nb represent number of the rollers. ξ=3/2 for the rolling ball bearing used in the intermediate support. θj is the rotation angle of the jth roller at time t. γ and ωc are bearing radial clearance and rotational speed of the bearing retainer respectively. uθj is the elastic radial deformation of the jth roller. r and R represent radius of the inner and outer ring. ωin and ωout are the rotational speeds of the inner and outer rings.

2.3. Equations of motion and numerical algorithm

Equation of motion of a nonlinear rotor system can be written as:

(15)
M u ¨ + G u ˙ + K u = F n o n l + F u n b ,

where Funb is the unbalance force vector. Fnonl is the nonlinear forces exerted by the SFDs and rolling bearing in this work. M, K, and G can be obtained quite readily with elements formulated in Section 2.1.

Mathematically, the model is a set of nonlinear second order differential equations. The computational efficiency of which depends on the numerical methods applied and the total DOFs of the rotor system. To combine accuracy of the finite element model and computational efficiency, the fixed interface modal synthesis is applied to reduce dimension of the mathematical model and thus the computational effort.

Eq. (15) can be rewritten as:

(16)
M I I M I J M J I M J J u ¨ I u ¨ J + G I I G I J G J I G J J u ˙ I u ˙ J + K I I K I J K J I K J J u I u J = F I u n b 0 + 0 F J n o n l ,

where uI can be interpreted as linear DOFs, with only unbalance forces acting on, while uJ represents the interface DOFs or nonlinear DOFs in this work, including all the nodes with nonlinear forces acting on.

According to fixed interface modal synthesis method, transformation between physical and modal coordinates is:

(17a)
u I u J = ϕ k ϕ c 0 I q k u J = T q k u J ,

where ϕk – the mass normalized normal mode matrix, ϕc – the mass normalized constrained mode, I – unity matrix, qk – the normal mode coordinates, T – the transformation matrix. ϕc is given as below:

(17b)
ϕ c = - K I I - 1 K I J .

In Eq. (16), let uJ= 0 and neglecting the gyroscopic effects gives:

(17c)
M I I u ¨ I + K I I u I = 0 .

Solving the eigenvalue problem corresponding to Eq. (17c) gives the mass normalized normal mode matrix ϕk and the modal frequency Ωk.

Substitute Eq. (17a) into Eq. (16) and premultiply with TT to obtain the reduced equations of motion:

(18a)
I M - I J M - J I M - J J q ¨ k u ¨ J + G - I I G - I J G - J I G - J J q ˙ k u ˙ J + K - I I 0 0 K - J J q k u J = ϕ k T F I u n b ϕ c T F I u n b + 0 F J n o n l ,

where:

(18b)
M - J J = M J J + M J I ϕ c + ϕ c T M I I ϕ c + M I J , M - I J = M - J I = ϕ k T M I I ϕ c + M I J , G - I I = ϕ k T G I I ϕ k , G - I I = G J J + G J I ϕ c + ϕ c T G I I ϕ c + G I J + c J J , G - I J = G - J I = ϕ k T G I I ϕ c + ϕ k T G I J , K - I I = d i a g Ω r 2         1 r n , K - J J = K I J + K J I ϕ c + k J J ,

where Ωr is the rth modal frequency obtained from Eq. (17c). FIunb is the unbalance force acting on linear DOFs. FJnonl is the nonlinear force acting on nonlinear DOFs.

And:

(19)
F J n o n l = - k J J + F J S F D F J B ,

where FJSFD=f1xsfdf1ysfdfNxsfdfNysfd is the nonlinear force vector exerted by SFDs which can be calculated with Eq. (13). FJB=fxbfyb-fxb-fyb is the nonlinear force vector exerted by rolling bearing which can be calculated with Eq. (14). kJJ and cJJ are given below:

(20)
k J J = k J B 0 0 0 ,           c J J = c J B 0 0 0 ,
(21)
k J B = d i a g k J i B ,         c J B = d i a g c J i B ,         1 i N ,
(22)
k J i B = k J i 0 0 k J i ,           c J i B = c J i 0 0 c J i ,

where kJi and cJi are the stiffness and damping coefficients of the elastic support, respectively. N is the number of elastic supports.

Rearranging Eq. (18) gives:

(23)
q ¨ k + G - I I q ˙ k + K - I I q k = ϕ k T F I u n b - M - I J u ¨ J - G - I J u ˙ J ,
(23b)
M - J J u ¨ J + G - J J u ˙ J + K - J J u J = F J n o n l + ϕ c T F I u n b - M - J I q ¨ k - G - J I q ˙ k .

In Eq. (23), the vector qk and uJ can be interpreted as linear and nonlinear DOFs of the reduced system respectively. Obviously, there is no nonlinear force acting on the linear DOFs. Thus, the explicit Newmark-β method applies to Eq.(23a) while implicit Newmark-β method applies to Eq. (23b).

According to assumptions of Newmark-β method, in time interval tntn+1:

(24a)
q ˙ k n + 1 u ˙ J n + 1 = q ˙ k n u ˙ J n + 1 - β q ¨ k n u ¨ J n + β q ¨ k n + 1 u ¨ J n + 1 Δ t ,
(24b)
q k n + 1 u J n + 1 = q k n u J n + q ˙ k n u ˙ J n Δ t + 0.5 - α q ¨ k n u ¨ J n + α q ¨ k n + 1 u ¨ J n + 1 Δ t 2 ,

where tn+1=tn+Δt and Δt is the time increment. The superscript n and n+1 denote tn and tn+1.

Following equations can be obtained from Eq. (24):

(25)
q ¨ k n + 1 = a q k n + 1 - A q n ,           q ˙ k n + 1 = b q k n + 1 - B q n ,
(26)
u ˙ J n + 1 = a u J n + 1 - A J n ,           u ˙ J n + 1 = b u J n + 1 - B J n .

With:

(27)
A q n = 1 α Δ t 2 q k n + 1 α Δ t q ˙ k n + 1 2 α - 1 q ¨ k n , A J n = 1 α Δ t 2 u k n + 1 α Δ t u ˙ k n + 1 2 α - 1 u ¨ k n , B q n = β α Δ t q k n + β α - 1 q ˙ k n + β α - 2 q ¨ k n , B J n = β α Δ t u k n + β α - 1 u ˙ k n + β α - 2 u ¨ k n , a = 1 α Δ t 2 ,         b = β α Δ t .

Substituting Eq. (24)-Eq. (27) into Eq. (23) yields:

(28a)
q k n + 1 = S q - 1 ϕ k T F I u n b   n + 1 - V q u J n + 1 + W q ,
(28b)
S J - V J S q - 1 V q u J n + 1 = F J n o n l   n + 1 + ϕ c T - V J S q - 1 ϕ k T F I u n b   n + 1 - V J S q - 1 W q + W J .

With:

(29a)
W q = M - I I A q n + G - I I B q n + M - I J A J n + G - I J B J n ,
(29b)
W J = M - J J A J n + G - J J B J n + M - J I A q n + G - J I B q n ,
(29c)
S q = a M - I I + b G - I I + K - I I ,
(29d)
S J = a M - J J + b G - J J + K - J J ,
(29e)
V q = a M - I J + b G - I J ,
(29f)
V J = a M - J I + b G - J I ,
(29g)
F J n o n l   n + 1 = F J n o n l   n u J n + 1 , u ˙ J n + 1 .

Using Eq. (26) to substitute for u˙Jn+1 in Eq. (29g) and substitute the new Eq. (29g) into Eq. (28b) to yield nonlinear equations of uJn+1. After solving Eq. (28b) for uJn+1 with numerical algorithms, qkn+1 can be obtained with uJn+1 substituted into Eq. (28a). Thus u˙Jn+1, u¨Jn+1, q˙kn+1, q¨kn+1 and can be solved from Eq. (25) and Eq. (26). The process continues to repeat and move forward to find uJn+2, and so on. When the iteration is over, the nonlinear response in physical coordinate system can be obtained from Eq. (17a).

Computational efficiency of solving Eq. (23) depends on solving Eq. (23b) while dimensions of Eq. (23b) depend on DOFs with nonlinear forces acting on. Thus, the computational efficiency can be greatly improved with the modeling technique and solving method described in this work.

To summarize, nonlinear model of the rotor system is established with finite element method and fixed interface modal synthesis method; subsequently an implicit time-domain method based on Newmark-β method is applied to solve the equations of motion of the reduced system thus dynamic characteristics can be obtained. The process of modeling and solving is described as follows: First, only the rotating components of the rotor assembly are modeled with the element described in Section 2.1 thus M, K and G in Eq. (15) can be obtained and partitioned as Eq. (16) for calculating ϕc and later use. Then, in accordance with Eq. (17c), modal analysis is conducted on the model with all supports pinned to obtain ϕk and Ωk. The left hand side of Eq. (18a) now is obtained with ϕk, ϕc and the partitioned M, K and G while the right hand side is formulated with Eq. (13), Eq. (14) and Eq. (19)-Eq. (22). Finally, Eq. (18a) is solved with the numerical method which is described by Eq. (23)-Eq. (29). Flowchart of the modeling and solving method is shown in Fig. 2.

Fig. 2. Flowchart of modeling and solving

 Flowchart of modeling and solving

2.4. Test rig description

An overall view of the dual-rotor test rig is given in Fig. 3. The test rig was designed to study the dynamic characteristics of co-and counter-rotating dual-rotor systems with 4 or 5 supports. In this work, the 4 supports counter-rotating dual-rotor system is studied. So, the additional support in Fig. 3 was removed. Thus, the studied rotor system is designed with 4 supports and 4 disks, two for each rotor. The test rig consists of two shafts disposed along the same axis, connected by an intermediate bearing. Each shaft is driven by its own motor thus their spin speeds can be different. Speed ratio of outer and inner rotor is –1.65. Rear support of the outer rotor is intermediate support in which a rolling ball bearing is adopted. Support I and II are front and rear supports of inner rotor respectively while Support III and IV are front and rear support of outer rotor. The squirrel cage + rolling bearing + SFD supporting scheme is adopted for support I, II and III. While support IV, or the intermediate bearing, is mounted without SFD. Four eddy current probes are placed to measure the displacement of disk 2 and 4, two for each disk. It is to be noted that in Fig. 3 the eddy current probes have not been placed in the right positions yet.

The test profiles include slow run-up and run-down through the speed range of the test rig. The angular acceleration for inner rotor during run-up and rundown is approximately 10 rpm/s. Besides, the steady state response are measured when inner rotor spins at 128 rad/s, 160 rad/s and 210 rad/s. The operational ranges for the inner and outer rotor are 0-232 rad/s and 0-382 rad/s respectively. The excitation is only due to the residual unbalance of the rotors.

Fig. 3. Dual-rotor test rig

 Dual-rotor test rig

3. Numerical results and experimental validation

3.1. Model introduction

The counter-rotating coaxial system studied in this paper is presented in Fig. 4.

The stationary coordinate system in Fig. 4 consists of three mutually perpendicular axes, ox, oy, oz, intersecting at the point o and axis of the rotor coincides with axis oz. The axes ox, oz are horizontal, oy is vertical.

Fig. 4. Structural diagram of the coaxial rotor system

 Structural diagram of the coaxial rotor system

Model parameters of the rotor system studied in this work are listed in Table 1-5. The Young’s modulus of the shaft is 196 GPa, mass density 7810 kg/m3, and shear modulus 75.5 GPa. To apply modal synthesis method, 40 modes are retained in the normal mode matrix ϕk. For the Newmark-β method, α= 0.25, β= 0.5.

Geometric dimensions and information of each element are listed in Table 1. Stiffness coefficients of elastic supports are listed in Table 2. Parameters of the intermediate bearing and SFDs are listed in Table 3 and 4. The unbalance configuration and inertia properties of each disk are listed in Table 5.

Table 1. Dimension and elements information of the rotor system

Node no.
Axial location (m)
Bearing/disk
Element no.
Outer diameter (m)
Inner diameter (m)
1
0
1
0.018
0.00
2
0.08143
2
0.018
0.00
3
0.16286
3
0.018
0.00
4
0.24429
4
0.018
0.00
5
0.24909
5
0.018
0.00
6
0.25479
6
0.018
0.00
7
0.28879
7
0.018
0.00
8
0.32279
8
0.018
0.00
9
0.35879
Disk no. 1
9
0.018
0.00
10
0.38369
10
0.018
0.00
11
0.40859
11
0.018
0.00
12
0.43349
12
0.018
0.00
13
0.43869
Bearing no. 1
13
0.018
0.00
14
0.44479
14
0.022
0.00
15
0.54752
15
0.022
0.00
16
0.65025
16
0.022
0.00
17
0.75298
17
0.022
0.00
18
0.85571
18
0.022
0.00
19
0.95844
19
0.022
0.00
20
1.06117
20
0.022
0.00
21
1.06517
Bearing no. 4
21
0.022
0.00
22
1.06867
22
0.022
0.00
23
1.08867
23
0.022
0.00
24
1.10867
Disk no. 2
24
0.022
0.00
25
1.14274
25
0.022
0.00
26
1.17681
26
0.022
0.00
27
1.21088
27
0.017
0.00
28
1.21488
Bearing no. 2
28
0.014
0.00
29
1.21838
29
0.014
0.00
30
1.23038
End of inner rotor
31
0.64200
30
0.035
0.03
32
0.66065
31
0.035
0.03
33
0.67930
Bearing no. 3
32
0.035
0.03
34
0.68650
33
0.038
0.03
35
0.71170
34
0.038
0.03
36
0.73690
35
0.038
0.03
37
0.76210
Disk no. 3
36
0.038
0.03
38
0.80784
37
0.038
0.03
39
0.85358
38
0.038
0.03
40
0.89932
39
0.038
0.03
41
0.94506
40
0.038
0.03
42
0.99080
Disk no. 4
41
0.038
0.03
43
1.01430
42
0.038
0.03
44
1.02030
43
0.070
0.03
45
1.03330
44
0.060
0.03
46
1.06030
45
0.060
0.03
47
1.06430
Bearing no. 4
46
0.060
0.03
48
1.07380

Table 2. Stiffness of elastic supports (squirrel cage)

Support I
Support II
Support III
Support IV
Stiffness (N/m)
1.45×106
2.21×105
9.29×105

Table 3. Parameters of the intermediate bearing

Radius of inner ring (mm)
Radius of outer ring (mm)
No. of rollers
Contact stiffness (N/m3/2)
Radial clearance (μm)
9.37
14.13
9
7.055×109
6

Table 4. Parameters of SFDs

Inner rotor
Outer rotor
SFD I
SFD II
SFD III
Radius R / mm
25
18
35
Length L / mm
15
15
20
Radial clearance c / mm
0.1
0.1
0.08
Dynamic viscosity μs / 10-2Pa·s
1.0752

Table 5. Unbalance configuration and inertia properties of disks

Inner rotor
Outer rotor
Disk 1
Disk 2
Disk 3
Disk 4
Unbalance (×10-5kg·m)
2
4
1
2
Mass (kg)
2.3386
2.3386
3.2590
1.6303
Polar moment of inertia (kg·m2)
0.00815
0.00815
0.01561
0.00661

Fig. 5. Spectrum cascade of the horizontal response of disk 2 and 4-numerical simulation

 Spectrum cascade of the horizontal response of disk 2 and 4-numerical simulation

a) Disk 2

 Spectrum cascade of the horizontal response of disk 2 and 4-numerical simulation

b) Disk 4

3.2. Unbalance response analysis

With the model established by method described in Section 2, steady-state unbalance response of the rotor system is obtained for rotational speed of the inner rotor varying from 4 rad/s to 1600 rad/s with a step length of 2 rad/s. Zero initial condition is adopted for the first step. For the rest steps, result of previous step is adopted as initial condition. Without loss of generality, unbalance response of disk 2 and disk 4 are analyzed to evaluate the coupling effect of the intermediate bearing.

Fig. 6. Spectrum cascade of the horizontal response of disk 2 and 4-experimental results

 Spectrum cascade of the horizontal response of disk 2 and 4-experimental results

a) Disk 2

 Spectrum cascade of the horizontal response of disk 2 and 4-experimental results

b) Disk 4

Spectrum cascades of the horizontal response of disk 2 and disk 4 with inner rotor operating in 4-400 rad/s are shown in Fig. 5 and Fig. 6. The simulation and experimental results are shown in Fig. 5 and Fig. 6 respectively. The numerically simulated spectrum cascades of the horizontal response of disk 2 and disk 4 with inner rotor operating in 400-1600 rad/s are shown in Figs. 7-9. Due to nonlinearities of SFD and the intermediate bearing as well as the multifrequency unbalance excitation of the inner and outer rotor, rich frequency components are observed in the response of disk 2 and disk 4. Following conclusions can be reached from Fig. 5-9:

(1) Responses of the inner and outer rotor are coupled because of the intermediate bearing. Frequency components, ω1 and ω2, corresponding to the unbalance excitation frequencies of the inner and outer rotor are dominant in the response of disk 2 and disk 4. The coupling effect also makes the spectrum cascades corresponding to disk 2 and disk 4 basically follow the same pattern with only slight differences in amplitude.

(2) Frequency components besides ω1 and ω2 in the responses of disk 2 and 4, such as 2ω1+ω2, ω1+2ω2, 3ω1+2ω2, 2ω1+3ω2, (ω2-ω1)/2, seen in Figs. 5-9, are mainly caused by the nonlinear forces of the SFD and the coupling effect of the intermediate bearing.

(3) By comparison between Fig. 5 and Fig. 7 it can be seen that the new combination frequency components, 2ω2-ω1,3ω2-2ω1emerge when inner rotor operating in 400-800 rad/s. However, 2ω2-ω1 and 3ω2-2ω1 are not observed for both disk 2 and 4 in Fig. 8 and Fig. 9 when the inner rotor operating in 800-1200 rad/s and 1200-1600 rad/s. Besides, 2ω1+ω2 is persistent in 0-1600 rad/s for disk 2 and 4, which can be seen from Fig. 5 and Figs. 7-9.

(4) Two critical speeds of the rotor system, 168 rad/s and 186 rad/s, are observed in rotational speed range 0-200 rad/s in Fig. 5. The first one, 168 rad/s, is excited by the outer rotor and the second one, 186 rad/s, is excited by the inner rotor. It can be obtained from the experimental data shown in Fig. 6 that the corresponding critical speeds of the rotor system are 173 rad/s and 189 rad/s. Discrepancies of critical speeds between numerical and experimental results are less than 3 %, which demonstrate the validity of the model established in this work.

Fig. 7. Spectrum cascade of the horizontal response of disk 2 and 4-numerical simulation

 Spectrum cascade of the horizontal response of disk 2 and 4-numerical simulation

a) Disk 2

 Spectrum cascade of the horizontal response of disk 2 and 4-numerical simulation

b) Disk 4

Fig. 8. Spectrum cascade of the horizontal response of disk 2 and 4-numerical simulation

 Spectrum cascade of the horizontal response of disk 2 and 4-numerical simulation

a) Disk 2

 Spectrum cascade of the horizontal response of disk 2 and 4-numerical simulation

b) Disk 4

Fig. 9. Spectrum cascade of the horizontal response of disk 2 and 4-numerical simulation

 Spectrum cascade of the horizontal response of disk 2 and 4-numerical simulation

a) Disk 2

 Spectrum cascade of the horizontal response of disk 2 and 4-numerical simulation

b) Disk 4

Fig. 10 shows bifurcation diagrams for horizontal displacement of disk 2 and disk 4 with rotor speed. Sampling period for the bifurcation diagram is 5×2π/ω1 because the rotational speed ratio of outer and inner rotor is –1.65 and 5 1.65×3. The horizontal axes of Fig. 10(a) and (b) are rotational speeds of inner and outer rotor respectively. And the vertical axes are horizontal displacement of disk 2 and disk 4.

Basically, disk 2 and 4 execute multiple period orbital motion with inner rotor operating in 4-390 rad/s, see Fig. 10(a), Fig. 10(b). However, with the increase of rotational speed, a transformation between periodic and chaotic motion can be observed. From 400 rad/s to 600 rad/s, the rotor execute chaotic motion which can be seen from Fig. 10(c), Fig. 10(d) and Fig. 11.

Fig. 10. Bifurcation diagrams of the horizontal response of disk 2 and disk 4 with rotor speed

 Bifurcation diagrams of the horizontal response of disk 2 and disk 4 with rotor speed

a) Disk 2

 Bifurcation diagrams of the horizontal response of disk 2 and disk 4 with rotor speed

b) Disk 4

 Bifurcation diagrams of the horizontal response of disk 2 and disk 4 with rotor speed

c) Disk 2

 Bifurcation diagrams of the horizontal response of disk 2 and disk 4 with rotor speed

d) Disk 4

 Bifurcation diagrams of the horizontal response of disk 2 and disk 4 with rotor speed

e) Disk 2

 Bifurcation diagrams of the horizontal response of disk 2 and disk 4 with rotor speed

f) Disk 4

 Bifurcation diagrams of the horizontal response of disk 2 and disk 4 with rotor speed

g) Disk 2

 Bifurcation diagrams of the horizontal response of disk 2 and disk 4 with rotor speed

h) Disk 4

Fig. 11. Disk 2: ω1= 500 rad/s, ω2= –825 rad/s

 Disk 2: ω1= 500 rad/s, ω2= –825 rad/s

a) Spectrogram

 Disk 2: ω1= 500 rad/s, ω2= –825 rad/s

b) Phase trajectory

 Disk 2: ω1= 500 rad/s, ω2= –825 rad/s

c) Poincare map

Fig. 12. Disk 2: ω1= 740 rad/s, ω2= –221 rad/s

 Disk 2: ω1= 740 rad/s, ω2= –221 rad/s

a) Spectrogram

 Disk 2: ω1= 740 rad/s, ω2= –221 rad/s

b) Phase trajectory

 Disk 2: ω1= 740 rad/s, ω2= –221 rad/s

c) Poincare map

Fig. 13. Disk 2: ω1= 820 rad/s, ω2= –1353 rad/s

 Disk 2: ω1= 820 rad/s, ω2= –1353 rad/s

a) Spectrogram

 Disk 2: ω1= 820 rad/s, ω2= –1353 rad/s

b) Phase trajectory

 Disk 2: ω1= 820 rad/s, ω2= –1353 rad/s

c) Poincare map

Fig. 14. Disk 2: ω1= 868 rad/s, ω2= –1432 rad/s

 Disk 2: ω1= 868 rad/s, ω2= –1432 rad/s

a) Spectrogram

 Disk 2: ω1= 868 rad/s, ω2= –1432 rad/s

b) Phase trajectory

 Disk 2: ω1= 868 rad/s, ω2= –1432 rad/s

c) Poincare map

Fig. 15. Disk 2: ω1= 950 rad/s, ω2= –1567.5 rad/s

 Disk 2: ω1= 950 rad/s, ω2= –1567.5 rad/s

a) Spectrogram

 Disk 2: ω1= 950 rad/s, ω2= –1567.5 rad/s

b) Phase trajectory

 Disk 2: ω1= 950 rad/s, ω2= –1567.5 rad/s

c) Poincare map

Fig. 16. Disk 2: ω1= 990 rad/s, ω2= –1633.5 rad/s

 Disk 2: ω1= 990 rad/s, ω2= –1633.5 rad/s

a) Spectrogram

 Disk 2: ω1= 990 rad/s, ω2= –1633.5 rad/s

b) Phase trajectory

 Disk 2: ω1= 990 rad/s, ω2= –1633.5 rad/s

c) Poincare map

Fig. 17. Disk 2: ω1= 1080 rad/s, ω2= –1782 rad/s

 Disk 2: ω1= 1080 rad/s, ω2= –1782 rad/s

a) Spectrogram

 Disk 2: ω1= 1080 rad/s, ω2= –1782 rad/s

b) Phase trajectory

 Disk 2: ω1= 1080 rad/s, ω2= –1782 rad/s

c) Poincare map

Fig. 18. Disk 2: ω1= 1300 rad/s, ω2= –2145 rad/s

 Disk 2: ω1= 1300 rad/s, ω2= –2145 rad/s

a) Spectrogram

 Disk 2: ω1= 1300 rad/s, ω2= –2145 rad/s

b) Phase trajectory

 Disk 2: ω1= 1300 rad/s, ω2= –2145 rad/s

c) Poincare map

In Fig. 11, continuous spectrum and the disordered discrete points on Poincare map all suggest chaotic motion. With the continuous increasing of rotational speed, from Fig. 12, 740 rad/s for inner rotor, it can be seen that although continuous spectrum still can be observed the discrete points on Poincare map has shown a pattern of surrounding four equilibrium points, which suggests a transformation to periodic motion is undergoing. With higher rotational speed, 820 rad/s, period 4 orbital motion is reached, which can be seen quite obviously from Fig. 10(e), Fig. 10(f) and Fig. 13. Then the process repeats with increasing rotational speeds, see Fig. 10(g), Fig. 10(h) and Fig. 14-Fig. 17.

Besides the pattern of transformation between different states of motion, it is found that the chaotic motion is closely related with the SFD. Fig. 11(a), Fig. 15(a) and Fig. 18(a) show the spectrogram of chaotic motion at different rotational speeds, the same low frequency component, approximately 30 Hz, can be seen quite obviously. This frequency component, 30 Hz or 188 rad/s, is almost the value of the first critical speed excited by the inner rotor. This phenomenon suggests that the SFD has an important influence on the state of motion of the rotor system.

3.3. Experimental validation

The comparison between numerical and experimental results has already been made in previous section, Fig. 5 and 6. Orbits of disk 2 and 4 under three different rotational speeds are compared in this section for further validation and analysis, see Fig. 19-Fig. 21.

Fig. 19. Orbit of disk 2 and disk 4 with ω1= 128 rad/s, ω2= –211 rad/s

 Orbit of disk 2 and disk 4 with ω1= 128 rad/s, ω2= –211 rad/s

a) Disk 2: numerical results

 Orbit of disk 2 and disk 4 with ω1= 128 rad/s, ω2= –211 rad/s

b) Disk 2: experimental results

 Orbit of disk 2 and disk 4 with ω1= 128 rad/s, ω2= –211 rad/s

c) Disk 4: numerical results

 Orbit of disk 2 and disk 4 with ω1= 128 rad/s, ω2= –211 rad/s

d) Disk 4: experimental results

Following conclusions can be reached by analysis of Fig. 19-Fig. 21:

(1) With dual unbalance excitation, disk 2 and 4 execute petal-shaped orbits. Especially when approaching critical speeds of the rotor system, Fig. 20.

(2) The rotor system is slightly anisotropic, which might be caused by the squirrel cage, see Fig. 19(b), Fig. 19(d), 20(b), Fig. 20(d), Fig. 21(b), Fig. 21(d).

(3) Under relatively low speed, the numerical results show great agreement with the experimental results.

Fig. 20. Orbit of disk 2 and disk 4 with ω1= 160 rad/s, ω2= –264 rad/s

 Orbit of disk 2 and disk 4 with ω1= 160 rad/s, ω2= –264 rad/s

a) Disk 2: numerical results

 Orbit of disk 2 and disk 4 with ω1= 160 rad/s, ω2= –264 rad/s

b) Disk 2: experimental results

 Orbit of disk 2 and disk 4 with ω1= 160 rad/s, ω2= –264 rad/s

c) Disk 4: numerical results

 Orbit of disk 2 and disk 4 with ω1= 160 rad/s, ω2= –264 rad/s

d) Disk 4: experimental results

Fig. 21. Orbit of disk 2 and disk 4 with ω1= 210 rad/s, ω2= –346 rad/s

 Orbit of disk 2 and disk 4 with ω1= 210 rad/s, ω2= –346 rad/s

a) Disk 2: numerical results

 Orbit of disk 2 and disk 4 with ω1= 210 rad/s, ω2= –346 rad/s

b) Disk 2: experimental results

 Orbit of disk 2 and disk 4 with ω1= 210 rad/s, ω2= –346 rad/s

c) Disk 4: numerical results

 Orbit of disk 2 and disk 4 with ω1= 210 rad/s, ω2= –346 rad/s

d) Disk 4: experimental results

4. Conclusions

A modeling method combining the finite element method and the fixed interface modal synthesis method has been developed in the work to establish the nonlinear model of a counter-rotating dual-rotor test rig for steady state nonlinear response analysis. First, a finite element model of the rotor system is established without considering supports, from which the mass, stiffness and gyroscopic matrices are obtained. Together with these matrices, the fixed interface modal synthesis method is applied to establish the nonlinear model of the rotor system in which the nonlinearities of SFD and intermediate bearing are considered. Subsequently, the Newmark-β method is improved to solve the nonlinear equations of motion. Then dynamic characteristics of the rotor system are investigated.

Conclusions are listed as following:

1) The modeling method developed in this work is fast and efficient in establishing nonlinear model of complex rotor systems. Combing with the improved Newmark-β method, nonlinear dynamic characteristics can be obtained accurately and efficiently.

2) Due to coupling effect of the intermediate bearing, dual unbalance excitation frequencies are dominant in the responses of inner and outer rotor. Besides, frequency components other than dual unbalance excitation frequencies in the responses of inner and outer rotor are mainly caused by coupling effect of the intermediate bearing.

3) Rotors execute stable periodic orbital motion in 4-390 rad/s. Transformations between periodic, quasi periodic and chaotic motion are reached with increasing rotational speed. A low frequency component which is approximately the first critical speed of the outer rotor is observed, which means the SFD has a great influence on the state of motion of the rotor system.

4) Orbit of the rotor system under dual unbalance excitation is petal shaped, especially when approaching critical speeds of the rotor system. Under relatively low speed, the validities of the modeling method and the solving method were demonstrated with the experimental results.

References

  1. Ji Lucheng Review and prospect of counter-rotating turbine. Aeroengine, Vol. 32, Issue 4, 2006, p. 49-53, (in Chinese). [Search CrossRef]
  2. Delbez A., Charlot G., Ferraris G., et al. Dynamic behavior of a counter-rotating multi rotor air turbine starter. American Society of Mechanical Engineers, 1993. [Search CrossRef]
  3. Ferraris G. Prediction of the dynamic behavior of non-symmetric coaxial co-or counter rotating rotors. Journal of Sound and Vibration, Vol. 195, Issue 4, 1996, p. 649-666. [Publisher]
  4. Huang Taiping The transfer matrix component mode synthesis for the eigensolutions of rotor systems. Journal of Nanjing Aeronautical Institute, Vol. 5, Issue 1, 1988, p. 143-155. [Search CrossRef]
  5. Huang Taiping, Luo Guihuo Optimization of rotor dynamics. Journal of Aerospace Power, Vol. 9, Issue 2, 1994, p. 113-116, (in Chinese). [Search CrossRef]
  6. Vance John M., Zeidan Fouad Y., Murphy Brian Machinery Vibration and Rotordynamics. Wiley, 2010. [Publisher]
  7. Hibner D. H. Dynamic response of viscous-damped multi-shaft jet engines. Journal of Aircraft, Vol. 12, Issue 4, 1975, p. 305-312. [Publisher]
  8. Glasgow D. A., Nelson H. D. Stability analysis of rotor-bearing systems using component mode synthesis. Journal of Mechanical Design, Transactions of the ASME, Vol. 102, Issue 2, 1980, p. 352-359. [Publisher]
  9. Li D. F., Gunter E. J. Component mode synthesis of large rotor systems. Journal of Engineering for Power, ASME, Vol. 104, Issue 3, 1982, p. 552-560. [Publisher]
  10. Li D. F., Gunter E. J. Study of the modal truncation error in the component mode analysis of a dual-rotor system. Journal of Engineering for Power, ASME, Vol. 104, Issue 3, 1982, p. 525-532. [Publisher]
  11. Gupta K., Gupta K. D., Athre K. Unbalance response of a dual rotor system: theory and experiment. Journal of Vibration and Acoustics, Transactions of the ASME, Vol. 115, Issue 4, 1993, p. 427-435. [Publisher]
  12. Gupta K., Kumar R., Tiwari M., et al. Effect of rotary inertia and gyroscopic moments on dynamics of two spool aero engine rotor. International Gas Turbine and Aeroengine Congress and Exposition, New York, 1993, p. 1-14. [Search CrossRef]
  13. Maharathi B. B., Dash P. R., Beheraf A. K. Dynamic behaviour analysis of a dual-rotor system using the transfer matrix method. International Journal of Acoustics and Vibration, Vol. 9, Issue 3, 2004, p. 115-128. [Publisher]
  14. Wang W., Kirkhope J. Component mode synthesis for multi-shaft rotors with flexible inter-shaft bearings. Journal of Sound and Vibration, Vol. 173, Issue 4, 1994, p. 537-555. [Publisher]
  15. Chiang H. W. D., Hsu C. N., Tu S. H. Rotor-Bearing analysis for turbo machinery single- and dual-rotor systems. Journal of Propulsion and Power, Vol. 20, Issue 6, 2004, p. 1096-1104. [Publisher]
  16. Michael Friswell I., Penny J. E. T., Garvey S. D., Lees A. W. Dynamics of Rotating Machines. Cambridge University Press, 2010. [Publisher]
  17. Hutchinson J. R. Shear coefficients for Timoshenko beam theory. Journal of Applied Mechanics, Vol. 68, 2001, p. 87-92. [Publisher]
  18. Narayanaswami R., Adelman H. M. Inclusion of transverse shear deformation in finite element displacement formulations. AIAA Journal, Vol. 12, Issue 11, 1974, p. 1613-1614. [Publisher]
  19. Chen Guo Nonlinear dynamic response analysis of unbalance rubbing coupling faults of rotor ball bearing stator coupling system. Journal of Aerospace Power, Vol. 22, Issue 10, 2007, p. 1771-1778, (in Chinese). [Search CrossRef]
  20. Yang Xi-guan, Luo Gui-huo, Tang Zhen-huan, et al. Modeling method and dynamic characteristics of high-dimensional counter-rotating dual rotor system. Journal of Aerospace Power, Vol. 29, Issue 3, 2014, p. 585-595, (in Chinese). [Search CrossRef]
  21. Zhou Hai-lun, Luo Gui-huo, Chen Guo, Wang Fei Analysis of the nonlinear dynamic response of a rotor supported on ball bearings with floating-ring squeeze film dampers. Mechanism and Machine Theory, Vol. 59, 2013, p. 65-77. [Publisher]

Articles Citing this One

Chinese Journal of Aeronautics
Xinxing Ma, Hui Ma, Haiqin Qin, Xumin Guo, Chenguang Zhao, Mingyue Yu
2021
Shock and Vibration
Dongxiong Wang, Nianxian Wang, Kuisheng Chen, Chun Ye
2019
Feng Zhi-zhuang, Chen Quan-long, Cheng Qi-you
2018