Analytical model for nonlinear vibration of flexible rotor system

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


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.

Model formulation
As shown in Fig. 1, the rotor system is considered as a single-span rotor supported by force , which can be divided into stiffness force and damping force. Noticing that, the force ′, which is perpendicular to 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.

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: 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 , namely the elastic force, is expressed as: where is the coefficient of linear elastic deformation, namely stiffness, is the coefficient of nonlinear elastic deformation. If = 0, Eq. (3) satisfies the linear Hooke's law, namely linear spring characteristics. If > 0, increases along with deformation increasing, namely hard spring characteristics. On the contrary, if < 0, decreases along with deformation increasing, namely soft spring characteristics [16].

Nonlinear damping
Generally, the damping dissipation process of non-conservative system [17] can be given as: Thus, the damping force is described as a generalized form: 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 = 2, Eq. (5) represents nonlinear viscous damping force, which is applied to high-speed movement [18]. Thus, the damping force is expressed as: where represents linear damping, is the coefficient of nonlinear damping.

Analytical model
Considering the single-span single-disc rotor system shown in Fig. 1, the parameters are as follows, for disc mass, for eccentricity, Ω 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 -axis direction can be described as: where: For briefness and clarity, some substitutions are given as: where is linear damping coefficient, is natural frequency, is nonlinear stiffness amplification coefficient, is nonlinear damping amplification coefficient, 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 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 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]: where is excitation force, is stiffness, is damping, is fluid circumferential average velocity ratio, is inertial mass, is rotational speed, is mass, is eccentricity, in Eq. (11)(12)(13). 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.

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, is introduced and defined as: Then, a nonlinear response can be written as a function of various time scales: where 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 multi-scale method, we suppose: where: Thus: Then, according to the power orders of , there are: : where ( = 0, 1) is the partial differential operator. The solution of Eq. (20) is: where is an unknown complex function, the overline represents the conjugate function. Take Eq. (22) into Eq. (21), we have: where represents the conjugate functions of each term before. The coefficient of first term must be zero to avoid duration term. At any point without acceleration, = 0, so = 1. There is: Assuming: where and are the functions of . According to the real part and the imaginary part respectively, can be obtained: where , are the integration constants. Take Eq. (26) and Eq. (22) into Eq. (23), can be solved, then the solution is: 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, , , and have an impact on the natural frequency but , does not.

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 multi-scale method, according to the power orders of , there are: : : The solution of Eq. (30) is: where is an unknown complex function, the overline represents the conjugate function. Take Eq. (32) into Eq. (31), we have: Similarly, to avoid duration term, the coefficient of first term must be zero. At any point without acceleration, = 0, so = 1. There is: Take Eq. (25) into Eq. (34), according to the real part and the imaginary part respectively, can be obtained.
where , are the integration constants. Take Eq. (35) and Eq. (32) into Eq. (33), can be solved, then the solution is: where is rotational speed ratio, is damping rotational speed ratio, represents amplitude, represents phase. Eq. (36) is the complete expression of solution , including steady-state and transient-state. When the time is large enough, decays to 0 rapidly, so the transient-state solution vanishes. Then Eq. (36) degenerates into the steady-state solution as follows: = cos(Ω + ) + cos(Ω + ) + cos(2Ω + ) + cos(3Ω + ) + . (38) 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 equals zero, = = = = 0, then Eq. (38) degenerates into the traditional linear response cos(Ω + ). The mechanism of nonlinear scale factor is fully clarified that the magnitude of reflects the degree of nonlinear influence, furthermore, it is also verifies the rationality and validity of the analytical model.

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 = , substituting and of Eq. (37), we have: To simplify: If = − + cos ≈ , 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. Table 1 shows the main parameters of the test rig.  [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×10 5 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

Experimental Facility and procedure
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 -horizontal and -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].

Results and discussion
The theoretical calculation is accomplished with the parameter assignment as follows, = 1 kg, = 1×10 4 N/m, = 20 Ns/m, = 1×10 -5 m, = 2.5×10 12 N/m 3 , = 8×10 4 Ns 2 /m 2 , thus = 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 is defined to describe the influence from transient-state response to steady-state response. As increases, the transient perturbation decays to zero quickly.  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 cos(Ω + ). 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 cos(Ω + ). But for the curve of working frequency FRE, there is only one peak point , 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. 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. Daren Yu and Jinfu Yang guide this work, organize the discussion and the analysis.

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