Full Text:   <2091>

CLC number: TH133.2

On-line Access: 2013-04-03

Revision Accepted: 2013-01-22

Crosschecked: 2013-03-06

Cited: 7

Clicked: 3771

Citations:  Bibtex RefMan EndNote GB/T7714

 Journal of Zhejiang University SCIENCE A 2013 Vol.14 No.4 P.268-280 http://doi.org/10.1631/jzus.A1200298

Investigation of the dynamic characteristics of a dual rotor system and its start-up simulation based on finite element method

 Author(s):  Zhong-xiu Fei1, Shui-guang Tong2, Chao Wei3 Affiliation(s):  1. Institute of Chemical Machinery Engineering, Zhejiang University, Hangzhou 310027, China; more Corresponding email(s):   feizhongxiu@zju.edu.cn Key Words:  Dual rotor system, Critical speed, Transient response, Finite element method (FEM) Share this article to： More <<< Previous Article|Next Article >>>

Zhong-xiu Fei, Shui-guang Tong, Chao Wei. Investigation of the dynamic characteristics of a dual rotor system and its start-up simulation based on finite element method[J]. Journal of Zhejiang University Science A, 2013, 14(4): 268-280.

@article{title="Investigation of the dynamic characteristics of a dual rotor system and its start-up simulation based on finite element method",
author="Zhong-xiu Fei, Shui-guang Tong, Chao Wei",
journal="Journal of Zhejiang University Science A",
volume="14",
number="4",
pages="268-280",
year="2013",
publisher="Zhejiang University Press & Springer",
doi="10.1631/jzus.A1200298"
}

%0 Journal Article
%T Investigation of the dynamic characteristics of a dual rotor system and its start-up simulation based on finite element method
%A Zhong-xiu Fei
%A Shui-guang Tong
%A Chao Wei
%J Journal of Zhejiang University SCIENCE A
%V 14
%N 4
%P 268-280
%@ 1673-565X
%D 2013
%I Zhejiang University Press & Springer
%DOI 10.1631/jzus.A1200298

TY - JOUR
T1 - Investigation of the dynamic characteristics of a dual rotor system and its start-up simulation based on finite element method
A1 - Zhong-xiu Fei
A1 - Shui-guang Tong
A1 - Chao Wei
J0 - Journal of Zhejiang University Science A
VL - 14
IS - 4
SP - 268
EP - 280
%@ 1673-565X
Y1 - 2013
PB - Zhejiang University Press & Springer
ER -
DOI - 10.1631/jzus.A1200298

Abstract:
Recently, the finite element method (FEM) has been commonly applied in the engineering analysis of rotor dynamics. Gyroscopic moments, rotary inertia, transverse shear deformation and gravity can be included in computational models of rotor-bearing systems. In this paper, a finite element model and its solution method are presented for the calculation of the dynamics of dual rotor systems. A typical structure with two rotor shafts is discussed and the procedure for obtaining the coupling motion equations of the subsystems is illustrated. A computer program is developed to solve critical speeds and to simulate the transient motion. The influence of gyroscopic moments on co-rotation and counter-rotation is analyzed, and the effect of the speed ratio on critical speed is studied. The dynamic characteristics under different conditions of increasing speed during start-up are demonstrated by comparison with transient nodal displacements. The presented model provides a complete foundation for further investigation of the dynamics of dual rotor systems.

1.  Introduction

Despite the extensive use of the transfer matrix method in rotor dynamics problems, the finite element technique has become more popular recently for its high efficiency and convenience of modeling. The finite element method (FEM) for rotors was first developed by Ruhl and Booker (1970). Nelson and Mc-Vaugh () extended it to include translational and rotary inertias, gyroscopic moments, and axial loads, and Nelson () considered the effect of shear deformation using the Timoshenko beam theory. With the application of this method, we are able to solve complex problems because the equations of motion are explicit and easier to couple with bearing flexibilities or other multiple excitations.

Dual rotor systems are widely used in aircraft engines because they can reduce the influence of gyroscopic moments generated by high and low pressure rotors when counter-rotating. Research on the dynamic characteristics of dual rotor structures is of great significance for lowering vibration and noise, and can enhance their overall performance, including their reliability and stability. In a dual rotor system, because of the inclusion of a special inter-bearing which connects the inter rotor with the outer rotor, the corresponding force becomes more complicated. The equations of the two subsystems should be coupled. This is the key difference between a dual and a single rotor system in terms of mathematical modeling. Critical speed prediction of a typical dual rotor system was first proposed by Rajan et al. (). Hsiao-Wei et al. (2002) carried out a systematic theoretical analysis and an analysis of the effects of the speed ratio and inter-bearing stiffness. However, little attention has been focused on the transient motion, which is becoming increasingly important in understanding the dynamic behavior when a rotor starts up, stops, or goes through a critical speed.

In this paper, based on the studies of Nelson () and Hsiao-Wei et al. (2002), a computational model of a dual rotor bearing system is established under an inertia coordinate. The procedure integrates the advantages of FEM and the Newmark method () to obtain natural frequencies, critical speed maps, and transient results. The basic process for extending this model to a more complex rotor-bearing system such as a triple rotor system and geared rotors is also presented.

2.  Component equations

A typical flexible rotor-bearing system with a single rotating shaft consists of a rotor, rotor segments and discrete bearings (Fig. 1) along with an inertia coordinate frame. The rotor is composed of discrete disks and the rotor segments have distributed mass and elasticity.

Fig.1
Coordinate system

Expressions of kinetic energy and strain energy are necessary to characterize the disk and shaft element. Lagrange’s formulation is applied to obtain the differential equations of motion for each component. The formulation is written in the form $$\frac{d}{{dt}}\left( {\frac{{\partial T}}{{\partial {{\dot u}_i}}}} \right) - \frac{{\partial T}}{{\partial {u_i}}} + \frac{{\partial U}}{{\partial {u_i}}} = {Q_i}$$, where T and U are the kinetic and strain energies, respectively, ui are generalized independent coordinates and Qi are generalized forces.

2.1.  Angular velocity and coordinates

OXYZ is an inertia coordinate, in which O is coincident with the center of the rigid disk and the X axis is collinear with the undeformed rotor center line. OX 3 Y 3 Z 3 is attached to the disk and rotates synchronously with the rotor. The relationship between them can be described by the following successive rotations: OXYZ rotates θz about OZ, yields OX 1 Y 1 Z 1; OX 1 Y 1 Z 1 rotates θy 1 about OY 1, yields OX 2 Y 2 Z 2; and OX 2 Y 2 Z 2 rotates θx 2 about OX 2, yields OX 3 Y 3 Z 3, where $${\theta _{y1}} \approx {\theta _y}$$, $${\theta _{x2}} = \phi$$, and $${\dot \theta _{x2}} = \dot \phi = \Omega$$; &phgr; is the angle of rotation about X and Ω is the spin speed, and the angular velocity of OX 3 Y 3 Z 3 is $${\left[ {\begin{array}{*{20}{c}} {{\omega _{x3}}}\quad{{\omega _{y3}}}\quad{{\omega _{z3}}} \end{array}} \right]^{T}} = Rot(x2,{\theta _{x2}}){\left[ {\begin{array}{*{20}{c}} {{{\dot \theta }_{x2}}}\quad0\quad0 \end{array}} \right]^T} + Rot(x2,{\theta _{x2}})Rot(y1,{\theta _{y1}}){\left[ {\begin{array}{*{20}{c}} 0\quad{{{\dot \theta }_{y1}}}\quad0 \end{array}} \right]^T} + Rot(x2,{\theta _{x2}})Rot(y1,{\theta _{y1}})Rot(z,{\theta _z}){\left[ {\begin{array}{*{20}{c}} 0\quad0\quad{{{\dot \theta }_z}} \end{array}} \right]^T}$$, where $$Rot(z,{\theta _z})$$, $$Rot(y1,{\theta _{y1}})$$ and $$Rot(x2,{\theta _{x2}})$$ are rotation matrices.

2.2.  Rigid disk

For a typical rotor disk, strain energy is neglected in view of the rigid body assumption, and it can be modeled as a four-degrees-of-freedom element. The vector of generalized coordinates $${u_D} = {\left[ {\begin{array}{*{20}{c}} V\quad W\quad B\quad \Gamma \end{array}} \right]^T}$$ includes the displacements V, W and the corresponding slopes B, Γ. With the application of the angular velocity of OX 3 Y 3 Z 3, the kinetic energy of the disk for lateral motion is given by $${T_D} = \frac{1}{2}{\left[ {\begin{array}{*{20}{c}} {\dot V} \\ {\dot W} \end{array}} \right]^T}\left[ {\begin{array}{*{20}{c}} {{m_D}}\quad0 \\ 0\quad{{m_D}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\dot V} \\ {\dot W} \end{array}} \right]\; + \frac{1}{2}{\left[ {\begin{array}{*{20}{c}} {\dot B} \\ {\dot \Gamma } \end{array}} \right]^T}\left[ {\begin{array}{*{20}{c}} {{J_d}}\quad0 \\ 0\quad{{J_d}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\dot B} \\ {\dot \Gamma } \end{array}} \right] - \Omega \dot \Gamma B{J_p} + \frac{1}{2}{J_p}{\Omega ^2}$$, where m D is the disk mass, J d denotes the diametral inertia and J p denotes the polar inertia. Following the Lagrangian approach, we can obtain: $${M_D}{\ddot u_D} - \Omega {G_D}{\dot u_D} = {Q_D}$$, where M D is the mass matrix of the disk element and its gyroscopic matrix is G D. The unbalance force (Fig. 2) caused by the eccentric mass is (Carrella, ) $$Q_D^U = {m_D}{\Omega ^2}{\left[ {\begin{array}{*{20}{c}} {{e_{y3}}}\quad{{e_{z3}}}\quad0\quad0 \end{array}} \right]^T}\cos (\Omega t) + {m_D}{\Omega ^2}{\left[ {\begin{array}{*{20}{c}} { - {e_{z3}}}\quad{{e_{y3}}}\quad0\quad0 \end{array}} \right]^T}\sin (\Omega t) = {m_D}{\Omega ^2}{\left[ {\begin{array}{*{20}{c}} {e\cos (\Omega t + \varphi )}\quad{e\sin (\Omega t + \varphi )}\quad0\quad0 \end{array}} \right]^T}$$, where e is the eccentricity, φ is the phase of unbalance, and the disk mass center is located at (e y3,e z3) in OX 3 Y 3 Z 3. When gravity needs to be considered, the item $$Q_D^G$$ is included (Appendix).

Fig.2
Unbalance force

2.3.  Finite shaft element

Fig. 3 shows a two-node shaft element with eight degrees of freedom. The time displacements of an arbitrary cross section in the element can be expressed by the product of its nodal displacements and shape functions (Michael et al., 2010).

Fig.3
Shaft element and generalized coordinates

According to the geometric analysis in Fig. 4, the slopes (B, Γ) can be given by the following relationships (Jei and Lee, ): $$\Gamma (s,t) = \frac{{\partial V(s,t)}}{{\partial s}}$$, $$B(s,t) = - \frac{{\partial W(s,t)}}{{\partial s}}$$.

Fig.4
Geometric relationship of generalized coordinates

For the shaft element the vector $${u_{\text{S}}} = {[{V_{\text{1}}}\;{W_{\text{1}}}\;{B_{\text{1}}}\;{\Gamma _{\text{1}}}\;{V_{\text{2}}}\;{W_{\text{2}}}\;{B_{\text{2}}}\;{\Gamma _{\text{2}}}]^T}$$ represents the generalized coordinates. The interpolation functions (N 1, N 2, N 3, and N 4) are used to calculate the nodal displacements, and the expressions are given by $$V = \left[ {\begin{array}{*{20}{c}} {{N_1}}\quad0\quad0\quad{{N_2}}\quad{{N_3}}\quad0\quad0\quad{{N_4}} \end{array}} \right]{u_S} = {N_V}{u_S}$$, $$W = \left[ {\begin{array}{*{20}{c}} 0\quad{{N_1}}\quad{ - {N_2}}\quad0\quad0\quad{{N_3}}\quad{ - {N_4}}\quad0 \end{array}} \right]{u_S} = {N_W}{u_S}$$, $$B = \left[ {\begin{array}{*{20}{c}} 0\quad{ - {N_1}^\prime }\quad{{N_2}^\prime }\quad0\quad0\quad{ - {N_3}^\prime }\quad{{N_4}^\prime }\quad0 \end{array}} \right]{u_S} = {N_B}{u_S}$$, $$\Gamma = \left[ {\begin{array}{*{20}{c}} {{N_1}^\prime }\quad0\quad0\quad{{N_2}^\prime }\quad{{N_3}^\prime }\quad0\quad0\quad{{N_4}^\prime } \end{array}} \right]{u_S} = {N_\Gamma }{u_S}$$.

The corresponding velocities are calculated by taking derivatives with respect to time. Similarly, the kinetic energy is obtained by computing an integral over the length of the element (Rao et al., ): ${\kern 0pt} \begin{matrix} {T_S} = \int_0^L {\frac{1}{2}\left( {\mu {{\left[ {\begin{array}{*{20}{c}} {\dot V} \\ {\dot W} \end{array}} \right]}^T}\left[ {\begin{array}{*{20}{c}} {\dot V} \\ {\dot W} \end{array}} \right] + {j_d}{{\left[ {\begin{array}{*{20}{c}} {\dot B} \\ {\dot \Gamma } \end{array}} \right]}^T}\left[ {\begin{array}{*{20}{c}} {\dot B} \\ {\dot \Gamma } \end{array}} \right]} \right)} ds \hfill \\ - \int_0^L {\Omega {j_p}} \dot \Gamma Bds + \int_0^L {\frac{1}{2}{j_p}} {\Omega ^2}ds \hfill \\ = \dot u_S^T\left( {\int_0^L {\frac{1}{2}\mu {{\left[ {\begin{array}{*{20}{c}} {{N_V}} \\ {{N_W}} \end{array}} \right]}^T}\left[ {\begin{array}{*{20}{c}} {{N_V}} \\ {{N_W}} \end{array}} \right]ds} } \right){{\dot u}_S} \hfill \\ + \dot u_S^T\left( {\int_0^L {\frac{1}{2}{j_d}{{\left[ {\begin{array}{*{20}{c}} {{N_B}} \\ {{N_\Gamma }} \end{array}} \right]}^T}\left[ {\begin{array}{*{20}{c}} {{N_B}} \\ {{N_\Gamma }} \end{array}} \right]ds} } \right){{\dot u}_S} \hfill \\ - \Omega \dot u_S^T\left( {\int_0^L {{j_p}N_\Gamma ^T{N_B}ds} } \right){{\dot u}_S} + \int_0^L {\frac{1}{2}{j_p}{\Omega ^2}ds} , \hfill \\ \end{matrix}$ where μ, j d and j p are the shaft element mass, diametral inertia and polar inertia per unit length, respectively.

For a differential element with a thickness of ds and located at s, the potential energy of elastic bending is (Rao, ) $$dU_S^B = \frac{1}{2}EI\left( {{{\left( {\frac{{{\partial ^2}y}}{{\partial {s^2}}}} \right)}^2} + {{\left( {\frac{{{\partial ^2}z}}{{\partial {s^2}}}} \right)}^2}} \right)ds = \frac{1}{2}EI{\left[ {\begin{array}{*{20}{c}} {V''} \\ {W''} \end{array}} \right]^T}\left[ {\begin{array}{*{20}{c}} {V''} \\ {W''} \end{array}} \right]ds$$, where E is the Young’s modulus, and I is the second moment of area.

Integrating the $$dU_S^B$$ and following the Lagrangian approach, we can obtain: $$\left( {M_S^T + M_S^R} \right){\ddot u_S} - \Omega {G_S}{\dot u_S} + {K_S}{u_S} = {Q_S}$$, where $$M_S^T$$, $$M_S^R$$, $${G_S}$$ and $${K_S}$$ are the translational mass, rotatory mass, gyroscopic and stiffness matrices of the shaft element, respectively.

Based on the principle of virtual displacement, the gravity force vector $$Q_S^G$$ can be deduced (Appendix). The virtual work of the gravity effect is $$\delta u_S^TQ_S^G = \int_0^L { - \mu g\delta {W^T}ds} = \int_0^L { - \mu g{{[{N_W}\delta {u_S}]}^T}ds} = \delta u_S^T\int_0^L { - \mu gN_W^Tds}$$.

2.4.  Bearing

With the assumption of little vibration, the inherent nonlinear characteristics of fluid-film bearings can be linearized at the static equilibrium position (Zhong, ; Abduljabbar et al., ). Then its governing equation can be expressed as $${M_b}{\ddot u_b} + {C_{b1}}{\dot u_b} + {K_{b1}}{u_b} = - {C_{b2}}\left( {{{\dot u}_b} - {{\dot \hat u}_s}} \right) - {K_{b2}}\left( {{u_b} - {{\hat u}_s}} \right)$$, where $${u_b} = {\left[ {{V_b},{W_b}} \right]^T}$$ are the displacements of the bearing node and $${\hat u_S} = {\left[ {{V_j},{W_j}} \right]^T}$$ are the displacements of the journal node. M b, C b1, K b1, C b2 and K b2 are included (Appendix), and the elements of the matrices are illustrated in Fig. 5. The generalized force acting on the journal is (Das et al., ) $${Q_b} = {C_{b2}}\left( {{{\dot u}_b} - {{\dot \hat u}_S}} \right) + {K_{b2}}\left( {{u_b} - {{\hat u}_S}} \right)$$.

Fig.5
Bearing stiffness and damping

In the case of a high stiffness foundation, V b=W b=0, and the equation is simplified to $${Q_b} = - {C_{b2}}{\dot \hat u_S} - {K_{b2}}{\hat u_S}$$.

3.  Assembly and system equations

3.1.  System of a single rotating shaft

Fig. 6 shows a rotor-bearing system with a single rotating shaft. The vibration equation of the rotor can be obtained by combining Eqs. (4) and (14) (Shiau et al., ): $${M^s}{\ddot u^s} - \Omega {G^s}{\dot u^s} + {K^s}{u^s} = Q_u^s + Q_d^s$$,

Typical rotor-bearing system

where M s, G s and K s are the mass, gyroscopic and stiffness matrices, respectively of the rotor system, and u s is the vector of its generalized coordinates. $$Q_u^s$$ is the unbalance force, and $$Q_d^s$$ denotes the interaction force between the rotor and bearings.

The vibration equation of the bearings is given by $$M_b^s\ddot u_b^s + C_b^s\dot u_b^s + K_b^su_b^s = - Q_b^s$$, where $$M_b^s$$, $$C_b^s$$ and $$K_b^s$$ are the mass, damping and stiffness matrices of the bearings system, respectively. $$Q_b^s$$ is the vector of the force and it is applied to the rotor by bearings. It can be calculated using the equation: $$Q_b^s = {\left[ {\begin{array}{*{20}{c}} {Q_b^i}\quad{Q_b^j}\quad{Q_b^k} \end{array}} \right]^T} = - {C_O}L{\dot u^s} - {K_O}L{u^s} + {C_O}{L_b}\dot u_b^s + {K_O}{L_b}u_b^s = \left[ {\begin{array}{*{20}{c}} { - {C_O}L}\quad{{C_O}{L_b}} \end{array}} \right]{[{\dot u^s}\;\dot u_b^s]^T} + \left[ {\begin{array}{*{20}{c}} { - {K_O}L}\quad{{K_O}{L_b}} \end{array}} \right]{[{u^s}\;u_b^s]^T}$$, with $$L{u^s} = {[\begin{array}{*{20}{c}} {{V_i}}\quad{{W_i}}\quad{{V_j}}\quad{{W_j}}\quad{{V_k}}\quad{{W_k}} \end{array}]^T}$$, $${L_b}u_b^s = {I_{6 \times 6}}{[\begin{array}{*{20}{c}} {{V_{bi}}}\quad{{W_{bi}}}\quad{{V_{bj}}}\quad{{W_{bj}}}\quad{{V_{bk}}}\quad{{W_{bk}}} \end{array}]^T}$$.

With the application of the location matrix L , $$Q_d^s$$ is given by the following relationship: $$Q_d^s = {L^T}Q_b^s$$.

Combining Eqs. (19) and (20) and substituting for $$Q_b^s$$, we obtain the final form used in programming for a single rotor system as follows: $$\left[ {\begin{array}{*{20}{c}} {{M^s}}\quad0 \\ 0\quad{M_b^s} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\ddot u}^s}} \\ {\ddot u_b^s} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} { - \Omega {G^s} + {L^T}{C_O}L}\quad{ - {L^T}{C_O}{L_b}} \\ { - {C_O}L}\quad{C_b^s + {C_O}{L_b}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\dot u}^s}} \\ {\dot u_b^s} \end{array}} \right]\; + \left[ {\begin{array}{*{20}{c}} {{K^s} + {L^T}{K_O}L}\quad{ - {L^T}{K_O}{L_b}} \\ { - {K_O}L}\quad{K_b^s + {K_{O}}{L_b}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{u^s}} \\ {u_b^s} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {Q_u^s} \\ 0 \end{array}} \right]$$.

3.2.  Dual rotor system

To acquire the vibration equation of a dual rotor system (Fig. 7), we couple the equations of the two subsystems: $$\left[ {\begin{array}{*{20}{c}} {{M^{s1}}}\quad0 \\ 0\quad{{M^{s2}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\ddot u}^{s1}}} \\ {{{\ddot u}^{s2}}} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{\Omega _1}{G^{s1}}}\quad0 \\ 0\quad{{\Omega _2}{G^{s2}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\dot u}^{s1}}} \\ {{{\dot u}^{s2}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{K^{s1}}}\quad0 \\ 0\quad{{K^{s2}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{u^{s1}}} \\ {{u^{s2}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {Q_u^{s1}} \\ {Q_u^{s2}} \end{array}} \right] + Q_{d1}^s + Q_{d2}^s$$,

Dual rotor system

where $$Q_{d1}^s$$ are the generalized forces derived from the base bearings, and $$Q_{d2}^s$$ are the generalized forces caused by the inter-bearing. The reaction force from the base bearings is $${Q_{B1}} = {\left[ {\begin{array}{*{20}{c}} {Q_1^b}\quad{Q_6^b}\quad{Q_7^b} \end{array}} \right]^T} = - {C_{B1}}{L_1}{[{\dot u^{s1}}\;{\dot u^{s2}}]^T} - {K_{B1}}{L_1}{[{u^{s1}}\;{u^{s2}}]^T}$$, where $$L_1^{12 \times 40}{[{u^{s1}}\;{u^{s2}}]^T} = {\left[ {\begin{array}{*{20}{c}} {{u_1}}\quad{{u_6}}\quad{{u_7}} \end{array}} \right]^T}$$.

The force from the inter-bearing is $${Q_{B2}} = - {C_{B2}}{L_2}{[{\dot u^{s1}}\;{\dot u^{s2}}]^T} - {K_{B2}}{L_2}{[{u^{s1}}\;{u^{s2}}]^T}$$, where $$L_2^{8 \times 40}{[{u^{s1}}\;{u^{s2}}]^T} = {\left[ {\begin{array}{*{20}{c}} {{u_4}}\quad{{u_{10}}} \end{array}} \right]^T}$$, $${C_{B2}} = \left[ {\begin{array}{*{20}{c}} {C_{4\_10}^b}\quad{ - C_{4\_10}^b} \\ { - C_{4\_10}^b}\quad{C_{4\_10}^b} \end{array}} \right]$$, $${K_{B2}} = \left[ {\begin{array}{*{20}{c}} {K_{4\_10}^{b}}\quad{ - K_{4\_10}^{b}} \\ { - K_{4\_10}^{b}}\quad{K_{4\_10}^{b}} \end{array}} \right]$$.

Utilizing the location matrices L 1 and L 2, $$Q_{d1}^s$$ and $$Q_{d2}^s$$ can be represented as Eq. (24). The final form used in programming for a dual rotor system is as follows: $$\left( {\left[ {\begin{array}{*{20}{c}} {{K^{s1}}}\quad0 \\ 0\quad{{K^{s1}}} \end{array}} \right] + L_1^T{K_{B1}}{L_1} + L_2^T{K_{B2}}{L_2}} \right)\left[ {\begin{array}{*{20}{c}} {{u^{s1}}} \\ {{u^{s2}}} \end{array}} \right]\; - \left( {\left[ {\begin{array}{*{20}{c}} {{\Omega _1}{G^{s1}}}\quad0 \\ 0\quad{{\Omega _2}{G^{s2}}} \end{array}} \right] - L_1^T{C_{B1}}{L_1} - L_2^T{C_{B2}}{L_2}} \right)\left[ {\begin{array}{*{20}{c}} {{{\dot u}^{s1}}} \\ {{{\dot u}^{s2}}} \end{array}} \right]\; + \left[ {\begin{array}{*{20}{c}} {{M^{s1}}}\quad0 \\ 0\quad{{M^{s2}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\ddot u}^{s1}}} \\ {{{\ddot u}^{s2}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {Q_u^{s1}} \\ {Q_u^{s2}} \end{array}} \right]$$.

4.  Solution of system equations

4.1.  Natural frequencies

The single rotor finite element (FE) model and the coupled FE model can both be simplified to a unified form as Eq. (34). The associated homogeneous equation needs to be solved to determine the natural frequencies of the system, which is an eigenvalue problem (Lee et al., ). $$M\ddot u + C\dot u + Ku = Q$$.

In the case of little damping and symmetrical matrices, usually in structural dynamics there are some standard procedures applicable for this large-scale eigenvalue problem. However, in rotor dynamics, due to the effects of anisotropic bearings and gyroscopic moments, the damping and stiffness matrices become asymmetric. Asymmetric linear dynamic systems can be transformed into equivalent symmetric systems by using non-singular linear transformations (Adhikari, ; ). Here, a special method is introduced that use state vectors, and the homogeneous equation can be written as follows (Garadin, ): $$\dot v = Av$$, where $$A = \left[ {\begin{array}{*{20}{c}} { - {M^{ - 1}}C}\quad{ - {M^{ - 1}}K} \\ I\quad0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\dot u} \\ u \end{array}} \right],v = \left[ {\begin{array}{*{20}{c}} {\dot u} \\ u \end{array}} \right]$$.

The eigenvalues and eigenvectors of matrix A are the corresponding natural frequencies and mode shapes of the rotor-bearing system. The natural frequencies depend on the rotor operating speed due to the gyroscopic forces, so the critical speeds are given by the intersections of the synchronous line and the natural frequency curves.

4.2.  Transient analysis

Mathematically, Eq. (34) represents a system of linear differential equations of second order, and $$M\ddot u$$, $$C\dot u$$, $$Ku$$ are the inertia, damping and elastic forces, respectively, all of which are time-dependent. In principle, the solution to the equation can be obtained by standard procedures such as the Runge-Kutta method, but it is very time-consuming for large scale matrices. In this study, the Newmark method is adopted for its high efficiency and numerical stability (Bathe, ).

This direct integration method is based on the assumption that the acceleration varies linearly between two time instances. The resulting expressions for the velocity and displacement vectors are written as $${\dot u_{t + \Delta t}} = {\dot u_t} + (1 - {\lambda _1}){\ddot u_t}\Delta t + {\lambda _1}{\ddot u_{t + \Delta t}}\Delta t$$, $${u_{t + \Delta t}} = {u_t} + {\dot u_t}\Delta t + \left( {\frac{1}{2} - {\lambda _2}} \right){\ddot u_t}\Delta {t^2} + {\lambda _2}{\ddot u_{t + \Delta t}}\Delta {t^2}$$, where &lgr; 1 and &lgr; 2 are the parameters used to obtain integration stability. Generally, $${\lambda _2} \geqslant {({\lambda _1} + 0.5)^2}/4$$ and $${\lambda _1} \geqslant 0.5$$. $${\ddot u_{t + \Delta t}}$$ and $${\dot u_{t + \Delta t}}$$ can be expressed in terms of $${u_t}$$, $${\dot u_t}$$, $${\ddot u_t}$$ and $${u_{t + \Delta t}}$$ and then they are substituted into Eq. (39) for solving $${u_{t + \Delta t}}$$. $$M{\ddot u_{t + \Delta t}} + C{\dot u_{t + \Delta t}} + K{u_{t + \Delta t}} = {Q_{t + \Delta t}}$$.

4.3.  Law of rising speed

In industrial applications, two laws of rising speed are of interest (Fig. 8), namely, the linear law and the exponential law (Lalanne and Ferraris, ).

Fig.8
Start-up of linear law (a) and exponential law (b)

For the linear law: ${\kern 0pt} \dot \phi (t) = \left\{ \begin{matrix} \alpha t,\;\;0 \leqslant t \leqslant {t_1}, \hfill \\ {\Omega _c},\;\;{t_1} \leqslant t, \hfill \\ \end{matrix} \right.$ ${\kern 0pt} \phi (t) = \left\{ \begin{matrix} \frac{1}{2}\alpha {t^2},\;\;0 \leqslant t \leqslant {t_1}, \hfill \\ \frac{1}{2}\alpha t_1^2 + {\Omega _c}(t - {t_1}),\;\;{t_1} \leqslant t, \hfill \\ \end{matrix} \right.$ ${\kern 0pt} \ddot \phi (t) = \left\{ \begin{matrix} \alpha ,\;\;0 \leqslant t \leqslant {t_1}, \hfill \\ 0,\;\;{t_1} \leqslant t, \hfill \\ \end{matrix} \right.$ where α is the angular acceleration.

For the exponential law: $$\dot \phi (t) = {\Omega _c}(1 - {e^{ - \gamma t}})$$, $$\phi (t) = \int_0^t {{\Omega _c}(1 - {e^{ - \gamma t'}})dt'}$$, $$\ddot \phi (t) = {\Omega _c}\gamma t{e^{ - \gamma t}}$$.

Due to the existence of angular acceleration, the excitation forces caused by unbalance mass become more complicated during transient motion, and the expression is $$\hat Q = {m_D}e\left[ {\begin{array}{*{20}{c}} {{\Omega ^2}\cos (\phi (t)) - \ddot \phi (t)\sin (\phi (t))} \\ {{\Omega ^2}\sin (\phi (t)) + \ddot \phi (t)\cos (\phi (t))} \\ 0 \\ 0 \end{array}} \right]$$.

5.  Numerical examples

To demonstrate the application and accuracy of the proposed program, simulation on a typical dual rotor system was carried out to determine the critical speed and transient response. Fig. 7 outlines the structure and the following Tables 13 list the detailed parameters.

Table 1

Shaft element data
 Shaft Element Length (cm) Inner radius (cm) Outer radius (cm) 1 1 10 0 1.1 2 20 0 1.1 3 15 0 1.1 4 5 0 1.1 5 10 0 1.1 2 6 7.5 2 2.6 7 15 2 2.6 8 7.5 2 2.6

• E=2.1×1011 Pa, μ=0.3, ρ=7800 kg/m3

• Table 2

Concentrated disk data
 Shaft Station Mass (kg) Polar inertia (kg·m2) Diametral inertia (kg·m2) 1 2 9.683 0.04900 0.02450 5 9.683 0.04900 0.02450 2 8 9.139 0.04878 0.02439 9 9.139 0.04878 0.02439

Table 3

Bearing data
 Bearing No. Station Stiffness (N/m) 1 1 2.60×107 2 6 1.75×107 3 7 1.75×107 4 4, 10 8.75×106

There are two operation modes in dual rotor systems: co-rotation and counter-rotation. The speed ratio is $${r_\omega } = {{{\omega _2}} \mathord{\left/ {\vphantom {{{\omega _2}} {{\omega _1}}}} \right. \kern-\nulldelimiterspace} {{\omega _1}}}$$ in co-rotation, and $${r_\omega } = - {{{\omega _2}} \mathord{\left/ {\vphantom {{{\omega _2}} {{\omega _1}}}} \right. \kern-\nulldelimiterspace} {{\omega _1}}}$$ in counter-rotation.

Gyroscopic moments are related to the direction of rotation, so research on the dynamic characteristics in both co-rotation and counter-rotation was conducted. Campbell diagrams of these two cases were obtained by calculation of natural frequencies. Figs. 9 and 10 show the variation in whirl frequency with inner rotor speed when rω =1.5 and rω =−1.5, respectively.

Fig.9
Campbell diagram of co-rotation in a dual rotor

Fig.10
Campbell diagram of counter-rotation in a dual rotor

Table 4 lists the data for the corresponding intersections. Because the influence of gyroscopic moments on dual rotor structures may decrease in counter-rotation, the critical speeds are smaller.

Table 4

Critical speeds
 Order Critical speeds (r/min) Inner rotor excitation Outer rotor excitation rω =1.5 rω =−1.5 rω =1.5 rω =−1.5 1 2429 2188 2309 2050 2 6026 5895 5914 5509 3 12 060 12 080 12 030 11 850

One common phenomenon that can be found in Figs. 9 and 10 is that the third order frequencies of both forward and backward whirl barely vary with rotation speed, while those of other orders vary significantly. In this specific case, this phenomenon indicates that the third order frequencies depend primarily on the bearing stiffness, while those of other orders depend more on the stiffness of the rotors.

In general, the stiffness of supports has a great effect on vibration performance; the inter-bearing stiffness, especially, is one of the key design parameters of a dual rotor system. The critical speeds of the first three orders were calculated with a set of such parameters, and the results are plotted in Fig. 11. The critical speeds of the first two orders are practically constant, while that of the third order increases. Furthermore, when counter-rotating, the critical speeds become more sensitive to the inter-bearing stiffness.

Fig.11
Effect of inter-bearing stiffness on critical speeds

A scanning calculation was carried out to find the resonant speeds of different rotating modes. Holding ω 1 at a constant value while varying ω 2, the corresponding critical speeds can be obtained from the Campbell diagrams (Table 5). Similarly, another group of data can be obtained with varying ω 1 and fixed ω 2. The critical speed map of the dual rotor system is drawn according to these two sets of results.

Table 5

Critical speeds for constant ω 1
 ω 1 Critical speed (r/min) 1st 2nd 3rd 0 2171 5720 12 029 2000 2352 5822 12 030 4000 2532 5917 12 030 6000 2707 6003 12 031 8000 2868 6082 12 031 10 000 3025 6155 12 031 12 000 3172 6222 12 032 14 000 3307 6282 12 032 16 000 3422 6338 12 032 18 000 3535 6389 12 032 20 000 3630 6436 12 032

In Fig. 12, the data for the intersections when rω =1.5 are approximately equal to those listed in Table 4. For the benefit of predicting resonance points, this kind of map is very helpful for choosing rω during the design process. It can be seen from the map or Table 5 that the speed ratio has almost no influence on the third order critical speed. The above analysis and the conclusion from Fig. 11 make it clear that the bearing stiffness has a significant impact on the natural frequencies of some specific orders.

Fig.12
Critical speed map of a dual rotor system

Transient simulation was carried out to examine the dynamic characteristics when a rotor system starts up and goes through the resonance region with different laws of rising speed. The parameters used were: simulation time: 0–6 s, time step: 0.0004 s, damping coefficient: β c=0.001, unbalance mass on node 5: f u=m u e=50 g·mm.

For the linear law:

For the exponential law:

The two sets of data in Fig. 13 are displacements of node 5 in the Y direction, corresponding to the linear and exponential laws, respectively. Similarly, the results in Fig. 14 are the transient amplitude. In practice, the exponential law is applied mostly for passing through the resonance region smoothly, as shown in the figure. Transient analysis can also be used to evaluate the stability of systems.

Fig.13
Transient displacements of node 5 for the linear law (a) and exponential law (b)

Fig.14
Transient amplitudes of node 5 for the linear law (a) and exponential law (b)

The whirl orbit is one of the most important features for representing the motion state and it is the essential theme in the time domain analysis of vibration signals. To investigate the influence of different damping coefficients, the laws of rising speed and the gravity effect, the following simulations were completed. The particular conditions are listed in Table 6. Four trajectory diagrams of node 5 are illustrated in Fig. 15. Figs. 15a and 15b present the whirl behavior with and without, respectively, the gravity considered. Figs. 15b and 15c indicate that β c affects the amplitude response explicitly. Fig. 15d shows the orbit of the entire start-up, showing that the stabilization occurs after passage through the critical speed.

Table 6

Simulation conditions
 Case Gravity β c Speed up Time plot (s) Fig. 15a Not 0.001 Linear 5.8–6 Fig. 15b Yes 0.001 Linear 5.8–6 Fig. 15c Yes 0.0002 Linear 5.8–6 Fig. 15d Yes 0.001 Exponential 0–6

Fig.15
Trajectory diagrams of node 5
The particular conditions of Figs .15a–15d are listed in Table 6

6.  Conclusions

The vibration equation of a dual rotor system was derived in detail, and by combining it with the Newmark method, a capable and efficient program was developed to solve its natural frequencies and simulate the transient response. FEM is practicable and effective for the engineering analysis of dual rotor systems, because the motion equations of subsystems can be coupled conveniently.

Gyroscopic moments will lead to different critical speeds in co-rotation and counter-rotation. During counter-rotation, the critical speeds decrease since the gyroscopic moments will be partially canceled.

The speed ratio rω will affect the dynamic characteristics of dual rotor system. A scanning calculation can be applied to draw a critical speed map, which is helpful for choosing this parameter during the design process.

From the analysis of the transient response, we found that the exponential law was better than the linear law because of the smaller amplitude. The stability of the rotor bearing system can also be estimated by the orbits.

References

[1] Abduljabbar, Z., ElMadany, M.M., AlAbdulwahab, A.A., 1996. Active vibration control of a flexible rotor. Computers and Structures, 58(3):499-511.

[2] Adhikari, S., 1999. Modal analysis of linear asymmetric non-conservative systems. Journal of Engineering Mechanics, 125(12):1372-1379.

[3] Adhikari, S., 2000. On symmetrizable systems of second kind. ASME Journal of Applied Mechanics, 67(4):797-802.

[4] Bathe, K.J., 1996.  Finite Element Procedures. Prentice-Hall Inc.,New Jersey :768-784.

[5] Carrella, A., Friswell, M.I., Zotov, A., Ewins, D.J., Tichonov, A., 2009. Using nonlinear springs to reduce the whirling of a rotating shaft. Mechanical Systems and Signal Processing, 23(7):2228-2235.

[6] Das, A.S., Nighil, M.C., Dutt, J.K., Irretier, H., 2008. Vibration control and stability analysis of rotor-shaft system with electromagnetic exciters. Mechanism and Machine Theory, 43(10):1295-1316.

[7] Friswell, M.I., Dutt, J.K., Adhikari, S., Lees, A.W., 2010. Time domain analysis of a viscoelastic rotor using internal variable models. International Journal of Mechanical Sciences, 52(10):1319-1324.

[8] Garadin, M., Kill, N., 1983. A new approach to finite element modeling of flexible rotors. Engineering Computations, 1(1):52-64.

[9] Chiang, H.W.D., Hsu, C.N., Jeng, W., 2002. Turbomachinery Dual Rotor-Bearing System Analysis. ASME Turbo Expo 2002: Power for Land, Sea, and Air, Amsterdam, the Netherlands :803-810.

[10] Jei, Y.G., Lee, C.W., 1988. Finite element model of asymmetrical rotor-bearing systems. Journal of Mechanical Science and Technology, 2(2):116-124.

[11] Lalanne, M., Ferraris, G., 1990.  Rotor Dynamics Prediction in Engineering. Wiley,New York :77-116.

[12] Lee, A.S., Ha, J.W., Choi, D.H., 2003. Coupled lateral and torsional vibration characteristics of a speed increasing geared rotor-bearing system. Journal of Sound and Vibration, 263(4):725-742.

[13] Nelson, H.D., 1980. A finite rotating shaft element using Timoshenko beam theory. ASME Journal of Mechanical Design, 102(4):793-803.

[14] Nelson, H.D., Mc-Vaugh, J.M., 1976. The dynamics of rotor-bearing systems using finite elements. ASME Journal of Engineering for Industry, 98(2):593-600.

[15] Newmark, N.M., 1959. A method of computation for structural dynamics. ASCE Journal of Engineering Mechanics, 85(3):67-94.

[16] Rajan, M., Nelson, H.D., Chen, W.J., 1986. Parameter sensitivity in the dynamics of rotor-bearing systems. ASME Journal of Vibration, Acoustics, Stress, and Reliability in Design, 108(2):197-206.

[17] Rao, J.S., 1996.  Rotor Dynamics (3rd Ed.). New Age International Publishers,London :69-106.

[18] Rao, J.S., Shiau, T.N., Chang, J.R., 1998. Theoretical analysis of lateral response due to torsional excitation of geared rotors. Mechanism and Machine Theory, 33(6):761-783.

[19] Shiau, T.N., Rao, J.S., Chang, J.R., Choi, S.T., 1999. Dynamic behavior of geared rotors. ASME Journal of Engineering for Gas Turbines and Power, 121(3):494-503.

[20] Zhong, Y.E., 1987.  Rotor Dynamics. (in Chinese), Tsinghua University Press,Beijing, China :176-195.