Numerical Analysis of Fluid‐Film‐Cavitation on Rotordynamic Vibration and Stability Behavior

The dynamic vibration and stability behavior of a rotor supported in journal bearings is investigated. In the lubricating film of a hydrodynamic bearing, different cavitation phenomena may occur. Cavitation in hydrodynamic bearings is most commonly caused by surrounding air, which is sucked into the bearing gap and also by outgassing of dissolved gases. Cavitation effects have to be considered in order to precisely analyze the performance of rotor/bearing systems. Numerous approaches have been proposed to treat cavitation using appropriate boundary conditions in the REYNOLDS fluid film equation, which is applied to calculate the hydrodynamic pressure distribution in the bearing gap. These approaches are generally divided into two types: non mass‐conserving and mass‐conserving cavitation models. In this paper, a mass‐conserving cavitation model based on a 2‐phase approach (lubricant/air mixture) is applied. For the 2‐phase model, a relationship has to be defined that describes the dependency between density of the oil/gas‐mixture and the pressure. Different density/pressure models are presented and discussed. Numerical run‐up simulations are performed to study the stability of rotor/bearing systems incorporating different cavitation models. The numerical model for the rotor/bearing system consists of a multibody model for the rotor and a finite element model for the fluid film. Both models are coupled by a co‐simulation approach. The different cavitation models are compared. The influence of the cavitation model and the cavitation parameters on the stability behavior of the rotor/bearing system are analyzed. The cavitation approaches are also compared with respect to accuracy and efficiency.


.1 Mass-conserving cavitation models
Different methods exist to physically and mathematically describe fluid film cavitation in hydrodynamic bearings [4]. Here, a 2-phase cavitation approach is used. Therefore, a mixture approach is applied to represent the mixture of lubricant and gas in the bearing gap. With the nonlinear lubricant fraction function ϑ according to ϑ(p) = ρ(p) ρ0 = η(p) η0 , ϑ (p) = ∂ϑ(p) ∂p , an averaged value for density ρ and viscosity η is calculated as a function of the local pressure p. Note that p 0 denotes the ambient pressure and ρ 0 = ρ(p 0 ), η 0 = η(p 0 ), where ρ 0 and η 0 denote the properties of the liquid lubricant, respectively. The function ϑ separates the full-film zone (p ≥ p 0 ) from the cavitation zone (p < p 0 ). In the full-film zone (ϑ ≈ 1, ϑ 1), the oil almost behaves like an incompressible fluid. In the cavitation zone (0 < ϑ < 1), ϑ decribes the mixture between oil and gas (air or dissolved lubricant).

(a) and (b). a) Regularized step function
In the first approach, a regularized step function is used, which is mathematically described by a third order polynomial for the transition zone (x 0 < x < 1, e.g. x 0 = 0.9) and linear functions with a small gradient for the high pressure zone (x ≥ 1) and the low pressure zone (x ≤ x 0 ) so that ϑ (x) is continuous, see Fig. 1(a). Using the regularized step function, numerical problems may occur. Due to the high gradient ϑ 1 in the transition zone, a large convection term is produced in the compressible REYNOLDS equation. In the high pressure zone and also in the low pressure zone, the small gradient ϑ 1 yields a very small coefficient for the ∂p ∂t -term in the time-dependent REYNOLDS equation. Both, very small and very large values of ϑ entail stiff differential equations [2].

and d) Polynomials of degree 2 and 3
The third and fourth ansatz are both polynomials in the region x ≤ 1: ϑ(x) = −x 2 + 2x + λx (degree 2) and ϑ(x) = x 2 (2x − 3) + λx (degree 3). The linear function λx (λ 1) is used to provide a small gradient ϑ (x) in the high pressure zone. e) Exponential function The lubricant fraction function can also be implemented by an exponential function according to ϑ(x) = 1 − e −k(x−x0) + λx. The following parameters have been chosen in this work: e1) k = 2, x 0 = 0; e2) k = 20, x 0 = 0.8; e3) k = 30, x 0 = 0.9; e4) k = 50, x 0 = 0.95. The linear function λx (λ 1) is used to provide a small gradient ϑ (x) in the high pressure zone. Using a 2-phase cavitation model, the compressible time-dependent REYNOLDS equation has to be used to calculate the pressure field p(φ,z). φ represents the circumferential coordinate,z represents the dimensionless axial coordinate of the bearing. H = 1 + D x sin(φ) − D y cos(φ) denotes the gap function, where D x and D y represent the dimensionless displacements of the journal center. The bearing has the radius R, the axial width B and the nominal radial gap size C. ω denotes the angular velocity of the rotor.
Here, a journal bearing with an axial groove is considered, see Fig. 1(c). The boundary conditions are shown in Fig. 1(d). At the axial groove, the DIRICHLET boundary conditions p = p in are applied, where p in denotes the oil inlet pressure. At the axial boundaries, special NEUMANN boudary conditions are used [2], namely R Hence, in the high pressure region (p ≥ p 0 ) ∂p ∂z 1, so that the pressure is almost the ambient pressure p 0 . In the cavitation region (p < p 0 ) ∂p ∂z 1, so that the axial flow is almost zero.

Non-mass-conserving cavitation models
To discuss the influence of fluid film cavitation on the rotordynamic behavior, we also consider non-mass conserving cavitation approaches. On the one hand, we use the well-known half-SOMMERFELD boundary conditions (GÜMBEL approach), where the incompressible REYNOLDS equation is used and the pressure below the ambient pressure is ignored. On the other hand, a penalty approach is applied. Therefore, the following REYNOLDS equation is considered with the penalty term Pen(p) = a e −b p−p 0 p 0 (a 1, b 1). This term may be interpreted as an artificial source term: fluid is provided into the cavitation zone so that p is always larger than p 0 . Note that the penalty approach produces smooth pressure functions at the cavitation boundaries, while the half-SOMMERFELD boundary conditions entire a non-smooth pressure profile.

Rotor/bearing co-simulation model
The rotor/bearing system is decomposed into several subsystems: one subsystem for the rotor (multibody system, MBS) and further subsystems for each fluid film (FE-systems). The MBS is coupled with the FE-subsystems by a co-simulation method [1] [2]. Here, an explicit sequential coupling approach (GAUSS-SEIDEL type) is used, where the MBS-subsystem is used as the master subsystem. A force/displacement decomposition technique is applied, i.e. output variables of the MBS subsystem are kinematical variables, whereas the output variables of the FE-subsystems are resultant forces/torques of the different fluid-films. The co-simulation is executed in the following way. Firstly, subsystem 1 (rotor, master) integrates from the macro-time point T n to T n+1 with the macro-time step H. For the numerical integration of the equations of motion of the MBS, the bearing forces/torques are required. Therefore, extrapolated bearing forces/torques are used. Secondly, the FE-subsystems (slaves) are integrated from T n to T n+1 . For the numerical integration of the FE-subsystem, kinematical quantities of subsystem 1 are required. Therefore, interpolated coupling variables are used.
Within the master-slave approach, a variable macro-time grid is used, which is defined by the step-size of the mastersubsystem. Here, the MBS is solved with a BDF solver (error = 1E − 6, H max = 1E − 6 s). The FE-subsystems are also solved with a BDF integrator (error = 1E − 5, H max = 1E − 6 s). The coupling variables are extrapolated and interpolated with quadratic polynomials.

Jeffcot rotor with large rotor mass
With three different cavitation models, run-up simulations have been performed for the case of lower external damping d a = 1 Ns/m (rotor stiffness c = 5000 N/mm) and higher external damping d a = 500 Ns/m (rotor stiffness c = 4000 N/mm). The rotor speed is increased linearly from 0 to 800 Hz in 2 s. Fig. 2(a) shows the dimensionless bearing eccentricities for the three bearing models (half-SOMMERFELD, penalty and 2-phase with ansatz function a)). As can be seen, the rotor systems become unstable (oil whirl/whip) at rotor speed of ≈ 200 Hz for the 2-phase model and ≈ 250 Hz for the half-SOMMERFELD and penalty model. The difference between the two non-mass conserving approaches is rather small. The model with the 2-phase approach becomes unstable earlier. Fig 2(b) shows the results of the simulations with higher external damping. The vibration of the rotor for all cavitation models can be divided into three phases. At the beginning, the rotor is stable and shows synchronous oscillations due to unbalance. With increasing rotor speed, the oil whirl/whip occurs with higher bearing eccentricities. Due to the higher external damping, the instability of the system can be passed through. At ≈ 1400 ms, the system becomes stable again and we only observe synchronous oscillations due to rotor unbalance. It can be detected that the 2-phase model becomes unstable earlier. The 2-phase model also shows higher bearing eccentricities in the unstable oil whirl/whip region.

Jeffcot rotor with small rotor mass
In order to investigate the influence of different 2-phase cavitation models, run-up simulations have been performed with the light JEFFCOT rotor. The rotor speed is increased linearly from 0 to 2400 Hz in 3 s. A comparison for eight different ansatz functions (see section 1.1) is shown in Fig. 2(c) and (d). In all simulations, an unstable oil whirl/whip region can be observed in the middle speed range. The basic oscillation behavior is quite similar. Two cavitation models, namely with ansatz b) and ansatz e1), predict however a significantly smaller oil whirl/whip region. These two models obviously entail a larger bearing damping. The simulation time with ansatz e1) is the most efficient one with a simulation time of ≈ 5 hours. The largest simulation time (≈ 187 hours) shows the model with ansatz a).

Turbocharger model
Finally, a turbocharger rotor is considered, which is supported with full-floating ring bearings. Run-up simulations have been carried out with three different bearing models: penalty approach, 2-phase model with ansatz a) and ansatz e1). The run-up speed is increased linearly from 0 to 1200 Hz in 3 s. Vibration of the compressor wheel in vertical direction is shown in www.gamm-proceedings.com