Dynamic behaviour of rotors supported by fluid-film bearings operated close to fluid-induced instability

This paper is focused on an analysis of rotating systems with fluid-film bearings, especially on their nonlinear behaviour in the course of developing fluid-induced instability. The studied system consists of Jeffcott rotor supported by a fluid-film bearing characterised by the Reynolds equation. The steady state response of the system is investigated by means of an approximate analytical solution of the Reynolds equation while the transient response of the system is investigated using a complex numerical solution. Results suggest that the rotor exercise a bounded chaotic motion if it becomes unstable. If the fluid-induced instability further develops, the motion actually becomes less chaotic and can be characterised as quasi-periodic.


Introduction
Fluid film bearings are widely used for supporting of various rotating machines because of their favourable properties, specifically low friction and wear and vibrationreducing capabilities.However, it has been known for more than 90 years that a fluid-induced instability can develop in a poorly designed bearing [1].The quality of the bearing design can be expressed by Sommerfeld number S 0 = ( p ψ 2 /η ω 0 , where p, ψ are a rated load relative clearance of the bearing, η is a dynamic viscosity of a fluid film and ω 0 is a rotor speed [2].In general, bearings with low S 0 -i.e.not enough loaded bearings, bearings with high-speed rotors or bearings with highly viscous filmsare prone to the the fluid-induced instability.
The instability, also known as oil whirl, was mathematically described at the turn of the 50s and 60s [3,4].Both authors also expressed a stability criterion for a linearised rotor-bearing system which is based on the Routh-Hurwitz criterion.The occurrence of the fluid-induced instability and the processes accompanying it were studied experimentally in the early 80s [5][6][7].
The oil whirl starts when hydrodynamical (HD) forces in the fluid film become higher than loading forces, which are predominantly caused by static loads.When this occurs, the HD forces start to push a journal -a part of a shaft, which rests in a bearing -in a forward circular motion around the bearing.The journal runs about the bearing with a period which in general lasts slightly longer than two revolutions.Therefore the oil whirl can be easily recognized by frequency which is 0.42-0.49to journal speed [6].
HD forces are determined by a pressure distribution of the fluid film, which is governed by a partial differential equation known as the Reynolds equation (RE) [8].The exact analytical solution of RE has not been known until 2012 [9], however, approximate analytical solutions for infinitely short [10,11] and infinitely long bearings [11] have been known since the 50's and later many researchers used perturbation methods in order to obtain approximate analytical solutions of full RE [12,13].As the computing power of personal computers had been rising, robust numerical methods including finite differential [14], finite element [14,15], finite volume [16] and element-based finite volume methods [17] were also employed.
Although all mentioned methods capture the phenomenon of oil whirl, there are relatively small number of studies which focus on nonlinear behaviours of rotors operated close to the instability threshold speed.A transient state which occurs in rotor-bearing system when a journal pass the threshold speed is studied in [13], where small amplitude perturbations are assumed in the bearing.The transient state of Jeffcott rotor supported by a gas bearing and accompanying bifurcations are studied in [18].
The main focus of this work is to perform an in-depth analysis of the transient state for Jeffcott rotor supported by journal bearing with consideration of all terms which are presented in RE.This may lead to better understanding of the oil whirl instability, and hence to understand how to design rotor systems which are operated in the unstable regime.
The rest of the article is outlined as follows.Section 2 presents a simple model on which computations of steadystate response were performed.Section 3 contains a brief description of an extended model for transient analysis.Model and simulation parameters, and results are discussed in sections 4 and 5, respectively.Finally, section 6 draws some conclusions.

Model for steady-state analysis
A general mathematical model of a simple rigid rotor which rotates at a constant angular velocity and is supported by a generic journal bearing is described in this section.

Equation of motion for a rigid rotor
We assume the rigid rotor of mass m which rotates at a constant angular velocity ω 0 and is subject to constant gravitational and out-of-balance forces f g and f u and to a vector of hydrodynamic forces f hd (see fig. 1).The equilibrium of acting and inertia forces can be written as [2] where M is a mass matrix and q(t) is an acceleration vector.If we further assume that the rotor does not perform yaw and pitch movements then its motion is given by horizontal displacement and vertical displacement x = x(t) and y = y(t) and Eq. ( 1) can be written in the form where g is the gravitational constant, U st is the static unbalance and F hd x and F hd y are horizontal and vertical components of vector f hd [2].
Alternatively, the motion of the rotor can be characterised by eccentricity e = e(t) and angle to the eccentricity φ = φ(t), which are related to x and y as follows Time derivatives ė and φ of e and φ are obviously A rigid rotor supported by a journal bearing.

Hydrodynamic forces in a journal bearing
Hydrodynamic forces are obtained by integrating hydrodynamic pressure in an oil film over a surface of the bearing.Hydrodynamic pressure p = p(s, z, t) is given by the Reynolds equation [8] ∂ ∂s where s and z are circumferential and axial coordinates, respectively, h = h(s, z, t) is a gap between the journal and the bearing and µ is dynamic viscosity of the oil film.The Reynolds equation ( 7) can be solved under various assumptions.Let us assume the half Sommerfeld conditions [10] and consider the following relations where l is the bearing length.Then radial and tangential components of hydrodynamic forces F hd,sb rad and F hd,sb tan shown in Fig. 1 can be expressed as [11] where ε = e /c r is the relative eccentricity.Equations ( 9) and ( 10) approximate the analytical solution of Reynolds equation ( 7) reasonably accurately if l < r holds [11].If l is from interval r, 4 r , the components of hydrodynamic forces have to be corrected.Bastani proposed [11] corrective functions in the form where f i , g i are cubic or quartic polynomial functions which are dependant on l /r and are expressed in [11].Hydrodynamic forces F hd rad and F hd tan are derived in a rotated coordinate system shown in fig. 1 and can be transformed to a fixed coordinate system from equation of motion (2) using relations

Solution strategy
Equation of motion (2) can be formally rewritten as where f = f( q, q, t) is a vector of all acting forces.Eq. ( 15) can be further transformed to the state-space; the resulting equation is then The equation of motion in this form can be integrated numerically by means of the Runge-Kutta method.

Model for transient analysis
An extension to the previous model allowing a transient analysis for ω 0 const.is described in this section.

Equation of motion
Let us assume that the angular speed of the rigid rotor from section 2.1 is time-dependant variable φ = φ(t).Since we want to be able to control φ, we introduce a controller and couple the controlle and the rotor with a rotational spring as shown in Fig. 2. The rotation of the rotor in such configuration is given by the equation where I is a moment of inertia of the rotor to the axis of rotation, d r , k r are rotational damping and stiffness of the rotational spring, respectively, Φ = Φ(t) and Ω = Ω(t) are arbitrary functions of controller's angular displacement and angular velocity and M f ric is a frictional moment [2].
Since frictional moment M f ric is of no concern because we do not analyse losses due friction nor thermodynamics of the oil film, M f ric is further neglected.With respect to the fact that angular velocity of the rotor φ is not constant, equations of motion are where the components of hydrodynamic forces FEM F hd x and FEM F hd y are obtained from finite-element-based numerical solution of the Reynolds equation [15].
In the first step of the solution process, hydrodynamic pressure p is obtained from Reynolds equation (7) under the assumptions tha i) the bearing is fully-flooded, ii) there is constant ambient pressure p a at edges of the bearing, and c) if hydrodynamic pressure p in any node of FE mesh drops below the predefined value p s , then p = p s is set in that particular node.
In the second step, p is integrated over the bearing surface.This integral yields FEM F hd x and FEM F hd y [9].Note that the approximate analytical solution of Eq. ( 7) which is described in section 2.2 is derived under the assumption that the angular speed of the journal is constant and therefore cannot be used for transient analysis.

Solution strategy
The system of eqs.(18) has to be transformed to the state-space in order to be integrable by standard numerical solvers.The transformation was briefly described in section 2.3.

Model and simulation parameters
Parameters of the system and their respective symbols for both steady-state and transient analysis can be found in tab. 1.The parameters are chosen so that the fluid-induced instability arises at ca. 2300 rpm.
The rotor was considered perfectly balanced in the case of the steady-state analysis.The simulations were performed in the speed interval n ∈ 2000, 2900 with the step 5 rpm.Time integration was performed in the interval t ∈ 0, 1.5 s with a variable time step.The results were then transformed to time series with the constant time step 10 −4 s.Only the time interval t ∈ 0.5, 1.5 s was subjected to subsequent post-processing.
The transient analysis was performed in the speed interval n ∈ 800, 4000 rpm with several values of constant angular acceleration and static unbalance U st .The selected values of U st correspond with balance quality grades G0, G6.3 and G16 per ISO 1940.Time integration was performed with a variable time step ∆t ∈ 5 × 10 −7 , 5 × 10 −5 s and the results were then transformed to time series with the constant time step 10 −4 s.

Results and discussion
The results of the aforementioned steady-state and transient analysis are presented and discussed in this section.

Steady-state analysis
Fig. 3a shows the results of the steady-state response analysis, more precisely there are plotted local maxima and minima of relative eccentricity ε.We assume that depicted trajectories are reasonably close to trajectories in respective limit cycles.It is apparent that until ca.2350 rpm the system is working at its equilibrium.If the rotor speed further grows, the trajectory will become elliptic as shown in fig.4a.This elliptic motion is happening at a frequency which is slightly lower than half of the rotor speed and might actually be chaotic because significant broadband noise is present [19]  as can be seen in fig.4c.It is safe to assume that around 2350 rpm a threshold speed of the fluid-induced instability can be found.After 2500 rpm mark is passed, the system comes into quasi-periodic motion with apparent response at the rotor speed and at even higher speed the motion is clearly limited by bearing clearance and the fluid-induced instability becomes fully developed.We can study the trajectories more rigorously, if we transform them to a coordinate system (χ, ψ) which rotates around z axis at the same speed as the rotor [19].
We further refer to the transformed trajectories as rotating trajectories.Fig. 4b shows that there are three distinct types of the rotating trajectories in the zone delimited by the threshold speed of the instability and by the speed at which the instability is fully-developed (we call this zone the transient zone in the following text).Firstly, there is almost circular trajectory which shows properties of a bounded chaotic motion [19].If the speed increases, the system comes into stable period doubling and tripling motion with less significant noise, possibly indicating quasiperiodic solution [13].Finally, when the instability fully develops, the rotating trajectory is circular and the rotor covers roughly 540 • in one cycle.All described phenomena are exhibited in the results of both steady-state and transient analysis with zero angular acceleration.However, all qualitative changes happen at higher speeds in the case of the transient analysis.

Transient analysis
Fig. 5a shows how the threshold speed of the instability is influenced by angular acceleration of the rotor.Apparently, if the angular acceleration is higher, the system will became unstable at higher speed.Fig. 5b shows that the influence of static unbalance on the threshold speed is rather small.Fig. 6 demonstrates there are three distinct types of rotating trajectories not only if the rotor is spinning at constant speed but also if its speed is increasing.Note that trajectories and spectra depicted in fig.6 are taken from the simulation of the perfectly balanced rotor with the lowest angular acceleration because this simulation offers the best frequency resolution of frequency analysis.
Similar trajectories are to be found in all simulated cases.However, if the rotor is unbalanced, there is an additional synchronous (1X) component present in the response.

Conclusions
This paper provided a brief yet coherent theory for dynamic behaviour of rigid rotors supported by fluid-film bearings occurring in the so called transient zone delimited by the threshold speed of the fluid-induced instability and by the speed at which the instability is fully-developed.
There are three distinct types of the rotor motion in the transient zone.The first type appears as a bounded chaotic motion.Interestingly, later two types are less chaotic and tend to be quasi-periodic.This behaviour is presumably caused by the fact that these motions are more limited by the geometry of the bearing.The later two types of motion are also characteristic by the response with synchronous (1X) components even if the rotor is perfectly balanced.
These three types of motion can be find not only if the rotor is analysed at constant speed but also if the angular speed is being constantly increased and if the rotor is unbalanced.
Preliminary studies suggest that the presented conclusions are valid also for journals of flexible rotors which are supported by multiple journal bearings.However, to prove that such phenomena happen in general is left to future work.

Figure 2 .
Figure2.A rigid rotor coupled to a controller.

Figure 3 .
Figure 3. Results of brute force bifurcation.Depicted are a) local extremes of steady-state responses and b) values of relative rotating displacement ψ if relative rotating displacement χ is 0. Relative rotating displacements are defined in Eqs.(21) and (22).Fig. 3b can be also perceived as the intersection of rotating trajectories with nψ plane.The intersection is highlighted in fig.4b.

Figure 4 .
Figure 4. Trajectories (obtained from steady-state analysis) of the rotor at selected speeds are depicted in a) fixed coordinates and b) rotating coordinates.Black dots show the position of the rotor at the beginning of each turn.There are also shown responses in the frequency domain in fig.4c.

Figure 5 .Figure 6 .
Figure 5.The results of transient analysis show that a) the threshold speed for the fluid-induced instability is dependent on the angular acceleration of the system and b) the static unbalance introduces new local extremes but almost does not influence the threshold speed.

Table 1 .
Parameters of analysed system.Items denoted with asterisk symbol ( * ) are used only during transient analysis.

Table 2 .
Summary of performed simulations.