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 multiscale method, the analytical solutions of steadystate and transientstate 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 transientstate influence on steadystate solution. The experimental result also verifies the rationality and validity of the analytical model and the analytical solutions.
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 highspeed miniature equipments, the primary perspective is to enhance rotational speed, system power, thrustweight 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, quasiperiod, 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 selfexcited 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 highspeed 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 highspeed 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 multiscale method, and then verified with experiment result.
2. Model formulation
As shown in Fig. 1, the rotor system is considered as a singlespan 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. 1Dynamic 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 onedimensional system [14] can be expressed as:
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:
The stiffness force ${f}_{k}$, namely the elastic force, is expressed as:
where $k$ is the coefficient of linear elastic deformation, namely stiffness, $\eta $ is the coefficient of nonlinear elastic deformation. If $\eta =$0, Eq. (3) satisfies the linear Hooke’s law, namely linear spring characteristics. If $\eta >$0, ${f}_{k}$ increases along with deformation increasing, namely hard spring characteristics. On the contrary, if $\eta <$0, ${f}_{k}$ decreases along with deformation increasing, namely soft spring characteristics [16].
2.2. Nonlinear damping
Generally, the damping dissipation process of nonconservative system [17] can be given as:
Thus, the damping force is described as a generalized form:
If $n=$0, Eq. (5) represents dry friction force. If $n=$1, Eq. (5) represents linear viscous damping force, which is applied to lowspeed movement in fluid. If $n=$ 2, Eq. (5) represents nonlinear viscous damping force, which is applied to highspeed movement [18]. Thus, the damping force ${f}_{c}$ is expressed as:
where $c$ represents linear damping, $\gamma $ is the coefficient of nonlinear damping.
2.3. Analytical model
Considering the singlespan singledisc 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:
where:
For briefness and clarity, some substitutions are given as:
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:
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 perioddoubling bifurcation, Hopf bifurcation, period3 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]:
On this basis, if ignore inertial force, and substitute coefficients for simple letters, Eq. (11) can be simplified as Alford model [6]:
Apparently, if the cross damping force is not considered, Eq. (11) can be simplified as Ravikovich model [12]:
$m\ddot{y}+2c\dot{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. (1113).
Fig. 2Schematic 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. Analytical solutions and analysis
3.1. Free vibration mode
Multiscale method can calculate both of the steadystate response and the transientstate response. In order to represent different time scales, ${T}_{n}$ is introduced and defined as:
Then, a nonlinear response can be written as a function of various time scales:
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:
Based on multiscale method, we suppose:
where:
Thus:
$\ddot{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:
where ${D}_{n}$ ($n=$0, 1) is the partial differential operator. The solution of Eq. (20) is:
where $C$ is an unknown complex function, the overline represents the conjugate function. Take Eq. (22) into Eq. (21), we have:
${\alpha}_{0}^{2}C\stackrel{}{C}{r}_{1}{\stackrel{}{r}}_{1}{e}^{2\mu {T}_{0}}{\alpha}_{0}^{2}{C}^{2}{r}_{1}^{2}{e}^{2{r}_{1}{T}_{0}}{\beta}_{0}^{2}{C}^{3}{e}^{3{r}_{1}{T}_{0}}+cc,$
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:
Assuming:
where $a$ and $\phi $ are the functions of ${T}_{1}$. According to the real part and the imaginary part respectively, $C$ can be obtained:
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:
$+\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)}{e}^{2\mu t+2j\sqrt{{\omega}_{0}^{2}{\mu}^{2}}t+j\frac{3\epsilon {\beta}_{0}^{2}{a}_{0}^{2}}{4\sqrt{{\omega}_{0}^{2}{\mu}^{2}}}t+2j{\phi}_{0}}$
$+\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)}{e}^{3\mu t+3j\sqrt{{\omega}_{0}^{2}{\mu}^{2}}t+j\frac{9\epsilon {\beta}_{0}^{2}{a}_{0}^{2}}{8\sqrt{{\omega}_{0}^{2}{\mu}^{2}}}t+3j{\phi}_{0}}+cc.$
Thus, the natural frequency is:
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:
Similarly, using multiscale method, according to the power orders of $\epsilon $, there are:
The solution of Eq. (30) is:
${r}_{1}=\mu +j\sqrt{{\omega}_{0}^{2}{\mu}^{2}},Q=\frac{1}{2}\frac{{\mathrm{\Omega}}^{2}r}{{\omega}_{0}^{2}{\mathrm{\Omega}}^{2}+2j\mu \mathrm{\Omega}},$
where $C$ is an unknown complex function, the overline represents the conjugate function. Take Eq. (32) into Eq. (31), we have:
${\alpha}_{0}^{2}Q\stackrel{}{Q}{\mathrm{\Omega}}^{2}{\alpha}_{0}^{2}C\stackrel{}{C}{r}_{1}{\stackrel{}{r}}_{1}{e}^{2\mu {T}_{0}}3{\beta}_{0}^{2}{Q}^{2}\stackrel{}{Q}{e}^{j\mathrm{\Omega}{T}_{0}}+{\alpha}_{0}^{2}{Q}^{2}{\mathrm{\Omega}}^{2}{e}^{2j\mathrm{\Omega}{T}_{0}}{\beta}_{0}^{2}{Q}^{3}{e}^{3j\mathrm{\Omega}{T}_{0}}$
$2j{\alpha}_{0}^{2}C{r}_{1}Q\mathrm{\Omega}{e}^{{r}_{1}{T}_{0}+j\mathrm{\Omega}{T}_{0}}+2j{\alpha}_{0}^{2}C{r}_{1}\stackrel{}{Q}\mathrm{\Omega}{e}^{{r}_{1}{T}_{0}j\mathrm{\Omega}{T}_{0}}6{\beta}_{0}^{2}C\stackrel{}{C}Q{e}^{2\mu {T}_{0}+j\mathrm{\Omega}{T}_{0}}$
$3{\beta}_{0}^{2}{C}^{2}Q{e}^{2{r}_{1}{T}_{0}+j\mathrm{\Omega}{T}_{0}}3{\beta}_{0}^{2}{C}^{2}\stackrel{}{Q}{e}^{2{r}_{1}{T}_{0}j\mathrm{\Omega}{T}_{0}}3{\beta}_{0}^{2}C{Q}^{2}{e}^{{r}_{1}{T}_{0}+2j\mathrm{\Omega}{T}_{0}}$
$3{\beta}_{0}^{2}C{\stackrel{}{Q}}^{2}{e}^{{r}_{1}{T}_{0}2j\mathrm{\Omega}{T}_{0}}{\alpha}_{0}^{2}{C}^{2}{r}_{1}^{2}{e}^{2{r}_{1}{T}_{0}}{\beta}_{0}^{2}{C}^{3}{e}^{3{r}_{1}{T}_{0}}+cc.$
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:
Take Eq. (25) into Eq. (34), according to the real part and the imaginary part respectively, $C$ can be obtained.
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:
$+\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}$
$\times \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)}$
$\times \mathrm{exp}\left(\mu t+j\sqrt{{\omega}_{0}^{2}{\mu}^{2}}tj\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}}tj\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}}$
$\times \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}}$
$\times \mathrm{exp}\left(\mu t+j\sqrt{{\omega}_{0}^{2}{\mu}^{2}}t2j\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)}$
$\times \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)$
To do further simplification, the solution $x$ is:
$+{a}_{0}{e}^{\mu t}\mathrm{c}\mathrm{o}\mathrm{s}\left(\phi t+{\phi}_{0}\right)+{A}_{01}{e}^{2\mu t}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\Omega}t+{A}_{40}{e}^{2\mu t}$
$+{A}_{11}{e}^{\mu t}\mathrm{c}\mathrm{o}\mathrm{s}\left(\left(\phi +\mathrm{\Omega}\right)t+{\theta}_{11}\right)+{A}_{11\text{'}}{e}^{\mu t}\mathrm{c}\mathrm{o}\mathrm{s}\left(\left(\phi \mathrm{\Omega}\right)t+{\theta}_{11\text{'}}\right)$
$+{A}_{21}{e}^{2\mu t}\mathrm{c}\mathrm{o}\mathrm{s}\left(\left(2\phi +\mathrm{\Omega}\right)t+{\theta}_{21}\right)+{A}_{21\text{'}}{e}^{2\mu t}\mathrm{c}\mathrm{o}\mathrm{s}\left(\left(2\phi \mathrm{\Omega}\right)t+{\theta}_{21\text{'}}\right)$
$+{A}_{12}{e}^{\mu t}\mathrm{c}\mathrm{o}\mathrm{s}\left(\left(\phi +2\mathrm{\Omega}\right)t+{\theta}_{12}\right)+{A}_{12\text{'}}{e}^{\mu t}\mathrm{c}\mathrm{o}\mathrm{s}\left(\left(\phi 2\mathrm{\Omega}\right)t+{\theta}_{12\text{'}}\right)$
$+{A}_{20}{e}^{2\mu t}\mathrm{c}\mathrm{o}\mathrm{s}\left(2\phi t+{\theta}_{20}\right)+{A}_{30}{e}^{3\mu t}\mathrm{c}\mathrm{o}\mathrm{s}\left(3\phi t+{\theta}_{30}\right),$
where:
${A}_{3}=\frac{\epsilon {\beta}_{0}^{2}}{4{\mathrm{\Omega}}^{2}\sqrt{{\left({\lambda}^{2}9\right)}^{2}+9{\delta}^{2}}}{A}_{0}^{3},{A}_{4}=\frac{\epsilon {\alpha}_{0}^{2}}{2{\lambda}^{2}}{A}_{0}^{2},{A}_{01}=\frac{3\epsilon {\beta}_{0}^{2}{a}_{0}^{2}}{2{\mathrm{\Omega}}^{2}r}{A}_{0}^{2},$
${A}_{11}=\frac{\epsilon {\alpha}_{0}^{2}{a}_{0}\lambda}{\sqrt{4{\lambda}^{2}{\delta}^{2}}+1}{A}_{0},{A}_{12}=\frac{3\epsilon {\beta}_{0}^{2}{a}_{0}}{8{\mathrm{\Omega}}^{2}\left(\sqrt{4{\lambda}^{2}{\delta}^{2}}+2\right)}{A}_{0}^{2},$
${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}_{40}=\frac{1}{2}\epsilon {\alpha}_{0}^{2}{a}_{0}^{2},{\theta}_{11}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{g}\frac{\delta}{\sqrt{4{\lambda}^{2}{\delta}^{2}}}+{\theta}_{0}+{\phi}_{0},{A}_{{11}^{\text{'}}}=\frac{\epsilon {\alpha}_{0}^{2}{a}_{0}\lambda}{\sqrt{4{\lambda}^{2}{\delta}^{2}}1}{A}_{0},$
${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},$
$\lambda =\frac{{\omega}_{0}}{\mathrm{\Omega}},\delta =\frac{2\mu}{\mathrm{\Omega}},\phi =\frac{1}{2}\mathrm{\Omega}\sqrt{4{\lambda}^{2}{\delta}^{2}}+\frac{3\epsilon {\beta}_{0}^{2}\left(2{A}_{0}^{2}+{a}_{0}^{2}\right)}{4\mathrm{\Omega}\sqrt{4{\lambda}^{2}{\delta}^{2}}},$
${\theta}_{0}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{g}\frac{\delta}{{\lambda}^{2}1},{\theta}_{1}=2{\theta}_{0},{\theta}_{2}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{g}\frac{2\delta}{{\lambda}^{2}4}+2{\theta}_{0},$
${\theta}_{3}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{g}\frac{3\delta}{{\lambda}^{2}9}+3{\theta}_{0},{\theta}_{21}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{g}\frac{\delta \sqrt{4{\lambda}^{2}{\delta}^{2}}+\delta}{3{\lambda}^{2}+1{\delta}^{2}+2\sqrt{4{\lambda}^{2}{\delta}^{2}}}+{\theta}_{0}+2{\phi}_{0},$
${{\theta}_{12}=2{\theta}_{0}+{\phi}_{0},\theta}_{20}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{g}\frac{\delta \sqrt{4{\lambda}^{2}{\delta}^{2}}}{2{\lambda}^{2}+{\delta}^{2}}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{g}\frac{\delta \sqrt{4{\lambda}^{2}{\delta}^{2}}}{3{\lambda}^{2}{\delta}^{2}}+2{\phi}_{0},$
${\theta}_{30}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{g}\frac{3\delta \sqrt{4{\lambda}^{2}{\delta}^{2}}}{8{\lambda}^{2}3{\delta}^{2}}+3{\phi}_{0},{\theta}_{{11}^{\text{'}}}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{g}\frac{\delta}{\sqrt{4{\lambda}^{2}{\delta}^{2}}}{\theta}_{0}+{\phi}_{0},$
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 steadystate and transientstate. When the time $t$ is large enough, ${e}^{\mu t}$ decays to 0 rapidly, so the transientstate solution vanishes. Then Eq. (36) degenerates into the steadystate solution as follows:
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 timeinvariant. Therefore, shorttime Fourier transform (STFT) can be used to obtain an approximate frequencydomain 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:
where:
Assuming ${A}_{a}=a$, substituting ${A}_{1}$ and ${\theta}_{1}$ of Eq. (37), we have:
To simplify:
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:
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 1Main 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 gasbearing rotor system test rig, gas supply system, as well as the vibration signal acquisition and analysis system. The test rig has a singlestage radial turbine coaxial with a shaft driven by highpressure air and maximum speed of 1.2×10^{5} rpm. A screwtype air compressor provides highpressure 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 Ktype 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 speedup process.
Fig. 3Layout 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) electropneumatic air regulator, 10) safety shutoff 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 010 kHz. The sensor signals through filters and amplifiers are input to the acquisition and analysis system for realtime monitoring, storage, and analysis [22].
5. Results and discussion
The theoretical calculation is accomplished with the parameter assignment as follows, $m=$1 kg, $k=$1×10^{4} N/m, $c=$20 Ns/m, $r=$1×10^{5} m, $\eta =$ 2.5×10^{12} N/m^{3}, $\gamma =$8×10^{4} Ns^{2}/m^{2}, thus ${\omega}_{0}=$ 100 rad/s, $\mu =$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 transientstate response to steadystate response. As ${t}_{1}$ increases, the transient perturbation decays to zero quickly.
Fig. 4 shows the analytical steadystate solution (${t}_{1}=$0 s) and transientstate solutions (${t}_{1}=$0.1, 0.2, 0.3, 0.4, 0.5 s) of Eq. (37). It is obvious from the partial enlargement that the transientstate solutions are more and more close to the steadystate solution as the transient time scale factor increases. As the transientsteady 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. 4The analytical steadystate solution (t1= 0 s) and transientstate solutions (t1= 0.1, 0.2, 0.3, 0.4, 0.5 s)
Fig. 5The spectrum transformed by FFT
Fig. 6The frequencyresponse curve of experimental result
If rotational speed is constant, the amplitudes and phases in Eq. (39) degenerate to constant, which satisfies the timeinvariant 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 speedup process shown in Fig. 6.
To further verify the nonlinear response characteristics, we put three frequencyresponse curves, including experimental result, steadystate 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 steadystate 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. 7Comparison of frequencyresponse curves
The timeamplitudefrequency spectrum of rotor speedup 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 fluidsolid 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. 8Timeamplitudefrequency spectrum of rotor speedup 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 multiscale 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 steadystate and transientstate are derived via multiscale 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 transientstate influence on steadystate, which can be obtained when transientsteady energy ratio decays to 10 %.
5) The steadystate 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.
References

Meng G. Retrospect and prospect to the research on rotordynamics. Journal of Vibration Engineering, Vol. 15, Issue 1, 2002, p. 19, (in Chinese).

Huang W., Wu X., Jiao Y., et al. Review of nonlinear rotor dynamics. Journal of Vibration Engineering, Vol. 13, Issue 4, 2000, p. 497509, (in Chinese).

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

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

Glienicke J. Spring Damping Constants of Turbine Bearings and the Effect on the Vibration Behavior of Single Rotor. Dissertation, TH, Karlsruhe, 1966, (in German).

Alford J. S. Protecting turbomachinery from selfexcited rotor whirl. Journal of Engineering for Power, Vol. 87, Issue 4, 1965, p. 333343.

Tanaka M., Hori Y. Stability characteristics of floating bush bearings. Journal of Lubrication Technology, Vol. 94, Issue 3, 1972, p. 248256.

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

Childs D. W. Dynamic analysis of turbulent annular seals based on hirs’ lubrication equation. Journal of Lubrication Technology, Vol. 105, Issue 3, 1983, p. 429436.

Bently D. E., Muszynska A. Role of circumferential flow in the stability of fluidhandling machine rotors. Rotordynamic Instability Problems in HighPerformance Turbomachinery (5th), 1988, p. 415430.

Muszynska A., Bently D. E. Antiswirl arrangements prevent rotor/seal instability. Journal of Vibration, Acoustics, Stress, and Reliability in Design, Vol. 111, Issue 2, 1989, p. 156162.

Ravikovich Yu. A. The Construction and Design of Sliding Bearing Units. Moscow, 1995, (in Russian).

Wu M., Yang J. The research on singledisc rotor nonlinear vibration model and mechanism. Mechatronics and Information Technology, PTS 1 and 2, Vols. 23, 2012, p. 954959.

Liu B., Peng J. Nonlinear Dynamics. Higher Education Press, Beijing, 2004, (in Chinese).

Khalil H. Nonlinear Systems. Electronic Industry Press, Beijing, 2005, (in Chinese).

Zhou J., Zhu Y. Nonlinear Vibration. Xi’an Jiaotong University Press, Xi’an, 1998, (in Chinese).

Hu H. Nonlinear Vibration. Aviation Industry Press, Beijing, 2000, (in Chinese).

Chen Y. Modern Analytical Method in Nonlinear Dynamics. Science Press, Beijing, 1992, (in Chinese).

Weng B., Li Y., Xu P. Engineering Nonlinear Vibration. Science Press, Beijing, 2007, (in Chinese).

Meng G., Xue Z. Study on nonlinear duffing characteristics of flexible rotor with SFDB. Journal of Aerospace Power, Vol. 4, Issue 2, 1989, p. 173197, (in Chinese).

Han D., Yang J., Chen C., et al. Experimental research on the dynamic characteristics of gashybrid bearingflexible rotor system. Journal of Vibroengineering, Vol. 16, Issue 5, 2014, p. 23632374.

Han D., Yang J., Zhao C., et al. An experimental study on vibration characteristics for gas bearingrotor system. Journal of Vibration Engineering, Vol. 25, Issue 6, 2012, p. 680685, (in Chinese).
About this article
The work is supported by National Science and Technology projects (Grant No. 2012BAA11B02), and the support is gratefully acknowledged.
Min Wu wrote this paper and take the charge the theoretical work and analysis. Shengbo Yang drew the pictures and a part of the analysis. Dongjiang Han take charge of the experimental work. Daren Yu and Jinfu Yang guide this work, organize the discussion and the analysis.