Vibration Properties of Dual-Rotor Systems under Base Excitation, Mass Unbalance and Gravity

Abstract

Rotor systems installed in a transportation system or under seismic excitations are considered to have a moving base. Although extensive research has been conducted on the dynamic behavior of the single-rotor system under base motions, few studies have dealt with the dynamics of dual-rotor systems, especially the counter-rotating dual-rotor systems used in airplane engines. Moreover, mass unbalance and gravity are unavoidable excitations for most rotor systems. Therefore, the vibration properties of a counter-rotating dual-rotor system with the coupled effects of base motions, mass unbalance and gravity are investigated in this paper for the first time. Using the Lagrange principle associated with the finite element method, a general model for dual-rotor systems under base motions was established by using Timoshenko beam elements, leading to a detailed analysis of the natural properties and harmonic responses of the system. The results revealed that different whirling modes (backward, forward or both) may be mutually excited. This research can be helpful for the design and vibration analysis of dual-rotor systems concerned with base motion.

1. Introduction

Rotating machines such as pumps, generators and ground gas turbines are fundamental equipment in the modern engineering applications. Gravity and mass unbalance from center-of-mass eccentricity along the rotor axis are unavoidable excitation sources for analyzing these rotor system problems. Moreover, base excitations should also be concerned for a rotor system mounted on a transportation system or under seismic excitation. These three excitations have mutual influences that induce intricate resonance problems. It is even more complicated for a counter-rotating dual-rotor system with inter-shaft bearings in airplane engines, compared with a single rotor system. Thus, it is important to reveal the vibration properties in detail for the dual-rotor system under base motions, mass unbalance and gravity.

As mentioned in the authors’ former work [1], many researchers have contributed to our understanding of the dynamics of the rotor system under base motions. The accurate description of base motions and simulated deformation of the rotor are the two priorities for rotor system modeling. For the former task, two groups of Euler angles were used to describe the relative rotations among the rotor, base and ground, and coordinate transformation theory made the model adaptable to arbitrary base motion as shown by, for example, Das [2,3], Dakel [4,5] and Chen [1]. For the latter task, the application of the finite element method made the model more realistic [4,5] compared with the Jeffcott rotor [6,7,8,9,10,11], a rigid disk with only four degrees of freedom [12,13,14,15,16,17,18,19], and the rigid disk in the Rayleigh–Ritz method for the shaft [20,21,22,23]. Especially compared with the Euler–Bernoulli beam element [24,25,26,27,28,29], the Timoshenko beam element [4,5,12,13,30,31,32,33,34,35,36,37] took the shear deformation of the shaft into account. Therefore, it was used to improve the model established by the author in [1].

Various types of base motion have been widely investigated to find out their effects on the dynamic behavior of the rotor systems. Random and shock excitations to bases have been studied by scholars to help maintain the stable operation of rotating machines under seismic excitations [6,12,13,14,30]. Transient or aperiodic base motion is involved mainly due to aircraft flight maneuvers [7,8,9,10,11,16,17,18]. Research into this aspect of the rotor system focuses on the effects of additional centrifugal force and gyroscopic momentum from the base motion on the nonlinear responses, which often refer to nonlinear factors. The periodic base motions are attracting more and more attention from scholars in recent years. The singular or combined effects of constant base rotation, harmonic base translation and harmonic base rotation are discussed in detail in [4,5,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37]. The specific stability problem induced by the nonlinearity from time-varying base rotations is a major topic of investigation [20,21,22,25,27,32,33]. In addition, some scholars [2,3,28,29] are devoted to minimizing the vibration induced by base motions by using an active control method or magnetic bearings. However, all the above studies are for single-rotor systems.

Extensive research has been conducted on the dynamic behavior of dual-rotor systems with a fixed base. Hibner [38] presented an efficient analytical technique to predict the vibratory response of an engine with nonlinear viscous damping. He pointed out that the aircraft engine structure must be analyzed in detail to provide all critical speeds and responses to imbalance. Glasgow et al. [39] performed a method of component mode synthesis for the analysis of multi-shaft rotor-bearing systems and maintained the essential dynamic characteristics of the lower modes after reducing the size of the overall problem. Li et al. [40] employed a component mode method to the dynamic analysis of large flexible rotor systems and a two-spool test rig was used to illustrate the capability and versatility of the component mode method. Gupta et al. [41] investigated the critical speeds, mode shape and unbalance responses both numerically and experimentally. The results are in reasonable agreement and show the appearance of a cross-excitation between the inner and outer shafts. Ferraris et al. [42] used a simple model of non-symmetric co-axial co-rotating or counter-rotating rotors to give analytical expression of the critical speeds, mass unbalance responses and sense of whirl. Chiang et al. [43] showed two methods to compute the critical speeds for a co-rotating dual-rotor system that has different advantages. In addition, the speed ratio and intershaft bearing stiffness were found out to be two important factors for the dynamic behavior of the system. Guskov et al. [44] examined the dynamics of a dual-rotor system both experimentally and numerically, and pointed out that taking into consideration the tilt stiffness in the system can improve the theoretical model. Zhang et al. [45] put forward a critical speed prediction method (similar to that of Chiang [43]) for a dual-rotor system based on the finite element method. Moreover, a vibration mode selection technique was presented to identify the backward and forward whirling modes.

Few studies have included base motions in the dynamic behavior analysis of a dual-rotor system. Chen et al. [46,47], though, did develop a dynamic model that included mass unbalance and base motions, using the Lagrange principle and finite element method. Steady-state characteristics during stable operation, and transient characteristics during the start-up of the co-rotating dual-rotor system were investigated adequately. Constant base rotation, harmonic translation and mass unbalance were the major forces involved in the investigation. However, the counter-rotating type for the dual-rotor system and excitations from gravity were not taken into account.

In a word, the synthesized influences of base motions, mass unbalance and gravity on the dynamic behavior of counter-rotating dual-rotor systems still lack in-depth investigation. Therefore, this article develops a general model for such a dynamic analysis under six degrees of base motions based on the Lagrange principle and finite element method and using Timoshenko beam elements to model the shaft. Constant base rotation, harmonic base translation, mass unbalance force and gravity were taken into account. The natural properties and harmonic responses of the counter-rotating dual-rotor system are analyzed in detail.

2. Modeling for a Dual-Rotor System under Base Motions

As shown in Figure 1, a counter-rotating dual-rotor system is mounted on a moving rigid base (the blue dashed box), which consists of a low-pressure rotor with a low rotating speed ${\mathsf{\Omega }}_{LP}$ and a high-pressure rotor with a higher rotating speed ${\mathsf{\Omega }}_{HP}$. The excitations due to the base motions, mass unbalance and gravity are taken into consideration. Three coordinate systems are defined in this paper: an inertial frame of reference $F_0\left(X_0,Y_0,Z_0\right)$, which is fixed to the ground; a non-inertial frame of reference $F_b\left(X_b,Y_b,Z_b\right)$, which moves with the rigid base; and a local frame of reference $F_l\left(X_l,Y_l,Z_l\right)$, which moves with the rotor. All three coordinate systems have relative movement to each other. For the dual-rotor model in this paper, longitudinal vibration and rotor torsion are ignored.

2.1. Modeling for the Disks, Unbalance Force and Linear Bearings

The modeling for the disks, unbalance force and linear bearings is similar to that in [1] and is quite simplified in this paper. Three base translations ($x_b$, $y_b$ and $z_b$) and three base rotations (pitching angle $\alpha$, rolling angle $\beta$ and yawing angle $\gamma$) are shown in Figure 2a. Two transverse displacements ($x_d$ and $z_d$) and two angular displacements (${\theta }_d$ and ${\psi }_d$) relative to the base for the disk are unknown and remain to be solved as shown in Figure 2b. Firstly, the kinetic energy $T_d$ and the gravitational potential energy $U_{dg}$ of the disk and the kinetic energy $T_u$ of the unbalanced mass are obtained as

$$T_d=\frac{1}{2}{m}_d{V}_{\mathrm{d}}{}^{\mathrm{T}}{V}_{\mathrm{d}}+\frac{1}{2}{\mathsf{\omega }}_{\mathrm{d}}{}^{\mathrm{T}}{I}_{\mathrm{d}}{\mathsf{\omega }}_{\mathrm{d}}$$

(1)

$$U_{dg}=m_{d}g{w}_{dg}$$

(2)

$$T_u=\frac{1}{2}{m}_u{V}_{\mathrm{u}}{}^{\mathrm{T}}{V}_{\mathrm{u}}$$

(3)

where $I_{\mathrm{d}}=\mathrm{diag}\left(\left[\begin{array}{ccc}{I}_{dx}& I_{dy}& I_{dz}\end{array}\right]\right)$; $m_d$ and $I_{\mathrm{d}}$ are the mass and mass-moment matrix of inertia of the disk; $V_{\mathrm{d}}$ and ${\mathsf{\omega }}_{\mathrm{d}}$ are the absolute translational velocity and absolute angular velocity of the disk, respectively; $g$ is the gravitational acceleration; $w_{dg}$ is the absolute vertical displacement of the disk; $m_u$ denotes the unbalanced mass; $V_{\mathrm{u}}$ is the absolute linear velocity of the unbalanced mass; and the superscript $\mathrm{T}$ is the matrix transposition symbol.

Secondly, the Lagrange principle is applied to achieve the final equations of motion (EOMs) of the disk as following

$$m_{\mathrm{d}}{\ddot{q}}_{\mathrm{d}}+\left({\mathsf{\Omega }}_r{g}_{\mathrm{d}0}+c_{\mathrm{d}}\right){\dot{q}}_{\mathrm{d}}+k_{\mathrm{d}}{q}_{\mathrm{d}}=f_{\mathrm{dtr}}+f_{\mathrm{dlb}}+f_{\mathrm{dan}}+f_{\mathrm{dg}}$$

(4)

where $q_{\mathrm{d}}$ is the general displacement vector of the disk; $m_{\mathrm{d}}$ and $g_{\mathrm{d}0}$ are the conventional mass and gyroscopic matrices of the disk, respectively; ${\mathsf{\Omega }}_r$ is ${\mathsf{\Omega }}_{LP}$ or ${\mathsf{\Omega }}_{HP}$ corresponding to the low-pressure or high-pressure rotor; $c_{\mathrm{d}}$ is the additional damping matrix of the disk due to base rotations; $f_{\mathrm{dtr}}$, $f_{\mathrm{dan}}$ and $f_{\mathrm{dg}}$ are external forces induced by base linear motions, base angular motions and gravity, respectively; and $f_{\mathrm{dlb}}$ is an additional external force due to both the base rotations and the position vector from the origin of the frame $F_b$ to the center of the left bearing.

$$c_{\mathrm{d}}=c_{\mathrm{d}1}{\omega }_{by}$$

(5)

$$k_{\mathrm{d}}\left(t\right)=k_{\mathrm{d}1}{\dot{\omega }}_{by}+k_{\mathrm{d}2}{\omega }_{bx}^2+k_{\mathrm{d}3}{\omega }_{by}^2+k_{\mathrm{d}4}{\omega }_{bz}^2+k_{\mathrm{d}5}{\omega }_{bx}{\omega }_{bz}+k_{\mathrm{d}6}{\mathsf{\Omega }}_r{\omega }_{by}$$

(6)

$$f_{\mathrm{dtr}}\left(t\right)=-p_{\mathrm{d}1}{a}_{bx}-p_{\mathrm{d}2}{a}_{bz}$$

(7)

$$f_{\mathrm{dlb}}\left(t\right)=p_{\mathrm{d}1}\left(x_{lb}\left({\omega }_{by}^2+{\omega }_{bz}^2\right)+y_{lb}\left({\dot{\omega }}_{bz}-{\omega }_{bx}{\omega }_{by}\right)-z_{lb}\left({\dot{\omega }}_{by}+{\omega }_{bx}{\omega }_{bz}\right)\right)$$

(8)

$$\begin{gathered} f_{\mathrm{dlb}}\left(t\right)=p_{\mathrm{d}1}\left(x_{lb}\left({\omega }_{by}^2+{\omega }_{bz}^2\right)+y_{lb}\left({\dot{\omega }}_{bz}-{\omega }_{bx}{\omega }_{by}\right)-z_{lb}\left({\dot{\omega }}_{by}+{\omega }_{bx}{\omega }_{bz}\right)\right) \\ +p_{\mathrm{d}2}\left(x_{lb}\left({\dot{\omega }}_{by}-{\omega }_{bx}{\omega }_{bz}\right)-y_{lb}\left({\dot{\omega }}_{bx}+{\omega }_{by}{\omega }_{bz}\right)+z_{lb}\left({\omega }_{bx}^2+{\omega }_{by}^2\right)\right) \end{gathered}$$

(9)

$$f_{\mathrm{dan}}\left(t\right)=s_{\mathrm{d}1}{\dot{\omega }}_{bx}+s_{\mathrm{d}2}{\dot{\omega }}_{bz}+s_{\mathrm{d}3}{\omega }_{bx}{\omega }_{by}+s_{\mathrm{d}4}{\omega }_{by}{\omega }_{bz}+s_{\mathrm{d}5}{\mathsf{\Omega }}_r{\omega }_{bx}+s_{\mathrm{d}6}{\mathsf{\Omega }}_r{\omega }_{bz}$$

(10)

$$f_{\mathrm{dg}}\left(t\right)=p_{\mathrm{d}1}g\cos \alpha \sin \beta -p_{\mathrm{d}2}g\cos \alpha \cos \beta$$

(11)

The final unbalancing forces are given by

$$f_{\mathrm{du}}=m_{u}e\left[\begin{array}{c}\left({\left({\mathsf{\Omega }}_r+{\omega }_{by}\right)}^2+{\omega }_{bz}^2\right)\cos \left({\mathsf{\Omega }}_{r}t+{\phi }_u\right)+\left({\dot{\omega }}_{by}+{\omega }_{bx}{\omega }_{bz}\right)\sin \left({\mathsf{\Omega }}_{r}t+{\phi }_u\right)\\ -\left({\left({\mathsf{\Omega }}_r+{\omega }_{by}\right)}^2+{\omega }_{bx}^2\right)\sin \left({\mathsf{\Omega }}_{r}t+{\phi }_u\right)+\left({\dot{\omega }}_{by}-{\omega }_{bx}{\omega }_{bz}\right)\cos \left({\mathsf{\Omega }}_{r}t+{\phi }_u\right)\\ 0\\ 0\end{array}\right]$$

(12)

where $e$ denotes the eccentricity of the unbalanced mass, and ${\phi }_u$ is the unbalanced angle with respect to the axis $O_b{X}_b$.

The linear bearings are concerned by placing the stiffness $k_b$ and damping $c_b$ of the bearings in the proper position in the stiffness matrices and damping matrices.

2.2. EOMs for the Shaft

A two-node Timoshenko beam replaced the Euler beam element from [1] to model the shear deformation effect on the shaft. The two nodes of the shaft element were $C_{s,j}$ and $C_{s,j+1}$, and their position vectors with respect to the base frame were defined as ${\left(O_b{C}_{s,j}\right)}_{F_b}={\left\{\begin{array}{ccc}{x}_{s,j}& y_{s,j}& z_{s,j}\end{array}\right\}}^{\mathrm{T}}$ and ${\left(O_b{C}_{s,j+1}\right)}_{F_b}={\left\{\begin{array}{ccc}{x}_{s,j+1}& y_{s,j+1}& z_{s,j+1}\end{array}\right\}}^{\mathrm{T}}$. The length of the shaft element was denoted as $l$. After the discretization of the shaft, the values of $y_{s,j}$, $y_{s,j+1}$ and $l$ remained constant. The general coordinate vector of the shaft element was

$$q_{\mathrm{s}}={\left[\begin{array}{cccccccc}{x}_{s,j}& z_{s,j}& {\theta }_{s,j}& {\psi }_{s,j}& x_{s,j+1}& z_{s,j+1}& {\theta }_{s,j+1}& {\psi }_{s,j+1}\end{array}\right]}^{\mathrm{T}}$$

(13)

After taking an arbitrary point $C_s$ with a distance of $y$ along $O_b{Y}_b$ from $C_{s,j}$ into consideration, its position vector was defined with respect to the base frame as ${\left(O_b{C}_s\right)}_{F_b}={\left\{\begin{array}{ccc}{x}_s& y_s& z_s\end{array}\right\}}^{\mathrm{T}}$, where $y_s=y_{s,j}+y$; $x_s$ and $z_s$ are unknowns; and ${\theta }_s$ and ${\psi }_s$ were attained by interposition:

$${\left[\begin{array}{cc}{x}_s\left(y\right)& z_s\left(y\right)\end{array}\right]}^{\mathrm{T}}=\mathsf{\psi }\left(y\right)q_{\mathrm{s}}$$

(14)

$${\left[\begin{array}{cc}{\theta }_s\left(y\right)& {\psi }_s\left(y\right)\end{array}\right]}^{\mathrm{T}}=\Phi \left(y\right)q_{\mathrm{s}}$$

(15)

For the Timoshenko beam element, the interposition function is given by

$$\mathsf{\psi }\left(y\right)=\left[\begin{array}{cccccccc}{\psi }_1& 0& 0& -{\psi }_2& {\psi }_3& 0& 0& -{\psi }_4\\ 0& {\psi }_1& {\psi }_2& 0& 0& {\psi }_3& {\psi }_4& 0\end{array}\right]$$

(16)

$$\Phi \left(y\right)=\left[\begin{array}{cccccccc}0& {\varphi }_1& {\varphi }_2& 0& 0& {\varphi }_3& {\varphi }_4& 0\\ -{\varphi }_1& 0& 0& {\varphi }_2& -{\varphi }_3& 0& 0& {\varphi }_4\end{array}\right]$$

(17)

where

$$s=\frac{y}{l},\Phi =\frac{12}{l^2}\frac{E{I}_s}{\kappa G{A}_s},G=\frac{E}{2\left(1+\mu \right)},A_s=\pi \left(r_{so}^2-r_{si}^2\right)$$

(18)

$$\kappa =\frac{6{\left(r_{si}^2+r_{so}^2\right)}^2{\left(1+\mu \right)}^2}{7r_{si}^4+34r_{si}^2{r}_{so}^2+7r_{so}^4+\mu \left(12r_{si}^4+48r_{si}^2{r}_{so}^2+12r_{so}^4\right)+{\mu }^2\left(4r_{si}^4+16r_{si}^2{r}_{so}^2+4r_{so}^4\right)}$$

(19)

$${\psi }_1=\frac{1}{1+\Phi }\left(1-3s^2+2s^3+\Phi \left(1-s\right)\right),{\psi }_2=\frac{ls}{1+\Phi }\left(1-2s+s^2+\frac{1}{2}\Phi \left(1-s\right)\right)$$

(20)

$${\psi }_3=\frac{s}{1+\Phi }\left(3s-2s^2+\Phi \right),{\psi }_4=\frac{ls}{1+\Phi }\left(-s+s^2-\frac{1}{2}\Phi \left(1-s\right)\right)$$

(21)

$${\varphi }_1=\frac{6s}{l\left(1+\Phi \right)}\left(-1+s\right),{\varphi }_2=\frac{1}{1+\Phi }\left(1-4s+3s^2+\Phi \left(1-s\right)\right)$$

(22)

$${\varphi }_3=\frac{6s}{l\left(1+\Phi \right)}\left(1-s\right),{\varphi }_4=\frac{s}{1+\Phi }\left(-2+3s+\Phi \right)$$

(23)

The absolute linear velocity $V_{\mathrm{s}}$ and absolute angular velocity ${\mathsf{\omega }}_{\mathrm{s}}$ of the point $C_s$ can be derived in similar steps as the disk. Then, the kinetic energy of a shaft element is given by

$$T_s=\frac{1}{2}{\rho }_s{A}_s\underset{0}{\overset{l}{\int }}\left(V_{\mathrm{s}}{}^{\mathrm{T}}{V}_{\mathrm{s}}\right)\mathrm{d}y+\frac{1}{2}\underset{0}{\overset{l}{\int }}\left({\mathsf{\omega }}_{\mathrm{s}}{}^{\mathrm{T}}{I}_{\mathrm{s}}{\mathsf{\omega }}_{\mathrm{s}}\right)\mathrm{d}y$$

(24)

where $I_{\mathrm{s}}=\mathrm{diag}\left(\left[\begin{array}{ccc}{i}_d& i_p& i_d\end{array}\right]\right)$; ${\rho }_s$, $A_s$, $i_d$ and $i_p$ are, respectively, the material density, area of shaft cross-section, the transverse and polar mass-moment of inertia per unit length of the shaft element. For a circular shaft, the $i_p=2i_d$ and $i_d=\rho I_s$ where $I_s=\pi \left(r_{so}^4-r_{si}^4\right)/4$ ($r_{so}$ is the shaft outer radius and $r_{si}$ is the shaft inner radius). The strain energy of a shaft element is given by

$$U_{sb}=\frac{1}{2}\underset{0}{\overset{l}{\int }}E{I}_s\left[{\left(\frac{\partial {\theta }_s}{\partial y}\right)}^2+{\left(\frac{\partial {\psi }_s}{\partial y}\right)}^2\right]\mathrm{d}y+\frac{1}{2}\underset{0}{\overset{l}{\int }}\kappa G{A}_s\left[{\left(\frac{\partial x_s}{\partial y}+{\psi }_s\right)}^2+{\left(\frac{\partial z_s}{\partial y}-{\theta }_s\right)}^2\right]\mathrm{d}y$$

(25)

where $E$ is the Young’s modulus of the shaft material.

Let $z_{s0}$ be the component of the absolute position vector ${\left(O_0{C}_s\right)}_{F_0}$ along the axis $O_0{Z}_0$. The gravitational potential energy of a shaft element is given by

$$U_{sg}={\rho }_s{A}_{s}g\underset{0}{\overset{l}{\int }}{z}_{s0}\mathrm{d}y$$

(26)

After substituting $T_s$, $U_{sb}$ and $U_{sg}$ into the Lagrange equation, simplifying ($\sin {\theta }_s\approx {\theta }_s$, $\sin {\psi }_s\approx {\psi }_s$, $\cos {\theta }_s\approx 1$, $\cos {\psi }_s\approx 1$) and ignoring the high-order terms, the EOMs for a shaft element was accomplished as

$$m_{\mathrm{s}}{\ddot{q}}_{\mathrm{s}}+\left({\mathsf{\Omega }}_r{g}_{\mathrm{s}0}+c_{\mathrm{s}}\right){\dot{q}}_{\mathrm{s}}+\left(k_{\mathrm{s}0}+k_{\mathrm{s}}\right)q_{\mathrm{s}}=f_{\mathrm{str}}+f_{\mathrm{slb}}+f_{\mathrm{san}}+f_{\mathrm{sg}}$$

(27)

where $m_{\mathrm{s}}$, $g_{\mathrm{s}0}$ and $k_{\mathrm{s}0}$ are the classical mass, gyroscopic and stiffness matrices of the Timoshenko beam element, respectively; $c_{\mathrm{s}}$ and $k_{\mathrm{s}}$ are additional time-varying damping and stiffness matrices due to the time-variable base rotations; $f_{\mathrm{str}}$, $f_{\mathrm{san}}$ and $f_{\mathrm{sg}}$ are external forces motivated by base translations, base rotations and gravity, respectively; and $f_{\mathrm{slb}}$ are additional external forces due to both the position vector $O_b{O}_{brg1}$ in Figure 2 and base rotations. The origin of the base frame is set right in the center of bearing 1; thus, $f_{\mathrm{slb}}=\mathbf{0}$. The precise coefficient matrices and force vectors are listed as follows:

$$m_{\mathrm{s}}=\underset{0}{\overset{l}{\int }}\rho A_s\mathsf{\psi }{\left(y\right)}^{\mathrm{T}}\mathsf{\psi }\left(y\right)\mathrm{d}y+i_d\underset{0}{\overset{l}{\int }}\Phi {\left(y\right)}^{\mathrm{T}}\Phi \left(y\right)\mathrm{d}y$$

(28)

$$g_{\mathrm{s}0}=i_p\underset{0}{\overset{l}{\int }}\Phi {\left(y\right)}^{\mathrm{T}}\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right]\Phi \left(y\right)\mathrm{d}y$$

(29)

$$\Gamma \left(y\right)={\psi }^{\text{'}\left(y\right)}+\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]\Phi \left(y\right)$$

(30)

$$k_{\mathrm{s}0}=\underset{0}{\overset{l}{\int }}E{I}_s\Phi {\left(y\right)}^{\mathrm{T}}\Phi \left(y\right)\mathrm{d}y+\underset{0}{\overset{l}{\int }}\kappa G{A}_s\Gamma {\left(y\right)}^T\Gamma \left(y\right)\mathrm{d}y$$

(31)

$$c_{\mathrm{s}}=\left(2\underset{0}{\overset{l}{\int }}\rho A_s\psi {\left(y\right)}^{\mathrm{T}}\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]\psi \left(y\right)\mathrm{d}y+\left(i_p-2i_d\right)\underset{0}{\overset{l}{\int }}\Phi {\left(y\right)}^{\mathrm{T}}\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right]\Phi \left(y\right)\mathrm{d}y\right){\omega }_{by}$$

(32)

$$k_{\mathrm{s}}=k_{\mathrm{s}1}{\dot{\omega }}_{by}+k_{\mathrm{s}2}{\omega }_{bx}^2+k_{\mathrm{s}3}{\omega }_{by}^2+k_{\mathrm{s}4}{\omega }_{bz}^2+k_{\mathrm{s}5}{\omega }_{bx}{\omega }_{bz}+k_{\mathrm{s}6}{\mathsf{\Omega }}_r{\omega }_{by}$$

(33)

$$k_{\mathrm{s}1}=\underset{0}{\overset{l}{\int }}\rho A_s\mathsf{\psi }{\left(y\right)}^{\mathrm{T}}\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]\mathsf{\psi }\left(y\right)\mathrm{d}y+\underset{0}{\overset{l}{\int }}\Phi {\left(y\right)}^{\mathrm{T}}\left[\begin{array}{cc}0& i_d\\ i_p-i_d& 0\end{array}\right]\Phi \left(y\right)\mathrm{d}y$$

(34)

$$k_{\mathrm{s}2}=\underset{0}{\overset{l}{\int }}\rho A_s\mathsf{\psi }{\left(y\right)}^{\mathrm{T}}\left[\begin{array}{cc}0& 0\\ 0& -1\end{array}\right]\mathsf{\psi }\left(y\right)\mathrm{d}y+\left(i_p-i_d\right)\underset{0}{\overset{l}{\int }}\Phi {\left(y\right)}^{\mathrm{T}}\left[\begin{array}{cc}0& 0\\ 0& -1\end{array}\right]\Phi \left(y\right)\mathrm{d}y$$

(35)

$$k_{\mathrm{s}3}=\underset{0}{\overset{l}{\int }}\rho A_s\mathsf{\psi }{\left(y\right)}^{\mathrm{T}}\left[\begin{array}{cc}-1& 0\\ 0& -1\end{array}\right]\mathsf{\psi }\left(y\right)\mathrm{d}y+\left(i_p-i_d\right)\underset{0}{\overset{l}{\int }}\Phi {\left(y\right)}^{\mathrm{T}}\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\Phi \left(y\right)\mathrm{d}y$$

(36)

$$k_{\mathrm{s}4}=\underset{0}{\overset{l}{\int }}\rho A_s\mathsf{\psi }{\left(y\right)}^{\mathrm{T}}\left[\begin{array}{cc}-1& 0\\ 0& 0\end{array}\right]\mathsf{\psi }\left(y\right)\mathrm{d}y+\left(i_p-i_d\right)\underset{0}{\overset{l}{\int }}\Phi {\left(y\right)}^{\mathrm{T}}\left[\begin{array}{cc}-1& 0\\ 0& 0\end{array}\right]\Phi \left(y\right)\mathrm{d}y$$

(37)

$$k_{\mathrm{s}5}=\underset{0}{\overset{l}{\int }}\rho A_s\mathsf{\psi }{\left(y\right)}^{\mathrm{T}}\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\mathsf{\psi }\left(y\right)\mathrm{d}y+\left(i_p-i_d\right)\underset{0}{\overset{l}{\int }}\Phi {\left(y\right)}^{\mathrm{T}}\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\Phi \left(y\right)\mathrm{d}y$$

(38)

$$k_{\mathrm{s}6}=i_p\underset{0}{\overset{l}{\int }}\Phi {\left(y\right)}^{\mathrm{T}}\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\Phi \left(y\right)\mathrm{d}y$$

(39)

$$f_{\mathrm{str}}=-p_{\mathrm{s}1}{a}_{bx}-p_{\mathrm{s}2}{a}_{bz}$$

(40)

$$p_{\mathrm{s}1}=\rho A_s\underset{0}{\overset{l}{\int }}\psi {\left(y\right)}^{\mathrm{T}}\left[\begin{array}{c}1\\ 0\end{array}\right]\mathrm{d}y,p_{\mathrm{s}2}=\rho A_s\underset{0}{\overset{l}{\int }}\psi {\left(y\right)}^{\mathrm{T}}\left[\begin{array}{c}0\\ 1\end{array}\right]\mathrm{d}y$$

(41)

$$\begin{gathered} f_{\mathrm{slb}}\left(t\right)=p_{\mathrm{s}1}\left(x_{lb}\left({\omega }_{by}^2+{\omega }_{bz}^2\right)+y_{lb}\left({\dot{\omega }}_{bz}-{\omega }_{bx}{\omega }_{by}\right)-z_{lb}\left({\dot{\omega }}_{by}+{\omega }_{bx}{\omega }_{bz}\right)\right) \\ +p_{\mathrm{s}2}\left(x_{lb}\left({\dot{\omega }}_{by}-{\omega }_{bx}{\omega }_{bz}\right)-y_{lb}\left({\dot{\omega }}_{bx}+{\omega }_{by}{\omega }_{bz}\right)+z_{lb}\left({\omega }_{bx}^2+{\omega }_{by}^2\right)\right) \end{gathered}$$

(42)

$$f_{\mathrm{san}}\left(t\right)=s_{\mathrm{s}1}{\dot{\omega }}_{bx}+s_{\mathrm{s}2}{\dot{\omega }}_{bz}+s_{\mathrm{s}3}{\omega }_{bx}{\omega }_{by}+s_{\mathrm{s}4}{\omega }_{by}{\omega }_{bz}+s_{\mathrm{s}5}{\mathsf{\Omega }}_r{\omega }_{bx}+s_{\mathrm{s}6}{\mathsf{\Omega }}_r{\omega }_{bz}$$

(43)

$$s_{\mathrm{s}1}=\underset{0}{\overset{l}{\int }}\rho A_s\left(y_{s,j}+y\frac{y_{s,j+1}-y_{s,j}}{l}\right)\mathsf{\psi }{\left(y\right)}^{\mathrm{T}}\left[\begin{array}{c}0\\ -1\end{array}\right]\mathrm{d}y+i_d\underset{0}{\overset{l}{\int }}\Phi {\left(y\right)}^{\mathrm{T}}\left[\begin{array}{c}-1\\ 0\end{array}\right]\mathrm{d}y$$

(44)

$$s_{\mathrm{s}2}=\underset{0}{\overset{l}{\int }}\rho A_s\left(y_{s,j}+y\frac{y_{s,j+1}-y_{s,j}}{l}\right)\mathsf{\psi }{\left(y\right)}^{\mathrm{T}}\left[\begin{array}{c}1\\ 0\end{array}\right]\mathrm{d}y+i_d\underset{0}{\overset{l}{\int }}\Phi {\left(y\right)}^{\mathrm{T}}\left[\begin{array}{c}0\\ -1\end{array}\right]\mathrm{d}y$$

(45)

$$s_{\mathrm{s}3}=\underset{0}{\overset{l}{\int }}\rho A_s\left(y_{s,j}+y\frac{y_{s,j+1}-y_{s,j}}{l}\right)\mathsf{\psi }{\left(y\right)}^{\mathrm{T}}\left[\begin{array}{c}-1\\ 0\end{array}\right]\mathrm{d}y+\left(i_p-i_d\right)\underset{0}{\overset{l}{\int }}\Phi {\left(y\right)}^{\mathrm{T}}\left[\begin{array}{c}0\\ -1\end{array}\right]\mathrm{d}y$$

(46)

$$s_{\mathrm{s}4}=\underset{0}{\overset{l}{\int }}\rho A_s\left(y_{s,j}+y\frac{y_{s,j+1}-y_{s,j}}{l}\right)\mathsf{\psi }{\left(y\right)}^{\mathrm{T}}\left[\begin{array}{c}0\\ -1\end{array}\right]\mathrm{d}y+\left(i_p-i_d\right)\underset{0}{\overset{l}{\int }}\Phi {\left(y\right)}^{\mathrm{T}}\left[\begin{array}{c}1\\ 0\end{array}\right]\mathrm{d}y$$

(47)

$$s_{\mathrm{s}5}=i_p\underset{0}{\overset{l}{\int }}\Phi {\left(y\right)}^{\mathrm{T}}\left[\begin{array}{c}0\\ -1\end{array}\right]\mathrm{d}y,s_{\mathrm{s}6}=i_p\underset{0}{\overset{l}{\int }}\Phi {\left(y\right)}^{\mathrm{T}}\left[\begin{array}{c}1\\ 0\end{array}\right]\mathrm{d}y$$

(48)

$$f_{\mathrm{sg}}\left(t\right)=p_{\mathrm{s}1}g\cos \alpha \sin \beta -p_{\mathrm{s}2}g\cos \alpha \cos \beta$$

(49)

2.3. FEM Equations of the Dual-Rotor System

2.3.1. EOMs for Dual-Rotor System under General Base Motions

The number of beam elements for the two shafts is denoted as $n_e$, and the number of nodes is $n_e+2$. After assembling equations for the disks, mass unbalance, linear bearings and shaft elements, the finite element equations for the dual-rotor system with base motions are obtained:

$$M\ddot{q}+\left({\mathsf{\Omega }}_{LP}{G}_{0,\mathrm{LP}}+{\mathsf{\Omega }}_{HP}{G}_{0,\mathrm{HP}}+C_0+C_{\mathrm{b}}\right)\dot{q}+\left(K_0+K_{\mathrm{b}}\right)q=F_{\mathrm{btr}}+F_{\mathrm{blb}}+F_{\mathrm{ban}}+F_{\mathrm{u}}+F_{\mathrm{g}}$$

(50)

where $q$ is the $4\left(n_e+2\right)\times 1$ dimension displacement vector. $M={\sum }^{ }\left(m_{\mathrm{d}}+m_{\mathrm{s}}\right)$; $G_0={\sum }^{ }\left(g_{\mathrm{d}0}+g_{\mathrm{s}0}\right)$; $C_0={\sum }^{ }\left(\xi M+\eta k_{\mathrm{s}0}+c_{\mathrm{brg}}\right);$ and $K_0={\sum }^{ }\left(k_{\mathrm{s}0}+k_{\mathrm{brg}}\right)$. These are, respectively, the conventional mass, gyroscopic, damping and stiffness matrices, and ${\sum }^{ }$ represents the assembly of EOMs of the sub-systems. $C_{\mathrm{b}}$ and $K_{\mathrm{b}}$ are additional time-variable damping and stiffness matrices due to base angular motions. The subscripts ‘tr’, ‘lb’, ‘an’, ‘u’ and ‘g’ on the right of the equation represent the sources of the external forces; namely, base translations, left bearing, base angular motions, mass unbalance and gravity. Part of these matrices and vectors are given by

$$C_{\mathrm{b}}=C_{\mathrm{b}1}{\omega }_{by}$$

(51)

$$K_{\mathrm{b}}=K_{\mathrm{b}1}{\dot{\omega }}_{by}+K_{\mathrm{b}2}{\omega }_{bx}^2+K_{\mathrm{b}3}{\omega }_{by}^2+K_{\mathrm{b}4}{\omega }_{bz}^2+K_{\mathrm{b}5}{\omega }_{bx}{\omega }_{bz}+\left(K_{\mathrm{b}6,\mathrm{LP}}{\mathsf{\Omega }}_{LP}+K_{\mathrm{b}6,\mathrm{HP}}{\mathsf{\Omega }}_{HP}\right){\omega }_{by}$$

(52)

$$F_{\mathrm{btr}}\left(t\right)=-p_{\mathrm{b}1}{a}_{bx}-p_{\mathrm{b}2}{a}_{bz}$$

(53)

$$\begin{gathered} F_{\mathrm{blb}}\left(t\right)=p_{\mathrm{b}1}\left(x_{lb}\left({\omega }_{by}^2+{\omega }_{bz}^2\right)+y_{lb}\left({\dot{\omega }}_{bz}-{\omega }_{bx}{\omega }_{by}\right)-z_{lb}\left({\dot{\omega }}_{by}+{\omega }_{bx}{\omega }_{bz}\right)\right) \\ +p_{\mathrm{b}2}\left(x_{lb}\left({\dot{\omega }}_{by}-{\omega }_{bx}{\omega }_{bz}\right)-y_{lb}\left({\dot{\omega }}_{bx}+{\omega }_{by}{\omega }_{bz}\right)+z_{lb}\left({\omega }_{bx}^2+{\omega }_{by}^2\right)\right) \end{gathered}$$

(54)

$$F_{\mathrm{ban}}\left(t\right)=s_{\mathrm{b}1}{\dot{\omega }}_{bx}+s_{\mathrm{b}2}{\dot{\omega }}_{bz}+s_{\mathrm{b}3}{\omega }_{bx}{\omega }_{by}+s_{\mathrm{b}4}{\omega }_{by}{\omega }_{bz}+\left(s_{\mathrm{b}5,\mathrm{LP}}{\mathsf{\Omega }}_{LP}+s_{\mathrm{b}5,\mathrm{HP}}{\mathsf{\Omega }}_{HP}\right){\omega }_{bx}+\left(s_{\mathrm{b}6,\mathrm{LP}}{\mathsf{\Omega }}_{LP}+s_{\mathrm{b}6,\mathrm{HP}}{\mathsf{\Omega }}_{HP}\right){\omega }_{bz}$$

(55)

$$F_{\mathrm{g}}\left(t\right)=p_{\mathrm{b}1}g\cos \alpha \sin \beta -p_{\mathrm{b}2}g\cos \alpha \cos \beta$$

(56)

where $C_{\mathrm{b}1}={\sum }^{ }\left(c_{\mathrm{d}1}+c_{\mathrm{s}1}\right)$; $K_{\mathrm{bi}}={\sum }^{ }\left(k_{\mathrm{di}}+k_{\mathrm{si}}\right)$; $p_{\mathrm{bj}}={\sum }^{ }\left(p_{\mathrm{dj}}+p_{\mathrm{sj}}\right)$; $s_{\mathrm{bi}}={\sum }^{ }\left(s_{\mathrm{di}}+s_{\mathrm{si}}\right)$; and ($\mathrm{i}=1,2,\cdots 6$, $\mathrm{j}=1,2$). Those vectors and matrices––including $C_{\mathrm{b}1}$, $K_{\mathrm{bi}}$, $p_{\mathrm{bj}}$ and $s_{\mathrm{bi}}$––are determined by geometric and material parameters and kept constant for a dual-rotor system.

2.3.2. EOMs for the Dual-Rotor System under Specific Base Motions

For the counter-rotating dual-rotor system under harmonic base translation along the $O_0{X}_0$ axis and base rotation around the $O_b{X}_b$ axis with constant angular velocity, the additional damping and stiffing matrices and the force vectors in Equation (50) are rewritten as

$$C_{\mathrm{b}}=\mathbf{0}$$

(57)

$$K_{\mathrm{b}}=K_{\mathrm{b}2}{\omega }_{bx}^2$$

(58)

$$F_{\mathrm{btr}}=p_{\mathrm{b}1}{x}_{b1}{\omega }_{xb}^2\cos \left({\omega }_{xb}t\right)$$

(59)

$$F_{\mathrm{blb}}\left(t\right)=\mathbf{0}$$

(60)

$$F_{\mathrm{ban}}=\left(s_{\mathrm{b}5,\mathrm{LP}}{\mathsf{\Omega }}_{LP}+s_{\mathrm{b}5,\mathrm{HP}}{\mathsf{\Omega }}_{HP}\right){\omega }_{bx}$$

(61)

$$F_{\mathrm{g}}=-p_{\mathrm{b}2}g\cos \left({\omega }_{bx}t\right)$$

(62)

$$F_{\mathrm{du}}=m_{u}e\left[\begin{array}{c}{\mathsf{\Omega }}_r^2\cos \left({\mathsf{\Omega }}_{r}t+{\phi }_u\right)\\ -\left({\mathsf{\Omega }}_r^2+{\omega }_{bx}^2\right)\sin \left({\mathsf{\Omega }}_{r}t+{\phi }_u\right)\\ 0\\ 0\end{array}\right]$$

(63)

For a counter-rotating dual-rotor system under harmonic base translation along the $O_0{X}_0$ axis and base rotation around the $O_b{Y}_b$ axis with constant angular velocity, the additional damping and stiffness matrices and the force vectors in Equation (50) are rewritten as

$$C_{\mathrm{b}}=C_{\mathrm{b}1}{\omega }_{by}$$

(64)

$$K_{\mathrm{b}}=K_{\mathrm{b}3}{\omega }_{by}^2+\left(K_{\mathrm{b}6,\mathrm{LP}}{\mathsf{\Omega }}_{LP}+K_{\mathrm{b}6,\mathrm{HP}}{\mathsf{\Omega }}_{HP}\right){\omega }_{by}$$

(65)

$$\begin{gathered} F_{\mathrm{btr}}\left(t\right)=p_{\mathrm{b}1}\frac{1}{2}{x}_{b1}{\omega }_{xb}^2\left(\cos \left(\left({\omega }_{xb}+{\omega }_{by}\right)t\right)+\cos \left(\left({\omega }_{xb}-{\omega }_{by}\right)t\right)\right) \\ +p_{\mathrm{b}2}\frac{1}{2}{x}_{b1}{\omega }_{xb}^2\left(\sin \left(\left({\omega }_{xb}+{\omega }_{by}\right)t\right)-\sin \left(\left({\omega }_{xb}-{\omega }_{by}\right)t\right)\right) \end{gathered}$$

(66)

$$F_{\mathrm{blb}}\left(t\right)=\mathbf{0}$$

(67)

$$F_{\mathrm{ban}}\left(t\right)=\mathbf{0}$$

(68)

$$F_{\mathrm{g}}\left(t\right)=p_{\mathrm{b}1}g\sin \left({\omega }_{yb}t\right)-p_{\mathrm{b}2}g\cos \left({\omega }_{yb}t\right)$$

(69)

$$F_{\mathrm{du}}=m_{u}e\left[\begin{array}{c}{\left({\mathsf{\Omega }}_r+{\omega }_{by}\right)}^2\cos \left({\mathsf{\Omega }}_{r}t+{\phi }_u\right)\\ -{\left({\mathsf{\Omega }}_r+{\omega }_{by}\right)}^2\sin \left({\mathsf{\Omega }}_{r}t+{\phi }_u\right)\\ 0\\ 0\end{array}\right]$$

(70)

3. Results and Discussions

Compared with a co-rotating dual-rotor system, a counter-rotating system possess more modern industrial applications and has more complex dynamic properties. Thus, a counter-rotating dual-rotor system was selected as the research target for this paper. The rotating direction of the low-pressure rotor is the same as the positive direction of the $O_b{Y}_b$ axis according to the right-hand screw rule, whereas that of the high-pressure rotor is the opposite. It is defined as forward whirling when the direction of the whirling is the same as the rotating direction of the low-pressure rotor, and backward whirling otherwise. The finite element division of the dual-rotor-bearing system is shown in Figure 3. The physical parameters for the shaft segments, disks and bearings are summarized in

Table 1. Shaft parameters for the dual-rotor-bearing system.

Rotors	Node–Node	$\mathit{l} \left(\mathbf{mm}\right)$	${\mathit{r}}_{\mathit{s}\mathit{i} }\left(\mathbf{mm}\right)$	${\mathit{r}}_{\mathit{s}\mathit{o}} \left(\mathbf{mm}\right)$
LP	1–2, 2–3	75	0	6
LP	3–4	180	0	6
LP	4–5, 5–6	120	0	6
LP	6–7, 7–8	50	0	6
HP	9–10, 10–11	120	19.9	27.5

,

Table 2. Disk parameters for the dual-rotor-bearing system.

Rotors	Node	${\mathit{m}}_{\mathit{d}} \left(\mathbf{kg}\right)$	${\mathit{I}}_{\mathit{d}}$$(\mathbf{k}\mathbf{g}\cdot {\mathbf{m}}^2)$	${\mathit{I}}_{\mathit{p}}$$(\mathbf{k}\mathbf{g}\cdot {\mathbf{m}}^2)$
LP	2, 7	0.7672	5.0924 × 10−4	9.9586 × 10−4
HP	10	0.1942	2.2109 × 10−4	4.3895 × 10−4

,

Table 3. Bearing parameters for the dual-rotor-bearing system.

Number	Node	${\mathit{k}}_{\mathit{b}}$$(\mathbf{k}\mathbf{g}\cdot {\mathbf{m}}^3)$	${\mathit{c}}_{\mathit{b}}$$(\mathbf{N}\cdot {\mathbf{m}}^{-1})$
1	1	5 × 107	100
2	2	5 × 107	100
3	6, 9	5 × 107	100
4	8	5 × 107	100
5	11	5 × 107	100

and

Table 4. Material parameters of the shaft.

Rotors	$\mathbf{Density}\mathit{\rho }$$(\mathbf{k}\mathbf{g}\cdot {\mathbf{m}}^{-3})$	$\mathbf{Modulus}\mathbf{of}\mathbf{Elasticity}\mathit{E} \left(\mathbf{Pa}\right)$	$\mathbf{Poisson}\text{’}\mathbf{s}\mathbf{Ratio}\mathit{\mu }$
LP	7800	2.06 × 1011	0.3
HP	2800	7 × 1010	0.34

.

3.1. Computation Method and Validation

From the results in Section 2.3.2, it is apparent that the mass, gyroscopic, damping and stiffing matrices on the left of the EOMs are constant and the forces on the right of the EOMs are constant or harmonic. Thus, the natural properties and responses of the linear dual-rotor system can be solved by direct matrix operations in the complex field.

Introducing a state-space vector $h=\left[\begin{array}{c}q\\ \dot{q}\end{array}\right]$ and assuming that $q=\mathsf{\varphi }{e}^{\lambda t}$, the second-order homogeneous EOMs for the dual-rotor system with constant base rotations can be transformed into first-order state space form

$$\left[\begin{array}{cc}{K}_0+K_{\mathrm{b}}& \mathbf{0}\\ \mathbf{0}& M\end{array}\right]\dot{h}+\left[\begin{array}{cc}\mathbf{0}& -\left(K_0+K_{\mathrm{b}}\right)\\ K_0+K_{\mathrm{b}}& {\mathsf{\Omega }}_{LP}{G}_{0,\mathrm{LP}}+{\mathsf{\Omega }}_{HP}{G}_{0,\mathrm{HP}}+C_{\mathrm{b}}\end{array}\right]h=\mathbf{0}$$

(71)

The equivalent eigenvalue problem is

$$-{\left[\begin{array}{cc}{K}_0+K_{\mathrm{b}}& \mathbf{0}\\ \mathbf{0}& M\end{array}\right]}^{\mathrm{T}}\left[\begin{array}{cc}\mathbf{0}& -\left(K_0+K_{\mathrm{b}}\right)\\ K_0+K_{\mathrm{b}}& {\mathsf{\Omega }}_{LP}{G}_{0,\mathrm{LP}}+{\mathsf{\Omega }}_{HP}{G}_{0,\mathrm{HP}}+C_{\mathrm{b}}\end{array}\right]\left[\begin{array}{c}\mathsf{\varphi }\\ \lambda \mathsf{\varphi }\end{array}\right]=\lambda \left[\begin{array}{c}\mathsf{\varphi }\\ \lambda \mathsf{\varphi }\end{array}\right]$$

(72)

The solved $\lambda$ and $\mathsf{\varphi }$ are the natural frequency and corresponding mode shape vector. For the constant forces assumed as $F_0$, the responses are solved by

$$q_0={\left[K_0+K_{\mathrm{b}}\right]}^T{F}_0$$

(73)

The harmonic forces assumed as $F_{\mathrm{v}}=V_{\mathrm{c}}\cos \left(\omega t\right)+V_{\mathrm{s}}\sin \left(\omega t\right)$ can be rewritten in a complex field as

$$F_{\mathrm{v}}=\left(V_{\mathrm{c}}-j{V}_{\mathrm{s}}\right)e^{j\omega t}$$

(74)

where $j$ is the imaginary unit. The steady-state solution is assumed to be

$$q_{\mathrm{v}}={\mathsf{\varphi }}_{\mathrm{v}}{e}^{j\omega t}$$

(75)

The harmonic responses are solved as

$${\mathsf{\varphi }}_{\mathrm{v}}={\left[-M{\omega }^2+j\omega \left({\mathsf{\Omega }}_{LP}{G}_{0,\mathrm{LP}}+{\mathsf{\Omega }}_{HP}{G}_{0,\mathrm{HP}}+C_0+C_{\mathrm{b}}\right)+K_0+K_{\mathrm{b}}\right]}^T\left(V_{\mathrm{c}}-j{V}_{\mathrm{s}}\right)$$

(76)

The results are the superposition of the constant and various harmonic responses.

The accuracy of the above computation method is validated by comparison of the results solved by the Runge–Kutta numerical integration method. The base motions are set as $x_b=x_{b1}\cos {\omega }_{xb}t$, $y_b=0$, $z_b=0$, $\alpha =\gamma =0$ and $\beta ={\omega }_{by}t$, where $x_{b1}=0.001 \mathrm{m}$, ${\omega }_{xb}=25 \mathrm{Hz}$, ${\omega }_{by}=10 \mathrm{Hz}$. The rotation speeds of the dual-rotor system are set as ${\mathsf{\Omega }}_{LP}=100\pi \mathrm{rad}/\mathrm{s}$, ${\mathsf{\Omega }}_{HP}=-1.2{\mathsf{\Omega }}_{LP}$. The mass unbalances located on the three disks are set as $m_{u1}{e}_1=m_{u2}{e}_2=m_{u3}{e}_3=1\times {10}^{-4}\mathrm{kg}\cdot \mathrm{m}$. The results of disk 1 in both time and spectra domain are shown in Figure 4. By comparison, one can see that the results of two methods are in very good agreement with each other, verifying the accuracy of solving by direct matrix operation. The natural properties and harmonic responses of the counter-rotating dual-rotor system under constant base rotations are analyzed in following sections using direct computation method.

3.2. Natural Properties Influenced by Base Rotations

The influence of base rotations on the natural properties of a dual-rotor system are analyzed in this section. For a counter-rotating system, the rotating speed relationship between the two rotors was set at ${\mathsf{\Omega }}_{HP}=-1.2{\mathsf{\Omega }}_{LP}$.

The Campbell diagrams with variable base rotations are summarized in Figure 5, where (a) and (c) represent the pitching base motion (${\omega }_{bx}\ne 0$), and (b) and (d) are for the rolling base motion (${\omega }_{by}\ne 0$). As for the pitching base motion (${\omega }_{bx}\ne 0$), the overall effect was not strong. The limited differences showed that both the backward and forward whirling frequencies decreased with increasing base angular velocity (see zoom on Ⅰ and Ⅱ in (a) and (c)). Thus, the critical rotating speeds excited by unbalance decreased a little (see point A, B, C, D, E and F in (c)). Compared to the pitching base motion, the effect of rolling base motion was more noticeable. With increasing base angular velocity, the backward whirling frequencies increased while forward whirling frequencies decreased. In consequence, the critical speeds excited by the low-pressure rotor unbalance decreased, while the critical speeds excited by high-pressure rotor unbalance increased.

3.3. Forced Vibration Due to Mass Unbalance with Base Rotations

The influences of the constant base rotations on the mass unbalance responses of the dual-rotor system are analyzed in this section. As mentioned above, the rotation speed relationship between the rotors is set at ${\mathsf{\Omega }}_{HP}=-1.2{\mathsf{\Omega }}_{LP}$.

For the mass unbalance on a low-pressure rotor, the mass unbalances on the three disks are set as $m_{u1}{e}_1=m_{u3}{e}_3=0$ and $m_{u2}{e}_2=1\times {10}^{-4}\mathrm{kg}\cdot \mathrm{m}$. The frequency responses are shown in Figure 6. Only the forward whirling modes are excited when the base is fixed (see the blue lines). For the pitching base motion shown in Figure 6a,b, the backward whirling modes were also excited when the pitching base motion is concerned. The resonant frequencies decreased a little with increasing angular velocity of the pitching base motion. For the rolling base motion shown in Figure 6c,d, only the forward whirling modes were excited. The resonant frequencies decreased remarkably, when the angular velocity of the rolling base motion increased.

For the mass unbalance on a high-pressure rotor, the mass unbalances are set at $m_{u1}{e}_1=m_{u2}{e}_2=0$ and $m_{u3}{e}_3=1\times {10}^{-4}\mathrm{kg}\cdot \mathrm{m}$. The frequency responses are shown in Figure 7. Only the backward whirling modes were excited when the base was fixed (see the blue lines). For the pitching base motion shown in Figure 7a,b, the forward whirling modes were also excited when the pitching base motion was concerned. The resonant frequencies decreased a little with increasing angular velocity of the pitching base motion. For the rolling base motion shown in Figure 7c,d, only the backward whirling modes were excited. The resonant frequencies decreased remarkably, when the angular velocity of the rolling base motion increased.

Here are some explanations for the above phenomenon. When the base was fixed, the system possessed perfect symmetry along the two linear displacement directions ($O_b{X}_b$ and $O_b{Z}_b$). Even the amplitudes of the two components of the unbalance force along these directions were equal. Thus, the unbalanced force was only in accordance with whirling modes along one direction: backward for the mass unbalance on the low-pressure rotor and forward for the mass unbalance on the high-pressure rotor. The rolling base motion did not break this symmetry, and the amplitudes of the two components of the unbalance force remained equal (see Equations (32), (36), (39), (64), (65) and (70)). However, when the pitching base motion was involved, the stiffness symmetry along $O_b{X}_b$ and $O_b{Z}_b$ was not maintained and the amplitudes of the two components of the unbalance force were not equal (see Equations (35), (58) and (63)).

3.4. Forced Vibration Due to Base Translations with Base Rotations

The influences of the base translations on the dual-rotor system are analyzed in this section. The counter-rotating dual-rotor system was taken as analysis target, and the rotating speeds of the LP and HP rotor were set at ${\mathsf{\Omega }}_{LP}=1000 \mathrm{rad}/\mathrm{s}$, ${\mathsf{\Omega }}_{HP}=-1.2{\mathsf{\Omega }}_{LP}$. The natural frequencies were calculated first and listed in

Table 5. Natural frequencies of the counter-rotating dual-rotor system (rad/s).

Orders	${\mathit{\omega }}_{\mathit{b}1}$	${\mathit{\omega }}_{\mathit{f}1}$	${\mathit{\omega }}_{\mathit{b}2}$	${\mathit{\omega }}_{\mathit{f}2}$	${\mathit{\omega }}_{\mathit{b}3}$	${\mathit{\omega }}_{\mathit{f}3}$
Without base rotations	461.3	542.4	1530.7	1554.5	2223.5	2230.8
With pitching base motion	459.3	540.4	1530.0	1553.9	2223.0	2230.4
With rolling base motion	483.3	520.6	1492.7	1592.5	2168.2	2286.0

.

The base translations were assumed to be time-periodic: $x_b=x_{b1}\cos {\omega }_{xb}t$, $y_b=y_{b1}\cos {\omega }_{yb}t$ and $z_b=z_{b1}\cos {\omega }_{zb}t$. The initial base angular angle was set to zero: $\alpha =\beta =\gamma =0$. For the base translation along the $O_0{X}_0\mathrm{axis}$, the base translation amplitude was set at $x_{b1}=0.001 \mathrm{m}$. The harmonic responses of the disk 1 with and without base rotations are shown in Figure 8, Figure 9 and Figure 10. The results revealed that not only the forward whirling but also the backward whirling resonances could be excited by harmonic base translations, while some differences still merit attention. For the base translations with pitching base motion, the response amplitudes changed remarkably; however, the resonant frequencies had little distinction compared with the results without base rotations because the effects of the pitching base motion were applied through changes in the stiffness matrices and unbalance force amplitude (see Equations (58) and (63)). For the base translations with rolling base motion, the unbalance forces consisted of two frequencies: combined translation and rolling motion (see Equation (66)). Backward whirling modes could only be excited by part of the frequency $\left|{\omega }_{xb}-{\omega }_{by}\right|$, while forward whirling modes could only be excited by another part $\left|{\omega }_{xb}+{\omega }_{by}\right|$.

Here are some explanations for the above phenomenon. For the case of fixed base and pitching base motion, all components of additional forces $F_{\mathrm{btr}}\left(t\right)$ changed as the cosine function of time $t$ (see Equation (59)). Thus, the excitation vector had components in accordance with all the whirling modes. However, in the case of rolling base motion, the two parts of the excitation vector only had components in accordance with backward whirling or forward whirling modes, respectively (see Equation (66)).

3.5. Forced Vibration Due to Gravity with Base Rotations

The influence of gravity on the dual-rotor system are analyzed in this section. For the counter-rotating system, the rotating speeds of the LP and HP rotors were set at ${\mathsf{\Omega }}_{LP}=1000 \mathrm{rad}/\mathrm{s}$ and ${\mathsf{\Omega }}_{HP}=-1.2{\mathsf{\Omega }}_{LP}$.

Gravity was constant and caused only constant displacements along its acting direction when the system was under a fixed base. It became a harmonic force with base angular frequency when the constant base rotations were concerned (see Equations (62) and (69)). The harmonic responses of the dual-rotor system due to gravity and constant base rotations are shown in Figure 11 and Figure 12. Both backward and forward whirling modes could be excited by gravity for the dual-rotor system with pitching base motion, but no resonance was excited with rolling base motion. The whirling frequencies changed with the value of base angular velocities as the additional stiffness and damping matrices were concerned with base angular velocities (see Equation (58), (64) and (65)).

4. Conclusions

The combined effects of base motions, mass unbalance and gravity on vibration properties of a counter-rotating dual-rotor system were studied in this paper for the first time. Some conclusions are summarized as follows:

(1)

Constant base rotations generated constant external stiffness and damping matrices, and so had some impact on its natural properties. The effects of rolling base motion are larger than that of pitching base motion.

(2)

When the base was fixed, mass unbalance on different rotors only excited whirling modes that were in the same direction with the rotor rotating. The constant pitching base motion caused other whirling modes to be excited by mass unbalance forces, while the constant rolling base motion did not.

(3)

Base translations can stimulate both backward and forward whirling modes of the counter-rotating dual-rotor system whether or not the constant pitching base motion was involved. It is interesting that when the constant rolling base motion and base translation along the horizontal direction were both taken into account, the combined frequencies of the base angular velocity and base translation frequency can induce different resonances.

(4)

Components along different directions of the gravity exhibited harmonic behavior when the constant base rotations were concerned. They can also excite different whirling modes, while one with rolling base motion did not. In addition, matrix transformation in a complex field was validated as a useful and efficient computing method.