# Analytical model for nonlinear vibration of flexible rotor system

## Min Wu1, Shengbo Yang2, Dongjiang Han3, Daren Yu4, Jinfu Yang5

1, 2, 4School of Energy Science and Engineering, Harbin Institute of Technology, Harbin, China

3, 5Institute of Engineering Thermophysics, Chinese Academy of Sciences, Beijing, China

1Corresponding author

Journal of Vibroengineering, Vol. 18, Issue 8, 2016, p. 4980-4994. https://doi.org/10.21595/jve.2016.17441
Received 22 July 2016; received in revised form 8 November 2016; accepted 18 November 2016; published 31 December 2016

Copyright © 2016 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Views 71
Abstract.

An analytical model is proposed to analyze a series of typical nonlinear behaviors of flexible rotor system, such as resonance, oscillation, whirl and whip. The model is constructed by introducing a defined nonlinear scale factor $\epsilon$, nonlinear stiffness and nonlinear damping. Based on multi-scale method, the analytical solutions of steady-state and transient-state are derived, and the nonlinear natural frequency and Frequency Response Equation (FRE) are obtained. A transient time scale factor ${t}_{1}$ is defined to reflect the transient-state influence on steady-state solution. The experimental result also verifies the rationality and validity of the analytical model and the analytical solutions.

Keywords: flexible rotor system, nonlinear vibration, analytical model, nonlinear scale factor, analytical solution.

#### 1. Introduction

Rotating machinery, as the core device of a dynamic system provides power and a series of necessary functions, which is widely used in aerospace industry, ship engineering and distributed energy resource in the modern time. Focus on the development of high-speed miniature equipments, the primary perspective is to enhance rotational speed, system power, thrust-weight ratio, and even to achieve integration, intelligence. In this context, the stability of a rotor system becomes vital to meet the requirements of high performance and long service. As the rotational speed increases, instability phenomenon caused by nonlinear vibrations can occur easily, which require more efficient stability analysis methods and engineering applications [1].

Generally, a rotor system is composed of rotor, bearing, shaft and sealing elements. The system responses to external excitations can be divided into three modes, linear dominant mode, nonlinear dominant mode and mixed mode. The linear vibration method is often used in traditional rotor dynamic researches while the nonlinear responses are not obvious. In conventional engineering, reasonable linearization can significantly reduce the difficulties of analysis and calculations, which makes linear rotor dynamic theory be fully developed and widely used [2]. On the contrary, when system presents obvious nonlinear behaviors, the magnitude of nonlinear responses will greater than the linear ones, which can lead to a series of nonlinear phenomenon, such as harmonic resonance, subharmonic resonance, ultraharmonic resonance, quasi-period, bifurcation, chaos for single degree of freedom (DOF) system and combination resonance additionally for multiple DOF system [3]. When the magnitudes of nonlinear and linear responses are the same, namely mixed mode, the nonlinear responses caused by interaction cannot be ignored, which may need other efficient methods, such as decoupling approach.

Many scholars made important research progresses. Lund model [4] was proposed in 1965 for the first time by solving Reynold’s equations using perturbation method. Considering both rotor and bearing, Lund presented an equation with 8 coefficients of linearized stiffness and damping force equivalent to the oil film force, which could be solved based on linear stability theory. According to the linearized oil film force model, the eigenvalue of the homogeneous equation can be used to calculate the speed threshold of system instability. Glienicke [5] presented a more systematic theoretical and experimental study of Lund model. A series of linearized equivalence of oil film force are established based on Lund’s theory. If the rotor whirl trajectory becomes ellipse, Lund model can be simplified into Alford model [6], in which the 8 coefficients can be simplified into principal terms and cross terms. Alford model presumed that cross stiffness and principal damping are the key factors to the instability problems of rotor system. The former’s increase can weaken system stability, but the latter’s is the opposite. Hori [7] studied the stability of rotor system. Ehrich [8] studied the effect of Alford force on the self-excited vibration. Based on short bearing theory, Childs [9] introduced added mass so as to describe the sealing fluid exciting force more accurately, and considered sealing fluid circumferential speed which Alford model ignored. In 1987, Muszynska and Bently [10, 11] defined fluid circumferential average velocity ratio to be the key parameter which could reflect oil film dynamic characteristics based on a large number of experiment data. Muszynska model can explain the instability mechanism of oil film for a high-speed rotor system, and especially describe the relationship of oil film cross stiffness (which causes instability), average velocity ratio, rotational speed and external damping. Ravikovich [12] studied whirl and whip caused by nonlinear factors of oil film force, and the influence on power and efficiency losses.

The problem of linear dominant mode can be dealt with the linear theory. However, when it comes to nonlinear dominant and mixed modes, the nonlinear responses usually cannot be represented accurately [13] by classic models based on linear theory, especially in response mechanism, rationality and reliability analysis. Moreover, both of the mechanism analysis and the reliability analysis cannot achieve the engineering requirements. Therefore, in this paper we propose an analytical model of a high-speed flexible rotor system, which can inherit and develop the traditional linear models, and analyze a series of typical nonlinear behaviors in the nonlinear dominant and mixed modes. The analytical solutions are also presented via multi-scale method, and then verified with experiment result.

#### 2. Model formulation

As shown in Fig. 1, the rotor system is considered as a single-span rotor supported by force $F$, which can be divided into stiffness force and damping force. Noticing that, the force $F\text{'}$, which is perpendicular to $F$ and caused by sealing force, is not considered in this paper, so that we can formulate this model with single DOF instead of the traditional crossing two DOF. By introducing nonlinear Hooke’s law and nonlinear damping force, we can construct nonlinear analytical model with single DOF.

Fig. 1. Dynamic analysis for a rotor

#### 2.1. Nonlinear stiffness

From the perspective of elastic theory, based on the hypothesis of small perturbation and deformation, traditional model for representing flexible rotor behaviors can be generally satisfied with Hooke’s law. However, under the circumstances of large perturbation or deformation, the elastic potential energy of one-dimensional system [14] can be expressed as:

(1)
$U\left(x\right)=\frac{1}{2}k{x}^{2}+\frac{1}{3}\xi {x}^{3}+\frac{1}{4}\eta {x}^{4}+\cdots .$

For engineering practice, the potential energy of many structures can be considered as symmetry [15]. Thus, the even order terms of Eq. (1) are retained:

(2)
$U\left(x\right)=\frac{1}{2}k{x}^{2}+\frac{1}{4}\eta {x}^{4}.$

The stiffness force ${f}_{k}$, namely the elastic force, is expressed as:

(3)
${f}_{k}=-\frac{dU}{dx}=-kx-\eta {x}^{3}.$

where $k$ is the coefficient of linear elastic deformation, namely stiffness, $\eta$ is the coefficient of nonlinear elastic deformation. If 0, Eq. (3) satisfies the linear Hooke’s law, namely linear spring characteristics. If 0, ${f}_{k}$ increases along with deformation increasing, namely hard spring characteristics. On the contrary, if 0, ${f}_{k}$ decreases along with deformation increasing, namely soft spring characteristics [16].

#### 2.2. Nonlinear damping

Generally, the damping dissipation process of non-conservative system [17] can be given as:

(4)
$m\stackrel{¨}{x}+q\left(\stackrel{˙}{x}\right)+kx=0.$

Thus, the damping force is described as a generalized form:

(5)

If 0, Eq. (5) represents dry friction force. If 1, Eq. (5) represents linear viscous damping force, which is applied to low-speed movement in fluid. If $n=$ 2, Eq. (5) represents nonlinear viscous damping force, which is applied to high-speed movement [18]. Thus, the damping force ${f}_{c}$ is expressed as:

(6)
${f}_{c}=-c\stackrel{˙}{x}-\gamma {\stackrel{˙}{x}}^{2},$

where $c$ represents linear damping, $\gamma$ is the coefficient of nonlinear damping.

#### 2.3. Analytical model

Considering the single-span single-disc rotor system shown in Fig. 1, the parameters are as follows, $m$ for disc mass, $r$ for eccentricity, $\mathrm{\Omega }$ for rotational frequency. And shaft mass is ignored, counterclockwise is positive. For general liquid or gas bearing, there is interaction of support and whirl, due to wedge flow characteristics. In this context, we consider only bearing supporting effect, so that based on Newton’s second law, equation of motion in the $x$-axis direction can be described as:

(7)
$m\frac{{d}^{2}{x}_{C}}{d{t}^{2}}=-kx-\eta {x}^{3}-c\stackrel{˙}{x}-\gamma {\stackrel{˙}{x}}^{2},$

where:

(8)
${x}_{C}=x+r\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Omega }t.$

For briefness and clarity, some substitutions are given as:

(9)
$\begin{array}{llll}\mu =\frac{c}{2m},& \epsilon {\beta }_{0}^{2}=\frac{\eta }{m},& {\omega }_{0}^{2}=\frac{k}{m},& \epsilon {\alpha }_{0}^{2}=\frac{\gamma }{m}\end{array},$

where $\mu$ is linear damping coefficient, ${\omega }_{0}$ is natural frequency, ${\beta }_{0}$ is nonlinear stiffness amplification coefficient, ${\alpha }_{0}$ is nonlinear damping amplification coefficient, $\epsilon$ is the nonlinear scale factor which denotes the impact of nonlinear responses of rotor system.

Then, the analytical model for nonlinear vibration of flexible rotor system can be written as:

(10)
$\stackrel{¨}{x}+2\mu \stackrel{˙}{x}+\epsilon {\alpha }_{0}^{2}{\stackrel{˙}{x}}^{2}+{\omega }_{0}^{2}x+\epsilon {\beta }_{0}^{2}{x}^{3}={\mathrm{\Omega }}^{2}r\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Omega }t.$

Noticing that, Eq. (10) has a similar form with Duffing equation [19], in which ${x}^{3}$ term is used to describe nonlinear elastic restoring force. Numerous researches show that the response characteristics of complex rotor system is consistent with the analysis of single DOF Duffing equation in many respects, such as period-doubling bifurcation, Hopf bifurcation, period-3 bifurcation, enter/exit chaos phenomenon [20].

This analytical model Eq. (10) deals with centrifugal force generated by unbalanced mass, essentially not as the same as Muszynska model, which aims to deal with the clearance force caused by sealing fluid (that is why the negative sign exists as a reaction force in Eq. (11), which is different from that in Eq. (10) as an action force). If the nonlinear scale factor $\epsilon$ is very small, the nonlinear terms can be neglected. Then considering linear interaction of stiffness and damping, i.e., cross stiffness and cross damping force, Eq. (10) can be transformed into Muszynska model [11]:

(11)
$-\left[\begin{array}{l}{F}_{x}\\ {F}_{y}\end{array}\right]=\left[\begin{array}{cc}k-{m}_{f}{\tau }^{2}{\omega }^{2}& \tau \omega c\\ -\tau \omega c& k-{m}_{f}{\tau }^{2}{\omega }^{2}\end{array}\right]\left[\begin{array}{l}x\\ y\end{array}\right]+\left[\begin{array}{cc}c& 2\tau \omega {m}_{f}\\ -2\tau \omega {m}_{f}& c\end{array}\right]\left[\begin{array}{l}\stackrel{˙}{x}\\ \stackrel{˙}{y}\end{array}\right]+\left[\begin{array}{cc}{m}_{f}& 0\\ 0& {m}_{f}\end{array}\right]\left[\begin{array}{l}\stackrel{¨}{x}\\ \stackrel{¨}{y}\end{array}\right].$

On this basis, if ignore inertial force, and substitute coefficients for simple letters, Eq. (11) can be simplified as Alford model [6]:

(12)
$-\left[\begin{array}{l}{F}_{x}\\ {F}_{y}\end{array}\right]=\left[\begin{array}{cc}{k}_{1}& {k}_{2}\\ -{k}_{2}& {k}_{1}\end{array}\right]\left[\begin{array}{l}x\\ y\end{array}\right]+\left[\begin{array}{cc}{c}_{1}& {c}_{2}\\ -{c}_{2}& {c}_{1}\end{array}\right]\left[\begin{array}{l}\stackrel{˙}{x}\\ \stackrel{˙}{y}\end{array}\right].$

Apparently, if the cross damping force is not considered, Eq. (11) can be simplified as Ravikovich model [12]:

(13)
$m\stackrel{¨}{x}+2c\stackrel{˙}{x}+2{k}_{xx}x+2{k}_{xy}y=-mr{\omega }^{2}\mathrm{s}\mathrm{i}\mathrm{n}\omega t,$
$m\stackrel{¨}{y}+2c\stackrel{˙}{y}+2{k}_{yy}y+2{k}_{yx}x=-mr{\omega }^{2}\mathrm{c}\mathrm{o}\mathrm{s}\omega t,$

where $F$ is excitation force, $k$ is stiffness, $c$ is damping, $\tau$ is fluid circumferential average velocity ratio, ${m}_{f}$ is inertial mass, $\omega$ is rotational speed, $m$ is mass, $r$ is eccentricity, in Eq. (11-13).

Fig. 2. Schematic of linear response and nonlinear response

Apparently, these models mentioned above used cross terms to describe nonlinear characteristics indirectly. In this paper, the analytical model introduces nonlinear stiffness force and nonlinear damping force instead of cross terms, i.e., Eq. (3) and Eq. (6), to directly describe the nonlinear effects with clear physical significance, and then to do further research on the interaction of linear and nonlinear forces, and the mechanism of nonlinear dynamic responses under the conditions of large displacement and high speed. Fig. 2 shows the linear and nonlinear response characteristics in frequency domain.

#### 3.1. Free vibration mode

Multi-scale method can calculate both of the steady-state response and the transient-state response. In order to represent different time scales, ${T}_{n}$ is introduced and defined as:

(14)

Then, a nonlinear response can be written as a function of various time scales:

(15)
$x=\sum _{n=0}^{m}{\epsilon }^{n}{x}_{n}\left({T}_{0},{T}_{1},{T}_{2},...,{T}_{m}\right),$

where $n$ depends on the calculation accuracy.

When an external excitation does not exist, the free vibration mode of the analytical model can be expressed as:

(16)
$\stackrel{¨}{x}+2\mu \stackrel{˙}{x}+\epsilon {\alpha }_{0}^{2}{\stackrel{˙}{x}}^{2}+{\omega }_{0}^{2}x+\epsilon {\beta }_{0}^{2}{x}^{3}=0.$

Based on multi-scale method, we suppose:

(17)
$x={x}_{0}\left({T}_{0},{T}_{1}\right)+\epsilon {x}_{1}\left({T}_{0},{T}_{1}\right),$

where:

(18)

Thus:

(19)
$\stackrel{˙}{x}=\frac{\partial {x}_{0}}{\partial {T}_{0}}+\epsilon \frac{\partial {x}_{0}}{\partial {T}_{1}}+\epsilon \frac{\partial {x}_{1}}{\partial {T}_{0}}+{\epsilon }^{2}\frac{\partial {x}_{1}}{\partial {T}_{1}},$
$\stackrel{¨}{x}=\frac{{\partial }^{2}{x}_{0}}{\partial {{T}_{0}}^{2}}+2\epsilon \frac{{\partial }^{2}{x}_{0}}{\partial {T}_{0}\partial {T}_{1}}+{\epsilon }^{2}\frac{{\partial }^{2}{x}_{0}}{\partial {{T}_{1}}^{2}}+\epsilon \frac{{\partial }^{2}{x}_{1}}{\partial {{T}_{0}}^{2}}+2{\epsilon }^{2}\frac{{\partial }^{2}{x}_{1}}{\partial {T}_{0}\partial {T}_{1}}+{\epsilon }^{3}\frac{{\partial }^{2}{x}_{1}}{\partial {{T}_{1}}^{2}}.$

Then, according to the power orders of $\epsilon$, there are:

(20)
${\epsilon }^{0}:{D}_{0}^{2}{x}_{0}+2\mu {D}_{0}{x}_{0}+{\omega }_{0}^{2}{x}_{0}=0,$
(21)
${\epsilon }^{1}:{D}_{0}^{2}{x}_{1}+2\mu {D}_{0}{x}_{1}+{\omega }_{0}^{2}{x}_{1}=-2{D}_{0}{D}_{1}{x}_{0}-2\mu {D}_{1}{x}_{0}-{\alpha }_{0}^{2}{\left({D}_{0}{x}_{0}\right)}^{2}-{\beta }_{0}^{2}{x}_{0}^{3},$

where ${D}_{n}$ (0, 1) is the partial differential operator. The solution of Eq. (20) is:

(22)

where $C$ is an unknown complex function, the overline represents the conjugate function. Take Eq. (22) into Eq. (21), we have:

(23)
${D}_{0}^{2}{x}_{1}+2\mu {D}_{0}{x}_{1}+{\omega }_{0}^{2}{x}_{1}=-\left[2{D}_{1}C\left({r}_{1}+\mu \right)+3{\beta }_{0}^{2}{C}^{2}\stackrel{-}{C}{e}^{-2\mu {T}_{0}}\right]{e}^{{r}_{1}{T}_{0}}$

where $cc$ represents the conjugate functions of each term before. The coefficient of first term ${e}^{{r}_{1}{T}_{0}}$ must be zero to avoid duration term. At any point without acceleration, ${T}_{0}=0$, so ${e}^{-2\mu {T}_{0}}=1$. There is:

(24)
$2{D}_{1}C\left(\mu +{r}_{1}\right)+3{\beta }_{0}^{2}{C}^{2}\stackrel{-}{C}=0.$

Assuming:

(25)
$C=\frac{1}{2}a{e}^{j\phi },$

where $a$ and $\phi$ are the functions of ${T}_{1}$. According to the real part and the imaginary part respectively, $C$ can be obtained:

(26)
$C=\frac{1}{2}{a}_{0}{e}^{j\frac{3{\beta }_{0}^{2}{a}_{0}^{2}}{8\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}}{T}_{1}+j{\phi }_{0}},$

where ${a}_{0}$, ${\phi }_{0}$ are the integration constants. Take Eq. (26) and Eq. (22) into Eq. (23), ${x}_{1}$ can be solved, then the solution $x$ is:

(27)
$x=\frac{1}{2}{a}_{0}{e}^{-\mu t+j\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}t+j\frac{3\epsilon {\beta }_{0}^{2}{a}_{0}^{2}}{8\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}}t+j{\phi }_{0}}-\frac{1}{4}\epsilon {\alpha }_{0}^{2}{a}_{0}^{2}{e}^{-2\mu t}$

Thus, the natural frequency is:

(28)
$\omega =\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}+\frac{3\epsilon {\beta }_{0}^{2}{a}_{0}^{2}}{8\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}}=\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}+\frac{3{a}_{0}^{2}}{8m\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}}\eta ,$

where the first term on the right side is linear natural frequency and the second is nonlinear natural frequency. Obviously, ${\omega }_{0}$, $\mu$, ${\beta }_{0}$ and $\eta$ have an impact on the natural frequency but ${\alpha }_{0}$, $\gamma$ does not.

#### 3.2. Forced vibration mode

When the external excitation exists, for instance centrifugal force caused by unbalanced mass of the rotor, the forced vibration mode of analytical model can be expressed as:

(29)
$\stackrel{¨}{x}+2\mu \stackrel{˙}{x}+\epsilon {\alpha }_{0}^{2}{\stackrel{˙}{x}}^{2}+{\omega }_{0}^{2}x+\epsilon {\beta }_{0}^{2}{x}^{3}={\mathrm{\Omega }}^{2}r\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Omega }t.$

Similarly, using multi-scale method, according to the power orders of $\epsilon$, there are:

(30)
${\epsilon }^{0}:{D}_{0}^{2}{x}_{0}+2\mu {D}_{0}{x}_{0}+{\omega }_{0}^{2}{x}_{0}=\frac{1}{2}{\mathrm{\Omega }}^{2}r\left({e}^{j\mathrm{\Omega }{T}_{0}}+{e}^{-j\mathrm{\Omega }{T}_{0}}\right),$
(31)
${\epsilon }^{1}:{D}_{0}^{2}{x}_{1}+2\mu {D}_{0}{x}_{1}+{\omega }_{0}^{2}{x}_{1}=-2{D}_{0}{D}_{1}{x}_{0}-2\mu {D}_{1}{x}_{0}-{\alpha }_{0}^{2}{\left({D}_{0}{x}_{0}\right)}^{2}-{\beta }_{0}^{2}{x}_{0}^{3}.$

The solution of Eq. (30) is:

(32)
${x}_{0}=C\left({T}_{1}\right){e}^{{r}_{1}{T}_{0}}+\stackrel{-}{C}\left({T}_{1}\right){e}^{{\stackrel{-}{r}}_{1}{T}_{0}}+Q{e}^{j\mathrm{\Omega }{T}_{0}}+\stackrel{-}{Q}{e}^{-j\mathrm{\Omega }{T}_{0}},$

where $C$ is an unknown complex function, the overline represents the conjugate function. Take Eq. (32) into Eq. (31), we have:

(33)
${D}_{0}^{2}{x}_{1}+2\mu {D}_{0}{x}_{1}+{\omega }_{0}^{2}{x}_{1}=-\left[2{D}_{1}C\left({r}_{1}+\mu \right)+6{\beta }_{0}^{2}CQQ+3{\beta }_{0}^{2}{C}^{2}\stackrel{-}{C}{e}^{-2\mu {T}_{0}}\right]{e}^{{r}_{1}{T}_{0}}$

Similarly, to avoid duration term, the coefficient of first term ${e}^{{r}_{1}{T}_{0}}$ must be zero. At any point without acceleration, ${T}_{0}=\text{0}$, so ${e}^{-2\mu {T}_{0}}=\text{1}$. There is:

(34)
$2{D}_{1}C\left({r}_{1}+\mu \right)+6{\beta }_{0}^{2}CQ\stackrel{-}{Q}+3{\beta }_{0}^{2}{C}^{2}\stackrel{-}{C}=0.$

Take Eq. (25) into Eq. (34), according to the real part and the imaginary part respectively, $C$ can be obtained.

(35)

where ${a}_{0}$, ${\phi }_{0}$ are the integration constants. Take Eq. (35) and Eq. (32) into Eq. (33), ${x}_{1}$ can be solved, then the solution $x$ is:

$x=\frac{1}{2}{a}_{0}\mathrm{exp}\left(-\mu t+j\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}t+j\frac{\frac{6\epsilon {\beta }_{0}^{2}{\mathrm{\Omega }}^{4}{r}^{2}}{{\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}\right)}^{2}+4{\mu }^{2}{\mathrm{\Omega }}^{2}}+3\epsilon {\beta }_{0}^{2}{a}_{0}^{2}}{8\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}}t+j{\phi }_{0}\right)$
$+\frac{1}{2}\frac{{\mathrm{\Omega }}^{2}r}{{\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}+2j\mu \mathrm{\Omega }}{e}^{j\mathrm{\Omega }t}-\frac{1}{4}\epsilon {\alpha }_{0}^{2}{a}_{0}^{2}{e}^{-2\mu t}$
$-\frac{\epsilon {\alpha }_{0}^{2}{\mathrm{\Omega }}^{6}{r}^{2}}{4{\omega }_{0}^{2}\left[{\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}\right)}^{2}+4{\mu }^{2}{\mathrm{\Omega }}^{2}\right]}-\frac{3\epsilon {\beta }_{0}^{2}{\mathrm{\Omega }}^{6}{r}^{3}}{8{\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}+2j\mu \mathrm{\Omega }\right)}^{3}\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}-2j\mu \mathrm{\Omega }\right)}{e}^{j\mathrm{\Omega }t}$
$+\frac{\epsilon {\alpha }_{0}^{2}{\mathrm{\Omega }}^{6}{r}^{2}}{4\left({\omega }_{0}^{2}-4{\mathrm{\Omega }}^{2}+4j\mu \mathrm{\Omega }\right){\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}+2j\mu \mathrm{\Omega }\right)}^{2}}{e}^{2j\mathrm{\Omega }t}-\frac{3\epsilon {\beta }_{0}^{2}{a}_{0}^{2}{\mathrm{\Omega }}^{2}r}{4\left[{\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}\right)}^{2}+4{\mu }^{2}{\mathrm{\Omega }}^{2}\right]}{e}^{-2\mu t+j\mathrm{\Omega }t}$
$-\frac{\epsilon {\beta }_{0}^{2}{\mathrm{\Omega }}^{6}{r}^{3}}{8\left({\omega }_{0}^{2}-9{\mathrm{\Omega }}^{2}+6j\mu \mathrm{\Omega }\right){\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}+2j\mu \mathrm{\Omega }\right)}^{3}}{e}^{3j\mathrm{\Omega }t}-\frac{\epsilon {\alpha }_{0}^{2}{a}_{0}{\mathrm{\Omega }}^{2}r\left(j\mu +\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}\right)}{2\left(2\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}+\mathrm{\Omega }\right)\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}+2j\mu \mathrm{\Omega }\right)}$
$×\mathrm{exp}\left(-\mu t+j\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}t+j\mathrm{\Omega }t+j\frac{\frac{6\epsilon {\beta }_{0}^{2}{\mathrm{\Omega }}^{4}{r}^{2}}{{\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}\right)}^{2}+4{\mu }^{2}{\mathrm{\Omega }}^{2}}+3\epsilon {\beta }_{0}^{2}{a}_{0}^{2}}{8\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}}t+j{\phi }_{0}\right)$
$-\frac{\epsilon {\alpha }_{0}^{2}{a}_{0}{\mathrm{\Omega }}^{2}r\left(j\mu +\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}\right)}{2\left(2\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}-\mathrm{\Omega }\right)\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}-2j\mu \mathrm{\Omega }\right)}$
$×\mathrm{exp}\left(-\mu t+j\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}t-j\mathrm{\Omega }t+j\frac{\frac{6\epsilon {\beta }_{0}^{2}{\mathrm{\Omega }}^{4}{r}^{2}}{{\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}\right)}^{2}+4{\mu }^{2}{\mathrm{\Omega }}^{2}}+3\epsilon {\beta }_{0}^{2}{a}_{0}^{2}}{8\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}}t+j{\phi }_{0}\right)$
$+\frac{3\epsilon {\beta }_{0}^{2}{a}_{0}^{2}{\mathrm{\Omega }}^{2}r\mathrm{e}\mathrm{x}\mathrm{p}\left(-2\mu t+2j\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}t+j\mathrm{\Omega }t+j\frac{\frac{6\epsilon {\beta }_{0}^{2}{\mathrm{\Omega }}^{4}{r}^{2}}{{\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}\right)}^{2}+4{\mu }^{2}{\mathrm{\Omega }}^{2}}+3\epsilon {\beta }_{0}^{2}{a}_{0}^{2}}{4\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}}t+2j{\phi }_{0}\right)}{8\left(-4{\mu }^{2}+3{\omega }_{0}^{2}+4j\mu \sqrt{{\omega }_{0}^{2}-{\mu }^{2}}+{\mathrm{\Omega }}^{2}+2j\mu \mathrm{\Omega }+4\mathrm{\Omega }\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}\right)\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}+2j\mu \mathrm{\Omega }\right)}$
$+\frac{3\epsilon {\beta }_{0}^{2}{a}_{0}^{2}{\mathrm{\Omega }}^{2}r\mathrm{e}\mathrm{x}\mathrm{p}\left(-2\mu t+2j\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}t-j\mathrm{\Omega }t+j\frac{\frac{6\epsilon {\beta }_{0}^{2}{\mathrm{\Omega }}^{4}{r}^{2}}{{\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}\right)}^{2}+4{\mu }^{2}{\mathrm{\Omega }}^{2}}+3\epsilon {\beta }_{0}^{2}{a}_{0}^{2}}{4\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}}t+2j{\phi }_{0}\right)}{8\left(-4{\mu }^{2}+3{\omega }_{0}^{2}+4j\mu \sqrt{{\omega }_{0}^{2}-{\mu }^{2}}+{\mathrm{\Omega }}^{2}-2j\mu \mathrm{\Omega }-4\mathrm{\Omega }\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}\right)\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}-2j\mu \mathrm{\Omega }\right)}$
$+\frac{3\epsilon {\beta }_{0}^{2}{a}_{0}{\mathrm{\Omega }}^{3}{r}^{2}}{32\left(\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}+\mathrm{\Omega }\right){\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}+2j\mu \mathrm{\Omega }\right)}^{2}}$
$×\mathrm{exp}\left(-\mu t+j\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}t+2j\mathrm{\Omega }t+j\frac{\frac{6\epsilon {\beta }_{0}^{2}{\mathrm{\Omega }}^{4}{r}^{2}}{{\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}\right)}^{2}+4{\mu }^{2}{\mathrm{\Omega }}^{2}}+3\epsilon {\beta }_{0}^{2}{a}_{0}^{2}}{8\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}}t+j{\phi }_{0}\right)$
$+\frac{3\epsilon {\beta }_{0}^{2}{a}_{0}{\mathrm{\Omega }}^{3}{r}^{2}}{32\left(-\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}+\mathrm{\Omega }\right){\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}-2j\mu \mathrm{\Omega }\right)}^{2}}$
$×\mathrm{exp}\left(-\mu t+j\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}t-2j\Omega t+j\frac{\frac{6\epsilon {\beta }_{0}^{2}{\Omega }^{4}{r}^{2}}{{\left({\omega }_{0}^{2}-{\Omega }^{2}\right)}^{2}+4{\mu }^{2}{\Omega }^{2}}+3\epsilon {\beta }_{0}^{2}{a}_{0}^{2}}{8\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}}t+j{\phi }_{0}\right)$
$+\frac{\epsilon {\alpha }_{0}^{2}{a}_{0}^{2}\left(2{\mu }^{2}-{\omega }_{0}^{2}-2j\mu \sqrt{{\omega }_{0}^{2}-{\mu }^{2}}\right)}{4\left(-4{\mu }^{2}+3{\omega }_{0}^{2}+4j\mu \sqrt{{\omega }_{0}^{2}-{\mu }^{2}}\right)}$
$×\mathrm{exp}\left(-2\mu t+2j\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}t+j\frac{\frac{6\epsilon {\beta }_{0}^{2}{\mathrm{\Omega }}^{4}{r}^{2}}{{\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}\right)}^{2}+4{\mu }^{2}{\mathrm{\Omega }}^{2}}+3\epsilon {\beta }_{0}^{2}{a}_{0}^{2}}{4\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}}t+2j{\phi }_{0}\right)$
$+\frac{\epsilon {\beta }_{0}^{2}{a}_{0}^{3}}{32\left(-3{\mu }^{2}+2{\omega }_{0}^{2}+3j\mu \sqrt{{\omega }_{0}^{2}-{\mu }^{2}}\right)}$
$×\mathrm{exp}\left(-3\mu t+3j\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}t+j\frac{\frac{18\epsilon {\beta }_{0}^{2}{\mathrm{\Omega }}^{4}{r}^{2}}{{\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}\right)}^{2}+4{\mu }^{2}{\mathrm{\Omega }}^{2}}+9\epsilon {\beta }_{0}^{2}{a}_{0}^{2}}{8\sqrt{{\omega }_{0}^{2}-{\mu }^{2}}}t+3j{\phi }_{0}\right)+cc.$

To do further simplification, the solution $x$ is:

(36)
$x={A}_{0}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\Omega }t+{\theta }_{0}\right)+{A}_{1}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\Omega }t+{\theta }_{1}\right)+{A}_{2}\mathrm{c}\mathrm{o}\mathrm{s}\left(2\mathrm{\Omega }t+{\theta }_{2}\right)+{A}_{3}\mathrm{c}\mathrm{o}\mathrm{s}\left(3\mathrm{\Omega }t+{\theta }_{3}\right)+{A}_{4}$

where:

(37)

${A}_{21}=\frac{3\epsilon {\beta }_{0}^{2}{a}_{0}^{2}}{4{\mathrm{\Omega }}^{2}\sqrt{{\left(3{\lambda }^{2}+1-{\delta }^{2}+2\sqrt{4{\lambda }^{2}-{\delta }^{2}}\right)}^{2}+{\left(\delta \sqrt{4{\lambda }^{2}-{\delta }^{2}}+\delta \right)}^{2}}}{A}_{0},$
${A}_{20}=\frac{\epsilon {\alpha }_{0}^{2}{a}_{0}^{2}\sqrt{{\left(-2{\lambda }^{2}+{\delta }^{2}\right)}^{2}+{\delta }^{2}\left(4{\lambda }^{2}-{\delta }^{2}\right)}}{4\sqrt{{\left(3{\lambda }^{2}-{\delta }^{2}\right)}^{2}+{\delta }^{2}\left(4{\lambda }^{2}-{\delta }^{2}\right)}},$
${A}_{30}=\frac{\epsilon {\beta }_{0}^{2}{a}_{0}^{3}}{4{\mathrm{\Omega }}^{2}\sqrt{{\left(8{\lambda }^{2}-3{\delta }^{2}\right)}^{2}+9{\delta }^{2}\left(4{\lambda }^{2}-{\delta }^{2}\right)}},$

${A}_{{21}^{\text{'}}}=\frac{3\epsilon {\beta }_{0}^{2}{a}_{0}^{2}}{4{\mathrm{\Omega }}^{2}\sqrt{{\left(3{\lambda }^{2}+1-{\delta }^{2}-2\sqrt{4{\lambda }^{2}-{\delta }^{2}}\right)}^{2}+{\left(\delta \sqrt{4{\lambda }^{2}-{\delta }^{2}}-\delta \right)}^{2}}}{A}_{0},$
${A}_{{12}^{\text{'}}}=\frac{3\epsilon {\beta }_{0}^{2}{a}_{0}}{8{\mathrm{\Omega }}^{2}\left(-\sqrt{4{\lambda }^{2}-{\delta }^{2}}+2\right)}{A}_{0}^{2},$

where $\lambda$ is rotational speed ratio, $\delta$ is damping rotational speed ratio, ${A}_{i}$ represents amplitude, ${\theta }_{i}$ represents phase. Eq. (36) is the complete expression of solution $x$, including steady-state and transient-state. When the time $t$ is large enough, ${e}^{-\mu t}$ decays to 0 rapidly, so the transient-state solution vanishes. Then Eq. (36) degenerates into the steady-state solution as follows:

(38)
$x={A}_{0}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\Omega }t+{\theta }_{0}\right)+{A}_{1}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\Omega }t+{\theta }_{1}\right)+{A}_{2}\mathrm{c}\mathrm{o}\mathrm{s}\left(2\mathrm{\Omega }t+{\theta }_{2}\right)+{A}_{3}\mathrm{c}\mathrm{o}\mathrm{s}\left(3\mathrm{\Omega }t+{\theta }_{3}\right)+{A}_{4}.$

Notice that, when rotational speed changes with time, Eq. (36) and Eq. (38) are not satisfied with the basic assumption of Fourier transform, which requires amplitude and phase to be time-invariant. Therefore, short-time Fourier transform (STFT) can be used to obtain an approximate frequency-domain solution.

If nonlinear scale factor $\epsilon$ equals zero, ${A}_{1}={A}_{2}={A}_{3}={A}_{4}=0$, then Eq. (38) degenerates into the traditional linear response ${A}_{0}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\Omega }t+{\theta }_{0}\right)$. The mechanism of nonlinear scale factor $\epsilon$ is fully clarified that the magnitude of $\epsilon$ reflects the degree of nonlinear influence, furthermore, it is also verifies the rationality and validity of the analytical model.

#### 3.3. Frequency response equation (FRE)

Generally, under an external excitation caused by an unbalanced force, the linear response belongs to single fundamental frequency response, but the nonlinear response may include a series of harmonic responses. For general engineering applications, working frequency FRE is often used to analyze instability problems. Considering the first two terms on the right side of Eq. (38), according to wave superposition principle, there is:

(39)
${A}_{0}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\Omega }t+{\theta }_{0}\right)+{A}_{1}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\Omega }t+{\theta }_{1}\right)={A}_{a}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\Omega }t+{\theta }_{a}\right),$

where:

(40)
${A}_{a}^{2}={A}_{0}^{2}+{A}_{1}^{2}+2{A}_{0}{A}_{1}\mathrm{c}\mathrm{o}\mathrm{s}\left({\theta }_{0}-{\theta }_{1}\right).$

Assuming ${A}_{a}=a$, substituting ${A}_{1}$ and ${\theta }_{1}$ of Eq. (37), we have:

(41)
${a}^{2}={A}_{0}^{2}+\frac{9{\epsilon }^{2}{\beta }_{0}^{4}}{16{\mathrm{\Omega }}^{4}{r}^{2}}{A}_{0}^{8}-\frac{3\epsilon {\beta }_{0}^{2}}{2{\mathrm{\Omega }}^{2}r}{A}_{0}^{5}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{0}.$

To simplify:

(42)
$-\frac{9}{16}{\epsilon }^{2}{\beta }_{0}^{4}{A}_{0}^{6}+\frac{3}{2}\epsilon {\beta }_{0}^{2}\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}\right){A}_{0}^{4}+\left[{\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}\right)}^{2}+4{\mu }^{2}{\mathrm{\Omega }}^{2}\right]{a}^{2}={\mathrm{\Omega }}^{4}{r}^{2}.$

If ${A}_{0}^{2}={a}^{2}-\frac{9{\epsilon }^{2}{\beta }_{0}^{4}}{16{\mathrm{\Omega }}^{4}{r}^{2}}{A}_{0}^{8}+\frac{3\epsilon {\beta }_{0}^{2}}{2{\mathrm{\Omega }}^{2}r}{A}_{0}^{5}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{0}\approx {a}^{2}$, the nonlinear FRE can be derived as:

(43)
$-\frac{9}{16}{\epsilon }^{2}{\beta }_{0}^{4}{a}^{6}+\frac{3}{2}\epsilon {\beta }_{0}^{2}\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}\right){a}^{4}+\left[{\left({\omega }_{0}^{2}-{\mathrm{\Omega }}^{2}\right)}^{2}+4{\mu }^{2}{\mathrm{\Omega }}^{2}\right]{a}^{2}={\mathrm{\Omega }}^{4}{r}^{2}.$

Thus, it can help to explain why the working frequency FRE cannot reflect the nonlinear characteristics accurately, or even analyze some nonlinear problems efficiently in engineering applications. The comparison between experimental result and analytical solution is given below.

#### 4. Experimental Facility and procedure

Table 1 shows the main parameters of the test rig.

Table 1. Main parameters of the test rig

 Relevant parameters Values Rotor material 40 Cr Rotor mass 10.36 kg Rotor length 698 mm Bearing inner diameter 50 mm Radial clearance 0.02 mm Clearance ratio 0.0016 Shaft diameter 50 mm Bearing length 60 mm Thrust disc diameter 90 mm Bearing span 525 mm

Fig. 3 [21] shows the schematic of the gas-bearing rotor system test rig, gas supply system, as well as the vibration signal acquisition and analysis system. The test rig has a single-stage radial turbine coaxial with a shaft driven by high-pressure air and maximum speed of 1.2×105 rpm. A screw-type air compressor provides high-pressure air with maximum pressure of 1.2 MPa and maximum mass flow of 3.90 kg/s. The airflow into the bearings is controlled by stabilizing valves and air regulators. A K-type thermocouple monitors the temperature of the supply gas. A calibrated mass flowmeter measures the mass flow rate into each bearing with an uncertainty of 0.195 kg/s. The regulators control the air pressure in proportion to the rotational speed in the speed-up process.

Fig. 3. Layout of the rotor system test rig and gas supply system: 1) air compressor, 2) gas tank, 3) filter, 4) dryer, 5) pressure gauge, 6) thermometer, 7) flowmeter, 8) pressure stabilizing valve, 9) electro-pneumatic air regulator, 10) safety shut-off valve, 11) computer, 12) data acquisition instrument, 13) eddy current sensor

Two pairs of eddy current displacement sensors, which are orthogonally positioned near test bearings, are used to measure the rotor lateral vibration amplitudes along the $X$-horizontal and $Y$-vertical axis. An eddy current key phase sensor mounted at the free end of the rotor is used to measure the phase signal. The displacement range is of 1 mm. The linearity is of 10.4 V/mm, and the frequency range is of 0-10 kHz. The sensor signals through filters and amplifiers are input to the acquisition and analysis system for real-time monitoring, storage, and analysis [22].

#### 5. Results and discussion

The theoretical calculation is accomplished with the parameter assignment as follows, 1 kg, 1×104 N/m, 20 Ns/m, 1×10-5 m, $\eta =$ 2.5×1012 N/m3, 8×104 Ns2/m2, thus ${\omega }_{0}=$ 100 rad/s, 10 rad/s. Using Matlab software, the curves of analytical solutions and working frequency FRE can be obtained. The transient time scale factor ${t}_{1}$ is defined to describe the influence from transient-state response to steady-state response. As ${t}_{1}$ increases, the transient perturbation decays to zero quickly.

Fig. 4 shows the analytical steady-state solution (0 s) and transient-state solutions (0.1, 0.2, 0.3, 0.4, 0.5 s) of Eq. (37). It is obvious from the partial enlargement that the transient-state solutions are more and more close to the steady-state solution as the transient time scale factor increases. As the transient-steady energy ratio decays to 10 %, the corresponding transient time scale factor ${t}_{1}$ is about 0.3 s. So, we take 0.3 s as the reference value. If ${t}_{1}$ is greater than 0.3 s, the transient energy is one order of magnitude lower than the steady energy, and thus can be negligible.

Fig. 4. The analytical steady-state solution (0 s) and transient-state solutions (0.1, 0.2, 0.3, 0.4, 0.5 s)

Fig. 5. The spectrum transformed by FFT

Fig. 6. The frequency-response curve of experimental result

If rotational speed is constant, the amplitudes and phases in Eq. (39) degenerate to constant, which satisfies the time-invariant condition of Fourier transform. As can be seen from Fig. 5 transformed by FFT, the magnitude of the third harmonic generation is three orders of magnitude lower than that of the fundamental frequency. The fundamental frequency is about 8000 rpm, which is a bit consistent with experimental result of speed-up process shown in Fig. 6.

To further verify the nonlinear response characteristics, we put three frequency-response curves, including experimental result, steady-state solution transformed by STFT, working frequency FRE, on the same graph respectively, and intercept the portion of fundamental frequency for comparing, which is shown in Fig. 7. Obviously, all of the curves take on resemblance in the linear vibration features caused by ${A}_{0}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\Omega }t+{\theta }_{0}\right)$. It is worth noting that, the AB segment of experimental result and the CD segment of steady-state solution both present nonlinear vibration characteristics caused by ${A}_{1}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\Omega }t+{\theta }_{1}\right)$. But for the curve of working frequency FRE, there is only one peak point $E$, which can be considered as an expression of the nonlinear feature, so that it is incomplete and inaccurate to analyze the nonlinear characteristics by just using the working frequency FRE.

Fig. 7. Comparison of frequency-response curves

The time-amplitude-frequency spectrum of rotor speed-up process is shown in Fig. 8, and especially in the area A, the whirl starts to appear at the rotational speed of 12487 rpm under the effect of gas film, which due to the fluid-solid coupling mechanism mainly caused by gas bearing used in this experiment. It needs to be noticed that, based on the large deformation and perturbation assumptions, the analytical model in this paper does not consider the subharmonic effect caused by gas film. Additional researches on whirl effect of bearing will be addressed in subsequent work.

Fig. 8. Time-amplitude-frequency spectrum of rotor speed-up process

#### 6. Conclusions

In this paper, with respect to reveal the mechanism of dynamic behaviors and nonlinear responses of a rotor system, an analytical model is constructed, and the analytical solutions are derived via multi-scale method. The conclusions are as follows.

1) By defining the nonlinear scale factor $\epsilon$, and introducing nonlinear stiffness and nonlinear damping, a nonlinear analytical model is presented, which aims to reflect the linear and nonlinear interaction.

2) The analytical solutions of steady-state and transient-state are derived via multi-scale method, and the nonlinear natural frequency and the working frequency FRE are obtained.

3) The mechanism of nonlinear scale factor $\epsilon$ is clarified, which can be used to reflect the degree of influence on nonlinear effects. The rationality and validity of the analytical model are also verified.

4) The transient time scale factor ${t}_{1}$ is defined to reflect transient-state influence on steady-state, which can be obtained when transient-steady energy ratio decays to 10 %.

5) The steady-state solution well reflects the nonlinear vibration characteristics of the fundamental frequency, which is consistent with the experimental result, rather than that of the working frequency FRE, in which the nonlinear characteristics is represented as one point. Therefore, the rationality of the analytical solutions is verified.

#### Acknowledgements

The work is supported by National Science and Technology projects (Grant No. 2012BAA11B02), and the support is gratefully acknowledged.

#### References

1. Meng G. Retrospect and prospect to the research on rotordynamics. Journal of Vibration Engineering, Vol. 15, Issue 1, 2002, p. 1-9, (in Chinese). [Search CrossRef]
2. Huang W., Wu X., Jiao Y., et al. Review of nonlinear rotor dynamics. Journal of Vibration Engineering, Vol. 13, Issue 4, 2000, p. 497-509, (in Chinese). [Search CrossRef]
3. Yang J., Yang S., Chen C., et al. Research on sliding bearings and rotor system stability. Journal of Aerospace Power, Vol. 23, Issue 8, 2008, p. 1420-1426, (in Chinese). [Search CrossRef]
4. Lund J. W. The stability of an elastic rotor in journal bearings with flexible, damped supports. Journal of Applied Mechanics, Vol. 32, Issue 4, 1965, p. 911-920. [Publisher]
5. Glienicke J. Spring Damping Constants of Turbine Bearings and the Effect on the Vibration Behavior of Single Rotor. Dissertation, TH, Karlsruhe, 1966, (in German). [Search CrossRef]
6. Alford J. S. Protecting turbomachinery from self-excited rotor whirl. Journal of Engineering for Power, Vol. 87, Issue 4, 1965, p. 333-343. [Publisher]
7. Tanaka M., Hori Y. Stability characteristics of floating bush bearings. Journal of Lubrication Technology, Vol. 94, Issue 3, 1972, p. 248-256. [Publisher]
8. Ehrich F. F. High order subharmonic response of high speed rotors in bearing clearance. Journal of Vibration, Acoustics, Stress, and Reliability in Design, Vol. 110, Issue 1, 1988, p. 9-16. [Publisher]
9. Childs D. W. Dynamic analysis of turbulent annular seals based on hirs’ lubrication equation. Journal of Lubrication Technology, Vol. 105, Issue 3, 1983, p. 429-436. [Publisher]
10. Bently D. E., Muszynska A. Role of circumferential flow in the stability of fluid-handling machine rotors. Rotordynamic Instability Problems in High-Performance Turbomachinery (5th), 1988, p. 415-430. [Search CrossRef]
11. Muszynska A., Bently D. E. Anti-swirl arrangements prevent rotor/seal instability. Journal of Vibration, Acoustics, Stress, and Reliability in Design, Vol. 111, Issue 2, 1989, p. 156-162. [Publisher]
12. Ravikovich Yu. A. The Construction and Design of Sliding Bearing Units. Moscow, 1995, (in Russian). [Search CrossRef]
13. Wu M., Yang J. The research on single-disc rotor nonlinear vibration model and mechanism. Mechatronics and Information Technology, PTS 1 and 2, Vols. 2-3, 2012, p. 954-959. [Search CrossRef]
14. Liu B., Peng J. Nonlinear Dynamics. Higher Education Press, Beijing, 2004, (in Chinese). [Search CrossRef]
15. Khalil H. Nonlinear Systems. Electronic Industry Press, Beijing, 2005, (in Chinese). [Search CrossRef]
16. Zhou J., Zhu Y. Nonlinear Vibration. Xi’an Jiaotong University Press, Xi’an, 1998, (in Chinese). [Search CrossRef]
17. Hu H. Nonlinear Vibration. Aviation Industry Press, Beijing, 2000, (in Chinese). [Search CrossRef]
18. Chen Y. Modern Analytical Method in Nonlinear Dynamics. Science Press, Beijing, 1992, (in Chinese). [Search CrossRef]
19. Weng B., Li Y., Xu P. Engineering Nonlinear Vibration. Science Press, Beijing, 2007, (in Chinese). [Search CrossRef]
20. Meng G., Xue Z. Study on nonlinear duffing characteristics of flexible rotor with SFDB. Journal of Aerospace Power, Vol. 4, Issue 2, 1989, p. 173-197, (in Chinese). [Search CrossRef]
21. Han D., Yang J., Chen C., et al. Experimental research on the dynamic characteristics of gas-hybrid bearing-flexible rotor system. Journal of Vibroengineering, Vol. 16, Issue 5, 2014, p. 2363-2374. [Search CrossRef]
22. Han D., Yang J., Zhao C., et al. An experimental study on vibration characteristics for gas bearing-rotor system. Journal of Vibration Engineering, Vol. 25, Issue 6, 2012, p. 680-685, (in Chinese). [Search CrossRef]

#### Cited By

 Journal of Vibration Engineering & Technologies Bin Wang, Yongtao Sun, Qian Ding 2020 Journal of Vibroengineering Min Wu, Shengbo Yang, Wen Bao, Jinfu Yang 2017