Nonlinear dynamic of flexible rotor supported by turbulent bearings with quadratic damping and temperature dependent viscosity

Abstract

A nonlinear dynamic model of a rotor supported by the turbulent journal bearing with quadratic damping under nonlinear suspension is presented in this study, and the lubricating oil is assumed to be temperature-dependent. The bifurcation diagrams, dynamic trajectories, power spectra, Poincaré map, Lyapunov exponent, and fractal dimension are used to verify the dynamic characteristics of the rotor-bearing system, and abundant subharmonic, and even chaotic motions are found. Temperature is one of the key parameters affecting the dynamic responses of the system, and the dynamic responses will perform non-periodic or even chaotic motions with higher working temperatures of lubricant fluid, as proved in this study. Nonlinear dynamic simulation analysis results will help engineers avoid unnecessary trial and error when designing and applying rotating machinery systems.

Keywords

Temperature-dependent viscosity quadratic damping turbulent bearing nonlinear dynamics

Introduction

Turbulent flow within lubrication systems in bearing assemblies can compromise stability at high rotational speeds and, in some cases, lead to failure. While rotational speed is a primary factor influencing Reynolds number, viscosity also plays a crucial role in initiating turbulent flow. Given that viscosity varies with temperature, relying on a constant viscosity value for bearing performance analysis can lead to significant inaccuracies. This issue is particularly relevant in industrial applications where precise bearing performance is critical. Consequently, some literature has concentrated on examining bearing systems in the context of turbulent effects to better understand and address these challenges. Soni et al.¹ modified Navier-Stokes equations with the turbulent effect and non-Newtonian lubrication by finite element method to solve. They found the combined effect of turbulent and non-Newtonian lubricant is to reduce the Sommerfield number as the Reynolds number and nonlinear factor increase. Betti et al.² applied the computational fluid dynamic method to solve the performance of the turbulent flow in tilting pad journal bearings, and the experiment also checked the numerical solutions. Xing et al.³ presented a turbulent lubrication model to predict the fluid pressure distribution in journal-axial coupled hydrodynamic bearing and the lubrication characteristics. Yang et al.⁴ proposed a numerical model to study the lubrication performances of coupled journal-thrust water-lubricated bearing with micro grooves. The results show that the turbulence effect can improve lubrication performances by enlarging the coupled journal-thrust bearing’s fluid pressure and load capacity. The laminar flow assumption no longer applies to the actual operating conditions. Zhu et al.⁵ studied the mixed lubrication performance of worn journal bearings supplied with axial end oil under turbulent effects. The simulation results show larger wear defects can significantly reduce the bearing’s maximum mean film pressure and friction power. Ebrahimi et al.⁶ performed a nonlinear dynamic analysis of a magnetically supported coaxial rotor in auxiliary bearings. The effects were investigated, including the inner shaft stiffness, inter-rotor bearing stiffness, auxiliary bearings stiffness, and disk position. Abundant periodic, subharmonic, quasi-periodic, and even chaotic motions were found. Ebrahimi⁷ also verifies an unbalanced flexible rotor model in active magnetic bearings through experiments. Using the Runge-Kutta method for motion integration, it incorporates gyroscopic moments, magnetic actuator coupling, and backup bearing contact forces. Nonlinear dynamic responses are analyzed via bifurcation plots, trajectories, Fast Fourier Transforms, Poincaré maps, and Lyapunov exponents across various speeds and eccentricities. Tests with a 5 g mass at different disk positions reveal complex phenomena like periodic, quasi-periodic, and chaotic motions. Results show that increasing rotational speed can stabilize the rotor, while higher unbalanced eccentricity requires more speed for stable operation, highlighting the need for careful design considerations.

As operating speed increases, the temperature of rotating machinery also rises. The viscosity of lubricating fluids in turbomachinery is highly sensitive to temperature variations. If temperatures fall too low, the machinery may fail to start, while excessively high temperatures can lead to lubrication failure. To gain a more comprehensive understanding of the system’s dynamic characteristics, it is essential to account for temperature effects in the analysis of turbomachinery dynamics. Taylor⁸ proposed an approximate model for short bearings including lubricant shear thinning, a non-Newtonian model. He also noted that modern automotive lubricants are multi-grade oils containing polymer additives at a few percent treatment rates, which will result in a significant dependence of lubricant viscosity on shear rate. Bukovnik et al.⁹ propose a modified thermal-elasto-hydrodynamic lubrication model for journal bearing, considering shear rate-dependent viscosity. Ighil et al.¹⁰ proposed the hydrodynamic performance analysis of the finite-length bearing and the short-bearing hypothesis, in which a global thermal study was also carried out to evaluate the effect of temperature. Maria-Eleni and Pantelis¹¹ presented the computational fluid dynamics (CFD) algorithm to investigate the growing pressure and temperature rise inside the bearing during operation. Baineni et al.¹² analyzed the performance of the hydrodynamic plain journal bearing combined with pressure, temperature, roughness, and cavitation effects. They prove that thermoelastohydrodynamic analysis of journal bearing is essential for better estimation of bearing performance. Sogo et al.¹³ found that the shear heat of the oil film causes a significant temperature distribution on the journal surface as rotational speeds increase. They also applied a simple unsteady model to efficiently analyze the oil film’s temperature, viscosity, and pressure distribution. Bagul et al.¹⁴ proposed journal bearing models developed for different speeds and eccentricity ratios to study the interaction between the fluid and the elastic behavior of the bearing. A thermohydrodynamic analysis of a circular journal bearing was also simulated using computational fluid dynamics methods. Kadam et al.¹⁵ discussed the decrease in viscosity with increasing temperature, explored through viscosity coupling. The dynamic water pressure field is obtained by solving the generalized Reynolds equation, and the finite element method is used to solve this equation numerically. The finite difference method is used in three-dimensional energy equations to predict the temperature distribution in fluid films. Torgal et al.¹⁶ derive an alternative method for finding the journal-bearing operating temperature. This article introduces some issues of hydrodynamic journal bearings based on thermohydrodynamic analysis, with a special focus on calculating bearing operating temperature based on a program developed for numerical simulation of bearing operation. Vaithianathan et al.¹⁷ thought different input parameters such as viscosity, speed, and shape (L/D ratio) effect on thermo-hydrodynamic lubrication model (i.e., Effect of heat on the viscosity is included in the computation) are carried out to arrive the performance parameters such as pressure developed in the bearing while in operation by using Ansys Fluent.

This study introduces a novel approach by analyzing the nonlinear dynamics of a flexible rotor supported by a turbulent journal bearing with quadratic damping, incorporating both nonlinear suspension and temperature-dependent viscosity. It uniquely addresses the combined effects of nonlinear suspension, nonlinear oil film force, and quadratic damping. Advanced analytical tools, such as bifurcation diagrams, dynamic trajectories, power spectra, Poincaré maps, fractal dimension, and Lyapunov exponents, are employed to comprehensively verify the dynamic behaviors of the rotor-bearing system. This innovative approach has significant applications in designing and optimizing rotor-bearing systems in various industrial settings, ensuring better performance and reliability.

Mathematical modeling

Temperature-dependent viscosity

It is well known that the viscosity coefficient varies with temperature. Assuming viscosity to be constant may be reasonable for low-speed applications or for a broader temperature range. However, for most modern rotating machinery, which operates at high speeds and over a wide temperature range, it is more appropriate to consider viscosity as temperature-dependent rather than constant. The temperature-dependent viscosity may be referred to by Taylor⁸ and could be performed as

μ = μ_{\infty} + \frac{μ_{0} - μ_{\infty}}{1 + \frac{γ}{γ_{c}}} = r μ_{0} + \frac{μ_{0} - r μ_{0}}{1 + \frac{γ}{γ_{c}}}

(1)

where $lo g_{10} γ_{c} = A + BT$ , $μ_{\infty} = r μ_{0}$ (the ratio r is defined in the region of $0.5 < r < 1.0$ ), $μ_{0} = ϑ \exp (\frac{T_{1}}{T_{2} + T})$ , $h = c (1 + ε \cos θ)$ , $γ = \frac{U}{h} = \frac{R ω}{c (1 + ε \cos θ)}$ . For SAE-15w/40, $r = 0.8$ , $A = 2.5$ , $B = 0.0225$ , $ϑ = 0.0388$ , $T_{1} = 1372.5 ° C$ , $T_{2} = 134.8 ° C$ . $lo g_{10} γ_{c} = 2.5 + 0.0225 T$ and $μ_{\infty} = 0.8 μ_{0}$ . Thus,

μ_{0} (T) = 0.0388 \exp (\frac{1372.5}{134.8 + T})

(2)

The temperature-dependent viscosity can be introduced by getting equation (2) into (1), as shown below

μ (T, θ) = 0.8 μ_{0} (T) + \frac{0.2 μ_{0} (T)}{1 + \frac{\frac{R ω}{c (1 + ε \cos θ)}}{10^{2.5 + 0.0225 T}}}

(3)

From equation (3), we know that viscosity is a function of both temperature (T) and the angle θ.

Journal bearings with turbulent effect

According to the hydrodynamic lubrication theory, the dynamic, turbulent modified Reynolds equation may be presented as the following type.

\frac{\partial}{\partial x} (\frac{h^{3} G_{θ}}{μ (T, θ)} \frac{\partial p}{\partial x}) + \frac{\partial}{\partial z} (\frac{h^{3} G_{z}}{μ (T, θ)} \frac{\partial p}{\partial z}) = \frac{U}{2} \frac{\partial h}{\partial x} + \frac{\partial h}{\partial t}

(4)

Where p is the pressure distribution, h is the film thickness, $h = c (1 + ε \cos θ)$ and $U = R ω$ , and $ω$ is the rotational speed. Then let’s operate the right-hand side terms of equation (4), that is, the first order derivative of the film thickness by location x and time t. Let $dx = Rd θ$

{\begin{matrix} \frac{\partial h}{\partial x} = \frac{\partial h}{R \partial θ} = - \frac{c ε}{R} \sin θ \\ \frac{\partial h}{\partial t} = c \overset{\cdot}{ε} \cos θ + c ε \overset{\cdot}{φ} \sin θ \end{matrix}

(5)

$G_{θ}$ and $G_{z}$ can be defined according to Constantinescu et al.¹⁸ indicated and the values are:

{\begin{matrix} G_{θ} = \frac{1}{12 + 0.0260 R {e^{*}}^{0.8265}} \\ G_{z} = \frac{1}{12 + 0.0198 R {e^{*}}^{0.741}} \end{matrix}

(6)

where $R e^{*}$ is local Reynolds number $R e^{*} = \frac{ρ Uh}{μ (T, θ)}$ and $μ (T, θ)$ is given in equation (3). In equation (6), when $G_{θ} = \frac{1}{12}$ and $G_{z} = \frac{1}{12}$ is supposed to be laminar conditions.

Introducing the rotating coordinate systems and taking into consideration the unsteady state of the journal rotation, that is, $\frac{\partial h}{\partial x} = - \frac{c ε}{R} \sin θ$ and $\frac{\partial h}{\partial t} = c \overset{\cdot}{ε} \cos θ + c ε \overset{\cdot}{φ} \sin θ$ , and then the modified turbulent Reynolds equation will be

\begin{matrix} \frac{\partial}{R^{2} \partial θ} (\frac{c^{3} {(1 + ε \cos θ)}^{3} G_{θ}}{μ (T, θ)} \frac{\partial p}{\partial θ}) + \frac{\partial}{\partial z} (\frac{c^{3} {(1 + ε \cos θ)}^{3} G_{z}}{μ (T, θ)} \frac{\partial p}{\partial z}) \\ = - 6 ω c ε \sin θ + 12 (c \overset{\cdot}{ε} \cos θ + c ε \overset{\cdot}{φ} \sin θ) \end{matrix}

(7)

Using the “short bearing approximation” (i.e., $\frac{L}{D} < 0.25$ , $\frac{\partial p}{\partial θ} < < \frac{\partial p}{\partial z}$ and $\frac{\partial p}{\partial θ} = 0$ )

Then, the pressure distribution in the circumferential direction $(\frac{\partial p}{\partial θ})$ could be neglected. Equation (7) could be modified as the following equation (8).

\frac{\partial^{2} p}{\partial z^{2}} = \frac{- 6 μ (T, θ) ω c ε \sin θ + 12 μ (T, θ) (c \overset{\cdot}{ε} \cos θ + c ε \overset{\cdot}{φ} \sin θ)}{c^{3} {(1 + ε \cos θ)}^{3} G_{z}}

(8)

Accompanying with the B. Cs. ${\begin{matrix} \frac{\partial p}{\partial z} = 0, z = 0 \\ p = 0, z = \pm \frac{L}{2} \end{matrix}$ and equation (3), that is, temperature-dependent viscosity introduced in the section “Temperature-dependent viscosity,” the pressure distribution will be

\begin{matrix} p (θ, z, t) = - \frac{3 μ (T, θ)}{c^{2} {(1 + ε \cos θ)}^{3} G_{z}} \\ [(ω - 2 \overset{\cdot}{φ}) ε \sin θ - 2 \overset{\cdot}{ε} \cos θ] (z^{2} - \frac{L^{2}}{4}) \end{matrix}

(9)

The resulting forces in the radial ( $f_{r}$ ) and tangential ( $f_{t}$ ) direction could be formulated as

{\begin{matrix} f_{r} = \int_{0}^{π} \int_{- \frac{L}{2}}^{\frac{L}{2}} p R \cos θ dzd θ \\ f_{t} = \int_{0}^{π} \int_{- \frac{L}{2}}^{\frac{L}{2}} p R \sin θ dzd θ \end{matrix}

(10)

and $f_{e} = - f_{r}, f_{φ} = - f_{t}$ , so

{\begin{matrix} f_{e} = \frac{L^{3} R}{2 c^{2}} \int_{0}^{π} \frac{μ (T, θ) [(ω - 2 \overset{\cdot}{φ}) ε \sin θ - 2 \overset{\cdot}{ε} \cos θ]}{G_{z} {(1 + ε \cos θ)}^{3}} \cos θ d θ \\ f_{φ} = - \frac{L^{3} R}{2 c^{2}} \int_{0}^{π} \frac{μ (T, θ) [(ω - 2 \overset{\cdot}{φ}) ε \sin θ - 2 \overset{\cdot}{ε} \cos θ]}{G_{z} {(1 + ε \cos θ)}^{3}} \sin θ d θ \end{matrix}

(11)

Where $G_{z} = \frac{1}{12 + 0.0198 R {e^{*}}^{0.741}} = \frac{1}{12 + 0.0198 {(\frac{ρ Uh}{μ (T, θ)})}^{0.741}}$ and it means that $G_{z}$ is also function of temperature T and angle $θ$ . The temperature-dependent lubricating forces in the vertical ( $F_{y}$ ) and horizontal ( $F_{x}$ ) direction could be summarized as

{\begin{matrix} F_{x} = f_{e} \cos φ - f_{φ} \sin φ \\ F_{y} = f_{e} \sin φ + f_{φ} \cos φ \end{matrix}

(12)

Dynamic equations for turbulent journal bearings with quadratic damping and temperature-dependent viscosity

A rigid rotor equipped with the turbulent journal bearing under nonlinear suspension as illustrated in Figure 1. Furthermore, the lubricating oil is assumed to be temperature dependent and quadratic damping is also considered for the turbulent journal bearing in this study. Firstly, the dynamic equations considering with nonlinear suspension for the bearing geometric center $O_{b}$ ( $X_{b}$ , $Y_{b}$ ) can be presented by the Newton’s second law.

{\begin{matrix} M_{b} {\overset{\cdot\cdot}{X}}_{b} + C_{b} \overset{\cdot}{X_{b}} + K_{1} X_{b} + K_{2} {X_{b}}^{3} - F_{x} = 0 \\ M_{b} \overset{\cdot\cdot}{Y_{b}} + C_{b} \overset{\cdot}{Y_{b}} + K_{1} Y_{b} + K_{2} {Y_{b}}^{3} - F_{y} = - M_{b} g \end{matrix}

(13)

where $M_{b}$ is the mass of the bearing and $C_{b}$ is the damping coefficient of the bearing. $K_{1}$ and $K_{2}$ are the first and third order coefficient of the stiffness respectively to represent nonlinear suspension effect.

Figure 1.

Rotor-bearing system equipped with nonlinear suspension and considering temperature dependent viscosity and quadratic damping.

Secondly, the dynamic equations considering quadratic damping for the geometric center of the rotor, $O_{r} (X_{r}$ , $Y_{r}$ ) can also be performed

{\begin{matrix} M_{r} {\overset{\cdot\cdot}{X}}_{r} + C_{r} {\overset{\cdot}{X}}_{r} | {\overset{\cdot}{X}}_{r} | + K_{s} (X_{r} - X_{j}) = M_{r} ρ ω^{2} \cos ϕ \\ M_{r} \overset{\cdot\cdot}{Y_{r}} + C_{r} \dot{Y_{r}} | \dot{Y_{r}} | + K_{s} (Y_{r} - Y_{j}) = M_{r} ρ ω^{2} \sin ϕ - M_{r} g \end{matrix}

(14)

where $X_{j} and Y_{j}$ are used to define the journal geometric center $O_{j}$ , $M_{r}$ is the mass of the rotor, $C_{r}$ is the damping coefficient of the shaft, $K_{s}$ is the stiffness of the shaft $ρ$ is the mass eccentricity of the rotor, $ω$ is the rotational speed of the rotor and $ϕ$ is rotational angle (ϕ = ωt). The term $C_{r} \dot{Y_{r}} | \dot{Y_{r}} |$ is used to simulate as quadratic damping.

Finally, the nonlinear oil film forces in the horizontal and vertical directions can be expressed from the force equilibrium in each direction with equation (12).

{\begin{matrix} F_{x} = f_{e} \cos φ - f_{φ} \sin φ = K_{s} (X_{r} - X_{j}) \\ F_{y} = f_{e} \sin φ + f_{φ} \cos φ = K_{s} (Y_{r} - Y_{j}) \end{matrix}

(15)

Rearranging equation (15) and we could get the following equations (See Appendix A).

{\begin{matrix} \overset{\cdot}{φ} = \frac{(Ψ_{7} - Ψ_{0} Ψ_{1}) Ψ_{0} Ψ_{6} - (Ψ_{8} - Ψ_{0} Ψ_{4}) Ψ_{0} Ψ_{3}}{{Ψ_{0}}^{2} Ψ_{2} Ψ_{6} - {Ψ_{0}}^{2} Ψ_{3} Ψ_{5}} \\ \overset{\cdot}{ε} = \frac{(Ψ_{8} - Ψ_{0} Ψ_{4}) Ψ_{0} Ψ_{2} - (Ψ_{7} - Ψ_{0} Ψ_{1}) Ψ_{0} Ψ_{5}}{{Ψ_{0}}^{2} Ψ_{2} Ψ_{6} - {Ψ_{0}}^{2} Ψ_{3} Ψ_{5}} \end{matrix}

(16)

Equations (13), (14), and (16) are the main governing equations used to simulate the nonlinear dynamic responses of the rotor-bearing system considering temperature-dependent viscosity and some linear effects, that is, nonlinear restoring force, nonlinear oil-film force quadratic damping force. Then, dimensionless equations could be rearranged.

{\begin{matrix} x_{b}^{,,} + \frac{2 c_{1}}{a_{1}} x_{b}^{,} + \frac{1}{{a_{1}}^{2}} x_{b} + \frac{α}{s^{2}} {x_{b}}^{3} - \frac{1}{2 r_{br} s^{2}} (x_{r} - x_{b} - ε \cos φ) = 0 \\ y_{b}^{,,} + \frac{2 c_{1}}{a_{1}} y_{b}^{,} + \frac{1}{{a_{1}}^{2}} y_{b} + \frac{α}{s^{2}} {y_{b}}^{3} - \frac{1}{2 r_{br} s^{2}} (y_{r} - y_{b} - ε \sin φ) + \frac{f_{b}}{s^{2}} = 0 \end{matrix}

(17)

{\begin{matrix} x_{r}^{,,} + \frac{2 c_{2}}{s} x_{r}^{,} | x_{r}^{,} | + \frac{1}{s^{2}} (x_{r} - x_{b} - ε \cos φ) = β \cos ϕ \\ y_{r}^{,,} + \frac{2 c_{2}}{s} y_{r}^{,} | y_{r}^{,} | + \frac{1}{s^{2}} (y_{r} - y_{b} - ε \sin φ) = β \sin ϕ - \frac{f_{r}}{s^{2}} = 0 \end{matrix}

(18)

{\begin{matrix} φ' = \frac{(Ψ_{7} - Ψ_{0} Ψ_{1}) Ψ_{0} Ψ_{6} - (Ψ_{8} - Ψ_{0} Ψ_{4}) Ψ_{0} Ψ_{3}}{({Ψ_{0}}^{2} Ψ_{2} Ψ_{6} - {Ψ_{0}}^{2} Ψ_{3} Ψ_{5}) ω} \\ ε' = \frac{(Ψ_{8} - Ψ_{0} Ψ_{4}) Ψ_{0} Ψ_{2} - (Ψ_{7} - Ψ_{0} Ψ_{1}) Ψ_{0} Ψ_{5}}{({Ψ_{0}}^{2} Ψ_{2} Ψ_{6} - {Ψ_{0}}^{2} Ψ_{3} Ψ_{5}) ω} \end{matrix}

(19)

where $x_{i} = \frac{X_{i}}{c} (i = b, r j),$ c is the clearance of oil film, $ϕ = ω t$ , $\frac{d}{dt} = ω \frac{d}{d ϕ}$ , $Ω^{2} = \frac{K_{s}}{M_{r}}, s^{2} = \frac{ω^{2}}{Ω^{2}}$ , $β = \frac{M_{r} ρ ω^{2}}{M_{r} c ω^{2}}$ , $f_{b} = \frac{M_{r} g}{K_{s} δ}$ , $c_{1} = \frac{C_{b}}{2 \sqrt{K_{1} M_{b}}}$ , $c_{2} = \frac{C_{r}}{2 \sqrt{K_{s} M_{r}}}$ , $R_{m} = \frac{M_{b}}{M_{r}}$ , $R_{k} = \frac{K_{s}}{K_{1}}$ , ${a_{1}}^{2} = R_{m} R_{k} a^{2}$ , $α = \frac{K_{2} c^{2}}{K_{s} R_{m}}$ .

Results and discussions

The dimensionless equations (17)–(19) describe the nonlinear dynamic characteristics of the rotor-bearing system under turbulent lubrication effect with the quadratic damping and temperature-dependent viscosity assumptions, and the fourth-order Runge-Kutta method is therefore used for solving the autonomous differential equations. The time step used for calculation in this study was pi/300, and the allowable error value was less than 0.0001. To ensure that the data used reaches a steady state, the time series data of the first 800 revolutions will not be used for dynamic behavior analysis. In addition, 60,000 data points from a time series of displacements of the rotor and bearing geometric centers will be used to generate attractors embedded in space. The range of the working temperature T is set to be 25°C ∼ 75 °C which is in the range of the Taylor,⁸ and the range of dimensionless rotating speed ratio s is set to be 0.01–2.76. The definition of local Reynolds number $R e^{*} = \frac{ρ Uh}{μ (T, θ)}$ is given according to Constantinescu et al.¹⁸ and $G_{θ} = \frac{1}{12}$ and $G_{z} = \frac{1}{12}$ is supposed to be laminar conditions. The values for the dimensionless parameters are set to be $c_{1} =$ 0.05, $R_{m} =$ 0.2, $R_{k} =$ 2.0, $β =$ 0.5, $f_{b} =$ 0.2, $c_{2} =$ 0.1, $α =$ 1.0.

In the beginning, the corresponding variations of the bearing geometric center ( $x_{b}$ , $y_{b}$ ) and rotor geometric center ( $x_{r}$ , $y_{r}$ ) coordinates of the return points in the Poincaré map are then plotted to form the bifurcation diagram. The bifurcation diagrams could provide useful information to observe the system’s overall dynamic behavior changes as the control parameters increase. To produce a bifurcation diagram for a dynamical system, the following steps are typically followed: First, identify the equations and parameters of the system, then select a control parameter. Next, numerical simulations are performed to obtain the stable states or periodic solutions under different parameter values. Finally, plot the results of the stable solutions as the parameter varies, thus creating the bifurcation diagram, which allows for a deeper exploration of the changes in the system’s behavior. Let’s use a simple example to illustrate how to generate a bifurcation diagram using mathematical formulas. Suppose we consider a simple dynamical system, such as the one-dimensional logistic equation, that is, $x_{n + 1} = r x_{n} (1 - x_{n})$ , and then we choose the parameter range, such as $1.2 < r < 5$ . We may choose an initial condition x₀ to perform multiple iterations (e.g., 1000 iterations) and plot the final stable state.

This study applies two control parameters, namely, the dimensionless rotational speed ratio, s, and the dimensionless unbalance coefficient β, to generate bifurcation diagrams. In each case, the bifurcation control parameter is varied with a constant step, and the state variables at the end of one integration step are taken as the initial values for the next step. In Figures 2 and 3, the bifurcation diagram is for the bearing geometric center and the rotor geometric center with dimensionless rotational speed ratio s (s = 0.01−2.76) at temperature T = 30°C, respectively. The dynamic responses perform almost the same behaviors as the control parameters, s increase for the rotor and the bearing. We found they all behave as periodic or 2T-period motions for the whole bifurcation parameters at temperature T = 30°C. However, as the working temperature rises to T = 75°C, the dynamic responses become more irregular, as shown in Figures 4 and 5. Comparing to the low-temperature case (T = 30°C), it is clear that the dynamic responses only perform periodic motions for low rotating speeds, that is, s < 0.87 and then they will be irregular dynamic motions or even chaotic motions for high speeds for the high-temperature case (T = 75°C). Temperature is apparently one of the key parameters affecting the dynamic responses of the system, and the dynamic responses will perform non-periodic or even chaotic motions with higher working temperatures of lubricant fluid. We selected six sets of data (s = 0.4, 0.6, 0.8, 1.6, 2.0, and 2.76) to describe the dynamic responses of the system, which exhibit periodic motions in the range 0.01 ≦ s ≦ 0.87, and display chaotic behavior when 0.87 < s < 2.76.

Figure 2.

Bifurcation diagram of the bearing geometric center using dimensionless rotational speed ratio s (s = 0.01–2.76) for temperature T = 30°C.

Figure 3.

Bifurcation diagram of the rotor geometric center using dimensionless rotational speed ratio s (s = 0.01–2.76) for temperature T = 30°C.

Figure 4.

Bifurcation diagram of the bearing geometric center using dimensionless rotational speed ratio s(s = 0.01–2.76) for temperature T = 75°C.

Figure 5.

Bifurcation diagram of the rotor geometric center using dimensionless rotational speed ratio s (s = 0.01–2.76) for temperature T = 75°C.

To more carefully define the dynamic characteristics under specific control parameters, the dynamic orbit, the Poincaré map, the power spectrum, the maximum Lyapunov exponent, and fractal dimension of the bearing geometric at dimensionless rotational speed ratio s = 1.8 and 2.0 for temperature T = 75°C as illustrated in Figures 6 and 7. We found that the dynamic response is not so regular. We can understand the system’s dynamic behavior through the distribution of mapping points in the Poincaré section. When the mapping point on the Poincaré section is a fixed point or a small number of discrete points, the system exhibits periodic motion. The system exhibits quasi-periodic motion when a closed curve appears in the Poincare cross-section. The system motion appears chaotic when dense points or fractal geometry appear in the Poincaré section. The dynamic responses are chaotic for s = 1.8 and 2.0 by finding some dense points in the Poincaré map. Then, numerous excitation frequencies could also be found in power spectra for s = 1.8 and 2.0, which means that the dynamic characteristics for the bearing geometric center show non-periodic vibrations at those rotating speeds. Further confirmation is carried out using the fractal dimension and maximum Lyapunov exponent to confirm whether the system enters the chaotic state. For an n-dimensional continuous dynamic system defined by a differential equation, whether it is an autonomous system or a non-autonomous system, for example, ${\overset{\cdot}{X}}_{i} = F_{i} (X_{1}, X_{2}, \dots, X_{n}, t)$ . Assuming that in the phase plane space trajectory of the initial state is $X_{0}$ and its nearby perturbation initial state is $X_{0} + Δ X_{0}$ . If the distance between the two initial states in the direction of the i-th dimension is $p_{i} (t)$ , then the Lyapunov exponent in the direction of the i-th dimension could be defined as $λ_{i} = lim_{t \to \infty} \frac{1}{t} (\ln \frac{p_{i} (t)}{p_{i} (0)})$ . It can be known from the above definition that it represents the degree of exponent convergence or divergence of the average distance between adjacent trajectories in the direction of the i-th dimension over a period of time t. Assume that the function F satisfies all the conditions for the existence of a unique solution to these two equations, and assume that its solution set in the phase space is $X (t) = Φ (t, t_{0}, X_{0})$ , where $Φ : R^{n + 2} \to R^{n}$ . By perturbating this dynamic trajectory in the state $X_{0}$ , the Jacobin matrix near state $X_{0}$ can be obtained as $M (t) = \frac{\partial Φ (t, t_{0}, X_{0})}{\partial X_{0}}$ and it satisfies $\frac{\partial M (t)}{\partial t} = \frac{\partial F}{\partial X} M (t)$ , in which $\frac{\partial F}{\partial X}$ is the $n \times n$ order Jacobin matrix of the vector function F. The change of perturbation quantity with time can be expressed as $Δ X (t) = M (t) Δ X_{0}$ . Let $τ$ to be the sampling time and then $t_{k} = k τ, (k = 0, 1, 2, \dots$ ). The Lyapunov exponent in the direction of the i-th dimension could be re-defined as $λ_{i} = lim_{n \to \infty} \frac{1}{n τ} \sum_{k = 1}^{n} \ln | M_{i} (t_{k}) |$ . For systems with chaotic motion, there is at least one positive Lyapunov exponent. Therefore, as long as it is known from the calculation that there is at least one positive Lyapunov exponent in the attractor, it can be determined that the system is a strange attractor and shows chaotic motion. For conservative or dissipative systems, as long as there is a positive Lyapunov exponent, it can be determined that the system’s motion exhibits the chaotic motion of a singular attractor. Therefore, in order to determine whether the system is in chaotic motion, we only need to find the maximum Lyapunov exponent of the system. It can be judged whether it is a chaotic motion by looking at whether it is positive. If the correlation integral is c(r), and $c (r) = lim_{N \to \infty} \frac{1}{N^{2}} \sum_{i, j = 1}^{N} H (r - | x_{i} - x_{j} |)$ , in which $H (x) = {\begin{matrix} 1, x > 0 \\ 0, x \leq 0 \end{matrix}$ is the Heavside function. When the distance, $| x_{i} - x_{j} |$ between the two states $x_{i}$ and $x_{j}$ in the phase space is less than the given value ( $r)$ , that is, $H (x) \neq 0$ , which means that the two points are related to each other and they contribute to the value of c(r); vice versa;, if the distance $| x_{i} - x_{j} |$ between $x_{i}$ and $x_{j}$ is greater than the given value of r, it means that the two points are unrelated to each other. The fractal dimension of the system can be determined by using the plot of the log value of the correlation function versus the log value of the radius of an N-dimensional hyper-sphere, that is, (log c(r)) versus (log r). As for the fractal dimension, the greatest integer M, which is less than the fractal dimension of the attracting set, and the number of points used to estimate the size should be less than 42^M. In short words, to calculate the Lyapunov exponent, you can typically follow these steps:

Initial condition: Select an initial condition x₀ and determine a small perturbation ε.

Iteration: As time progresses, iteratively compute the evolution of the system, recording the states x_n and the perturbed states x_n + ε.

Calculate perturbation growth: Compute the growth factor of the perturbation and calculate the limit of the Lyapunov exponent for sufficiently large n, that is, $λ_{n} = lim_{n \to \infty} \frac{1}{n} (\ln \frac{[x_{n} + ε - x_{n}]}{[ε]})$ .

Figure 6.

Dynamic orbit, Poincaré map, power spectrum, maximum Lyapunov exponent, and fractal dimension of the bearing geometric center at dimensionless rotational speed ratio s = 1.8 for temperature T = 75°C: (a) Dynamic orbit, (b) Poincaré map, (c) power spectrum, (d) maximum Lyapunov exponent and (e) fractal dimension.

Figure 7.

Dynamic orbit, Poincaré map, power spectrum, maximum Lyapunov exponent, and fractal dimension of the bearing geometric center at dimensionless rotational speed ratio s = 2.0 for temperature T = 75°C: (a) Dynamic orbit, (b) Poincaré map, (c) power spectrum, (d) maximum Lyapunov exponent and (e) fractal dimension.

It is observed that the maximum Lyapunov exponents are all positive for the bearing and rotor geometric center s = 1.8 and 2.0 for temperature T = 75°C, as shown in Figures 6 –9, so we could say chaotic motions existing at the above parameters. The fractal dimension is 1.38, 1.41, 1.22, and 1.37 sequentially for Figures 6 –9. Therefore, we prove again that the system motion appears chaotic at those parameters.

Figure 8.

Dynamic orbit, Poincaré map, power spectrum, maximum Lyapunov exponent, and fractal dimension of the rotor geometric center at dimensionless rotational speed ratio s = 1.8 for temperature T = 75°C: (a) Dynamic orbit, (b) Poincaré map, (c) power spectrum, (d) maximum Lyapunov exponent and (e) fractal dimension.

Figure 9.

Dynamic orbit, Poincaré map, power spectrum, maximum Lyapunov exponent, and fractal dimension of the rotor geometric center at dimensionless rotational speed ratio s = 2.0 for temperature T = 75°C: (a) Dynamic orbit, (b) Poincaré map, (c) power spectrum, (d) maximum Lyapunov exponent and (e) fractal dimension.

The unbalanced parameter, β, is one of the control parameters used to generate bifurcation diagrams for observing dynamic responses in rotating machinery. Figure 10 presents the bifurcation diagram of the bearing’s geometric center using the dimensionless unbalanced parameter β (ranging from 0.001 to 0.345) at a temperature of 25°C. For clarity, dynamic trajectories and the Poincaré map are also provided for key parameters: β = 0.1, β = 0.25, β = 0.26, and β = 0.34. As β increases, the amplitude of the system generally increases, but the dynamic responses become less regular. Most of the dynamic responses exhibit nT-periodic motions for various values of β. However, exceptions occur at β = 0.25 and β = 0.34, where the dynamic trajectories reveal aperiodic motions, and the Poincaré map indicates chaotic behavior. Figure 11 is the bifurcation diagram of the rotor geometric center using dimensionless unbalanced parameter $β$ ( $β = 0.001 - ~ 0.345$ ) for temperature T = 25°C. The dynamic responses are almost the same as the case of the bearing, except the amplitudes are a little larger than the bearing case. In the higher temperature case, that is, T = 75°C, various non-periodic motions are found for the bearing and rotor geometric center in Figures 12 and 13. The dynamic responses will return regularly subharmonic notions with 8T-period. Still, the amplitude will enlarge as $β$ increases. We take some key parameters into the discussion, for example, $β = 0.04$ , $β = 0.06$ , $β = 0.10$ , $β = 0.12$ , $β = 0.14$ are for chaotic motions and $β = 0.16$ , $β = 0.20$ , $β = 0.30$ are for nT-periodic motions. Once again, we found that temperature significantly impacts the dynamic stability of the system. Therefore, overlooking the effects of temperature could lead to misjudgments about the system’s motion state and result in unexpected vibrations or failures.

Figure 10.

Bifurcation diagram of the bearing geometric center using dimensionless unbalance parameter $β$ ( $β = 0.001 - ~ 0.345$ ) for temperature T = 25°C.

Figure 11.

Bifurcation diagram of the rotor geometric center using dimensionless unbalance parameter $β$ ( $β = 0.001 - ~ 0.345$ ) for temperature T = 25°C.

Figure 12.

Bifurcation diagram of the bearing geometric center using dimensionless unbalance parameter $β$ ( $β = 0.001 - ~ 0.345$ ) for temperature T = 75°C.

Figure 13.

Bifurcation diagram of the rotor geometric center using dimensionless unbalance parameter $β$ ( $β = 0.001 - ~ 0.345$ ) for temperature T = 75°C.

To compare the difference of nonlinear dynamic responses with and without quadratic damping effect, the bifurcation diagrams of the bearing and rotor geometric centers without quadratic damping effect using dimensionless rotational speed ratio s (s = 0.01–2.76) for temperature T = 75°C are presented in Figures 14 and 15. It is obviously that non-periodic responses only exist at the range 1.43 ≦ s ≦ 1.65, and they behave periodic or sub-harmonic for most of the rotating speed ratios. The results show that the ignorance of the quadratic damping effect for rotor-bearing systems with turbulent flow may cause misjudgment of the dynamic response in those systems and further lead to the failure of the system in design.

Figure 14.

Bifurcation diagram of the bearing geometric center without quadratic damping effect using dimensionless rotational speed ratio s (s = 0.01–2.76) for temperature T = 75°C.

Figure 15.

Bifurcation diagram of the rotor geometric center without quadratic damping effect using dimensionless rotational speed ratio s (s = 0.01–2.76) for temperature T = 75°C.

The findings of this study underscore the importance of considering nonlinear dynamics in rotor-bearing systems, particularly under varying operational conditions. The nonlinear dynamic analysis conducted reveals that factors such as quadratic damping and temperature-dependent viscosity significantly influence the behavior of the rotor system.

In line with these findings, Żywica et al.¹⁹ highlighted the critical role of dynamic interactions in the rotor-bearings-support structure system of a multi-stage organic Rankine cycle (ORC) microturbine. Their detailed analysis demonstrates how the integration of various components, including the high-speed rotor, bearings, and support structures, can reveal unfavorable dynamical phenomena that may compromise the reliability and performance of microturbines in combined heat and power (CHP) applications. This emphasizes the necessity of a comprehensive understanding of dynamic characteristics to prevent potential failures during operation.

Additionally, the work by Chen et al.²⁰ emphasizes the significance of nonlinear oil film behavior in sliding bearings under heavy loads. Their study provides a quadratic model for the nonlinear oil film force, which complements the present analysis by illustrating how nonlinearity impacts dynamic characteristic testing. They found that while nonlinearity is negligible at lower loads, it becomes pronounced under heavy loads, leading to substantial identification errors in stiffness and damping coefficients. This insight aligns with our results, where higher rotational speeds and increased lubrication temperatures also contribute to non-periodic or chaotic dynamic responses.

Collectively, these studies reinforce the notion that robust diagnostic systems are essential for monitoring dynamic behaviors in turbomachinery. By integrating the findings from Żywica et al.¹⁹ and Chen et al.²⁰ with our analysis, this research highlights the importance of incorporating nonlinear effects and dynamic interactions in the design and operational strategies of rotor machinery. This approach will ultimately enhance the reliability and efficiency of such systems in real-world applications.

Conclusions

This study presents a nonlinear dynamic analysis of a flexible rotor supported by bearings with turbulent lubrication under nonlinear suspension. The analysis incorporates the effects of quadratic damping and temperature-dependent viscosity. As the rotational speed increases to a critical value, the Reynolds number rises, resulting in turbulent flow and higher lubrication temperatures. Therefore, accounting for quadratic damping and temperature-dependent viscosity is crucial for a comprehensive nonlinear dynamic analysis of turbomachinery.

Bifurcation diagrams, dynamic orbits, Poincaré maps, power spectra, maximum Lyapunov exponents, and fractal dimensions are utilized to analyze the dynamic responses of the rotor-bearing system, revealing strong nonlinear effects, including nonlinear suspension, quadratic damping, and nonlinear oil film force. These factors contribute to various non-periodic motions and even chaotic behaviors. Notably, the temperature-dependent viscosity shows that it should not be assumed constant; higher lubricant temperatures can lead to non-periodic or chaotic dynamic responses.

The findings of this study have significant practical implications for engineers and designers in the field of rotor machinery. By understanding the impact of temperature and lubrication on dynamic behavior, professionals can better anticipate and mitigate undesirable vibrations in real-world applications. The simulation results provide valuable insights that can inform the design and operational strategies for turbomachinery, ultimately leading to enhanced reliability and performance in industrial applications. This study aims to guide future analyses and designs in rotor machinery and related fields, ensuring more stable and efficient operations.

Footnotes

Appendix A

Let $Ψ_{0} = - \frac{L^{3} R}{2 c^{2}}$ , $Ψ_{1} = - \int_{0}^{π} \frac{μ (T, θ) ω ε \sin θ \cos θ}{G_{z} {(1 + ε \cos θ)}^{3}} d θ$ , $Ψ_{2} = \int_{0}^{π} \frac{2 μ (T, θ) ε \sin θ \cos θ}{G_{z} {(1 + ε \cos θ)}^{3}} d θ$ , $Ψ_{3} = \int_{0}^{π} \frac{2 μ (T, θ) c o s θ c o s θ}{G_{z} {(1 + ε c o s θ)}^{3}} d θ$ , $Ψ_{4} = \int_{0}^{π} \frac{μ (T, θ) ω ε \sin θ \sin θ}{G_{z} {(1 + ε \cos θ)}^{3}} d θ$ , $Ψ_{5} = \int_{0}^{π} \frac{- 2 μ (T, θ) ε \sin θ \sin θ}{G_{z} {(1 + ε \cos θ)}^{3}} d θ$ , $Ψ_{6} = \int_{0}^{π} \frac{- 2 μ (T, θ) \cos θ \sin θ}{G_{z} {(1 + ε \cos θ)}^{3}} d θ$ and $f_{e}$ and $f_{φ}$ could be simplified as follows

(A1)

{\begin{matrix} f_{e} = Ψ_{0} (Ψ_{1} + Ψ_{2} \overset{\cdot}{φ} + Ψ_{3} \overset{\cdot}{ε}) \\ f_{φ} = Ψ_{0} (Ψ_{4} + Ψ_{5} \overset{\cdot}{φ} + Ψ_{6} \overset{\cdot}{ε}) \end{matrix}

Substituting equation (A1) into equation (15) and rearranging equations for the oil film forces considering with temperature effect in the horizontal and vertical direction. Meanwhile, Let $Ψ_{7} = K_{s} (X_{r} - X_{b} - e \cos φ) \cos φ + K_{s} (Y_{r} - Y_{b} - e \sin φ) \sin φ Ψ_{8} = K_{s} (Y_{r} - Y_{b} - e \sin φ) \cos φ - K_{s} (X_{r} - X_{b} - e \cos φ) \sin φ$ and then

(A2)

{\begin{matrix} Ψ_{0} Ψ_{2} \overset{\cdot}{φ} + Ψ_{0} Ψ_{3} \overset{\cdot}{ε} = Ψ_{7} - Ψ_{0} Ψ_{1} \\ Ψ_{0} Ψ_{5} \overset{\cdot}{φ} + Ψ_{0} Ψ_{6} \overset{\cdot}{ε} = Ψ_{8} - Ψ_{0} Ψ_{4} \end{matrix}

(A3)

{\begin{matrix} \overset{\cdot}{φ} = \frac{(Ψ_{7} - Ψ_{0} Ψ_{1}) Ψ_{0} Ψ_{6} - (Ψ_{8} - Ψ_{0} Ψ_{4}) Ψ_{0} Ψ_{3}}{{Ψ_{0}}^{2} Ψ_{2} Ψ_{6} - {Ψ_{0}}^{2} Ψ_{3} Ψ_{5}} \\ \overset{\cdot}{ε} = \frac{(Ψ_{8} - Ψ_{0} Ψ_{4}) Ψ_{0} Ψ_{2} - (Ψ_{7} - Ψ_{0} Ψ_{1}) Ψ_{0} Ψ_{5}}{{Ψ_{0}}^{2} Ψ_{2} Ψ_{6} - {Ψ_{0}}^{2} Ψ_{3} Ψ_{5}} \end{matrix}

Handling Editor: Sharmili Pandian

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Cai-Wan Chang-Jian

Data availability

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

References

Soni

Vakharia

DP.

Static analysis of finite hydrodynamic journal bearing in turbulent regime with non-Newtonian lubricant. Tribol Online 2015; 10: 246–161

Betti

Forte

Ciulli

Turbulence effects in tilting pad journal bearings: a review. Lubricants 2022; 10: 171.

Xing

Xue

Han

, et al. Turbulent lubrication model for journal-axial coupled hydrodynamic bearings without grooves. Tribol Int 2023; 190: 109036.

Yang

Cai

Wang

, et al. Numerical analysis of turbulence effect for coupled journal-thrust water-lubricated bearing with micro grooves. J Tribol 2023; 145(8): 084101.

Zhu

Zhang

Wang

, et al. Mixed turbulent lubrication analysis of worn journal bearing with end oil supply. Tribol Int 2024; 193: 109380.

Ebrahimi

Ghayour

Khanlo

HM.

Effects of some design parameters on bifurcation behavior of a magnetically supported coaxial rotor in auxiliary bearings. Eng Comput 2017; 34(7): 2379–2395.

Ebrahimi

Experimental and theoretical investigations of unbalance flexible rotor in active magnetic bearings considering backup bearings contacts. Noise Vibr Worldwide 2022; 53(4–5): 161–171.

Taylor

RI.

The inclusion of lubricant shear thinning in the short bearing approximation. Proc Inst Mech Eng 1999; 213: 35–46.

Bukovnik

Offner

Čaika

, et al. Thermo-elasto-hydrodynamic lubrication model for journal bearing including shear rate-dependent viscosity. Lubr Sci 2007; 19(4): 231–245.

10.

Ighil

Bounif

Maspeyrot

Thermo-hydrodynamic study of the journal bearing under static load. Proc IMechE C J Mech Eng Sci 2008; 222(9): 1801–1809.

11.

Maria-Eleni

Pantelis

GN.

Simulation of thermo-hydrodynamic behavior of journal bearings, lubricating with commercial oils of different performance. Simul Modell Pract Theory 2020; 104: 102128.

12.

Baineni

Borse

Chippa

Analysis of thermoelastohydrodynamic lubrication of journal bearing including the effect of surface roughness and cavitation. J Inst Eng India Ser D 2022; 3: 1. DOI: 10.1007/s40033-022-00335-z.

13.

Sogo

Inoue

Hatakenaka

Thermo-hydrodynamic lubrication analysis of journal bearings using a simplified THL model and prediction of onset speed of instability. Trans JSME (in Japanese) 2023; 89(924): 23–00068.

14.

Bagul

Patil

RY.

Thermal analysis of journal bearing using CFD software for performance enhancement. Int J Sci Technol Eng 2017; 3.

15.

Kadam

Banwait

Laroiya

SC.

Thermo hydrodynamic analysis of plain journal bearing with modified viscosity temperature equation. Int J Mech Eng Technol (IJMET) 2014; 5: 31.

16.

Torgal

Saini

Thermo-hydrodynamic analysis of journal bearing to find out equivalent temperature. Int J Sci Technol Eng 2015; 1(12): 201–207.

17.

Vaithianathan

Annamalai

. Ansys fluent analysis of a thermo-hydrodynamically lubricated journal bearing. In: Natarajan

Prakash

Sankaranarayanasamy

(eds.) Recent advances in manufacturing, automation, design and energy technologies. Lecture notes in mechanical engineering. Singapore: Springer, 2021, pp. 615–620.

18.

Constantinescu

Eng

Tech

M.S.

Theory of turbulent lubrication. Institution of Mechanical Engineers, 1959,173, pp 615–620.

19.

Żywica1

Breńkacz1

Bagiński

Interactions in the rotor-bearings-support structure system of the multi-stage ORC microturbine. J Vibr Eng Technol 2018; 6: 369–377.

20.

Chen

Tang

, et al. Influence of oil film nonlinearity on identification accuracy of dynamic characteristic coefficient of heavy-duty sliding bearing. J Brazil Soc Mec Sci Eng 2023; 45: 233.