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 quasiperiodic motion, chaotic motion and periodic motion again. The quasiperiodic 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.
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, counterrotating coaxial system has been applied to many aeroengines, such as RB211, GE120, F119 and so on. The counterrotating technology is beneficial to improve the fuel consumption rate and thrustweight 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 counterrotating dualrotor 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 [715]. 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 dualrotor 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 multirotor system. In 2004, HsiaoWei [15] established the mathematical model of a dualrotor 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 counterrotating 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 counterrotating 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$\beta $ 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 $\theta $ and $\phi $ for angular displacements in $X$ and $Y$ directions respectively.
The nodal displacement vector can be described as:
The displacement vector within the element can be interpolated as:
The shape functions in Eq. (2) are:
Fig. 1Diagram of beam element
In Eq. (3), ${l}_{e}$ is the element length and $g$ is formulated as given below:
where ${E}_{e}$ the Young’s Modulus, ${I}_{e}$ the area moment of inertia, ${G}_{e}$ the shear modulus, ${A}_{e}$ the crosssectional area. $\kappa $ is the shear coefficient, and for a hollow circular shaft section [16, 17]:
where $\mu $ is the Poisson’s ratio, $\lambda $ is the ratio of the inner shaft radius to the outer shaft radius. Hence for a solid shaft $\lambda \to 0$.
With rotary inertia and shear effects considered, the kinetic energy and strain energy for a shaft element is:
where ${\rho}_{e}$ denotes the density of the material and $\gamma $ is the transverse shear strain [18].
And thus the element mass matrix:
The element stiffness matrix:
where:
The element gyroscopic matrix:
where:
For disk element, the mass and gyroscopic matrices are:
where ${m}_{d}^{i}$, ${I}_{d}^{i}$ and ${I}_{p}^{i}$ 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]:
where $x$ and $y$ are the horizontal and vertical displacements of the journal. ${I}_{j}$ ($j=$1, 2, 3) are Sommerfeld integrals. $R$ and $L$ are radius and length of the SFD respectively. ${\mu}_{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:
where the superscript $ir$ and $or$ in Eq.(14c) denote inner ring and outer ring of the bearing. ${k}_{n}$ in Eq.(14a) is the contact stiffness between rollers and rings. ${N}_{b}$ represent number of the rollers. $\xi =3/2$ for the rolling ball bearing used in the intermediate support. ${\theta}_{j}$ is the rotation angle of the$j$th roller at time $t$. $\gamma $ and ${\omega}_{c}$ are bearing radial clearance and rotational speed of the bearing retainer respectively. ${u}_{\theta j}$ is the elastic radial deformation of the $j$th roller. $r$and $R$ represent radius of the inner and outer ring. ${\omega}_{in}$ and ${\omega}_{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:
where ${\mathbf{F}}^{unb}$ is the unbalance force vector. ${\mathbf{F}}^{nonl}$ is the nonlinear forces exerted by the SFDs and rolling bearing in this work. $\mathbf{M}$, $\mathbf{K}$, and $\mathbf{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:
where ${\mathbf{u}}_{I}$ can be interpreted as linear DOFs, with only unbalance forces acting on, while ${\mathbf{u}}_{J}$ 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:
where ${\mathbf{\varphi}}_{k}$ – the mass normalized normal mode matrix, ${\mathbf{\varphi}}_{c}$ – the mass normalized constrained mode, $\mathbf{I}$ – unity matrix, ${\mathbf{q}}_{k}$ – the normal mode coordinates, $\mathbf{T}$ – the transformation matrix. ${\mathbf{\varphi}}_{c}$ is given as below:
In Eq. (16), let ${\mathbf{u}}_{J}=\mathbf{}$0 and neglecting the gyroscopic effects gives:
Solving the eigenvalue problem corresponding to Eq. (17c) gives the mass normalized normal mode matrix ${\mathbf{\varphi}}_{k}$ and the modal frequency ${\mathrm{\Omega}}_{k}$.
Substitute Eq. (17a) into Eq. (16) and premultiply with ${\mathbf{T}}^{T}$ to obtain the reduced equations of motion:
where:
where ${\mathrm{\Omega}}_{r}$ is the $r$th modal frequency obtained from Eq. (17c). ${\mathbf{F}}_{I}^{unb}$ is the unbalance force acting on linear DOFs. ${\mathbf{F}}_{J}^{nonl}$ is the nonlinear force acting on nonlinear DOFs.
And:
where ${\mathbf{F}}_{J}^{SFD}=\left[\begin{array}{ccc}{f}_{1x}^{sfd}& {f}_{1y}^{sfd}& \cdots \begin{array}{cc}{f}_{Nx}^{sfd}& {f}_{Ny}^{sfd}\end{array}\end{array}\right]$ is the nonlinear force vector exerted by SFDs which can be calculated with Eq. (13). ${\mathbf{F}}_{J}^{B}=\left[\begin{array}{ccc}{f}_{x}^{b}& {f}_{y}^{b}& \begin{array}{cc}{f}_{x}^{b}& {f}_{y}^{b}\end{array}\end{array}\right]$ is the nonlinear force vector exerted by rolling bearing which can be calculated with Eq. (14). ${\mathbf{k}}_{JJ}$ and ${\mathbf{c}}_{JJ}$ are given below:
where ${k}_{Ji}$ and ${c}_{Ji}$ are the stiffness and damping coefficients of the elastic support, respectively. $N$ is the number of elastic supports.
Rearranging Eq. (18) gives:
In Eq. (23), the vector ${\mathbf{q}}_{k}$ and ${\mathbf{u}}_{J}$ 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$\beta $ method applies to Eq.(23a) while implicit Newmark$\beta $ method applies to Eq. (23b).
According to assumptions of Newmark$\beta $ method, in time interval $\left[\begin{array}{cc}{t}_{n}& {t}_{n+1}\end{array}\right]$:
where ${t}_{n+1}={t}_{n}+\mathrm{\Delta}t$ and $\mathrm{\Delta}t$ is the time increment. The superscript $n$ and $n+1$ denote ${t}_{n}$ and ${t}_{n+1}$.
Following equations can be obtained from Eq. (24):
With:
Substituting Eq. (24)Eq. (27) into Eq. (23) yields:
With:
Using Eq. (26) to substitute for ${\dot{\mathbf{u}}}_{J}^{\mathrm{n}+1}$ in Eq. (29g) and substitute the new Eq. (29g) into Eq. (28b) to yield nonlinear equations of ${\mathbf{u}}_{J}^{n+1}$. After solving Eq. (28b) for ${\mathbf{u}}_{J}^{n+1}$ with numerical algorithms, ${\mathbf{q}}_{k}^{n+1}$ can be obtained with ${\mathbf{u}}_{J}^{n+1}$ substituted into Eq. (28a). Thus ${\dot{\mathbf{u}}}_{J}^{n+1}$, ${\ddot{\mathbf{u}}}_{J}^{n+1}$, ${\dot{\mathbf{q}}}_{k}^{n+1}$, ${\ddot{\mathbf{q}}}_{k}^{n+1}$ and can be solved from Eq. (25) and Eq. (26). The process continues to repeat and move forward to find ${\mathbf{u}}_{J}^{n+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 timedomain method based on Newmark$\beta $ 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 $\mathbf{M}$, $\mathbf{K}$ and $\mathbf{G}$ in Eq. (15) can be obtained and partitioned as Eq. (16) for calculating ${\mathbf{\varphi}}_{c}$ and later use. Then, in accordance with Eq. (17c), modal analysis is conducted on the model with all supports pinned to obtain ${\mathbf{\varphi}}_{k}$ and ${\mathrm{\Omega}}_{k}$. The left hand side of Eq. (18a) now is obtained with ${\mathbf{\varphi}}_{k}$, ${\mathbf{\varphi}}_{c}$ and the partitioned $\mathbf{M}$, $\mathbf{K}$ and $\mathbf{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. 2Flowchart of modeling and solving
2.4. Test rig description
An overall view of the dualrotor test rig is given in Fig. 3. The test rig was designed to study the dynamic characteristics of coand counterrotating dualrotor systems with 4 or 5 supports. In this work, the 4 supports counterrotating dualrotor 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 runup and rundown through the speed range of the test rig. The angular acceleration for inner rotor during runup 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 0232 rad/s and 0382 rad/s respectively. The excitation is only due to the residual unbalance of the rotors.
Fig. 3Dualrotor test rig
3. Numerical results and experimental validation
3.1. Model introduction
The counterrotating 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. 4Structural diagram of the coaxial rotor system
Model parameters of the rotor system studied in this work are listed in Table 15. The Young’s modulus of the shaft is 196 GPa, mass density 7810 kg/m^{3}, and shear modulus 75.5 GPa. To apply modal synthesis method, 40 modes are retained in the normal mode matrix ${\mathbf{\varphi}}_{k}$. For the Newmark$\beta $ method, $\alpha =$ 0.25, $\beta =$ 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 1Dimension 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 2Stiffness of elastic supports (squirrel cage)
Support I  Support II  Support III  Support IV  
Stiffness (N/m)  1.45×10^{6}  2.21×10^{5}  9.29×10^{5} 
Table 3Parameters of the intermediate bearing
Radius of inner ring (mm)  Radius of outer ring (mm)  No. of rollers  Contact stiffness (N/m^{3/2})  Radial clearance (μm) 
9.37  14.13  9  7.055×109  6 
Table 4Parameters 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 ${\mu}_{s}$ / 10^{2}^{}Pa·s  1.0752 
Table 5Unbalance configuration and inertia properties of disks
Inner rotor  Outer rotor  
Disk 1  Disk 2  Disk 3  Disk 4  
Unbalance (×10^{5}^{}kg·m)  2  4  1  2 
Mass (kg)  2.3386  2.3386  3.2590  1.6303 
Polar moment of inertia (kg·m^{2})  0.00815  0.00815  0.01561  0.00661 
Fig. 5Spectrum cascade of the horizontal response of disk 2 and 4numerical simulation
a) Disk 2
b) Disk 4
3.2. Unbalance response analysis
With the model established by method described in Section 2, steadystate 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. 6Spectrum cascade of the horizontal response of disk 2 and 4experimental results
a) Disk 2
b) Disk 4
Spectrum cascades of the horizontal response of disk 2 and disk 4 with inner rotor operating in 4400 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 4001600 rad/s are shown in Figs. 79. 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. 59:
(1) Responses of the inner and outer rotor are coupled because of the intermediate bearing. Frequency components, ${\omega}_{1}$ and ${\omega}_{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 ${\omega}_{1}$ and ${\omega}_{2}$ in the responses of disk 2 and 4, such as ${2\omega}_{1}+{\omega}_{2}$, ${\omega}_{1}+{2\omega}_{2}$, ${3\omega}_{1}+2{\omega}_{2}$, ${2\omega}_{1}+3{\omega}_{2}$, $({\omega}_{2}{\omega}_{1})/2$, seen in Figs. 59, 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\omega}_{2}{\omega}_{1}\text{,}$${3\omega}_{2}2{\omega}_{1}$emerge when inner rotor operating in 400800 rad/s. However, ${2\omega}_{2}{\omega}_{1}$ and ${3\omega}_{2}2{\omega}_{1}$ are not observed for both disk 2 and 4 in Fig. 8 and Fig. 9 when the inner rotor operating in 8001200 rad/s and 12001600 rad/s. Besides, ${2\omega}_{1}+{\omega}_{2}$ is persistent in 01600 rad/s for disk 2 and 4, which can be seen from Fig. 5 and Figs. 79.
(4) Two critical speeds of the rotor system, 168 rad/s and 186 rad/s, are observed in rotational speed range 0200 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. 7Spectrum cascade of the horizontal response of disk 2 and 4numerical simulation
a) Disk 2
b) Disk 4
Fig. 8Spectrum cascade of the horizontal response of disk 2 and 4numerical simulation
a) Disk 2
b) Disk 4
Fig. 9Spectrum cascade of the horizontal response of disk 2 and 4numerical simulation
a) Disk 2
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\times 2\pi /{\omega}_{1}$ because the rotational speed ratio of outer and inner rotor is –1.65 and $5\approx $ 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 4390 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. 10Bifurcation diagrams of the horizontal response of disk 2 and disk 4 with rotor speed
a) Disk 2
b) Disk 4
c) Disk 2
d) Disk 4
e) Disk 2
f) Disk 4
g) Disk 2
h) Disk 4
Fig. 11Disk 2: ω1= 500 rad/s, ω2= –825 rad/s
a) Spectrogram
b) Phase trajectory
c) Poincare map
Fig. 12Disk 2: ω1= 740 rad/s, ω2= –221 rad/s
a) Spectrogram
b) Phase trajectory
c) Poincare map
Fig. 13Disk 2: ω1= 820 rad/s, ω2= –1353 rad/s
a) Spectrogram
b) Phase trajectory
c) Poincare map
Fig. 14Disk 2: ω1= 868 rad/s, ω2= –1432 rad/s
a) Spectrogram
b) Phase trajectory
c) Poincare map
Fig. 15Disk 2: ω1= 950 rad/s, ω2= –1567.5 rad/s
a) Spectrogram
b) Phase trajectory
c) Poincare map
Fig. 16Disk 2: ω1= 990 rad/s, ω2= –1633.5 rad/s
a) Spectrogram
b) Phase trajectory
c) Poincare map
Fig. 17Disk 2: ω1= 1080 rad/s, ω2= –1782 rad/s
a) Spectrogram
b) Phase trajectory
c) Poincare map
Fig. 18Disk 2: ω1= 1300 rad/s, ω2= –2145 rad/s
a) Spectrogram
b) Phase trajectory
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. 14Fig. 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. 19Fig. 21.
Fig. 19Orbit of disk 2 and disk 4 with ω1= 128 rad/s, ω2= –211 rad/s
a) Disk 2: numerical results
b) Disk 2: experimental results
c) Disk 4: numerical results
d) Disk 4: experimental results
Following conclusions can be reached by analysis of Fig. 19Fig. 21:
(1) With dual unbalance excitation, disk 2 and 4 execute petalshaped 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. 20Orbit of disk 2 and disk 4 with ω1= 160 rad/s, ω2= –264 rad/s
a) Disk 2: numerical results
b) Disk 2: experimental results
c) Disk 4: numerical results
d) Disk 4: experimental results
Fig. 21Orbit of disk 2 and disk 4 with ω1= 210 rad/s, ω2= –346 rad/s
a) Disk 2: numerical results
b) Disk 2: experimental results
c) Disk 4: numerical results
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 counterrotating dualrotor 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$\beta $ 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$\beta $ 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 4390 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

Ji Lucheng Review and prospect of counterrotating turbine. Aeroengine, Vol. 32, Issue 4, 2006, p. 4953, (in Chinese).

Delbez A., Charlot G., Ferraris G., et al. Dynamic behavior of a counterrotating multi rotor air turbine starter. American Society of Mechanical Engineers, 1993.

Ferraris G. Prediction of the dynamic behavior of nonsymmetric coaxial coor counter rotating rotors. Journal of Sound and Vibration, Vol. 195, Issue 4, 1996, p. 649666.

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. 143155.

Huang Taiping, Luo Guihuo Optimization of rotor dynamics. Journal of Aerospace Power, Vol. 9, Issue 2, 1994, p. 113116, (in Chinese).

Vance John M., Zeidan Fouad Y., Murphy Brian Machinery Vibration and Rotordynamics. Wiley, 2010.

Hibner D. H. Dynamic response of viscousdamped multishaft jet engines. Journal of Aircraft, Vol. 12, Issue 4, 1975, p. 305312.

Glasgow D. A., Nelson H. D. Stability analysis of rotorbearing systems using component mode synthesis. Journal of Mechanical Design, Transactions of the ASME, Vol. 102, Issue 2, 1980, p. 352359.

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. 552560.

Li D. F., Gunter E. J. Study of the modal truncation error in the component mode analysis of a dualrotor system. Journal of Engineering for Power, ASME, Vol. 104, Issue 3, 1982, p. 525532.

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. 427435.

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. 114.

Maharathi B. B., Dash P. R., Beheraf A. K. Dynamic behaviour analysis of a dualrotor system using the transfer matrix method. International Journal of Acoustics and Vibration, Vol. 9, Issue 3, 2004, p. 115128.

Wang W., Kirkhope J. Component mode synthesis for multishaft rotors with flexible intershaft bearings. Journal of Sound and Vibration, Vol. 173, Issue 4, 1994, p. 537555.

Chiang H. W. D., Hsu C. N., Tu S. H. RotorBearing analysis for turbo machinery single and dualrotor systems. Journal of Propulsion and Power, Vol. 20, Issue 6, 2004, p. 10961104.

Michael Friswell I., Penny J. E. T., Garvey S. D., Lees A. W. Dynamics of Rotating Machines. Cambridge University Press, 2010.

Hutchinson J. R. Shear coefficients for Timoshenko beam theory. Journal of Applied Mechanics, Vol. 68, 2001, p. 8792.

Narayanaswami R., Adelman H. M. Inclusion of transverse shear deformation in finite element displacement formulations. AIAA Journal, Vol. 12, Issue 11, 1974, p. 16131614.

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. 17711778, (in Chinese).

Yang Xiguan, Luo Guihuo, Tang Zhenhuan, et al. Modeling method and dynamic characteristics of highdimensional counterrotating dual rotor system. Journal of Aerospace Power, Vol. 29, Issue 3, 2014, p. 585595, (in Chinese).

Zhou Hailun, Luo Guihuo, Chen Guo, Wang Fei Analysis of the nonlinear dynamic response of a rotor supported on ball bearings with floatingring squeeze film dampers. Mechanism and Machine Theory, Vol. 59, 2013, p. 6577.