Nonlinear analysis of rotor-stator dynamic system with gear transmission under water lubrication bearing

Abstract

This study presents a comprehensive nonlinear dynamic analysis of a water-lubricated bearing–gear–rotor system by incorporating the effects of turbulent mixed lubrication and system damping. A dimensionless model was developed and numerically solved using the fourth-order Runge–Kutta method. The dynamic responses under various dimensionless rotational speed ratios (s) and damping ratios ( $\xi$ ) were systematically investigated through bifurcation diagrams, dynamic orbits, Poincaré maps, power spectra, fractal dimension, and Lyapunov exponent analysis. The results reveal that the system exhibits chaotic behavior at low speed ratios (s<1.2) and low damping values ( $\xi <0.14$ ), characterized by non-periodic motion, scattered Poincaré return points, and a broad excitation spectrum. As speed or damping increases, the system transitions into periodic states, though intermittent chaos reappears within specific damping intervals. At $\xi =0.1,$ all components, that is, bearing, gear, and rotor, demonstrated chaotic dynamics, as confirmed by non-integer fractal dimensions (2.05, 1.60, and 1.97) and positive Lyapunov exponents (0.22, 0.12, and 0.16). These findings highlight the complex influence of speed and damping on the stability and performance of water-lubricated systems. The results provide valuable insights for the design and control of high-performance, low-viscosity fluid-lubricated machinery operating under nonlinear and potentially chaotic conditions.

Introduction

Since water is an environmentally friendly and non-toxic lubricant, water-lubricated bearings are widely used in fields with high environmental requirements, avoiding the pollution issues associated with oil-based lubricants. Water-lubricated bearings are extensively employed in marine applications, particularly for ship bearings in seawater, as water is a natural lubrication source. This reduces the need for oil-based lubricants and is well-suited for marine environments. Water-lubricated systems use water as the lubricant, eliminating the risk of contamination from oils and complying with safety standards in the food and pharmaceutical industries. Water-lubricated systems can also serve a dual purpose of lubrication and cooling, especially in high-temperature working environments, such as in high-performance motors and power generation equipment that require efficient heat dissipation. Accordingly, water-lubricated bearing systems are primarily applied in fields that require environmental protection, cleanliness, and special working conditions, owing to their advantages of being environmentally friendly, safe, low-cost, and reducing pollution.

Qiao et al. ¹ develop a turbulent mixed lubrication model for water-lubricated rubber bearings (WLRBs) based on an extended average Reynolds equation incorporating turbulence. Their experimental validation confirms the model’s effectiveness. The study reveals how turbulence affects mixed lubrication under varying load, rotational speed, axial misalignment, and radial clearance. Xie et al. ² present a theoretical and experimental investigation into the micro asperity contact load ratios and lubrication regime transitions in water-lubricated stern tube bearings. Experimental validation supports the accuracy of their revised model, providing valuable insights for defining contact load ratios and predicting lubrication behavior in such bearings. Zhu et al. ³ propose a generalized turbulent lubrication model for journal bearings with isotropic rough surfaces by integrating the Hashimoto turbulent model and the Patir–Cheng model for laminar flow. Using the finite difference method to solve both lubrication and energy equations, they analyze how surface roughness and journal misalignment affect bearing performance. Mallya et al. ⁴ analyze the static performance of misaligned multiple axial groove water-lubricated bearings under turbulent flow conditions. Their study shows that misalignment significantly influences pressure distribution and load-carrying capacity. The numerical model accounts for turbulence effects and geometric asymmetry, offering insights into bearing behavior under off-design operating conditions. Chen et al. ⁵ investigate the transient tribo-dynamic behavior of water-lubricated bearing–rotor systems under varying design parameters. Their study considers dynamic loading and fluid–structure interaction, showing how different bearing designs impact system stability, vibration characteristics, and lubrication performance during transient operations. Xiang et al. ⁶ explore the nonlinear dynamic behavior of water-lubricated bearings under mixed thermoelastohydrodynamic (TEHD) conditions. Their study integrates thermal effects, elastic deformation, and mixed lubrication characteristics to evaluate system stability. Chen et al. ⁷ perform a tribo-dynamic-wear coupling analysis on water-lubricated bearings with journal surface imperfections subjected to repeated start-stop cycles. Their study reveals how surface defects and transient operating conditions jointly affect wear evolution, lubrication film stability, and dynamic responses, providing a comprehensive understanding of bearing durability in real-world applications. Wang et al. ⁸ investigate the mixed thermal-visco-hyperelastic hydrodynamic lubrication behavior of water-lubricated rubber bearings operating in deep-sea conditions. Their work considers complex interactions of thermal effects, viscoelastic material properties, and hydrodynamic lubrication, providing insights into performance and reliability under extreme environmental pressures. Cheng et al. ⁹ examine the viscosity changes of water glycol under deep-sea conditions characterized by high pressure and low temperature. Their experimental study reveals how these extreme conditions influence the fluid’s rheological behavior, which is critical for lubrication performance in subsea applications. Xie et al. ¹⁰ perform a fluid-structure interaction study on a novel water-lubricated bearing considering particle-dynamic effects. Their combined theoretical and experimental approach reveals how suspended particles influence lubrication performance and structural responses, improving understanding of bearing behavior in contaminated water environments. Litwin et al. ¹¹ investigate the effects of polymer bearing materials and lubricating groove configurations on the wear behavior of journal bearings operating with contaminated water lubrication. Their experimental results highlight how material choice and groove layout significantly influence wear resistance and lubrication effectiveness under harsh conditions. Litwin and Kropp ¹² conduct theoretical and experimental investigations on sliding bearings with sintered bronze bushes lubricated by contaminated water containing solid particles. Their study examines how contamination affects lubrication film stability and bearing wear, providing insights for bearing design in polluted environments. Wodtke and Litwin ¹³ investigate the thermal effects on water-lubricated stern tube bearings through combined experimental and theoretical approaches. Their results highlight the impact of temperature variations on lubrication performance and bearing stability, providing guidance for design improvements under thermal loading. Song et al. ¹⁴ propose an entropy-based model to predict the fluid–solid–thermal coupled wear behavior of journal bearings subjected to repeated start-stop cycles. Their approach integrates thermodynamic principles with tribological and fluid dynamics analyses to enhance the accuracy of wear predictions under transient operating conditions.

In summary, previous studies have significantly advanced the understanding of water-lubricated bearing performance by incorporating turbulence,^1,3,4 mixed lubrication and asperity contact,^2,6 thermo-elasto-hydrodynamic (TEHD) and thermal effects,^6,8,13 fluid–structure interactions,^5,10 and wear mechanisms under transient or contaminated conditions.^7,11,12,14 These contributions have provided valuable insights into lubrication regimes, stability boundaries, and material effects. However, most of the existing models primarily focus on static or quasi-static performance and local effects, such as pressure distribution, load capacity, or wear evolution, while relatively few works have systematically examined the global nonlinear dynamics of coupled bearing–gear–rotor systems under turbulent mixed lubrication. In particular, comprehensive bifurcation analysis, characterization of chaotic responses, and parametric studies on damping and speed ratios remain insufficiently addressed in the literature. To bridge these gaps, the present study establishes a dimensionless nonlinear dynamic model of a water-lubricated bearing–gear–rotor system that incorporates turbulence and damping effects. By employing bifurcation diagrams, orbit trajectories, Poincaré sections, power spectra, fractal dimensions, and Lyapunov exponents, the study provides a systematic framework to analyze stability and complex dynamic behaviors. This approach not only complements but also extends previous research by offering a holistic understanding of the interplay between turbulence, damping, and system dynamics. In this study, the rotor is flexible, and it is supported by two nonlinear journal bearings (water-lubricated), with the motion of each bearing and the rotor influenced by forces from the bearings and gears. The rotor’s deflection, due to its flexibility, is modeled within the equations, and the nonlinear nature of the bearing suspension adds complexity to the dynamics of the system, reflecting real-world bearing behaviors under lubrication conditions. Analytical tools, including bifurcation diagrams, dynamic trajectories, power spectra, Poincaré map, fractal dimension, and Lyapunov exponent, are used to verify the dynamic characteristics of the rotor-bearing-gear system.

Mathematical modeling

In the present study, the model is developed under several assumptions. The lubricant is treated as an isothermal, incompressible, and inertia-less fluid, while mixed lubrication is modeled using a modified average Reynolds equation that includes turbulence and roughness corrections. The bearing is considered linearly elastic, with deformation calculated via the influence coefficient method. High-frequency structural vibrations and gear surface wear are neglected, and the bearing ends are assumed fixed, with rotor-bearing-gear coupling restricted to the horizontal plane. The proposed model also has certain limitations. High-frequency structural vibrations and nonlinear effects of the bearing material are not considered. Turbulence corrections approximate the average flow behavior, so instantaneous high Reynolds number effects are not captured. Effects of fluid contamination or particle dynamics are ignored, and asynchronous behavior between the rotor, bearings, and gears is not addressed, which could be explored in future studies.

The modified average Reynolds equation is used to model the lubrication in the bearing when it operates under mixed lubrication conditions, where the load is primarily supported by contact between rough surfaces. This occurs at high loads and low speeds. In this state, the lubrication is influenced significantly by the surface roughness of the bearing, which alters the behavior of the water film. For the analysis, an isothermal, incompressible, and inertialess flow model is employed through the use of the average Reynolds equation ¹⁵

$\frac{1}{R^2}\frac{\partial }{\partial \theta }\left({\varphi }_{\theta }{h}^3\frac{\partial p}{\partial \theta }\right)+\frac{\partial }{\partial z}\left({\varphi }_z{h}^3\frac{\partial p}{\partial z}\right)=6\mathrm{\Omega }{\mu \varphi }_c\frac{\partial h}{\partial \theta }+6\Omega \mu \sigma \frac{\partial {\varphi }_s}{\partial \theta }+12\mu {\varphi }_c\frac{\partial h}{\partial t}$

(1)

To incorporate turbulence into the average Reynolds equation, turbulence correction factors were introduced based on the bulk-flow theory of Hirs ¹⁶ and the Ng–Pan model. ¹⁷ In these models, the effect of turbulence is included by modifying the pressure and shear flow factors ( ${\varphi }_{\theta }$ , ${\varphi }_z$ , ${\varphi }_s$ ) as functions of the local Reynolds number. Accordingly, in equation (1), ${\varphi }_{\theta }$ and ${\varphi }_z$ represent not only the roughness-induced flow factors (Patir & Cheng model) but also include the turbulent correction to account for the inertia effect of high Reynolds number flows. When the Reynolds number is small, the flow remains laminar and the factors tend to unity, that is, ${\varphi }_{\theta }$ , ${\varphi }_z$ → 1. However, as the Reynolds number increases, the turbulence correction reduces the effective pressure and shear flow contributions, thereby capturing the transition to turbulent lubrication. Unlike the conventional Reynolds equation, which assumes laminar flow between smooth surfaces, the present formulation extends the average Reynolds equation in three important aspects. First, turbulence corrections are introduced through the flow factors ( ${\varphi }_{\theta }$ , ${\varphi }_z$ , ${\varphi }_s$ ) based on the bulk-flow theory of Hirs ¹⁶ and the Ng–Pan model. ¹⁷ Second, the model incorporates mixed lubrication effects by including a contact factor ( ${\varphi }_c$ ) that accounts for asperity interactions between rough surfaces. Third, elastic deformation of the bearing bush is considered using the influence coefficient method, allowing the film thickness to couple fluid pressure with structural compliance. This integration transforms the average Reynolds equation into a mixed lubrication framework with turbulence and deformation effects, which is further recast into a dimensionless nonlinear form for system-level dynamic analysis. The novelty of the proposed flow model therefore lies in combining turbulence, surface roughness, and structural deformation in a unified formulation suitable for analyzing the nonlinear dynamics of the bearing-gear-rotor system. Where θ and z are the bearing circumferential and axial coordinates, respectively, R is the bearing radius, ${\varphi }_{\theta }$ and ${\varphi }_z$ are the pressure flow factors, ${\varphi }_s=\frac{2E}{R\mathrm{\Omega }}\left(-\frac{h^3}{12\mu }\frac{\partial p}{\partial x}\right)$ is the shear flow factor ( ${\varphi }_{\theta }$ , ${\varphi }_z\to 1\hspace{0.17em}\&\hspace{0.5em}{\varphi }_s\to 0\hspace{0.17em}as\hspace{0.17em}\frac{h}{\sigma }\to \infty$ ), ${\varphi }_c$ is the contact factor, p is the water film pressure, ℎ is the water film thickness, μ is the viscosity, ω is the journal rotational angular speed, and $\sigma =\sqrt{{{\sigma }_1}^2+{{\sigma }_2}^2}$ is the surface composite roughness, ${\sigma }_1$ and ${\sigma }_2$ are standard deviations of ${\delta }_1$ and ${\delta }_2$ (random roughness amplitude of surfaces). The water film thickness equation is divided into two components: the nominal water film thickness generated by the rotor’s motion is $h_r=c\left(1+\epsilon \hspace{0.17em}\mathit{cos}\hspace{0.17em}\theta \right)$ , and the elastic deformation of the bearing bush is defined as $h_e$ . This relationship can be expressed as

$h\left(\theta ,z,t\right)=h_r\left(\theta ,z\right)+h_e\left(\theta ,z,t\right)$

(2)

Because the bearing bush typically has a small elastic modulus, deformation during operation is significant. The deformation of the bearing bush is calculated by the influence coefficient method ¹⁸

$h_e\left(\theta ,z,t\right)={\Sigma }_{\varsigma }{\Sigma }_{\kappa }G\left({\theta }_I,z_j,{\theta }_{\zeta },z_{\kappa }\right)(p\left({\theta }_{\zeta },z_{\kappa }\right)+p_a\left({\theta }_{\zeta },z_{\kappa }\right)$

(3)

where $p_a\left({\theta }_{\zeta },z_{\kappa }\right)=k_c{E}_c{F}_{25}\left(H_d\right)$ is the asperity contact pressure, ¹⁹ $k_c=\frac{16\sqrt{2}}{15\pi {\left(\sigma {\beta }_0\mu \right)}^2\sqrt{\frac{\sigma }{{\beta }_0}}}$ is the coefficient of the contact stiffness, $E_c=\frac{1}{\left(\frac{\left(1-{{\upsilon }_j}^2\right)}{E_j}+\frac{\left(1-{{\upsilon }_b}^2\right)}{E_b}\right)}$ is the composite elastic modulus between the journal and the bearing, $F_{25}\left(H_d\right)=\{\begin{array}{cc}4.4086\times {10}^{-5}\times {\left(4-H_d\right)}^{6.804}& for\hspace{0.5em}{H}_d\le 4\\ 0& for\hspace{0.5em}{H}_d>4\end{array}$ is a function that depends on the dimensionless film thickness ratio $H_d$ and which is also defined by. ¹⁹ In previous studies, bearing deformation was often neglected or simplified by assuming a rigid bearing bush, which may be reasonable for metallic bearings with high stiffness. However, for water-lubricated polymer bearings, the relatively low elastic modulus makes the elastic deformation comparable to the nominal film thickness, thus significantly influencing pressure distribution and load-carrying capacity. In the present study, the bearing deformation is calculated by the influence coefficient method, which allows the local film thickness to be coupled directly with both hydrodynamic pressure and asperity contact pressure. This approach ensures that the interaction between fluid film pressure and structural compliance is realistically represented. The originality of this formulation lies in extending the mixed lubrication framework to include the structural compliance of water-lubricated bearings, which is essential for accurately predicting system stability and dynamic responses under high load conditions.

From the pressure distribution results and asperity contact pressure distribution with different L/D ratios (i.e., L/D = 0.4, 0.8, 1.2, 1.6, and 2.0) shown in Figure 1(a) and (b), it can be observed that short bearings (with a smaller L/D ratio) may exhibit higher load capacity. This is because short bearings typically have a smaller axial length, which, under the same size constraints, results in a relatively larger contact area. This larger contact area allows the bearing to better withstand the applied load. This effect influences the pressure distribution in the bearing, particularly under operational conditions. In short bearings, the film thickness and dynamic behavior of the fluid may lead to higher pressures, which could enhance load capacity. Moreover, the asperity contact pressure, which arises from the microscopic roughness of the bearing surfaces, also plays a crucial role in the overall pressure distribution. The presence of asperities can increase localized contact pressure in areas where the lubricant film is thinner or where surface roughness causes direct contact between the bearing surfaces. This asperity contact pressure is influenced by the roughness height and distribution, and its interaction with the fluid film pressure can further elevate the total pressure in certain regions. The combination of fluid film pressure and asperity contact pressure determines the overall performance and load capacity of the bearing. Shorter bearings with higher surface roughness might experience more significant asperity contact pressures, leading to variations in the pressure distribution compared to longer bearings.

Figure 1.

(a) The pressure distribution and (b) The asperity contact pressure distribution.

The total load capacity of a water-lubricated bearing can be considered as the sum of the contributions from the hydrodynamic pressure in the lubricating film and the asperity contact between the bearing and the journal. The modified Reynolds equation will divided into two parts with the short bearing approximation (i.e., L/D < 1.0、∂p/∂θ << ∂p/∂z, and then we can set ∂p/∂θ = 0.) and the long bearing approximation (i.e., L/D > 1.0、∂p/∂z << ∂p/∂θ, and then we can set ∂p/∂z = 0.) could be performed as

$\frac{{\partial }^{2}p}{\partial z^2}=-\frac{6\mu c}{h^3}\left[\left(\mathrm{\Omega }-2\dot{\phi }\right)\epsilon \frac{{\varphi }_c}{{\varphi }_z}\mathit{sin}\hspace{0.17em}\theta -2\dot{\epsilon }\frac{{\varphi }_c}{{\varphi }_z}\mathit{cos}\hspace{0.17em}\theta \right]\text{for}\hspace{0.17em}\text{short}\hspace{0.17em}\text{bearing}$

(4)

$\frac{{\partial }^{2}p}{\partial {\theta }^2}=-\frac{6\mu c}{\left(\frac{{\varphi }_{\theta }}{R^2}+\frac{E}{6\mu R^2\mathrm{\Omega }\sigma }\right)h^3}\left[\left(\mathrm{\Omega }-2\dot{\phi }\right)\epsilon \frac{{\varphi }_c}{{\varphi }_z}\mathit{sin}\hspace{0.17em}\theta -2\dot{\epsilon }\frac{{\varphi }_c}{{\varphi }_z}\mathit{cos}\hspace{0.17em}\theta \right]\text{for}\hspace{0.17em}\text{long}\hspace{0.17em}\text{bearing}$

(5)

The pressure distribution for short bearing with the boundary conditions $\{\begin{array}{l}\frac{\partial p}{\partial z}=0,z=0\\ p=0,z=\pm \frac{L}{2}\end{array}$ will be

$p_{short}=\frac{6\mu c}{h^3}\left[\left(\mathrm{\Omega }-2\dot{\phi }\right)\epsilon \frac{{\varphi }_c}{{\varphi }_z}\mathit{sin}\hspace{0.17em}\theta -2\dot{\epsilon }\frac{{\varphi }_c}{{\varphi }_z}\mathit{cos}\hspace{0.17em}\theta \right]\left(z^2-\frac{L^2}{4}\right)$

(6)

The pressure distribution for long bearing

$p_{long}={\int }_0^{\pi }{\int }_0^{\pi }\frac{6\mu c}{\left(\frac{{\varphi }_{\theta }}{R^2}+\frac{E}{6\mu R^2\mathrm{\Omega }\sigma }\right){\left(c\left(1+\epsilon \hspace{0.17em}\mathit{cos}\hspace{0.17em}\theta \right)+h_e\left(\theta ,z,t\right)\right)}^3}\left[\left(\mathrm{\Omega }-2\dot{\phi }\right)\epsilon \frac{{\varphi }_c}{{\varphi }_z}\mathit{sin}\hspace{0.17em}\theta -2\dot{\epsilon }\frac{{\varphi }_c}{{\varphi }_z}\mathit{cos}\hspace{0.17em}\theta \right]d\theta d\theta$

(7)

The water film force in the radial and tangential direction for short and long bearing approximation could be presented by integrating the function $p_{short}$ and $p_{long}$ both with $\theta$ and z. That is

$\{\begin{array}{l}{f}_{e,short}=\frac{L^{3}R\mu }{{2c}^2}{\int }_0^{\pi }\frac{\left[\left(\mathrm{\Omega }-2\dot{\phi }\right)\epsilon \frac{{\varphi }_c}{{\varphi }_z}\mathit{sin}\hspace{0.17em}\theta -2\dot{\epsilon }\frac{{\varphi }_c}{{\varphi }_z}\mathit{cos}\hspace{0.17em}\theta \right]}{{\left(1+\epsilon \hspace{0.17em}\mathit{cos}\hspace{0.17em}\theta \right)}^3}\mathit{cos}\hspace{0.17em}\theta d\theta \\ \hspace{0.17em}{f}_{\phi ,short}=-\frac{L^{3}R\mu }{{2c}^2}{\int }_0^{\pi }\frac{\left[\left(\mathrm{\Omega }-2\dot{\phi }\right)\epsilon \frac{{\varphi }_c}{{\varphi }_z}\mathit{sin}\hspace{0.17em}\theta -2\dot{\epsilon }\frac{{\varphi }_c}{{\varphi }_z}\mathit{cos}\hspace{0.17em}\theta \right]}{{\left(1+\epsilon \hspace{0.17em}\mathit{cos}\hspace{0.17em}\theta \right)}^3}\mathit{sin}\hspace{0.17em}\theta d\theta \end{array}$

(8)

$\{\begin{array}{l}{f}_{e,long}=6\mu cL{\int }_0^{\pi }{\int }_0^{\pi }{\int }_0^{\pi }\begin{array}{l}\frac{6\mu c}{\left(\frac{{\varphi }_{\theta }}{R^2}+\frac{E}{6\mu R^2\mathrm{\Omega }\sigma }\right){\left(c\left(1+\epsilon \hspace{0.17em}\mathit{cos}\hspace{0.17em}\theta \right)+h_e\left(\theta ,z,t\right)\right)}^3}[\left(\mathrm{\Omega }-2\dot{\phi }\right)\epsilon \frac{{\varphi }_c}{{\varphi }_z}\mathit{sin}\hspace{0.17em}\theta -\\ 2\dot{\epsilon }\frac{{\varphi }_c}{{\varphi }_z}\mathit{cos}\hspace{0.17em}\theta ]d\theta d\theta \hspace{0.17em}\mathit{cos}\hspace{0.17em}\theta d\theta \end{array}\\ \hspace{0.17em}{f}_{\phi ,long}=6\mu cL{\int }_0^{\pi }{\int }_0^{\pi }{\int }_0^{\pi }\begin{array}{l}\frac{6\mu c}{\left(\frac{{\varphi }_{\theta }}{R^2}+\frac{E}{6\mu R^2\mathrm{\Omega }\sigma }\right){\left(c\left(1+\epsilon \hspace{0.17em}\mathit{cos}\hspace{0.17em}\theta \right)+h_e\left(\theta ,z,t\right)\right)}^3}[\left(\mathrm{\Omega }-2\dot{\phi }\right)\epsilon \frac{{\varphi }_c}{{\varphi }_z}\mathit{sin}\hspace{0.17em}\theta -\\ 2\dot{\epsilon }\frac{{\varphi }_c}{{\varphi }_z}\mathit{cos}\hspace{0.17em}\theta ]d\theta d\theta \hspace{0.17em}\mathit{sin}\hspace{0.17em}\theta d\theta \end{array}\end{array}$

(9)

Then the water film force in the horizontal and vertical direction for short and long bearing approximation could be

$\{\begin{array}{l}{F}_x=f_e\hspace{0.17em}\mathit{cos}\hspace{0.17em}\phi -f_{\phi }\hspace{0.17em}\mathit{sin}\hspace{0.17em}\phi \\ F_y=f_e\hspace{0.17em}\mathit{sin}\hspace{0.17em}\phi +f_{\phi }\hspace{0.17em}\mathit{cos}\hspace{0.17em}\phi \end{array}$

(10)

Figure 2 shows a flexible rotor supported by two water-lubricated journal bearings under nonlinear suspension. The rotor’s geometric center is denoted by $\left(x_r,y_r\right)$ , and the bearing geometric centers are at $\left(x_b,y_b\right)$ . The nonlinear bearing forces depend on the relative displacements between the rotor and the bearings, which are modeled as nonlinear functions of both displacement and velocity. The water-lubricated journal bearings support the rotor and exert forces that are dependent on the rotor’s motion. These forces are nonlinear due to the nature of the bearing suspension, which typically involves complex fluid dynamics and elastic properties of the bearing material. The motion equations for the bearings are governed by the following nonlinear equations

$\{\begin{array}{l}{m}_b\ddot{x_b}+c_b\dot{x_b}+k_{b1}\left(x_b-x_r\right)+k_{b2}\left({x_b}^3-{x_r}^3\right)-F_x=0\hspace{0.17em}\\ m_b\ddot{y_b}+c_b\dot{y_b}+k_{b1}\left(y_b-y_r\right)+k_{b2}\left({y_b}^3-{y_r}^3\right)-F_y=-m_{b}g\end{array}$

(11)

where $m_b$ is the mass of the bearing and $c_b$ is the damping coefficient of the bearing. $k_{b1}$ and $k_{b2}$ are the first and third order stiffness, respectively, to present nonlinear suspension effect. $F_x$ and $F_y$ are the nonlinear bearing support forces in the x-direction and y-direction, respectively, which depend on bearing displacement and velocity. These forces are defined in equation (10) and their dynamic responses will be compared with those predicted by short and long bearing approximations, respectively. Equation (11) represents a two-degree-of-freedom Duffing-type nonlinear oscillator for the bearing-rotor system, commonly used in rotordynamics studies.^20,21 In the present study, the equation is extended by including the fluid-film forces $F_x$ and $F_y$ derived from the turbulence- and roughness-modified mixed lubrication model, thereby coupling the nonlinear structural dynamics with hydrodynamic and asperity contact effects. This enables the investigation of complex dynamic behaviors such as periodic, quasi-periodic, and chaotic motions under realistic operating conditions.

Figure 2.

Flexible rotor-gear system supported by water-lubricated bearings with nonlinear suspension.

The flexible rotor is influenced by the forces from the bearings and the gears, as well as any external forces applied to it. The rotor’s dynamic behavior is significantly influenced by the flexibility of the rotor itself, which leads to deflections and vibrations

$\{\begin{array}{l}{m}_r\ddot{x_r}+c_r\dot{x_r}+k_r\left(x_r-x_b\right)-F_{g,x}=m_r\rho {\mathrm{\Omega }}^2\hspace{0.17em}\mathit{cos}\hspace{0.17em}\varphi \hspace{0.17em}\\ m_r\ddot{y_r}+c_r\dot{y_r}+k_r\left(y_r-y_b\right)-F_{g,y}=m_r\rho {\mathrm{\Omega }}^2\hspace{0.17em}\mathit{sin}\hspace{0.17em}\varphi -m_{r}g\end{array}$

(12)

where $m_r$ is the mass of the rotor, $c_r$ is the damping coefficient of the rotor, $k_r$ is the spring constant for the rotor. $\rho$ is the mass eccentricity of the rotor, $\omega$ is the rotational speed of the rotor and $\varphi$ is rotational angle (ϕ = ωt). $F_{g,x}=-C_m\left[{\dot{X}}_p-{\dot{X}}_g-E_p\mathrm{\Omega }\hspace{0.17em}\sin \left(\mathrm{\Omega }t\right)\right]-K_m\left[X_p-X_g-E_p\hspace{0.17em}\cos \left(\mathrm{\Omega }t\right)\right]$ , and $F_{g,y}=M_g{E}_g{{\mathrm{\Omega }}_g}^2\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_2-M_g{E}_g{\ddot{\theta }}_2\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_2+C_m\left[{\dot{X}}_p-{\dot{X}}_g-E_p\mathrm{\Omega }\hspace{0.17em}\sin \left(\mathrm{\Omega }t\right)\right]+K_m\left[X_p-X_g-E_p\hspace{0.17em}\cos \left(\mathrm{\Omega }t\right)\right]$ .

The gears, including the pinion and the large gear, interact with the rotor and each other (see Figure 3). These interactions involve forces transmitted through the contact points of the gears. The pinion and the gear have their own moments of inertia, damping, and spring constants, influencing the overall dynamics of the system. For the pinion gear, its motion in the x-direction and y-direction is governed by

$\{\begin{array}{l}{J}_p\ddot{x_p}+c_p\dot{x_p}+k_p\left(x_p-x_r\right)+F_{p,x}=0\hspace{0.17em}\\ J_p\ddot{y_p}+c_p\dot{y_p}+k_p\left(y_p-y_r\right)+F_{p,y}=0\end{array}$

(13)

Figure 3.

Interaction between the gear and the pinion.

For the large gear, its motion in the x-direction and y-direction is given by

$\{\begin{array}{l}{J}_g\ddot{x_g}+c_g\dot{x_g}+k_g\left(x_g-x_r\right)+F_{g,x}=0\hspace{0.17em}\\ J_g\ddot{y_g}+c_g\dot{y_g}+k_g\left(y_g-y_r\right)+F_{g,x}=0\end{array}$

(14)

where $F_{p,x}=C_m\left[{\dot{X}}_p-{\dot{X}}_g-E_p\mathrm{\Omega }\hspace{0.17em}\sin \left(\mathrm{\Omega }t\right)\right]+K_m\left[X_p-X_g-E_p\hspace{0.17em}\cos \left(\mathrm{\Omega }t\right)\right]$ , and $F_{p,y}=M_p{E}_p{{\mathrm{\Omega }}_p}^2\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_1-M_p{E}_p{\ddot{\theta }}_1\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_1-C_m\left[{\dot{Y}}_p-{\dot{Y}}_g-E_p\mathrm{\Omega }\hspace{0.17em}\cos \left(\mathrm{\Omega }t\right)\right]+K_m\left[Y_p-Y_g-E_p\hspace{0.17em}\sin \left(\mathrm{\Omega }t\right)\right]$ .

Equations (1), (2), and (12) are the nonlinear dynamic equations used to simulate the dynamic responses of the rotor-bearing system, considering temperature-dependent viscosity in this study, and then integrating them to be the following dimensionless equations

$\{\begin{array}{l}{B}_x^{,,}+\frac{2\xi }{s_1}{B}_x^{,}+\frac{1}{{s_1}^2}\left(B_x-R_x\right)+\frac{\alpha }{s^2}\left({B_x}^3-{R_x}^3\right)-\frac{F_x}{m_{r}c{\mathrm{\Omega }}^2}=0\hspace{0.17em}\\ B_y^{,,}+\frac{2\xi }{s_1}{B}_y^{,}+\frac{1}{{s_1}^2}\left(B_y-R_y\right)+\frac{\alpha }{s^2}\left({B_y}^3-{R_y}^3\right)-\frac{F_y}{m_{r}c{\mathrm{\Omega }}^2}+\frac{f_b}{{s_1}^2}=0\end{array}$

(15)

$\{\begin{array}{l}{R}_x^{,,}+\frac{2{\xi }_2}{s}{R}_x^{,}+\frac{1}{s^2}\left(R_x-B_x\right)-\frac{F_{g,x}}{m_{r}c{\mathrm{\Omega }}^2}=\beta \hspace{0.17em}\mathit{cos}\hspace{0.17em}\varphi \hspace{0.17em}\\ R_y^{,,}+\frac{2{\xi }_2}{s}{R}_y^{,}+\frac{1}{s^2}\left(R_y-B_y\right)-\frac{F_{g,y}}{m_{r}c{\mathrm{\Omega }}^2}=\beta \hspace{0.17em}\mathit{sin}\hspace{0.17em}\varphi -\frac{f_r}{s^2}\end{array}$

(16)

$\{\begin{array}{l}{P}_x^{,,}+\frac{2{\xi }_3}{s_2}{P}_x^{,}+\frac{1}{{s_2}^2}\left(P_x-R_x\right)+\frac{F_{p,x}}{m_{r}c{\mathrm{\Omega }}^2}=0\hspace{0.17em}\\ P_y^{,,}+\frac{2{\xi }_3}{s_2}{P}_y^{,}+\frac{1}{{s_2}^2}\left(P_y-R_y\right)+\frac{F_{p,y}}{m_{r}c{\mathrm{\Omega }}^2}=0\end{array}$

(17)

$\{\begin{array}{l}{G}_x^{,,}+\frac{2{\xi }_4}{s_3}{G}_x^{,}+\frac{1}{{s_3}^2}\left(G_x-R_x\right)+\frac{F_{g,x}}{m_{r}c{\mathrm{\Omega }}^2}=0\hspace{0.17em}\\ G_y^{,,}+\frac{2{\xi }_4}{s_3}{G}_y^{,}+\frac{1}{{s_3}^2}\left(G_y-R_y\right)+\frac{F_{g,y}}{m_{r}c{\mathrm{\Omega }}^2}=0\end{array}$

(18)

where c is the clearance of water film, $\varphi =\omega t$ , $\frac{d}{dt}=\omega \frac{d}{d\varphi }$ , ${\mathrm{\Omega }}^2=\frac{k_r}{m_r},s^2=\frac{{\omega }^2}{{\mathrm{\Omega }}^2}$ , $\beta =\frac{m_r\rho {\omega }^2}{m_{r}c{\omega }^2}$ , $f_b=\frac{m_{r}g}{k_{r}c}$ , $\xi =\frac{c_b}{2\sqrt{k_b{m}_b}}$ , ${\xi }_2=\frac{c_r}{2\sqrt{k_r{m}_r}}$ , ${\xi }_3=\frac{c_p}{2\sqrt{k_p{m}_p}}$ , ${\xi }_4=\frac{c_g}{2\sqrt{k_g{m}_g}}$ , $R_{mb}=\frac{m_b}{m_r}$ , $R_{kb}=\frac{k_b}{k_r}$ , ${s_1}^2=R_{mb}{R}_{kb}{s}^2$ , $R_{mp}=\frac{m_p}{m_r}$ , $R_{kp}=\frac{k_p}{k_r}$ , ${s_2}^2=R_{mp}{R}_{kp}{s}^2$ , $R_{mg}=\frac{m_g}{m_r}$ , $R_{kp}=\frac{k_g}{k_r}$ , ${s_3}^2=R_{mp}{R}_{kp}{s}^2$ .

Compared with previous models, the present formulation emphasizes system-level nonlinear dynamics rather than localized lubrication effects. For example, Qiao et al. ¹ extended the average Reynolds equation with turbulence correction factors to predict static pressure distribution in WLRBs, while Zhu et al. ³ integrated the Hashimoto turbulence model with the Patir & Cheng roughness model and applied the finite difference method to analyze thermal and roughness influences on journal bearings. These approaches are valuable but remain limited to static or quasi-static performance and local lubrication phenomena. In contrast, the present study develops a comprehensive dimensionless nonlinear model for the coupled bearing-gear-rotor system. The proposed formulation incorporates both turbulence and mixed lubrication under dynamic operating conditions, thereby enabling the investigation of global nonlinear responses and stability boundaries. Furthermore, while earlier studies often employed finite difference or finite element techniques for localized lubrication analysis, the present work combines the Runge–Kutta method with nonlinear dynamical system tools including bifurcation diagrams, Poincaré maps, power spectra, fractal dimensions, and Lyapunov exponents to provide a systematic characterization of periodic, quasi-periodic, and chaotic motions. This integration highlights the novelty of our approach and its contribution to advancing the understanding of water-lubricated bearing–gear–rotor dynamics.

Results and discussions

In this study, the fourth-order Runge–Kutta (RK4) method is employed to numerically integrate the nonlinear governing equations of the coupled bearing–gear–rotor system. Although more advanced techniques such as adaptive time-stepping methods or implicit solvers are available, RK4 is chosen due to its simplicity, robustness, and high accuracy for problems with moderate stiffness. The step size is carefully selected to ensure numerical stability and convergence of the solution. Moreover, RK4 allows straightforward implementation of the Poincaré map, bifurcation diagram, and Lyapunov exponent calculations, which are essential for the nonlinear dynamic analysis conducted in this work. Specifically, it is used to solve the dimensionless nonlinear dynamic equations governing the rotor–bearing system. A fixed time step of π/300 is selected to ensure sufficient temporal resolution, with the local truncation error controlled to be less than 1×10 ^-4 . To ensure that only steady-state responses are analyzed, the initial 800 rotor revolutions are discarded from the time series data. Subsequently, 60,000 data points representing the relative displacement between the bearing and rotor geometric centers are collected and utilized to reconstruct the system’s attractor in phase space using delay-coordinate embedding. From these embedded dynamics, several nonlinear diagnostic tools are constructed, including the orbital trajectories, power spectral densities, Poincaré sections, and bifurcation diagrams. In addition, two key nonlinear invariants, that is, the fractal (correlation) dimension and the maximum Lyapunov exponent (MLE), are computed to further characterize the system’s dynamic regime (e.g., periodic, quasi-periodic, or chaotic).

The optimal time delay for phase space reconstruction is determined via the autocorrelation function, as suggested by Nayfeh, ²² yielding a delay time approximately equal to one-third of a rotor revolution. The correlation dimension is then calculated using the Grassberger–Procaccia algorithm, which involves plotting the logarithm of the correlation integral log[C(r)] versus the logarithm of the radius log(r) of an N-dimensional hypersphere. According to Smith, ²³ the embedding dimension M should be the greatest integer less than the actual fractal dimension of the attractor, and the number of data points used for estimation must satisfy the condition N<42 ^M to ensure statistical reliability. In the bifurcation analysis, two dimensionless control parameters are varied: the rotational speed ratio s and the dimensionless damping coefficient $\xi$ . These parameters govern the system’s transition across different dynamic regimes and are critical in identifying stability boundaries and nonlinear phenomena such as bifurcations and chaos.

Since the supporting forces generated by the water film are inherently nonlinear and result from the integration of pressure distributions, it is essential to first examine the numerical results of both the fluid film pressure distribution and the asperity contact pressure distribution under various aspect ratios (L/D), as illustrated in Figure 1. These results serve as a foundation for the subsequent analysis of the system’s nonlinear dynamic characteristics. According to the numerical analysis, short bearings (i.e., those with smaller L/D ratios) may exhibit enhanced load-carrying capacity due to several important factors. First, for a fixed outer diameter and comparable overall size, short bearings provide a relatively larger effective contact area. This allows for a higher load per unit area, thereby improving the bearing’s load-supporting capability. Second, short bearings typically operate with thinner lubricant films under dynamic conditions. This leads to an increase in hydrodynamic pressure within the film, which further contributes to improved load capacity. Third, the fluid dynamic behavior within short bearings differs from that in longer bearings. These differences can result in localized pressure intensification due to more abrupt changes in film geometry, thereby enhancing the bearing’s ability to support dynamic loads. However, it is also important to consider system stability. Longer bearings tend to offer better dynamic stability because of their extended axial length, which allows for smoother pressure distribution and improved damping characteristics. Therefore, selecting an appropriate bearing design requires balancing load capacity and dynamic stability based on the specific requirements of the application. If load-carrying capacity is the primary objective and other system parameters can be adjusted to maintain sufficient stability, then a short bearing may be the preferred choice. On the other hand, if operational stability and durability are of greater concern, a longer bearing may be more suitable. In this study, a short bearing configuration is chosen in order to investigate its nonlinear dynamic behavior. This choice is motivated by the fact that short bearings are more prone to complex dynamic responses, making them ideal for studying phenomena such as bifurcation, quasi-periodicity, and chaos within water-lubricated bearing systems.

Figure 4 presents the bifurcation diagrams of the geometric centers of the bearing, gear, and rotor in the horizontal direction, plotted with respect to the dimensionless rotational speed ratio s ranging from 0.10 to 6.10. These diagrams reveal the global dynamic characteristics of the system under varying speed ratios. Specifically, the coordinates of the return points in the Poincaré map, namely the bearing geometric center (B _x , B _y ), the gear geometric center (G _x , G _y ), and the rotor geometric center (R _x , R _y ) are extracted and plotted to construct the bifurcation diagrams. This analysis provides insight into the transitions between different dynamic states, such as periodic, quasi-periodic, and chaotic motions across a wide speed range. From the bifurcation diagram results, it can be observed that the dynamic responses of the bearing, gear, and rotor exhibit non-periodic vibrations at lower rotational speeds (i.e., when s <1.2). As the system transitions into higher speed regimes, the responses evolve into periodic motions. Furthermore, it is evident that at higher speeds, the vibrational response of the gear becomes desynchronized from that of the bearing and rotor. This behavior deviates significantly from what is typically observed in traditional oil-lubricated bearing systems, where components often exhibit more synchronized and smoother transitions in dynamic behavior. This discrepancy may be attributed to the distinct fluid-structure interaction characteristics and nonlinear properties of water-lubricated bearings, including lower viscosity, higher compressibility effects, and the presence of mixed lubrication or asperity contact. These factors can introduce additional complexity into the system dynamics, making the response more sensitive to speed variations and structural coupling effects, particularly under high-speed conditions. To gain a clearer understanding of the system’s dynamic characteristics, we refer to Figures 5 and 6, which illustrate the dynamic orbits and Poincaré maps of the bearing geometric center (B _x ) and the gear geometric center (G _x ) for various dimensionless speed ratios s = 0.1, 0.5, 1.0, 1.2, 2.0, and 6.0 (arranged from left to right and top to bottom). At lower speed ratios (s = 0.1, 0.5, 1.0, and 1.2), the dynamic orbits exhibit highly irregular, non-periodic oscillations, while the corresponding Poincaré maps display a dense set of scattered points. These characteristics are indicative of chaotic motion, suggesting that both the bearing and gear geometric centers are operating within a chaotic regime at these speeds. As the speed ratio increases to s = 2.0 and s = 6.0, the bearing’s dynamic response transitions into a periodic vibration state, reflecting greater stability. Interestingly, under these high-speed conditions, the gear displays a distinct 8T-periodic response, suggesting the emergence of a high-order subharmonic motion or frequency locking between meshing dynamics and system rotation. This behavior may lead to gear-bearing dynamic decoupling, potentially causing mismatched excitation frequencies between components. Such a mismatch can amplify vibrational energy transfer, increase wear, and compromise the long-term operational stability of the system, particularly under high-speed or misaligned conditions. Further investigation into coupling stiffness, damping nonlinearity, and the influence of tooth contact dynamics is necessary to fully understand the root cause of this asynchronous behavior.

Figure 4.

Bifurcation diagram of the bearing geometric center, gear geometric center, and rotor geometric center in the horizontal direction using dimensionless rotational speed ratio s (s = 0.10∼6.10).

Figure 5.

Dynamic orbits and Poincaré maps of the bearing geometric center (B _x ) for s = 0.1, 0.5, 1.0, 1.2, 2.0, and 6.0. (Arranged from left to right and top to bottom).

Figure 6.

Dynamic orbits and Poincaré maps of the gear geometric center (G _x ) for s = 0.1, 0.5, 1.0, 1.2, 2.0, and 6.0. (Arranged from left to right and top to bottom).

To provide a detailed explanation of the observed chaotic behavior, we refer to Figure 7, Figure 8, and Figure 9, which present the dynamic orbit, Poincaré map, power spectrum, fractal dimension, and maximum Lyapunov exponent of the bearing geometric center, gear center, and rotor center at a dimensionless rotational speed ratio of s = 0.33. The dynamic orbits exhibit highly irregular and unstable non-periodic oscillations, while the corresponding Poincaré maps show dense clusters of discrete scattered points, which are typical characteristics of a chaotic system. In addition, the power spectra reveal the presence of numerous excitation frequencies, further indicating that the system does not exhibit simple periodic behavior under this speed condition. To further confirm the presence of chaos, we computed the fractal dimension and maximum Lyapunov exponent for each component. The fractal dimensions of the bearing geometric center, gear center, and rotor center are 1.76, 1.38, and 1.68, respectively, all of which are non-integer values, suggesting the presence of a strange attractor in each case. Furthermore, the maximum Lyapunov exponents for the bearing, gear, and rotor are 0.17, 0.07, and 0.16, respectively, all of which are positive values. These results provide strong evidence that at s = 0.33, the dynamic responses of the bearing, gear, and rotor all fall into a chaotic vibration state.

Figure 7.

Dynamic orbit, Poincaré map, power spectrum, fractal dimension, and maximum Lyapunov exponent of the bearing geometric center at dimensionless rotational speed ratio s = 0.33.

Figure 8.

Dynamic orbit, Poincaré map, power spectrum, fractal dimension, and maximum Lyapunov exponent of the gear geometric center at dimensionless rotational speed ratio s = 0.33.

Figure 9.

Dynamic orbit, Poincaré map, power spectrum, fractal dimension, and maximum Lyapunov exponent of the rotor geometric center at dimensionless rotational speed ratio s = 0.33.

The dimensionless damping ratio, $\xi$ is also an important parameter for evaluating the nonlinear dynamic characteristics of the system. This is because damping significantly influences the system’s energy dissipation capacity, stability, and response amplitude under periodic and non-periodic excitations. Variations in the damping ratio can lead to bifurcations, transitions between periodic and chaotic motions, and changes in resonance behavior. To illustrate this, we refer to Figure 10, which presents the bifurcation diagrams of the bearing geometric center, gear geometric center, and rotor geometric center in the horizontal direction under a fixed dimensionless rotational speed ratio of s = 0.33, while the dimensionless damping ratio $\xi$ is varied from 0.001 to 0.397. From the bifurcation diagrams, it can be observed that at low damping levels ( $\xi <0.14$ ), the system exhibits simultaneously unstable and non-periodic oscillations across all components, with relatively large vibration amplitudes. This behavior is indicative of chaotic motion and results from insufficient damping to suppress the energy transferred through nonlinear interactions within the system. As the damping ratio increases beyond $\xi >0.14$ , the system responses gradually stabilize and transition into periodic vibrations, suggesting that the increased damping effectively suppresses chaotic dynamics and enhances overall system stability. However, within the intermediate range of 0.214< $\xi <0.27$ , the system once again enters a regime of non-periodic and unstable oscillations, although the vibration amplitudes become significantly smaller compared to the low-damping case. This may be attributed to the complex interplay between partial energy dissipation and residual nonlinear excitation, resulting in the re-emergence of chaotic motion at a lower energy level. Finally, when the damping ratio exceeds $\xi >0.27$ , the system responses across the bearing, gear, and rotor consistently return to periodic vibration states. In this high-damping regime, the system achieves better synchronization and dynamic suppression, indicating that sufficient damping is critical to maintaining stable operation and avoiding chaotic instabilities. To gain a clearer understanding of the system’s dynamic characteristics, we refer to Figures 11 and 12, which illustrate the dynamic orbits and Poincaré maps of the bearing geometric center (B _x ) and the gear geometric center (G _x ) at various dimensionless damping ratios $\xi =0.001,0.01,0.1,0.16,0.24,\text{and}\hspace{0.17em}0.3$ (arranged from left to right and top to bottom). The dynamic orbit diagrams reveal unstable oscillatory behavior, while the corresponding Poincaré maps show multiple clusters of return points, indicating that at $\xi =0.001,0.01,0.1,\text{and}\hspace{0.17em}0.24$ , the system enters a state of chaotic vibration. In contrast, at $\xi =0.16\hspace{0.17em}\text{and}\hspace{0.17em}0.3,$ both the bearing and gear complete only one loop per cycle, and the Poincaré maps display a single return point. This clearly suggests that under these damping conditions, the system exhibits a stable periodic motion. To further clarify the presence of chaotic behavior, Figure 13, Figure 14, and Figure 15 present the dynamic orbit, Poincaré map, power spectrum, fractal dimension, and maximum Lyapunov exponent for the bearing geometric center, gear center, and rotor center at a dimensionless damping ratio of $\xi$ = 0.1. The dynamic orbits display highly irregular and non-periodic trajectories, while the corresponding Poincaré maps exhibit dense clusters of scattered return points. These features are characteristic of chaotic motion. The power spectra also reveal a broad range of excitation frequencies, indicating that the system does not follow simple periodic behavior at this damping level. To quantitatively confirm the existence of chaos, the fractal dimensions of the bearing, gear, and rotor centers were calculated as 2.05, 1.60, and 1.97, respectively. These non-integer values suggest the presence of strange attractors in the system. Furthermore, the maximum Lyapunov exponents for the same components were found to be 0.22, 0.12, and 0.16. Since all values are positive, they indicate sensitive dependence on initial conditions, which is a key indicator of chaotic dynamics. These findings clearly demonstrate that at $\xi$ = 0.1, the bearing, gear, and rotor all exhibit chaotic vibration behavior.

Figure 10.

Bifurcation diagram of the bearing geometric center, gear geometric center, and rotor geometric center in the horizontal direction using dimensionless damping ratio $\xi$ ( $\xi$ = 0.001∼0.397) for s = 0.33.

Figure 11.

Dynamic orbits and Poincaré maps of the bearing geometric center at dimensionless damping ratios $\xi$ = 0.001, 0.01, 0.1, 0.16, 0.24, and 0.3. (Arranged from left to right and top to bottom).

Figure 12.

Dynamic orbits and Poincaré maps of the gear geometric center at dimensionless damping ratios $\xi$ = 0.001, 0.01, 0.1, 0.16, 0.24, and 0.3. (Arranged from left to right and top to bottom).

Figure 13.

Dynamic orbit, Poincaré map, power spectrum, fractal dimension, and maximum Lyapunov exponent of the bearing geometric center at dimensionless damping ratio $\xi$ = 0.1.

Figure 14.

Dynamic orbit, Poincaré map, power spectrum, fractal dimension, and maximum Lyapunov exponent of the gear geometric center at dimensionless damping ratio $\xi$ = 0.1.

Figure 15.

Dynamic orbit, Poincaré map, power spectrum, fractal dimension, and maximum Lyapunov exponent of the rotor geometric center at dimensionless damping ratio $\xi$ = 0.1.

Conclusions

This study investigated the nonlinear dynamic behavior of a water-lubricated bearing-gear-rotor system under varying dimensionless parameters, with a focus on the effects of rotational speed ratio and damping ratio. Through bifurcation analysis, dynamic orbits, Poincaré maps, power spectra, fractal dimensions, and maximum Lyapunov exponents, the system’s transition between periodic and chaotic states was thoroughly examined. The results indicate that at lower rotational speed ratios (s<1.2), the system consistently exhibits non-periodic and unstable oscillations. As the speed increases, the system transitions into a periodic state; however, asynchronous behavior is observed between the gear and the bearing–rotor components at high speeds, deviating from classical oil-lubricated system behavior. This distinction may be attributed to the unique turbulence and mixed lubrication mechanisms inherent in water-lubricated systems.

Moreover, the dimensionless damping ratio $\xi$ was shown to significantly influence system stability. At low damping levels ( $\xi <0.14$ ), the system enters chaotic motion with large amplitude responses. As damping increases beyond $\xi =0.14$ , the system transitions into periodic motion, though a narrow range (0.214< $\xi <0.27$ ) reveals a reemergence of chaotic behavior with reduced amplitudes. For $\xi >0.27$ , stable periodic oscillation is regained. Detailed analysis at $\xi =0.1$ confirms the presence of chaos, supported by non-integer fractal dimensions and positive Lyapunov exponents across all components. These findings underscore the critical role of damping and speed in determining the nonlinear dynamic state of water-lubricated systems. In conclusion, the dynamic behavior of water-lubricated bearing systems is highly sensitive to operating parameters. Understanding these nonlinear responses is essential for the design, control, and fault diagnosis of rotating machinery operating under similar lubrication conditions.

In this study, the nonlinear dynamic responses of a water-lubricated bearing-gear-rotor system were investigated under different speeds and damping ratios. However, the asynchronous behavior among the rotor, bearings, and gears at high rotational speeds was not explicitly addressed. Such behavior may arise from complex fluid-structure interactions, time delays in contact dynamics, or localized instabilities. A more in-depth study focusing on this asynchronous mechanism is recommended for future research, which would help clarify its root causes and provide further guidance for high-speed marine transmission design.