The Dynamic Analysis of Two-Rotor Three-Bearing System

A finite element model considering the shear effect and gyroscopic effect is developed to study the linear and nonlinear dynamic behavior of two-rotor three-bearing system named N+1 configuration with rub-impact in this paper. The influence of rotational speed, eccentric condition, and the stiffness of coupling on the dynamic behavior of N+1 configuration and the propagation of motion are discussed in detail. The linear rotordynamic analysis included an evaluation of rotor critical speed and unbalance response. The results show that the critical speed and unbalance response of rotors are sensitive to coupling stiffness in N+1 configuration. In the nonlinear analysis, bifurcation diagram, shaft-center trajectory, amplitude spectrum, and Poincaré map are used to analyze the dynamic behavior of the system.The results of the research transpire that these parameters have the great effects on the dynamic behavior of the system. The response of the system with rub-impact shows abundant nonlinear phenomena. The system will exhibit synchronous periodic motion, multiperiodic motion, quasiperiodic motion, and chaotic motion patterns under rotor-stator rub interaction conditions.The dynamic response is more complicated for flexible coupling and twomass eccentricities than that of system with rigid coupling and one mass eccentricity.


Introduction
Ultra supercritical steam turbine technology can meet notably the requirements for high efficiencies to reduce both fuel costs and emissions as well as for a reliable supply of electric energy at low cost [1,2].A typical ultra supercritical turbo-set comprises three separate turbine modules operating at different pressure and temperature levels.These modules are the high pressure turbine (HP), intermediate pressure turbine (IP), and, depending on the cooling water conditions, one, two, or three low pressure turbines (LP).The generator is directly coupled to the last LP turbine.There are two configurations for the shafting supporting structure of ultra supercritical steam turbine unit.In the first configuration named 2N configuration, the HP rotor, IP rotor, and LP rotors are supported by two journal bearings, respectively.In the second configuration, the HP rotor is supported by two journal bearings, and IP rotor and LP rotors are supported by only one journal bearing, respectively.Rotors are coupled through mechanical couplings.The configuration where  rotors are supported by +1 bearings is known as N+1 configuration.N+1 configuration has many advantages.For example, compact shafting structure and small bearing span exist in N+1 configuration to minimize foundation deformation for the bearing load and the alignment of shafting.As a consequence, the ultra supercritical unit with N+1 configuration is easy to install and maintain, and the investment of ultra supercritical power plant can significantly be reduced.Due to the specific supporting structure of N+1 configuration, however, rotor vibration is more greatly affected by adjacent rotor in the ultra supercritical turbo-set with N+1 configuration than with 2N configuration.It is very difficult to analyze and diagnose the vibration fault due to the decrease in the numbers of vibration monitoring locations for the ultra supercritical turbo-set with N+1 configuration.In addition, unstable rotor vibration was caused easily in N+1 configuration.The research about rotor vibration of the ultra supercritical turbo-set with N+1 configuration is still in the preliminary and exploratory stage at present.
Most large rotating machinery consists of multiple rotors, which are coupled through mechanical couplings.There are many types of industrial coupling some of which are rigid, gear, and flexible type.The couplings' function is to transmit torque from the driver to the driven rotor.Coupling plays an important role in turbomachinery.Pennacchi et al. [3] propose a complete and original method to simulate 2 International Journal of Rotating Machinery the behavior of real shaft line with rigid coupling and analyze nonlinear character of the system.A rigid mechanical coupling model and a flexible coupling model are built to examine the lateral and torsional responses of two rotating shafts, respectively, with theoretical and numerical analysis in [4,5].Redmond [6] presented a model which enables dynamic analysis of flexibly coupled misaligned shafts and investigated the influence of coupling-stiffness anisotropy, bearing nonlinearity, mass unbalance, and static torquetransmission effects.
In turbomachinery, such as ultra supercritical turbo-set with N+1 configuration, the requirement of high efficiency has led to reduced clearances between rotating and stator components.With tighter clearances, rotor-stator rub interaction faults are prominent in high performance turbomachinery operations.Rotor physical contact with a stationary element of rotating machine and the subsequent rubbing at the contact area, which involves several major physical phenomena such as sliding friction, impact, and modification of the system stiffness, contribute to the potentially serious malfunctions in rotating machinery that may lead to machine failure.The malfunction may be mild under light rub or complete destruction of the machine under severe rub [7].Light partial rubs and full annular rubs would cause major changes in the rotor normally synchronous vibrations.Heavy rubs may cause significant impacts as well as subsynchronous and supersynchronous vibrations of the shaft.The system vibration often contains some very complicated phenomena including not only periodic motion but also chaotic motion.
The nonlinear dynamic behavior and stability of rotorbearing system subjected to rub-impact have been studied by scholars and engineers for some time.Different types of fault symptoms have been reported by Muszynska [8] and Ahmad [9] who give a detailed literature review about this complex physical phenomenon.The majority of published works have been focused on the global dynamic performance of rotor-bearing system by using different mathematical models.Khanlo et al. [10] studied the chaotic vibration of a rotating flexible continuous shaft-disk system with rubimpact including the Coriolis and centrifugal effects.The results show that the effect of centrifugal force was greater than that of the Coriolis force on the occurrence of the rub-impact, and the rub-impact would occur at lower speed ratios under the influence of the Coriolis and centrifugal forces.Torkhani et al. [11] presented a model suitable to be employed for rotor system under light, medium, and heavy partial rubs and the corresponding predictions are compared with the experimental results based on the big test rig supported by elliptical bearing.Chu and Lu [12] investigated experimentally the nonlinear dynamic vibration of some different rotor-bearing systems including one and two rotors with single-and multi-disks.The experiment results show that the system motion generally contains the multiple harmonic components and fractional harmonic components.Shen et al. [13] studied the dynamic behavior of a rub-impact rotor-bearing system with initial permanent bow.
Based on the short bearing theory, the nonlinear oil film forces from the journal bearing are obtained.The influence of rotating speeds, initial permanent bow lengths, and phase angles on the response of system is investigated.Wang et al. [14] investigated the nonlinear dynamic characteristics of a single disk rotor system supported by oil film journal bearings using lumped-mass model.The effects of rotating speed and damping ratio on dynamic behavior of system are discussed.The results reveal a complex dynamic behavior comprising periodic, multiperiodic, chaotic, and quasiperiodic response.Zhang and Wang [15] investigated the effects of phase difference between disks on rotor dynamic response, highlighting field balancing of the rotor as a possible solution in rub-impact mitigation.Roques et al. [16] introduced a rotor-stator model of a turbo generator in order to investigate speed transients with rotor-diaphragm rubbing caused by an accidental blade-off imbalance.The highly nonlinear equations due to contact conditions are solved through an explicit prediction, correction time-marching procedure combined with the Lagrange multiplier approach dealing with a node-to-line contact strategy.The results highlight the significant role of the friction coefficient together with the diaphragm modeling, from rigid to fully flexible, in the interaction phenomenon.Based on finite element theory with contact analysis, Behzad et al. [17] developed an algorithm to investigate the rotor-to-stator partial rubbing vibration and obtained the responses of the partial rubbing for different rotational speed.In addition, the sensitivity of the initial clearance, the stator stiffness, the damping parameter, and the coefficient of friction is analyzed.Ma et al. [18] investigated numerically and experimentally the fault characteristics of a two-disk rotor system with a point-point contact rubbing model between a disc and an elastic rod.
From the investigations cited above, it is known that much work has been done, respectively, on the research of the single span rotor system with fewer disks and rubimpact.Based on authors' best knowledge, few publications have considered the nonlinear dynamic behavior of rotorbearing system with N+1 configuration subjected to rubimpact.In this paper, the dynamic model of two-rotor threebearing system with rub-impact is established using finite element method.The effect of bearings is expressed in form of stiffness and damping coefficients which are varying with rotating speed.Attention is paid to the influence of rotational speed, eccentric condition, and the stiffness of coupling on the dynamic behavior of two-rotor three-bearing system and the propagation of motion.Bifurcation diagram, shaft-centre trajectory, amplitude spectrum, and Poincaré map are used to analyze the dynamic behavior of the system.

Physical Model and Governing Equations
In the right side of the equation, the first one is the sliding velocity term that gives rise to the bearing stiffness coefficients and the second one is the squeeze-film velocity term that gives rise to the bearing damping coefficients. = (, ), ℎ = ℎ(, ), 0 ≤  ≤ 2, −/2 ≤  ≤ /2, and  denotes viscosity.Also, (, ) and ℎ(, ) are the film pressure distribution and the film thickness distribution, respectively.And  is the hydrodynamic-active axial length of the journal bearing.Here only (, ) is the "unknown" pressure and all other parameters are specified.By solving this equation and integrating (, ) over the journal cylindrical surface, the bearing forces in  and  direction can be obtained as given in ( The corresponding eight stiffness and damping coefficients are compactly expressed using subscript notation, as given in (3).Here,  and  denote coordinate directions: In most of the previous papers, (1) was solved by neglecting either the axial pressure flow term ("long bearing" solution) or the circumferential pressure flow term ("short bearing" solution).With either approximation, the RLE is reduced to an ordinary differential equation.These two approximate solutions provide an upper bound and lower bound on bearing forces, respectively.In this paper, the bearing stiffness is based on a full 2D numerical solution algorithm.The rotor system is supported by elliptical bearings.Parameters of bearing used in this paper are listed in Table 1 and the eight coefficients of the type of bearing were evaluated and plotted versus operating speed in Figure 1.

Rub-Impact Force.
The rub-impact is assumed to take place at the large disk.There is an initial clearance of   between the disk and stator.The rubbing model is shown in Figure 2. Without considering thermal effect from friction, the relationship between radial rotor displacement and contact force at disk 2 can be defined as where  = √ 2 +  2 is the radial displacement. and  are the displacement components of the geometric centre of the disk in  and  directions. is the friction coefficient between disk and stator.  is the radial stiffness of the stator.The components of the rub-impact forces in  and  directions can be expressed as (5)

Dynamic Governing Equation.
The finite element method is used to establish the dynamic governing equation of rotors system based on the Timoshenko beam theory including the shear effect, gyroscopic effect, and transverse torsion in this paper.A two-rotor three-bearing system is considered.The rotor-bearing system that includes two shafts and four disks is supported by three journal bearings.The left rotor including two disks is supported by two elliptical journal bearings, and the right rotor including two disks is supported by one elliptical journal bearing.Two shafts are coupled through a mechanical coupling.Very often coupling is modeled as an elastic component with isotropic translational stiffness and rotational stiffness between two stations.The system is discretized by use of beam elements as shown in Figure 3.The left rotor is 970 mm in length and the right rotor is 835 mm.The two-rotor diameter is 50 mm with the two big disks having a diameter of 270 mm and a length of 55 mm and the two small disks having a diameter of 123 mm and a length of 70 mm.The length is broken down into 41 Timoshenko beam elements in the finite element model.Bearings are located at node 6, node 20, and node 39, respectively.Every node has 4 degrees of freedom, two translational and two rotational:  and  are the translation freedoms in  and  directions, and   and   denote the rotation angles of the cross section around axes  and , respectively.The whole nodal displacements can be expressed as   = (        )  ,  = 1, 2, . . ., 42.By integrating all the local matrices of shaft elements, disks, bearings, and unbalances, the dynamic behavior of the whole rotors system is described as the following second order differential equation:   where , , and  are the mass matrix, stiffness matrix, and damping matrix due to structure and bearings, respectively. is the gyroscopic matrix of the rotors system with the angular velocity . is the load vector, which is composed of three parts, that is, the unbalance forces, rotor-stator contact forces, and gravity force.The damping matrix includes Rayleigh damping matrix and gyroscopic matrix.The Newmark method is used to solve the transient differential equations.The chosen time step is 2/( * 300) s, so the convergence of the quantities of interest with respect to the time step is achieved.

Numerical Calculation and Discussion
In this paper, the main object is to evaluate the effects of rotational speed, eccentric condition, and the stiffness of coupling on the dynamic behavior of rotor-bearing system with N+1 configuration and the propagation of motion as it related to a rub-impact condition.The radial stiffness of stator is 4.57 N/m.The clearance and friction coefficient between disk and stator are 0.06 mm and 0.1, respectively.The speed range considered was from 300 to 15,000 rpm.

Critical Speed and Mode.
To evaluate the effects of the coupling parameters on the dynamic response of the rotorbearing system with N+1 configuration, critical speed and mode for different values of coupling stiffness were analyzed   in the rotor dynamic model.To begin the analysis, a critical speed map of the system was developed as shown in Figure 4.The first, third, 4th, and 5th critical speeds are more sensitive to coupling parameter when the value of coupling stiffness is less than 1.07 N/m and elevated with the increasing of coupling stiffness.While the value of coupling stiffness is greater than 1.07 N/m, the change about the first, third, 4th, and 5th critical speeds is not obvious, and the damped critical speed tends to be stable.With the change of coupling stiffness, however, the value of the second critical speed of the system almost was constant.The second critical speed of the system was not sensitive to coupling property in that it was the first critical speed of the left rotor supported by two bearings.
Figure 5 shows the mode shapes for the first five synchronous critical speeds of the system when the value of coupling stiffness is 1.06 N/m and 1.09 N/m, respectively.The first mode shape of the system is corresponding to the first critical speed of right rotor.The second, the third, the fourth, and the fifth mode shapes of the system are corresponding to the first critical speed of left rotor, the second critical speed of the right rotor, the second critical speed of the left rotor, and the third critical speed of the right rotor, respectively.

Effect of Rotor Unbalance on the Dynamic Response.
The rotors always have some amount of residual unbalance no matter how well they are balanced.These residual unbalance forces are discrete and located at different planes with different amplitudes.18 cases were analyzed in the rotor dynamic model of rotor-bearing system with N+1 configuration (Table 2).The parameters chosen for evaluation include the location of unbalance, the amplitude of unbalance, and the stiffness of coupling.Unbalance forces with   different amplitude were applied at the two big disk locations, respectively, and both two big disk locations to excite the rotor in a manner that would produce a rub-impact condition.The linear rotor dynamic analysis was then performed for the 18 cases in order to evaluate the effect of the coupling characteristics on the dynamic response of the system.
As shown in Figures 6-8, two big disks radial displacements of a set of cases are considered.The displacements rapidly increase around 4500-5500 rpm as the rotor passes through the first critical speed.The displacement of big disk with unbalance is larger than that of another big disk without unbalance.The greater the amplitude of unbalance is, the greater the displacement of disk is.The unbalance response of two disks, shown in Figures 9-11, was sensitive to coupling stiffness when the amplitude of unbalance was constant.For the left big disk with unbalance, the unbalance response of left big disk was not sensitive to coupling stiffness in that the left rotor is supported by two bearings.When the right big disk is unbalanced, the response of right big disk was very sensitive to coupling property.The first critical speed excited by unbalance force increases with the increase of coupling stiffness for the right big disk.

Effect of Rotational Speed on the Dynamic Response.
Firstly, the effect of rotational speed on the dynamic response of the system is discussed in this section because the speed is one of the most important parameters of a rotor-bearing system.Maintaining other parameters unchanged and taking rotational speed as the control parameter, the nonlinear dynamic behavior of the rotor supported by three elliptical bearings at different rotating speeds was investigated.Two cases, one in which left rotor has mass eccentricity and the other in which both these two rotors have mass eccentricities, are considered.The value of the mass eccentricity is 0.001 kg⋅m.Figure 12 shows bifurcation diagrams of  displacement of the right big disk at different rotating speeds when the radial stiffness of coupling   is 1.05, 1.07, and 1.09 N/m for two mass eccentricities.For one mass eccentricity, the bifurcation diagrams of displacement of the right big disk are plotted in Figure 13 when   are 1.05, 1.07, and 1.09 N/m.The maximum radial displacements of the right big disk at different rotating speeds and the radial   stiffness of coupling for two cases are obtained, as shown in Figure 14.
In nonlinear dynamics, the bifurcation diagram often is used to provide a summary of essential dynamics of the system with varying of governing parameter.As shown in Figures 12 and 13, the most complicated case is Figure 12(a).When rotational speed is small and less than 283 rad/s, the motions of the system are stable with synchronous periodic motion.The system will show double and multiperiodic motion when the speed is between 293 and 346 rad/s.Then the system shows the period-one motion again for a large range from 356 to 775 rad/s.As the rotational speed is changed from 785 to 995 rad/s, the multiperiodic motion appears again.When the rotational speed locates between 1026 and 1340 rad/s, the system will become irregular and enter into the chaotic motion.In a narrow range from 1351 to 1382 rad/s of rotational speed, the system shows the multiperiod motion again.Once the speed is higher than 1393 rad/s, the rotor vibrations are chaotic motion.
For Figure 12(b), with the increasing of rotating speed, the system will experience period-one motion, double periodic motion, period-one motion, double periodic motion, quasiperiodic motion, period-one motion, quasiperiodic motion, chaotic motion, double and multiperiodic motion, and chaotic motion in sequence.When   = 1.09N/m and considering two mass eccentricities, the range of rotational speed for chaotic motion is smaller than those when   = 1.05 and 1.07 N/m.But the motion is still complicated including period-one motion, double and multiperiodic motion, and chaotic motion.When only left rotor has mass eccentricity, the motions are more complicated for flexible coupling than those for rigid coupling in which the periodone motions and multiperiodic motion are prominent.
From Figure 14, it can be seen that when both two rotors have mass eccentricities, the rub-impact is easier to take place even though the rotational speed is small.For one mass eccentricity or two mass eccentricities, the lower the value of the radial stiffness of coupling is, the more complicated the motion of system is.In order to make the system more stable, it is better to choose the rigid coupling to connect two rotors and avoid large and more mass eccentricities.
Besides the bifurcation diagram, shaft-centre trajectory, amplitude spectrum, and Poincaré map are used to further illustrate the dynamic behaviors of the presented system when given parameters.The results are shown in Figures 15  and 16.When  = 869 rad/s and   = 1.07N/m for two mass eccentricities, the system exhibits a quasiperiodic motion.The Poincaré section map is a closed curve.The shaftcentre trajectory has a regular pattern.In Fourier amplitude spectrum, there is two-peak amplitude.Nearly all of the vibration energy is distributed between the subsynchronous and synchronous components, but especially in the synchronous component.For one mass eccentricity and keeping the other parameters unchanged, from Figure 16 it can be known that there is only a single point in the Poincaré map and the shaft-centre trajectory is a closed circle.Furthermore, Fourier amplitude spectrum also has one-peak amplitude at 1X component.All these phenomena mean that the system is synchronous period-one motion.

Effect of Stiffness of Coupling on the Dynamic Response.
For multirotor system, the coupling will play an important role in the dynamic behavior of system.In this section, the effect of radial stiffness of the coupling on the dynamic response of the system is discussed in detail.Figure 17 shows the bifurcation diagrams of  displacement of the right big disk at different radial stiffness of coupling   when the rotational speed  is 314, 995, and 1445 rad/s for one mass eccentricity.When the rotational speed is 314 rad/s, Figure 17(a) indicates that the system only appears synchronous periodic motion.For  = 995 rad/s, the system mainly exhibits double periodic motion and period-one motion as shown in Figure 17 are shown in Figure 18.There are two discrete points in the Poincaré map and the shaft-centre trajectory is two closed circles.When the coupling is flexible and  = 1445 rad/s, the system shows chaotic motion in Figure 17(c).The response characteristics of a chaotic motion are plotted in Figure 19.The shaft-centre trajectory is disorder, and the Poincaré map shows many discrete points.Correspondingly, Fourier amplitude spectrum is continuous and the vibration energy is mainly distributed in the subsynchronous component.But for rigid coupling, the system is period-one motion.All the results reveal that the motion of system is complicated when the coupling is very flexible.

Conclusions
The linear and nonlinear dynamic behavior of a two-rotor system with four disks supported by three elliptical bearings are investigated by using finite element method including the shear effect and gyroscopic effect in this paper.The effect of bearings is expressed in form of stiffness and damping coefficients which are varying with rotating speed.The effects of two important parameters, the rotational speed and the radial stiffness of coupling, on the dynamic response of the system are examined in detail.The linear analysis results show that the critical speed and unbalance response of rotors are sensitive to coupling stiffness in N+1 configuration.Furthermore, two cases, one in which only the left rotor has mass eccentricity and the other in which both two rotors have mass eccentricities, are considered in the nonlinear analysis.For the different combination of system parameters, numerical results can be obtained and plotted by bifurcation diagrams, shaft-centre trajectory, Poincaré map, and amplitude spectrum.The results indicate that the dynamic response of the presented system is very abundant.When the coupling is flexible, the response of system will comprise synchronous motion with period-one, double, and multiperiodic motion, quasiperiodic motion, and chaotic motion with varying of rotational speed.Generally, the dynamic behavior of the system for flexible coupling is more complicated than that for rigid coupling and the system motion generally contains the multiple harmonic components and fractional harmonic components.The system mainly is synchronous periodic motion as the increasing of radial stiffness of coupling.On the other hand, the results also show that when every rotor has mass eccentricity the response appears more nonlinear phenomena than that when only left rotor has mass eccentricity, especially for rigid coupling.According to the above analysis, it can be concluded that, in order to make the system more stable, it is better to use the rigid coupling to connect two rotors and avoid large and more mass eccentricities.The results developed in this study are useful in identifying, understanding, and avoiding some undesirable vibration behaviors in these types of rotating machinery.

REL: Reynolds Lubrication Equation 𝑝:
The film pressure ℎ: The film thickness : Lubricating oil viscosity : The hydrodynamic-active axial length of the journal bearing : Angle variable : The stiffness matrix of the whole rotors system : The damping matrix of the whole rotors system : The gyroscopic matrix of the whole rotors system {}: The displacement vector {}: The load vector.

Figure 1 :
Figure 1: Dynamic stiffness (a) and damping coefficients (b) of bearing with the speed increase.

Figure 2 :
Figure2: The mechanic model of rotor system with radial rubimpact.

Figure 3 :
Figure 3: The finite element model of rotors system.

Figure 4 :
Figure 4: Critical speed map of the rotor-bearing system with N+1 configuration.

Figure 10 :
Figure 10: Disk displacements for the 1.06 N/m, 1.07 N/m, 1.08 N/m, and 1.09 N/m coupling stiffness and 1.0 kg⋅mm ∠0 ∘ mass unbalance at the left big disk.

Figure 15 :
Figure 15: The response characteristics of system at  = 869 rad/s and   = 1.07N/m for two mass eccentricities.

Figure 16 :
Figure16: The response characteristics of system at  = 869 rad/s and   = 1.07N/m for one mass eccentricity.

Figure 17 :
Figure 17: The bifurcation diagrams versus the radial stiffness of coupling for one mass eccentricity.

Figure 18 :
Figure 18: The response characteristics of system at  = 995 rad/s and   = 3.65 N/m for one mass eccentricity.

Figure 19 :
Figure 19: The response characteristics of system at  = 1445 rad/s and   = 8.05 N/m for one mass eccentricity.

Table 1 :
Parameters of bearings.

Table 2 :
Cases with different unbalance parameters and different values of coupling stiffness.