The Influence of Coordinate Systems on the Stability Analysis of Lateral–Torsional Coupled Vibration

Abstract

Stability analysis of lateral–torsional coupled vibration is obligatory for rotating machinery, such as aero-engines. However, the state-of-the-art method may lead to stability misjudgment under different coordinate systems. The cause of this misjudgment has not yet been well explored. The purpose of this paper is to clarify the error source of the stability analysis in a more comprehensive manner. A vertical Jeffcott rotor model including torsion vibration is built, and the Lagrange approach is applied to establish the motion equations. The coordinate transformation matrix is used to transfer the motion equations into the rotating coordinate system, making the coefficients of the motion equation constants. The differences in the unstable speed regions in the two coordinate systems are captured. The limitations of the Floquet theory and Hill’s determinant analysis in the stability estimation of the lateral–torsional coupled vibration are explained. It is found that, for Hill’s method, increasing the number of the harmonic truncation cannot correct the misjudgment, and the matrix truncation is the fundamental error source. The above research provides more accurate theoretical support for the analysis of the lateral–torsional coupling instability of rotors.

1. Introduction

Present day operating gaps between rotating and stationary components are kept to a minimum to improve the aerodynamic performance of aero-engines. Moreover, the unbalance and eccentricity may cause undesirable lateral vibration, which are produced by manufacturing errors, the fan blade being off, etc. [1,2,3]. The consequence would be a rotating element coming into contact with the stationary part, namely a rotor stator rub [4].

The rub can significantly change the dynamic behavior of rotor systems. There has been extensive research on the rotor–stator contact phenomenon [5,6]. As a result, the rotor exhibits non-smooth nonlinearities and a complicated vibration motion, such as periodic, quasi-periodic [7,8], and chaotic vibrations with sub-synchronous and super-synchronous frequency components [9,10,11]. Chu et al. [12] and Sun et al. [13] used a nonlinear model with piecewise linear stiffness and demonstrated that the rub impact between the rotor and the stator exhibited periodic, quasi-periodic, and chaotic vibrations. Qin et al. [14] observed grazing bifurcation and chaos in the response of a rubbing rotor. Sinha [15] presented an analytical method to illustrate the nonlinear dynamic effect of blades rubbing against the rigid outer case in rotating machinery. Using the complex nonlinear modes, Hong et al. [16] analyzed the nonlinear dynamic characteristics of a rotor system with an additional constraint due to the rub impact.

Researchers have also found that the torsional vibration caused by the rub could not be neglected [17]. The dynamic characteristics of the rotor system are seriously affected by the coupling of the vibrations in lateral and torsional directions under conditions such as gear meshing [18], unbalanced excitation [19,20], and cracking [21]. The subsequent lateral–torsional coupled vibration may lead to response amplification and even instability, resulting in the catastrophic failure of machines. Al-Bedoor [22] carried out a transient analysis with torsional motion. He observed a super-harmonic of order 3, with possible instability, in the case of torsional vibrations. Huang [23] studied the characteristics of the torsional vibration of a rotor that was unbalanced, and found the torsional vibration of the shaft could result in lateral vibration with a bisynchronous frequency. From the modal and response perspectives, Hong et al. [24] discussed the lateral–torsional coupling effects in detail. The results showed that veering and lock-in phenomena occurred between the frequencies of the forward whirl mode and the torsional mode, and lock-in could lead to a kind of principal instability.

Different modeling and simulation methods have been developed to meet the various analysis requirements of lateral–torsional coupled vibration. In terms of modeling, the lumped parameter model [25] preserves the dynamic characteristics of the target modes of the system with just a few degrees of freedom (DOFs), and they are widely used to qualitatively reveal the instability and chaos of the rotor caused by bending–torsional coupling or nonlinear friction–impact force, with high computational efficiency. Shi et al. [26] proposed a model with six DOFs to investigate the vibrational stability of the unbalanced rotor bending–torsional coupling vibration of a hydro-turbine generator unit. Recently, researchers have gradually adopted the finite element model [27] due to the improvement in computing ability. Although the computational burden significantly increases, it can quantitatively predict the dynamic characteristics of the rotor in engineering. Zheng et al. [28] built an FE model of a rub-impact rotor during hovering flight; the disk was five DOFs. To study coupled lateral–torsional dynamic problems, researchers have employed different numerical integration techniques to investigate this nonlinear phenomenon. Mokhtar et al. [29] adopted a forward increment Lagrange multiplier method with a central difference integration scheme to calculate the transient response of the rotor-stator contact. Cao et al. [30] used a fourth-order Runge–Kutta method to solve problems, including the lateral and torsional analysis of a flexible rotor-bearing system, which comprised gyroscopic effects, nonlinear short journal bearings, nonlinear short squeeze film dampers, and external nonlinear forces/torques. Transient analysis with the typical time integration method can reveal the change process of the response, especially the instability amplification, but it needs more computational efforts. Li et al. [31] proposed a general algorithm named the linear and nonlinear nodes separation method to improve the efficiency of the transient numerical simulation of the coupled lateral–torsional dynamic behavior. Additionally, experimental investigations were performed to validate the aforementioned numerical methods in the analysis of the lateral–torsional coupled vibration [32,33].

An appropriate coordinate system should be chosen for the dynamic analysis of rotor systems. Most studies adopt the inertial coordinate system, which is fixed to the ground. The inertial frame is suitable for axisymmetric and non-axisymmetric rotors with strong applicability [34]. However, the time-varying governing motion equations in the inertia frame are complex [35]. The Floquet–Hill formulation has become an essential and strong tool for computing these equations and performing stability analysis of lateral–torsional coupled vibrations [36]. Because of the small harmonic truncation order in the method, the eigenvalues are redundant, which could cause stability misjudgment. Currently, the misjudgment can be partially corrected by filtering out the spectrum of the solutions based on the eigenvalues [37] and eigenvectors [38,39] of Hill’s matrix. Nonetheless, in some cases, the redundant roots have been estimated as correct ones. The effectiveness of the filtering strategy might be related to the case. This leads to a lack of an extensive understanding of the misjudgment, which regrettably still occurs in many rotor stability studies [24,35,40]. On the other hand, some researchers have started to use the rotating coordinate system; the frame rotates at the same speed as the rotor system [41]. The time-varying coefficients of the motion equations can be converted to constants by the rotational transformation matrix, and then the difficulty of analysis is significantly reduced. The above research shows the advantages of rotor stability analyses in the rotating frame.

Although the motion equations of the rotor system have different forms in the two coordinate systems, the results should theoretically be consistent. However, discrepancies may arise resulting in differences in the stability judgment for some motions [42]. The inconsistency of the stability results between the two frames which has not been fully revealed deserves further investigation. Therefore, in this paper, we study the influence of the choice of the coordinate system on the analysis of lateral–torsional coupled vibration, and discuss the reasons for the differences between the corresponding stability results. A three DOFs vertical Jeffcott rotor system including torsion vibration is built (Section 2). The analysis of lateral–torsional coupled vibration in both inertial and rotating coordinate systems is performed (Section 3). The influence of the coordinate system choice is revealed by comparing the results of time-domain response and modal analysis in both coordinate systems, and the cause of misjudgment is further explained from the perspective of the truncation of motion equations (Section 4).

2. Model Formulation

A vertical Jeffcott rotor system considering torsional vibration is shown in Figure 1. The rotor has a mass $m$ imbalance with eccentricity $e$. The rotor’s mass moment of inertia around the $z$-axis is denoted as $J_{\mathrm{p}}$. Assuming rotor symmetry, the lateral stiffness values are set to $k_{\mathrm{x}}=k_{\mathrm{y}}=k$, and the equivalent damping coefficients are defined as $c_{\mathrm{x}}=c_{\mathrm{y}}=c$. The torsional stiffness and the damping coefficient are represented as $k_{\mathrm{r}}$ and $c_{\mathrm{r}}$, respectively. The rotation speed of the rotor is denoted as $\Omega$ and the torsional displacement as $\theta$. The arrow in Figure 1 represents the rotational direction of the rotor.

The 3-DOFs rotor model in the inertial coordinate includes two lateral motions ($x$, $y$) and one torsional motion ($\theta$). By assuming small displacements around the static equilibrium position and constant rotational speed, a periodically time-varying system written in the inertial coordinate system is obtained (detailed in Appendix A) [24]:

$$\left\{\begin{array}{l}m\ddot{x}-me\ddot{\theta }\sin \Omega t+c\dot{x}-2me\Omega \dot{\theta }\cos \Omega t+kx+me{\Omega }^2\sin \Omega t=me{\Omega }^2\cos \Omega t\\ m\ddot{y}+me\ddot{\theta }\cos \Omega t+c\dot{y}-2me\Omega \dot{\theta }\sin \Omega t+ky-me{\Omega }^2\cos \Omega t=me{\Omega }^2\sin \Omega t\\ (m{e}^2+J_{\mathrm{p}})\ddot{\theta }-me\ddot{x}\sin \Omega t+me\ddot{y}\cos \Omega t+c_{\mathrm{r}}\dot{\theta }+k_{\mathrm{r}}\theta =0\end{array}\right.$$

(1)

In Equation (1), the first and second equations are force-balance equations separately in the $x$-direction and $y$-direction, and the third equation is a torque-balance equation. Equation (1) shows that the torsional vibration couples with the lateral vibration by the inertia terms of equations, namely $-me\ddot{\theta }\sin \Omega t$ and $me\ddot{\theta }\cos \Omega t$ in the two lateral motion equations and ($-me\ddot{x}\sin \Omega t+me\ddot{y}\cos \Omega t$) in the torsional motion equation, which are caused by the mass of the unbalance $e$. Furthermore, the presence of time-varying coefficients $\sin \Omega t$ and $\cos \Omega t$, caused by the coupling between lateral and torsional vibrations, makes the calculation of the motion equations more complex. System parameters considered in the present study are listed in

Table 1. Rotor’s physical and simulation parameters [24].

Parameter	Variable	Value	Unit
Mass	$m$	1	$\mathrm{kg}$
Lateral damping coefficient	$c$	0	$\mathrm{N}·\mathrm{m}/\mathrm{s}$
Lateral stiffness	$k$	104	$\mathrm{N}/\mathrm{s}$
Mass moment of inertial	$J_{\mathrm{p}}$	0.5	$\mathrm{kg}·{\mathrm{m}}^2$
Torsional damping coefficient	$c_{\mathrm{r}}$	0	$\mathrm{N}·\mathrm{m}·\mathrm{s}/\mathrm{rad}$
Torsional stiffness	$k_{\mathrm{r}}$	1250	$\mathrm{kg}·{\mathrm{m}/\mathrm{s}}^2$
Eccentricity	$e$	100	$\mathrm{mm}$
Rotation speed	$\Omega$	varying	$\mathrm{rpm}$

.

3. Analysis in Two Coordinate Systems

In this section, analytical methods for rotor stability analysis in both inertial and rotating coordinate systems are derived, and modal analysis is carried out to predict unstable regions of the rotor. Subsequently, a comparison is made between the results obtained in the two coordinate systems.

3.1. Analysis in Inertial Coordinate System

3.1.1. Analytical Method

In Section 2, the governing equations of the unbalanced rotor system under the inertial frame are obtained by Lagrange’s approach. The coefficients $\sin \Omega t$ and $\cos \Omega t$ with a period $T=2\pi /\Omega$ are periodic time-varying, leading to solutions to the motion equations periodic with the same period $T$. One can solve the aforementioned equations by the well-known Floquet theory [43] and Hill’s method [44].

Let us transform the periodic coefficients $\sin \Omega t$ and $\cos \Omega t$ into exponential form. Consider Euler’s formula:

$$\cos \Omega t=\frac{e^{\mathrm{i}\Omega t}+e^{-\mathrm{i}\Omega t}}{2},\sin \Omega t=\frac{-\mathrm{i}{e}^{\mathrm{i}\Omega t}+\mathrm{i}{e}^{-\mathrm{i}\Omega t}}{2}$$

(2)

where $\mathrm{i}=\sqrt{-1}$. Then, replace Equation (2) in Equation (1) and write it as matrix form. To analyze the modal characteristics of the unbalanced rotor system, set the external force vector equal to zero, and we obtain the following equation:

$$(M+M_t{e}^{\mathrm{i}\Omega t}+{\tilde{M}}_t{e}^{-\mathrm{i}\Omega t})\ddot{q}+(C+C_t{e}^{\mathrm{i}\Omega t}+{\tilde{C}}_t{e}^{-\mathrm{i}\Omega t})\dot{q}+(K+K_t{e}^{\mathrm{i}\Omega t}+{\tilde{K}}_t{e}^{-\mathrm{i}\Omega t})q=0$$

(3)

The matrices in Equation (3) are included in Appendix B.

The time-varying coefficients $\sin \Omega t$ and $\cos \Omega t$ are now substituted by exponential terms $e^{-\mathrm{i}\Omega t}$ and $e^{\mathrm{i}\Omega t}$. The solution of Equation (3) is periodic with the same period $T=2\pi /\Omega$. Thus, the solution of Equation (3) can be assumed in the following form according to Floquet’s theory:

$$\left\{\begin{array}{l}q=\mathit{\varphi }{e}^{\lambda t}\\ \mathit{\varphi }(t)=\mathit{\varphi }(t+T_{\mathrm{k}}),\hspace{0.33em}{T}_{\mathrm{k}}=2\pi /\Omega \end{array}\right.$$

(4)

where $\varphi$ is the quasi-mode vector, which is periodic with the same period $T$, and $\lambda$ is the quasi-eigenvalue [45]. The period function $\varphi$ can be expressed by the general Fourier series:

$$\mathit{\varphi }={\displaystyle \sum _{j=-\infty }^{+\infty }{I}_j{e}^{\mathrm{i}j\Omega t}}$$

(5)

where $I_j$ is the harmonic contribution of the $j$th order harmonic. By substituting Equations (4) and (5) into Equation (3), we obtain the following equation:

$$\begin{array}{l}(M+M_t{e}^{\mathrm{i}\Omega t}+{\tilde{M}}_t{e}^{-\mathrm{i}\Omega t}){\displaystyle \sum _{j=-\infty }^{+\infty }(-j{\Omega }^2+2\mathrm{i}\lambda j\Omega +{\lambda }^2)I_j{e}^{\mathrm{i}j\Omega t}}+\\ (C+C_t{e}^{\mathrm{i}\Omega t}+{\tilde{C}}_t{e}^{-\mathrm{i}\Omega t}){\displaystyle \sum _{j=-\infty }^{+\infty }(\mathrm{i}j\Omega +\lambda )I_j{e}^{\mathrm{i}j\Omega t}}+\\ (K+K_t{e}^{\mathrm{i}\Omega t}+{\tilde{K}}_t{e}^{-\mathrm{i}\Omega t}){\displaystyle \sum _{j=-\infty }^{+\infty }{I}_j{e}^{\mathrm{i}j\Omega t}}=0\end{array}$$

(6)

By separately setting the coefficients of each exponential term $e^{\mathrm{i}j\Omega t}$ to zero in Equation (6), an infinite set of algebraic equations are then obtained:

$$\left\{\begin{array}{c}⋮\\ \begin{array}{l}\{[{(\mathrm{i}(j-2)\Omega +\lambda )}^2{M}_t+(\mathrm{i}(j-2)\Omega +\lambda )C_t+K_t]I_{j-2}\\ +[{(\mathrm{i}(j-1)\Omega +\lambda )}^{2}M+(\mathrm{i}(j-1)\Omega +\lambda )C+K]I_{j-1}\\ +[{(\mathrm{i}j\Omega +\lambda )}^2{\tilde{M}}_t+(\mathrm{i}j\Omega +\lambda ){\tilde{C}}_t+{\tilde{K}}_t]I_j\}{e}^{\mathrm{i}(j-1)\Omega t}=0\\ \{[{(\mathrm{i}(j-1)\Omega +\lambda )}^2{M}_t+(\mathrm{i}(j-1)\Omega +\lambda )C_t+K_t]I_{j-1}\\ +[{(\mathrm{i}j\Omega +\lambda )}^{2}M+(\mathrm{i}j\Omega +\lambda )C+K]I_j\\ +[{(\mathrm{i}(j+1)\Omega +\lambda )}^2{\tilde{M}}_t+(\mathrm{i}(j+1)\Omega +\lambda ){\tilde{C}}_t+{\tilde{K}}_t]I_{j+1}\}{e}^{\mathrm{i}j\Omega t}=0\\ \{[{(\mathrm{i}j\Omega +\lambda )}^2{M}_t+(\mathrm{i}j\Omega +\lambda )C_t+K_t]I_j\\ +[{(\mathrm{i}(j+1)\Omega +\lambda )}^{2}M+(\mathrm{i}(j+1)\Omega +\lambda )C+K]I_{j+1}\\ +[{(\mathrm{i}(j+2)\Omega +\lambda )}^2{\tilde{M}}_t+(\mathrm{i}(j+2)\Omega +\lambda ){\tilde{C}}_t+{\tilde{K}}_t]I_{j+2}\}{e}^{\mathrm{i}(j+1)\Omega t}=0\end{array}\\ ⋮\end{array}\right.$$

(7)

Let us write Equation (7) as the infinite determinant form:

$$\left[\begin{array}{ccccccc}\ddots & ⋮& ⋮& ⋮& ⋮& ⋮& ⋰\\ \cdots & A_{-2}& {\Delta }_{+1}& 0& 0& 0& \cdots \\ \cdots & {\Delta }_{-1}& A_{-1}& {\Delta }_{+1}& 0& 0& \cdots \\ \cdots & 0& {\Delta }_{-1}& A_0& {\Delta }_{+1}& 0& \cdots \\ \cdots & 0& 0& {\Delta }_{-1}& A_{+1}& {\Delta }_{+1}& \cdots \\ \cdots & 0& 0& 0& {\Delta }_{-1}& A_{+2}& \cdots \\ ⋰& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right]\left[\begin{array}{c}⋮\\ I_{-2}\\ I_{-1}\\ I_{0\begin{array}{c}\end{array}}\\ I_{+1}\\ I_{+2}\\ ⋮\end{array}\right]=0$$

(8)

where:

$$\begin{array}{l}{A}_j=(-j^2{\Omega }^2+2\mathrm{i}\lambda j\Omega +{\lambda }^2)M+(\mathrm{i}j\Omega +\lambda )C+K;\\ {\Delta }_{+1}=(-{(j+1)}^2{\Omega }^2+2\mathrm{i}\lambda (j+1)\Omega +{\lambda }^2){\tilde{M}}_t+(\mathrm{i}(j+1)\Omega +\lambda ){\tilde{C}}_t+{\tilde{K}}_t;\\ {\Delta }_{-1}=(-{(j-1)}^2{\Omega }^2+2\mathrm{i}\lambda (j-1)\Omega +{\lambda }^2)M_t+(\mathrm{i}(j-1)\Omega +\lambda )C_t+K_t.\end{array}$$

(9)

Due to the infinite expansion of Equation (5), the determinant associated with Equation (8) is called the Hill infinite determinant [44].

One reconstructs Equation (8) into Hill’s eigenvalue problem:

$$({\lambda }^2\stackrel{⌢}{M}+\lambda \stackrel{⌢}{C}+\stackrel{⌢}{K})\mathit{\varphi }=0$$

(10)

and the matrices in Equation (10) are also included in Appendix B.

Equation (10) is derived from the original motion Equation (3). Thus, the modal solutions of Equation (10) are equivalent to the solutions of the unbalanced rotor system. Because of the inclusion of the damping matrix $\stackrel{⌢}{C}$, Equation (10) can be expressed in the state space:

$$A\Phi =\lambda B\Phi$$

(11)

where:

$$A=\left[\begin{array}{cc}-\stackrel{⌢}{C}& -\stackrel{⌢}{K}\\ \stackrel{⌢}{K}& 0\end{array}\right],B=\left[\begin{array}{cc}\stackrel{⌢}{M}& 0\\ 0& \stackrel{⌢}{K}\end{array}\right],\Phi =\left[\begin{array}{c}\lambda \mathit{\varphi }\\ \mathit{\varphi }\end{array}\right]$$

(12)

According to Hill’s determinant convergence, the number of harmonic significantly contributing to the quasi-modes $\varphi$, namely the Floquet solution, is finite [45]. A finite truncated harmonic order $j_{\max }$ then needs to be defined. By solving the eigenproblem of Equation (11) in the state space truncated to the $j_{\max }$ harmonic order for a given $\Omega$, one can obtain $n\times (2j_{\max }+1)$ eigenvalues ${\lambda }_l$, where $n=3$ is the number of degrees of freedom of the unbalanced rotor system. According to Equations (4) and (5), each eigenvalue is associated with a complex eigenvector ${\varphi }_l$, expressed in the time domain as:

$${\mathit{\varphi }}_l(t)={\displaystyle \sum _{-j_{\max }}^{j_{\max }}{I}_{jl}{e}^{(\mathrm{i}j\Omega +{\lambda }_l)t}}$$

(13)

where $I_{jl}$ is the $j$th harmonic of the $l$th eigenvector of the unbalanced rotor system. Then, the modal analysis of the rotor system in the inertial frame can be implemented. The stability of the $l$th quasi-mode is determined by the nonzero real parts of the corresponding eigenvalues obtained from Equation (11).

3.1.2. Modal Frequency and Stability

Modal analysis of the unbalanced rotor in the inertial frame is performed in this section, and the parameter values are listed in Table 1. The truncated harmonic order is $j_{\max }=2$, and the step of the rotation speed is 5 rpm. Modal angular frequencies and the real part of eigenvalues varying with the rotation speed are demonstrated in Figure 2. ${\omega }_{\mathrm{b}+,j}$,${\omega }_{\mathrm{b}-,j}$ and ${\omega }_{\mathrm{t},j}$ ($j=-1,0,1$) in Figure 2a represent the $j$th-order harmonic frequency of the forward whirl, backward whirl, and torsional modes, respectively. There exist three notable coupling regions in Figure 2a, and eigenvalues with positive real parts form three regions in the corresponding speed range in Figure 2b, which are demonstrated by the blue, red and green regions, respectively. These three regions are unstable regions. Region A1 corresponds to the region of rotation speed where the -1st-order harmonic frequency of the forward whirl mode is equal to zero. This kind of instability region, caused only by harmonic frequency, is named the ‘secondary instability region’. Region B1 and C1 are caused by the coupling between harmonic and fundamental frequencies of forward whirl modes and torsional modes, termed the ‘principal instability region’ [45]. To conclude, when the truncated harmonic order is $j_{\max }=2$, the number of the instability regions caused by the coupling between lateral modes and torsional modes is two, and the two regions are mostly overlapped.

3.2. Analysis in Rotating Coordinate System

Now, we transform Equation (1) into rotating coordinates $\eta (t)$, $\varsigma (t)$, which are attached to the disk. The coordinate transformation is as follows:

$$\mathbf{X}=\mathbf{T}\mathbf{Y}$$

(14)

where

$$\mathbf{X}=\left\{\begin{array}{c}x\\ y\\ \theta \end{array}\right\},\mathbf{Y}=\left\{\begin{array}{c}\eta \\ \varsigma \\ \theta \end{array}\right\},\mathbf{T}=\left[\begin{array}{ccc}\cos \Omega t& -\sin \Omega t& 0\\ \sin \Omega t& \cos \Omega t& 0\\ 0& 0& 1\end{array}\right]$$

(15)

Then, Equation (1) expressed in rotating coordinates has constant coefficients:

$$\left\{\begin{array}{l}m\ddot{\eta }+c\dot{\eta }-2m\Omega \dot{\varsigma }-2me\Omega \dot{\theta }+(k-m{\Omega }^2)\eta -c\Omega \varsigma =0\\ m\ddot{\varsigma }+me\ddot{\theta }+c\dot{\varsigma }+2m\Omega \dot{\eta }+(k-m{\Omega }^2)\varsigma +c\Omega \eta -me{\Omega }^2\theta =0\\ (m{e}^2+J_p)\ddot{\theta }+me\ddot{\varsigma }+c_{\mathrm{r}}\dot{\theta }+2me\Omega \dot{\eta }+k_{\mathrm{r}}\theta -me{\Omega }^2\varsigma =0\end{array}\right.$$

(16)

and the matrix form of Equation (16) can be written as:

$$\mathbf{M}\ddot{Y}+\mathbf{C}\dot{\mathbf{Y}}+\mathbf{K}\mathbf{Y}=0$$

(17)

The matrices in Equation (17) are given in Appendix B in detail.

The modal characteristics of the rotor system in the rotating frame are obtained. The instability prediction results differ from those in the inertial frame. There exist two unstable regions, as shown in Figure 3. Similar to Region A1 in Figure 2, Region A2 in Figure 3 is also caused by the frequency of the whirl mode being equal to zero, namely, the secondary instability region. The coupling between the whirl mode and torsional mode results in an additional unstable region, referred to as Region B2. Upon comparing the unstable regions depicted in Figure 2 with those in Figure 3, it can be observed that Figure 3a has one less coupling region compared to Figure 2a, and Figure 3b has one less unstable region compared to Figure 2b.

3.3. Comparison of Results in Two Coordinate Systems

The differences between the modal analysis results in two frames are shown in Figure 4 in detail. Because the unstable regions are symmetric about the $x$-axis, Figure 4b displays only the top half of the unstable regions. In both the rotating and inertial frames, there exist Region A3 and Region C3, while Region B3 is exclusively observed in the modal analysis results of the rotating frame. From Figure 4b, unstable speed ranges can be determined. If we only consider the principal instability region, the range of unstable speeds in the inertial frame is from $\Omega$ = 1305 rpm to $\Omega$ = 1624 rpm. However, in the rotational frame, this range is shifted and becomes from $\Omega$ = 1356 rpm to $\Omega$ = 1624 rpm. In other words, the speed range from $\Omega$ = 1305 rpm to $\Omega$ = 1356 rpm is the unstable speed range in the inertial frame but the stable speed range in the rotating frame. It is unreasonable for one rotor system to have different modal calculation results in two coordinate systems.

4. Verification of the Results

In Section 3, we found that the unstable ranges obtained in the rotating frame are not consistent with those obtained in the inertial frame. A time-domain analysis is carried out to verify which one is the true unstable range in this section. The reasons for the disparities in stability prediction are also elucidated.

4.1. Time Domain Results

The fourth-order Runge–Kutta method is performed in both coordinate systems to investigate the unstable speed ranges. Figure 5 and Figure 6 display the sequential displacements along the $x$-axis and the corresponding rotor orbits at the following spin speeds: $\Omega$ = 1355 rpm, $\Omega$ = 1356 rpm, $\Omega$ = 1357 rpm, $\Omega$ = 1624 rpm, and $\Omega$ = 1625 rpm in the inertial frame. Similarly, Figure 7 and Figure 8 demonstrate the results in the rotating frame. An initial displacement is given in the time domain simulation. Given that the rotor system is undamped, the response will exhibit periodic behavior (limit cycle) when it is stable and diverge when it is unstable. Accordingly, the system is stable at rotation speeds $\Omega$ = 1355 rpm, $\Omega$ = 1356 rpm, and $\Omega$ = 1625 rpm, but unstable at $\Omega$ = 1357 rpm and $\Omega$ = 1624 rpm. Especially, the system is stable when $\Omega$ = 1357 rpm both in two frames from the time domain analyses.

However, Hill’s method in the inertial frame gives inconsistent results, i.e., the system is judged as unstable at $\Omega$ = 1357 rpm, which is different from the time domain results. Therefore, Region B3 in Figure 4 is not a real unstable speed range. The coupling region B2 in Figure 2a and unstable region B2 in Figure 2b obtained by Floquet theory and Hill’s method in the inertial coordinate system are false.

4.2. Tracing to the Source of Misjudgment

4.2.1. On the Truncation Order in Hill’s Method

The influence of the truncated harmonic order in Hill’s method is studied in this section. As the number of truncated harmonics increases, the number of modal frequency coupling regions rises from two to six, shown in Figure 9. Meanwhile, the unstable regions obtained by the real part of eigenvalues do not change when the number of the harmonic truncation order increases. Therefore, we can conclude that the number of harmonic truncation in Hill’s method does not impact the prediction of instability for an unbalanced rotor in the inertial frame.

Furthermore, one can notice that, as the number of the harmonic truncation order $j_{\max }$ increases, all the newly identified coupling regions correspond to the unstable region 3. The false unstable region 2 cannot be eliminated by increasing the truncation order and is always at the top of the coupling regions.

4.2.2. On the Correspondence of Solutions in Two Frames

The coordinate transformation matrix $T$ in Equation (14) has time-varying terms, which may affect the solutions of motion equations and lead to misjudgment. The cause of the false unstable region is studied from the perspective of the solutions of the motion equations in this section.

The modal vectors ${\varphi }_i$ and eigenvalues ${\lambda }_i$ of the rotor system can be obtained by modal analysis at a certain speed in both two frames. Then, the solutions of the motion equations $X_i$ can be derived by Equation (18), wherein $i$ stands for the mode number. To verify the effect of the coordinate transformation matrix $T$ on the solutions, Equations (19) and (20) are used to transfer the solution from the inertial (rotating) frame to that in the rotating (inertial) frame, respectively.

$${\mathbf{X}}_i={\mathit{\varphi }}_i{e}^{{\lambda }_{i}t}$$

(18)

$${\mathbf{X}}_i^{’}=\mathbf{T}{\mathbf{Y}}_i$$

(19)

$${\mathbf{Y}}_i^{’}={\mathbf{T}}^{-1}{\mathbf{X}}_i$$

(20)

Taking one rotational speed ${\Omega }_i$ in the false unstable region, we calculate the solutions ${\mathit{x}}_i$ and ${\mathit{y}}_i$ of the motion equations given in Equations (11) and (17) in both frames. Equations (19) and (20) are used to transfer ${\mathit{x}}_i$ and ${\mathit{y}}_i$ into ${\mathit{y}}_i^{’}$ and ${\mathit{x}}_i^{’}$ respectively, wherein ${\mathit{x}}_i^{’}$ is the transformed solution in the inertial frame and ${\mathit{y}}_i^{’}$ is that in the rotating frame. The convergence of Hill’s method is achieved by $j_{\max }=2$ so that the number of modes is nine in the inertial frame. The number of modes of the system is three in the rotating frame. By substituting solutions ${\mathit{x}}_i$, ${\mathit{x}}_i^{’}$, ${\mathit{y}}_i$ and ${\mathit{y}}_i^{’}$ into Equations (1) and (16) in two coordinate systems, whether the solution satisfies the motion equations can be judged. The check mark (✓) means that the solution of the mode satisfies the equations, and the cross mark (✗) means that it does not, as demonstrated in

Table 2. Rotor’s physical and simulation parameters.

Solutions	${\mathit{M}}_1$	${\mathit{M}}_2$	${\mathit{M}}_3$	${\mathit{M}}_4$	${\mathit{M}}_5$	${\mathit{M}}_6$	${\mathit{M}}_7$	${\mathit{M}}_8$	${\mathit{M}}_9$
obtained in inertial frame (${\mathit{x}}_i$)	✓	✗	✗	✗	✗	✓	✓	✗	✗
transferred to rotating frame (${\mathit{y}}_i^{’}$)	✓	✗	✗	✗	✗	✓	✓	✗	✗
obtained in rotating frame (${\mathit{y}}_i$)	✓	-	-	-	-	✓	✓	-	-
transferred to inertial frame (${\mathit{x}}_i^{’}$)	✓	-	-	-	-	✓	✓	-	-

, wherein $M_k$ is the $k$th mode.

From Table 2, it can be observed that out of nine solutions $x_i$ obtained in the inertial frame, six of them do not satisfy the motion equations. Similarly, the same issue arises for the solutions $y_i^{’}$ obtained in the rotating frame through coordinate transformation. On the contrary, all the solutions $y_i$ obtained in the rotating frame satisfy the motion equations, as well as the solutions $x_i^{’}$ in the inertial frame after coordinate transformation. Additionally, one can notice that even though the convergence of Hill’s method is achieved by $j_{\max }=3$, six out of the fifteen solutions still fail to satisfy the motion equations in both two frames. Therefore, it is reasonable to conclude that the solutions that do not satisfy the motion equations in the inertial frame are responsible for the false modal frequency coupling region and unstable region observed in the results of the modal analysis. Moreover, increasing the truncated order will not eliminate this error.

4.2.3. On the Matrix Truncation in the Solving Process

Since the error in Hill’s method in the inertial frame is independent of the harmonic truncation, it is necessary to explore the calculation process of the motion equations to determine the origin of the error. If the truncated harmonic number is 2, Equation (6) will become Equation (21). Equation (22) is the matrix form of Equation (21), where ${\Delta }_{-1}$ and ${\Delta }_{+1}$ are given in detail in Equation (8).

The coefficient matrix in Equation (22) is not a square matrix; therefore, the modal analysis cannot be performed directly. To obtain a square coefficient matrix, we need to remove ${\Delta }_{-1}$ in the first row and ${\Delta }_{+1}$ in the last row of the matrix, thus Equation (23) is obtained. Correspondingly, the term of $I_{-2}$ and $I_{+2}$ in Equation (21) should be removed as well; therefore, the first and the last equation in Equation (21) do not hold any more. This phenomenon provides an explanation for the consistent occurrence of six sets of solutions that fail to satisfy the motion equations in the inertial frame.

In summary, the truncation of motion equations causes the misjudgments of stability. The false unstable region in the inertial frame cannot be eliminated by only increasing the number of the harmonic truncation order in Hill’s method. The absence of several terms by matrix truncation means that the modal solutions do not satisfy the motion equations. The false modal frequency coupling region always appears at the top of the coupling regions, where one needs to pay more attention to the stability results.

$$\left\{\begin{array}{l}\{[{(-2\mathrm{i}\Omega +\lambda )}^2{M}_t+(-2\mathrm{i}\Omega +\lambda )C_t+K_t]I_{-2}\\ +[{(-\mathrm{i}\Omega +\lambda )}^{2}M+(-\mathrm{i}\Omega +\lambda )C+K]I_{-1}\\ +[{\lambda }^2{\tilde{M}}_t+\lambda {\tilde{C}}_t+{\tilde{K}}_t]I_0\}{e}^{-\mathrm{i}\Omega t}=0\cdots \cdots \cdots \cdots \cdots \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot (1)\\ \{[{(-\mathrm{i}\Omega +\lambda )}^2{M}_t+(-\mathrm{i}\Omega +\lambda )C_t+K_t]I_{-1}\\ +[{\lambda }^{2}M+\lambda C+K]I_0\\ +[{(\mathrm{i}\Omega +\lambda )}^2{\tilde{M}}_t+(\mathrm{i}\Omega +\lambda ){\tilde{C}}_t+{\tilde{K}}_t]I_{+1}\}{e}^{0\mathrm{i}\Omega t}=0\cdots \cdot \cdot \cdot (2)\\ \{[{\lambda }^2{M}_t+\lambda C_t+K_t]I_0\\ +[{(\mathrm{i}\Omega +\lambda )}^{2}M+(\mathrm{i}\Omega +\lambda )C+K]I_{+1}\\ +[{(2\mathrm{i}\Omega +\lambda )}^2{\tilde{M}}_t+(2\mathrm{i}\Omega +\lambda ){\tilde{C}}_t+{\tilde{K}}_t]I_{+2}\}{e}^{\mathrm{i}\Omega t}=0\cdot \cdot \cdot \cdot (3)\end{array}\right.$$

(21)

$$\left[\begin{array}{ccccc}{\Delta }_{-1}& A_{-1}& {\Delta }_{+1}& 0& 0\\ 0& {\Delta }_{-1}& A_0& {\Delta }_{+1}& 0\\ 0& 0& {\Delta }_{-1}& A_{+1}& {\Delta }_{+1}\end{array}\right]\left[\begin{array}{c}{I}_{-2}\\ I_{-1}\\ I_{0\begin{array}{c}\end{array}}\\ I_{+1}\\ I_{+2}\end{array}\right]=0$$

(22)

$$\left[\begin{array}{ccc}{A}_{-1}& {\Delta }_{+1}& 0\\ {\Delta }_{-1}& A_0& {\Delta }_{+1}\\ 0& {\Delta }_{-1}& A_{+1}\end{array}\right]\left[\begin{array}{c}{I}_{-1}\\ I_{0\begin{array}{c}\end{array}}\\ I_{+1}\end{array}\right]=0$$

(23)

5. Conclusions

In this paper, we study the influence of coordinate system choice on stability estimation of lateral–torsional coupled vibration in an unbalanced rotor and try to explain the reason for the stability misjudgment of Hill’s method in inertial coordinates. The findings of this research can offer enhanced theoretical supports for the stability analysis of the lateral–torsional coupling vibration of rotors.

The unstable regions obtained by modal analyses are different in the inertial and rotating frames. The unstable region exhibits a wider range of speeds in the inertial frame compared to that in the rotating frame. It is confirmed by the time integral that the redundant unstable region obtained in the inertial frame is false.

The misprediction of unstable regions under the inertial frame can be attributed to the two truncations in Floquet’s theory and Hill’s method, namely the harmonic truncation and the matrix truncation. Harmonic truncation makes Hill’s method infinite determinant finite, leading to a non-square matrix. Subsequently, false solutions are generated, which finally give rise to a misjudgment of the rotor stability.

We suggest utilizing the rotating frame to analyze the lateral–torsional coupled vibration if the rotor can be assumed as axially symmetric. This approach ensures the constant coefficients of motion equations, avoids truncations in the calculation process, and enhances the precision of rotor system instability prediction. If the rotor system is not axially symmetric due to complex supports or cracks, converting all coefficients of motion equations to constants in the rotating frame becomes infeasible. In such cases, the errors caused by the harmonic truncation and the matrix truncation should be taken into account when using the Floquet–Hill formulation. Particularly, when two close unstable regions appear, a false unstable region caused by erroneous eigenvalues may exist. To avoid the stability misjudgment in the analysis of lateral–torsional coupled vibration, the eigenvalues should be filtered to extract the ones satisfying the motion equations of the rotor system by either time domain calculation or eigenvalue filtering method.

Further works will focus on unbalanced rotors in engineering. Experimental studies will be also carried out to verify the mispredicted unstable results.