1 Introduction

The field of rotor dynamics is concerned with the research of dynamic and stability characteristics of rotating machinery, and it plays an important role in improving the safety and performance of the entire systems. The modeling of rotor systems and dynamic characteristics analysis are the fundamental research content in the field of rotor dynamics. The rotor system is usually supported by bearings and influenced by internal phenomena that rotor rotates around a single axis. Recently the multi-disk and multi-span rotor system is becoming an important field for rotor dynamics research (Chen, 2009).

Many studies have been conducted on rotor modeling and dynamic characteristics analysis. The Muszynska model highlights the seal force nonlinear characteristics with clear physical meaning (Muszynska and Bently, 1990). Al-Nahwi et al. (2003) analyzed the principle and interaction of steam excitation on the Jeffcott rotor system combining it with the Moore-Greitzer flow field model. Luo et al. (2007) built a periodical time-variable high-dimensional dynamic rotor system based on the rotor’s finite element model and investigated the stability of the system. Cheng et al. (2008) studied the nonlinear dynamic behaviors of a rotor/bearing/seal coupled system with Muszynska’s seal forces and Capone’s oil-film forces. The influence of the rotation speed, seal clearance, and eccentricity of the rotor were analyzed. de Castro et al. (2008) modeled a flexible rotor with a central disk under unbalanced excitation and validated a complete nonlinear solution to simulate the fluid-induced instability during run-up and run-down. Wang et al. (2009) established a nonlinear mathematical model for orbital motion of the rotor under the influence of leakage flow through an interlocking seal and used the fourth-order Runge-Kutta method to solve it. Particular attention was placed on the serpentine flow path by spatially separating the aerodynamic force on the rotor surface into two parts, e.g., the seal clearance and the cavity volume (Wang et al., 2009). Okabe and Cavalca (2009) developed an analytical model of a tilting pad bearing based on the short bearing assumption with the turbulence effect included. They found that the bearing model with turbulent flow effect generated higher hydrodynamic forces when compared to the one without this effect, which highlights the importance of considering such phenomenon during the analysis of high speed hydrodynamic bearings. Li et al. (2011a; 2011b) applied the Hamilton principle and the finite element method (FEM) to construct a novel nonlinear model of a rotor/bearing/seal system. The dynamic behavior of the system is illustrated by bifurcation diagrams, large Lyapunov exponents, phase trajectory diagrams, and Poincare maps. Li et al. (2012) presented the effects of journal misalignment on the transient flow of a finite grooved journal bearing using a new 3D computational fluid dynamics analysis method. Based on the FEM and the Lagrange equation, Zhou et al. (2014) proposed a novel nonlinear model of a double disc rotor-seal system, including the coupled effects of the gravity force of the discs, Muszynska’s nonlinear seal fluid dynamic force, and the mass eccentricity of the discs. Other researchers focused on the dynamic analysis of rotor systems supported by gas bearings. Půst and Kozánek (2007) calculated the dynamic characteristics of bearings at different revolutions, which took into account the inertia properties of tilting pads. Ertas et al. (2010) tested a rotor system using a 70-mm diameter damped gas bearing reaching ultra-high speeds of 50000 r/min and experimentally evaluated the ability of the damped gas bearing to withstand large rotor excursions. Rashidi et al. (2010a; 2010b) studied the preload effect on the behavior of a rigid rotor supported by gas-lubricated noncircular journal bearings. Results of this study revealed how the complex dynamic behavior of two types of noncircular bearing systems, comprising periodic, KT-periodic, and quasi-periodic responses of the rotor center, varies with changes in the preload value.

While the above studies were very significant for the rotor design and improvement for academic researches, they simplified the whole rotor/bearing/seal system as a rotor supported by two bearings without taking the concept of multi-rotor into consideration. The simplification of the traditional rotor/bearing/ seal system may lead to some errors in numerical simulation compared with the multi-rotor system, but it would be more difficult to analyze the dynamic orbits of all the geometric centers of rotors.

Recently, the multi-disk and multi-span rotor systems have been the subject of study in the domain of rotor dynamics. Ding and Leung (2005) built a test rotor rig consisting of two flexibly coupled shafts, each supported at the ends by two hydrodynamic bearings. A mathematical model for the test rig was developed and the non-stationary processes of the system were analyzed numerically. Wu and Jing (2008) established a dynamic model for a two-span bearing-rotor system with nonlinear oil film that was close to a practical situation. The influence of the seal force in the model was not considered. The bifurcation and chaos behaviors of the system under different work conditions were studied. Chang-Jian (2010) investigated the vibration characteristics of two rotors equipped with long journal bearings placed at both ends. Different cross sections, shaft lengths, bearing masses, and different bearing approximations were presented to analyze and discuss the differences of the dynamic responses. Luo et al. (2012) set up a nonlinear dynamic model of a two-span rotor bearing system with two cracks in the shafts and analyzed the characteristics of this system. Fei et al. (2013) discussed a typical structure with two rotor shafts and illustrated the procedure for obtaining the coupling motion equations of the subsystems using the FEM.

In this paper, a dynamic model of a two-disk and two-span rotor/bearing/seal system is created based on the Hamilton principle and the FEM. The rotor system consists of two flexible coupled shafts, each of which is supported by two short journal bearings and attached with one disk. The motion equations of the rotor system were solved with the fourth-order Runge-Kutta method. The nonlinear dynamic behavior of this two-span flexible rotor system was analyzed by the bifurcation diagrams, time-history diagrams, phase trajectories, and Poincare maps. The numerical results indicate that rotational speed, the nonlinear seal force, and the oil-film force have an effect on the stability of the rotor system.

2 A nonlinear model of a rotor system based on the Hamilton principle and FEM

The model of a two-span rotor system supported by four journal bearings is shown in Fig. 1. Two identical Jeffcott rotors are connected with a flexible coupler, each of which is supported by journal bearings on both ends of the rotor. In the following description, the oil-film force model is established by the short bearing theory and the nonlinear seal force is described by the Muszynska model. The expressions of the kinetic energy and strain energy are derived by a projection angle method. In Fig. 1, d is the diameter of the shaft, b d and d d are the breadth and diameter of the disks, and L and l c are the length of one span rotor and coupling, respectively.

Fig. 1
figure 1

The physical model of a two-disk and two-span rotor system

Full size image

2.1 Nonlinear oil-film force model

The short bearing unsteady oil-film force model is

((1))

where F ox and F oy are the bearing oil-film force in the x direction and y direction, respectively; is the Sommerfeld coefficient, μ is the sectional shear correction factor, Ω is the rotational speed of the rotor, X=x/c o and Y=y/c o are defined as the dimensionless translation displacements at the bearing, c o, R o, and L o are the average clearance, radius, and length of the bearing, respectively. is dimensionless time. K o and C o can be written as

((2))

where

2.2 Nonlinear seal force model

The Muszynska model can be used to describe the nonlinear characteristic of a steam excitation force very well. The expression of this model can be described as follows (Muszynska and Bently, 1990):

((3))

where F sx and F sy are the seal force in the x direction and y direction, respectively; K f, D f, and τ f are the nonlinear functions of the translation displacements x and y at the disk. They can be written as

((4))

where K o, D o, and m f can be obtained from the Childs equation (Childs, 1983).

2.3 Finite element model of a rotor system

The finite element model of a rotor system is established by applying the FEM as Fig. 2, where m d is the disk mass. There are six nodes in total, which are distributed on the bearings, disks, and five shaft elements. The displacement vector of the shaft element between nodes 1 and 2 can be written as

((5))

The translation displacements and rotation displacements of a random point in the shaft element can be approximated as

((6))
((7))

where D is the rotation displacement shape function matrix, and N is the translation displacement shape function matrix. The expressions of N i and D i (i=1, 2, 3, 4) can be found in (Zeng, 2004).

Fig. 2
figure 2

The finite element model of a rotor system

Full size image

2.4 Hamilton principle

A general dynamic problem with N degrees of freedom should be expressed by N differential equations. By using the Hamilton principle, there is only one equation. This means that the Hamilton principle is highly recapitulative for dynamics problems and independent of the choice of coordinate systems, while the Newton equation changes along with the coordinate systems. Therefore, the Hamilton principle is more convenient to use for analyzing the dynamic systems.

The extended Hamilton principle can be expressed as

((8))

where T and U are kinetic energy and potential energy of the rotor system, respectively, and W is the work done by non-conservative forces and any forces not accounted for in the potential energy function.

The kinetic energy of the shaft element (T s) and the disk (T d) under the fixed coordinate system can be written as

((9))
((10))

where l is the length of the shaft element; ρ is the material density of shaft unit; I d and I p are the diametric inertia moment and polar inertia moment of the shaft element, respectively; and J d and J p are the diametric inertia moment and polar inertia moment of the disk, respectively.

Considering the bending and shearing distortion of the shaft element, the strain energy of the shaft element can be obtained as

((11))

where E is the Young’s modulus of elasticity, I is the area rotary inertia of shaft unit, k is the shape factor, G is the gravity vector of the rotor system, A is the sectional area of shaft unit, , and .

By substituting Eqs. (6) and (7) into Eqs. (9) and (11), the following energy equation of the shaft element can be obtained:

((12))
((13))

where

and

are the mass matrix and the gyroscope matrix of the shaft element, respectively;

is the stiffness matrix of the shaft element.

2.5 Motion equation

The shaft element motion equation can be achieved by substituting Eqs. (12) and (13) into Eq. (8) as follows:

((14))
((15))

where

C e and C d are the damping matrixes of shaft unit and disk, respectively; F e and F d are the vectors of forces acting on shaft unit and disk, respectively; and G e and G d are the gravity vectors of the shaft unit and disk, respectively.

The motion equation of the rotor system can be achieved by assembling all the shaft elements and the disks while ignoring the rotation displacement, which is given as follows:

((16))

where q=[x 1 y 1 x 2 y 2 … x 12 y 12]T is the nodal displacement vector of the rotor system, K is the stiffness matrix of the rotor system, G is the gravity vector of the rotor system, and

is the vector of the forces acting on the rotor system, where F dx =m d r d Ω 2cos(Ω t) and F dy =m d r d Ω 2sin(Ω t) are the eccentric force acting on the disk in the x and y directions, respectively.

We can define the dimensionless time and dimensionless displacement as follows:

((17))

Therefore, Eq. (16) can be transformed as

((18))

3 Numerical analysis and discussion

In this study, time-history diagrams, bifurcation diagrams, phase trajectory diagrams, and Poincare maps are used to illustrate the nonlinear dynamics of the rotor system. In the dynamic system, a bifurcation diagram shows the possible long-term values (equilibria/fixed points or periodic orbits) of a system as a function of a bifurcation parameter in the system. A Poincare map can be interpreted as a discrete dynamical system with a state space that is one dimension smaller than the original continuous dynamical system. By numerical integration of the equations of the rotor motion, the y coordinate, which denotes the dimensionless speed in the x or y direction, can be plotted versus the x coordinate, which stands for the dimensionless displacement in the x or y direction as the time increases. This produces a phase trajectory diagram.

The numerical analysis was performed by using the fourth-order Runge-Kutta method and implemented in MATLAB. The main partial parameters of the numerical calculation are listed in Table 1. The bifurcation diagrams at bearings and disks are plotted as shown in Fig. 3. Trajectory diagrams, time-history diagrams, and Poincare maps of disks and bearings at different rotational speeds are also listed to analyze the dynamic characteristics of the rotor system in this study. For simplicity, all the diagrams are only illustrated in the x direction. In the numerical calculation, the structural damping of the bearing is 2300 N∙s/m, and the disk’s is 5300 N∙s/m. These dampings are described as vectors and added to the material damping of the rotor system. The stiffness of the flexible coupler is very small, which can be deemed as 1/100 of its nearby shaft element. The dimensionless time ranges from 0 to 800π, and the step is 0.01π. Thus, 80001 data can be achieved. The latter 40001 data are plotted, while the first 40000 data are omitted considering the iteration veracity. In addition, the seal length, seal pressure drop, structural damping, and the stiffness of the rotor system are varied for the dynamic characteristic analysis of the rotor system’s stability.

Table 1 Partial test parameters
Full size table
Fig. 3
figure 3

Bifurcation diagram at bearing 1 (a), disk (b), and bearing 2 (c) with rotational speed increasing

Full size image
Fig. 4
figure 4

Numerical analysis results at Ω=300 rad/s

(a) Time-history diagram at bearing 1; (b) Trajectory diagram at bearing 1; (c) Poincare map at bearing 1; (d) Time-history diagram at the disk; (e) Trajectory diagram at the disk; (f) Poincare map at the disk; (g) Time-history diagram at bearing 2; (h) Trajectory diagram at bearing 2; (i) Poincare map at bearing 2; (j) Frequency spectrum at bearing 2

Full size image
Fig. 5
figure 5

Numerical analysis results at Ω=764 rad/s

(a) Time-history diagram at bearing 1; (b) Trajectory diagram at bearing 1; (c) Poincare map at bearing 1; (d) Time-history diagram at the disk; (e) Trajectory diagram at the disk; (f) Poincare map at the disk; (g) Time-history diagram at bearing 2; (h) Trajectory diagram at bearing 2; (i) Poincare map at bearing 2; (j) Frequency spectrum at bearing 2

Full size image

Fig. 3 shows the bifurcation diagrams at bearing 1, bearing 2, and the disks. Figs. 4–8 illustrate the time-history diagrams, the trajectory diagrams, and the Poincare maps with different rotational speeds at bearing 1, bearing 2, and the disks. In Fig. 3, it can be determined that at a lower speed, Ω≤700 rad/s, the system motion is maintained steady, which is a periodic motion and the amplitude is limited. When Ω>700 rad/s, the system turns into a triple periodic motion, which is shown as three isolated points in the corresponding Poincare map and three obvious frequency parts in the frequency spectrum. The trajectory diagram also reflects this characteristic. Meanwhile, the amplitude is greatly increased. When the rotational speed increases to 720–760 rad/s, the rotor system represents quasi-periodic motion, which is shown as a whole cycle in the Poincare map. The motion of the system becomes triple periodic motion in the range of 764 rad/s≤ Ω≤792 rad/s. With increasing rotational speed, the response of the rotor system becomes more complicated, and the quasi-periodic motion and multi-periodic motion happen alternatively. For example, the quasi-periodic motion happens at 900 rad/s (Fig. 6), the septuple periodic motion happens at 1216 rad/s (Fig. 7), and the quasi-periodic motion occurs at 1400 rad/s (Fig. 8), which is shown as a whole cycle in the Poincare map and some frequency parts are unable to be commonly divided in the spectrum of the frequency.

Fig. 6
figure 6

Numerical analysis results at Ω=900 rad/s

(a) Time-history diagram at bearing 1; (b) Trajectory diagram at bearing 1; (c) Poincare map at bearing 1; (d) Time-history diagram at the disk; (e) Trajectory diagram at the disk; (f) Poincare map at the disk; (g) Time-history diagram at bearing 2; (h) Trajectory diagram at bearing 2; (i) Poincare map at bearing 2; (j) Frequency spectrum at bearing 2

Full size image
Fig. 7
figure 7

Numerical analysis results at Ω=1216 rad/s

(a) Time-history diagram at bearing 1; (b) Trajectory diagram at bearing 1; (c) Poincare map at bearing 1; (d) Time-history diagram at the disk; (e) Trajectory diagram at the disk; (f) Poincare map at the disk; (g) Time-history diagram at bearing 2; (h) Trajectory diagram at bearing 2; (i) Poincare map at bearing 2; (j) Frequency spectrum at bearing 2

Full size image
Fig. 8
figure 8

Numerical analysis results at Ω=1400 rad/s

(a) Time-history diagram at bearing 1; (b) Trajectory diagram at bearing 1; (c) Poincare map at bearing 1; (d) Time-history diagram at the disk; (e) Trajectory diagram at the disk; (f) Poincare map at the disk; (g) Time-history diagram at bearing 2; (h) Trajectory diagram at bearing 2; (i) Poincare map at bearing 2; (j) Frequency spectrum at bearing 2

Full size image

4 Conclusions

In this paper, a model of a two-disk and two-span rotor/bearing/seal system is created based on the Hamilton principle, which is more accurate and easier for a numerical solution. The nonlinear coupling vibration of the rotor system is investigated by the fourth-order Runge-Kutta method. The figures achieved with various rotational speeds show the complexity of the nonlinear vibration and the bifurcation behavior of the rotor system. The numerical results indicate that the damping reduction and the increase of the seal pressure drop will make the bifurcation of the rotor system happen earlier, while the increase of the stiffness will prevent the bifurcation points from appearing too early and also greatly reduce the amplitude of the system motion. Based on the numerical analysis, it can be concluded that the rotational speed, the nonlinear seal force, the oil-film force, and the stiffness of the coupling have a great effect on the stability of a rotor system.

This study can enhance understanding of the nonlinear dynamics of rotor systems and be helpful in choosing some designing parameters which affect the stability of the rotor systems. Further work will be done to analyze the nonlinear dynamics of the multi-disk and multi-span rotor systems.