Nonlinear dynamic analysis of gear-rotor-bearing system equipped with HSFD under hydraulic actuator active control

Abstract

In this study, the nonlinear dynamic analysis of the hybrid squeeze film damper mounted on a gear-rotor-bearing system is presented. The strong nonlinear oil film force, nonlinear rub-impact force between the rotor and stator, and the gear meshing force are considered in this study. The bifurcation diagrams are plotted using a dimensionless imbalance factor, damping coefficient, and speed ratio as control parameters in this study. The beginning of chaotic motion is identified by phase diagrams, power spectrum, Poincaré maps, bifurcation diagrams, maximum Lyapunov exponent, and the fractal dimension of the system. The numerical results provide useful insights into the design and development of a gear-rotor-bearing system for rotating machinery that operates at high speeds and highly nonlinear conditions. Hybrid squeeze film damper is used to suppress irregular dynamic behaviors and reduce greater amplitude, even if the rotational speed of the rotating machine is higher, which is also confirmed in this study.

Keywords

Hybrid squeeze film damper gear-rotor-bearing hydraulic active control

Introduction

Gears play an important role in the transmission system, especially for rotor bearing systems and turbomachinery, therefore, the consideration of gears transmission in rotor-bearing systems is needed. The rub-impact effect is a common and important phenomenon that occurs between the rotor and the stator. These effects may cause the rotating system to become instable. The dynamic analysis of the gear-rotor-bearing system with rub-impact effect is therefore important for analyzing rotating machinery. In order to control the dynamic stability of the gear-rotor-bearing system with the shock effect, we take the hybrid squeeze film damper (HSFD) to control this kind of system. The HSFD is actually a useful tool, and has been well developed and used in many fields. A squeeze-film damper bearing is actually a special type of journal bearing with its journal mechanically prevented from rotating, but it is free to vibrate within the clearance space. The squeezing action produces hydrodynamic forces in the fluid film. A squeeze-film damper bearing can be designed such that the journal can statically find its own position within the clearance, or be held centrally within the clearance, by retaining the springs. If the retaining springs are not used, the influence of contact and wear at zero speed will occur and it will add to the complexity of analysis of the squeeze-film damper bearing. In a squeeze-film damper bearing-rotor system, the fluid support pressure is generated entirely by the motion of the journal and depends on the viscosity of the lubricating fluid. However, the hydrodynamic pressure around the bearing is nonlinear, as it may cause fairly large vibrations of the rotor complicating the analysis of this system.

Some literatures have been studying on the subjects about squeeze-film dampers, and have also found many interesting and outstanding results. Morgan and Cameron¹ were the first researchers to focus on the analysis of porous metal bearings under hydrodynamic conditions. Many theoretical and experimental studies have performed on the analysis of the performance. Savoulides and Breuer² analyzed the low-order model of a very short hybrid gas bearing. The analysis shows that the best way to run bearings in a hybrid model is by using the bearing-dependent low-speed hydrostatics and high-speed hydrodynamics. Adiletta³ proposed a study of dynamics of an unbalanced rigid rotor suppored by squeeze film dampers with two-lobe wave bearings. Simulation results showed that two-lobe wave geometry influences the lubricating behaviors of the system and improves some unstable branches and the whirling phenomenon. Faris et al.⁴ showed their observation of a rotor-bearing system containing composite and non-composite squeeze film dampers. They found that the performance of squeeze film damper improved with increasing oil pressure and length/diameter ratio within the range tested. Zhang et al.⁵ presented the cell mapping method and experimental observations of multiple-objective design optimization of squirrel cage for squeeze film damper. Both the theoretical results and experimental results proved that the multi-objective optimization method offers a novel and promising tool for squirrel cage design in real industries. Hsu et al.⁶ discussed the nonlinear dynamic with effects of flywheel eccentricity in a turbine generator with a squeeze film damper. The effect of flywheel eccentricity was proved to be an important factor. Arghir et al.⁷ performed experiments on a flexure pivot tilting pad hybrid gas bearings for load-on-pad (LOP) and load-between-pad (LBP) configuration under various imbalanced conditions. San Andrés and Ryu⁸ experimented with a curved pivot tilt pad hybrid gas bearing to perform a pad loading and pad load configuration under various unbalanced conditions. San Andrés and Ryu⁹ introduced a strategy to employ an inexpensive air pressure regulator to control the supply pressure into the hybrid bearings, and both the theoretical results and experimental results showed that it reduced or even eliminated high amplitudes of rotor motion while crossing the system’s critical speeds. El-Shafei¹⁰ proposed a concept of actively controlling high speed rotating machinery through an experimental and analytical method. The control mechanism composed by a HSFD and found the usage of a HSFD will be able to reduce the amplitude of vibration of the rotating machinery. Li and Chen¹¹ showed that a HSFD with piezoelectric crystal electrohydraulic active control not only overcomes the bi-stable problem, but also restrains the shock in the rotor supporting system. Chen et al.¹² performed a theoretical analysis of a hybrid squeeze-film damper-mounted rigid rotor with active control and found an abundant amount of HSFD can suppress nonlinear dynamic responses. Amina et al.¹³ applied finite difference method to simulate the performance of three- and four-lobe hybrid journal bearings. They proved that the usage of a hydrostatic squeeze film damper loaded between pads is significantly influenced by the flow regimes, and the numerical results presented in this work can be useful to the bearing designers. Recently, Chang-Jian^14–16 presented gear dynamic analysis mounted with HSFD combined with turbulent flow effect and rub-impact effect in the bearings, and proved that HSFD can be used to suppress those irregular or undesired vibrations occurring in the rotating machines.

To avoid the undesirable non-synchronous vibrations, hybrid squeeze-film dampers (a combination of hydrostatic and hydrodynamic ranges) with active controls are designed and simulated in this study. Numerical examples are given for a gear-rotor-bearing system supported by HSFD with two pairs of hydraulic actuators with PD controllers. To help the gear-rotor-bearing system to avoid an undesirable non-synchronous vibration, a suitable increased proportional gain is also applied to this system. The reminder of this paper is organized as follows. The Mathematical modeling section derives dynamic models for the gear-rotor-bearing system with the nonlinear suspension effect, nonlinear rub-impact effect, strongly nonlinear gear mesh force, and strongly nonlinear oil-film force under the suspension of HSFD. The Numerical results and discussions section presents the numerical analysis results obtained for the behavior of the gear-rotor-bearing system under various operational conditions. Finally, the Conclusions section presents some brief conclusions.

Mathematical modeling

Figure 1 shows a gear-rotor-bearing system equipped with HSFD under the assumptions of nonlinear rub-impact force, nonlinear fluid film force, nonlinear suspension effect, and nonlinear gear meshing force. Nonlinear rub-impact force can be found between the rotor and the stator. The bearing consists of four hydrostatic chambers and four hydrodynamic regions. The structure of this type of bearings should be popularized to consist of 2 N (N = 2, 3, 4…) hydrostatic chambers, and 2 N hydrodynamic regions. In this study, the oil pressure distribution model in HSFD was proposed to integrate the pressure distribution in the dynamic and static pressure regions, as described in Section 2.1.

Figure 1.

Schematic illustration of hybrid squeeze-film damper mounted on the gear-rotor-bearing system, cross section of HSFD and model of force diagram for pinion and gear.

The instant oil film supporting force for HSFD

In order to analyze the pressure distribution, the Reynolds equation of constant lubricant properties and incompressibility should be assumed, and then the Reynolds equation is as follows

\frac{1}{R^{2}} \frac{\partial}{\partial θ} (\frac{h^{3}}{μ} \frac{\partial p}{\partial θ}) + \frac{\partial}{\partial z} (\frac{h^{3}}{μ} \frac{\partial p}{\partial z}) = 6 \frac{U}{R} \frac{\partial h}{\partial θ} + 12 \frac{\partial h}{\partial t}

(1)

The support area of the HSFD should be divided into three areas: static pressure zone, rotational direction dynamic pressure zone, and axial dynamic pressure zone, as shown in Figure 1. In the part of HSFD with $- a \leq z \leq a$ , the long bearing theory (Appendix 1) is assumed and the Reynolds equation is solved with the boundary condition of static pressure region $p_{c, i}$ acquiring the pressure distribution $p_{0} (θ)$ . In the part of HSFD with $a \leq | z | \leq \frac{L}{2}$ , the short bearing theory (Appendix 1) is assumed and solves the Reynolds equation with the boundary condition of $p (z, θ) |_{z = \pm a} = p_{0} (θ)$ and $p (z, θ) |_{z = \pm L / 2} = 0$ , yielding the pressure distribution in the axis direction dynamic pressure region $p (z, θ)$ . Finally, a formula of pressure distribution in the whole supporting region is obtained.

According to the aforementioned conditions, the instantaneous oil film pressure distribution is as follows.

The instant pressure in rotating direction within the range of $- a \leq z \leq a$ is

p_{0} (θ) = {\begin{matrix} \begin{array}{l} p_{c, i}; \frac{π}{2} (i - 1) - \frac{β}{2} - φ_{b} \leq θ \leq \frac{π}{2} (i - 1) + \frac{β}{2} - φ_{b} \\ p_{i} (θ); \frac{π}{2} (i - 1) + \frac{β}{2} - φ_{b} \leq θ \leq \frac{π}{2} i - \frac{β}{2} - φ_{b} \end{array} & , i = 1, 2, 3, 4 \end{matrix}

(2)

where

\begin{array}{l} p_{i} (θ) = p_{c, i} + \frac{6 μ R^{2}}{δ^{2}} \frac{\dot{ε}}{ε} [\frac{1}{{(1 + ε \cos θ)}^{2}} - \frac{1}{{(1 + ε \cos θ_{i 1})}^{2}}] \\ + c_{1} \int_{θ_{i 1}}^{θ} \frac{1}{δ^{3} {(1 + ε \cos θ)}^{3}} d θ - \int_{θ_{i 1}}^{θ} \frac{12 {\dot{φ}}_{b} μ δ ε R^{2} \cos θ}{δ^{3} {(1 + ε \cos θ)}^{3}} d θ \\ c_{1} = \frac{p_{c, i + 1} - p_{c, i} - \frac{6 μ R^{2}}{δ^{2}} \frac{\dot{ε}}{ε} [\frac{1}{{(1 + ε \cos θ_{i 2})}^{2}} - \frac{1}{{(1 + ε \cos θ_{i 1})}^{2}}] + \int_{θ_{i 1}}^{θ_{i 2}} \frac{12 {\dot{φ}}_{b} μ δ ε R^{2} \cos θ}{δ {(1 + ε \cos θ)}^{3}} d θ}{\int_{θ_{i 1}}^{θ_{i 2}} \frac{1}{δ^{3} {(1 + ε \cos θ)}^{3}} d θ} \\ θ_{i 1} = (i - 1) \frac{π}{2} + \frac{β}{2} - φ_{b}, θ_{i 2} = i \frac{π}{2} - \frac{β}{2} - φ_{b}, i = 1, 2, 3, 4 \end{array}

(3)

The instant pressure in the axis direction within the range of $a \leq | z | \leq \frac{L}{2}$ is

p (θ, z) = (\frac{L}{2} - | z |) {[A_{1} (θ) {\dot{φ}}_{b} ε + A_{2} (θ) \dot{ε}] (a - | z |) + p_{0} (θ) \frac{1}{L / 2 - a}}

(4)

where

A_{1} (θ) = \frac{6 μ δ \sin θ}{δ^{3} {(1 + ε \cos θ)}^{3}}; A_{2} (θ) = \frac{6 μ δ \cos θ}{δ^{3} {(1 + ε \cos θ)}^{3}}

The instantaneous oil film force of the different elements is determined by integrating equations (2) and (4) onto the area of the journal sleeve. In the static pressure zone, the force is shown

F_{r s} = \sum_{i = 1}^{4} p_{c i} 2 a R [\sin (\frac{π}{2} (i - 1) + \frac{β}{2} - φ_{b}) - \sin (\frac{π}{2} (i - 1) - \frac{β}{2} - φ_{b})]

(5)

F_{τ s} = \sum_{i = 1}^{4} p_{c i} 2 a R [\cos (\frac{π}{2} (i - 1) - \frac{β}{2} - φ_{b}) - \cos (\frac{π}{2} (i - 1) + \frac{β}{2} - φ_{b})]

(6)

In the dynamic pressure region in the direction of rotation, the force is shown

F_{r c} = \sum_{i = 1}^{4} \int_{θ_{i 1}}^{θ_{i 2}} p_{i} (θ) R 2 a \cos θ d θ

(7)

F_{τ c} = \sum_{i = 1}^{4} \int_{θ_{i 1}}^{θ_{i 2}} p_{i} (θ) R 2 a \sin θ d θ

(8)

In the axial dynamic pressure region, the force is shown.

F_{r a} = \int_{- L / 2}^{- a} d z \int_{0}^{2 π} p (θ, z) R \cos θ d θ + \int_{L / 2}^{a} d z \int_{0}^{2 π} p (θ, z) R \cos θ d θ

(9)

F_{τ a} = \int_{- L / 2}^{- a} d z \int_{0}^{2 π} p (θ, z) R \sin θ d θ + \int_{L / 2}^{a} d z \int_{0}^{2 π} p (θ, z) R \sin θ d θ

(10)

The damping forces generated in the radial and tangential directions are determined by summing the above supporting forces. It is shown below.

F_{r} = F_{r s} + F_{r c} + F_{r a}

(11)

F_{τ} = F_{τ s} + F_{τ c} + F_{τ a}

(12)

The gear meshing force

$L_{p y}$ and $L_{g y}$ are the centrifugal forces of the pinion and gear vertical gear meshing direction, $G_{p y}$ and $G_{g y}$ are the inertial forces in the direction of the vertical gear meshing of the pinion and gear, $W_{c x}$ is the dynamic gear meshing force in the horizontal direction, and $W_{c y}$ is the dynamic gear mesh force in the vertical direction. $L_{p y}$ , $L_{g y}$ , $G_{p y}$ , $G_{g y}$ , $W_{c x}$ and $W_{c y}$ can be performed as

L_{p y} = m_{p} e_{p} ω_{p}^{2} \sin θ_{1}

(13)

L_{g y} = m_{g} e_{g} ω_{g}^{2} \sin θ_{2}

(14)

G_{p y} = m_{p} e_{p} {\ddot{θ}}_{1} c o s θ_{1}

(15)

G_{g y} = m_{g} e_{g} {\ddot{θ}}_{2} c o s θ_{2}

(16)

W_{c x} = C_{m} ({\dot{X}}_{p} - {\dot{X}}_{g} - e_{p} ω \sin (ω t)) + K_{m} (X_{p} - X_{g} - e_{p} \cos (ω t))

(17)

W_{c y} = C_{m} ({\dot{Y}}_{p} - {\dot{Y}}_{g} - e_{p} ω \cos (ω t)) + K_{m} (Y_{p} - Y_{g} - e_{p} \sin (ω t))

(18)

The rub-impact force

From Figure 1, the radial impact force $f_{n}$ and the tangential rub force $f_{t}$ could be expressed as

f_{n} = (e - δ) k_{c}

(19)

f_{t} = (f + b v) f_{n}, if e \geq δ

(20)

Then we could acquire the rub-impact forces in the horizontal and vertical directions.

R_{x} = \frac{(e - δ) k_{c}}{e} [X - (f + b v) Y]

(21)

R_{y} = \frac{(e - δ) k_{c}}{e} [(f + b v) X + Y]

(22)

Dynamics equation

Applying the principle of force balance, the force acting on the center of journal 1, i.e. O_j1 (X_j1, Y_j1), center of journal 2, i.e. O_j2 (X_j2, Y_j2), equations of motion of O_g(X_g, Y_g) and O_p(X_p, Y_p) in Cartesian coordinate form and equations of motion of the center of bearing 1 (X₁, Y₁), the center of bearing 2 (X₂, Y₂) and the center of the rotor under the assumption of nonlinear suspension are given by

m_{p} {\ddot{X}}_{p} + C {\dot{X}}_{p} + K X_{p} = W_{c x} + F_{x 1} + K X_{p 0}

(23)

m_{p} {\ddot{Y}}_{p} + C {\dot{Y}}_{p} + K Y_{p} = L_{p y} - G_{p y} - W_{c y} - m_{p} g + F_{y 1} + K Y_{p 0}

(24)

m_{g} {\ddot{X}}_{g} + C {\dot{X}}_{g} + K X_{g} = - W_{c x} + F_{x 2} + K X_{g 0}

(25)

m_{g} {\ddot{Y}}_{g} + C {\dot{Y}}_{g} + K Y_{g} = L_{p y} - G_{p y} + W_{c y} - m_{g} g + F_{y 2} + K Y_{g 0}

(26)

m_{r} {\ddot{X}}_{r} + C_{r} {\dot{X}}_{r} + K_{r} X_{r} = m_{r} ρ ω^{2} \cos φ + R_{x} + F_{x} + K_{r} X_{0}

(27)

m_{r} {\ddot{Y}}_{r} + C_{r} {\dot{Y}}_{r} + K_{r} Y_{r} = m_{r} ρ ω^{2} \sin φ + R_{y} + F_{y} + K_{r} Y_{0} - m_{r} g

(28)

where F_x1, F_y1, F_x2, F_y2, F_x and F_y are given as follows.

F_{x 1} = f_{e 1} cos ϕ_{1} + f_{ϕ 1} sin ϕ_{1} = K (X_{p} - X_{j 1}) / 2

(29)

F_{y 1} = f_{e 1} sin ϕ_{1} - f_{ϕ 1} cos ϕ_{1} = K (Y_{p} - Y_{j 1}) / 2

(30)

F_{x 2} = f_{e 2} cos ϕ_{2} + f_{ϕ 2} sin ϕ_{2} = K (X_{g} - X_{j 2}) / 2

(31)

F_{y 2} = f_{e 2} sin ϕ_{2} - f_{ϕ 2} cos ϕ_{2} = K (Y_{g} - Y_{j 2}) / 2

(32)

F_{x} = K_{p 2} (X_{r} - X_{g})

(33)

F_{y} = K_{p 2} (Y_{r} - Y_{g})

(34)

In which f_e1 and f_φ1 are the viscous damping forces in the radial and tangential directions for the center of journal 1, respectively, and f_e2 and f_φ2 are the viscous damping forces in the radial and tangential directions for the center of journal 2, respectively, and f_e and f_φ are the viscous damping forces in the radial and tangential directions, respectively.

The origin of the o-xyz-coordinate system is taken to be the bearing center O _b . Dividing these two equations by $m c ω^{2}$ and defining a nondimensional time $ϕ = ω t$ and a speed parameter $s = \frac{ω}{ω_{n}}$ , one obtains the following non-dimensionalized equations of motion

\begin{array}{l} x_{p}^{″} = - \frac{2 ξ_{2}}{s} {x^{'}}_{p} - \frac{1}{s^{2}} (x_{p} - x_{1} - ε_{1} \cos φ_{1}) + β \cos (ϕ / 4) - \frac{2 ξ_{3}}{s} ({x^{'}}_{p} - {x^{'}}_{g} - E_{p} \sin ϕ) \\ - \frac{Λ}{s^{2}} (x_{p} - x_{g} - E_{p} \cos ϕ) + \frac{B}{s} (\frac{x_{p} F_{r} - y_{p} F_{τ}}{ε}) + \frac{x_{p 0}}{s^{2}} \end{array}

(35)

\begin{array}{l} y_{p}^{″ ​} = - \frac{2 ξ_{2}}{s} {y^{'}}_{p} - \frac{1}{s^{2}} (y_{p} - y_{1} - ε_{1} \sin φ_{1}) + β \sin (ϕ / 4) - \frac{2 ξ_{3}}{s} ({y^{'}}_{p} - {y^{'}}_{g} - E_{p} \cos ϕ) \\ - \frac{Λ}{s^{2}} (y_{p} - y_{g} - E_{p} \sin ϕ) - \frac{f}{s^{2}} + \frac{B}{s} (\frac{y_{p} F_{r} + x_{p} F_{τ}}{ε}) + \frac{y_{p 0}}{s^{2}} \end{array}

(36)

\begin{array}{l} x_{g}^{″ ​} = - \frac{2 ξ_{4}}{s} {x^{'}}_{g} - \frac{1}{s^{2}} (x_{g} - x_{2} - ε_{2} \cos φ_{2}) + β_{g} \cos (ϕ / 8) + \frac{2 ξ_{5}}{s} ({x^{'}}_{p} - {x^{'}}_{g} - E_{p} \sin ϕ) \\ - \frac{Λ_{g}}{s^{2}} (x_{p} - x_{g} - E_{p} \cos ϕ) + \frac{B}{s} (\frac{x_{g} F_{r} - y_{g} F_{τ}}{ε}) + \frac{x_{g 0}}{s^{2}} \end{array}

(37)

\begin{array}{l} y_{g}^{​ ″ ​} = - \frac{2 ξ_{4}}{s} {y^{'}}_{g} - \frac{1}{s^{2}} (y_{g} - y_{2} - ε_{2} \sin φ_{2}) + β_{g} \sin (ϕ / 8) + \frac{2 ξ_{5}}{s} ({y^{'}}_{p} - {y^{'}}_{g} - E_{p} \cos ϕ) \\ + \frac{Λ_{g}}{s^{2}} (y_{p} - y_{g} - E_{p} \sin ϕ) - \frac{f_{g}}{s^{2}} + \frac{B}{s} (\frac{y_{g} F_{r} + x_{g} F_{τ}}{ε}) + \frac{y_{g 0}}{s^{2}} \end{array}

(38)

x_{r}^{″} + \frac{2 ξ_{2}}{s_{3}} x_{r}^{'} + \frac{1}{s_{3}^{2}} x_{r} = β_{r} \cos φ + \frac{R_{x}}{p} + \frac{B}{s} (\frac{x_{r} F_{r} - y_{r} F_{τ}}{ε}) + \frac{x_{0}}{s^{2}}

(39)

y_{r}^{″} + \frac{2 ξ_{2}}{s_{3}} y_{r}^{′} + \frac{1}{s_{3}^{2}} y_{r} = β_{r} \sin φ + \frac{R_{y}}{p} - \frac{f_{r}}{s^{2}} + \frac{B}{s} (\frac{y_{r} F_{r} + x_{r} F_{τ}}{ε}) + \frac{y_{0}}{s^{2}}

(40)

Equations (35) to (40) describe the nonlinear dynamic of the gear-rotor-bearing system. In the current study, the approximate solutions of these coupled nonlinear differential equations are obtained using a fourth-order Runge–Kutta numerical scheme. In order to control the gear-rotor-bearing system, four hydraulic actuators are designed to control the oil flow in each oil chamber to change the pressure of the four oil chambers (for i = 1, 2, 3, 4), as shown in Figure 2 and Figure 3. The kinetic equations in the four oil chambers are discussed and well presented in Chen et al.¹²

Figure 2.

The flow rate control structure of HSFD (only shows the 2nd oil chamber).

Figure 3.

The change of flow rate in ith oil chamber.

Numerical results and discussions

Without hydraulic active control

In the simulations, the system parameters are assigned the following values: $B = 0.3$ , $β = 0.35$ , $β_{g} = 0.35$ , $f = 0.2$ , $f_{g} = 0.2$ , $ξ_{1} = 0.01$ , $ξ_{2} = 0.02$ , $ξ_{3} = 0.015$ , $ξ_{4} = 0.02$ , $ξ_{5} = 0.015$ , $C_{1 p} = 2.0$ , $Λ = 0.04$ , $Λ_{g} = 0.04$ .

The fourth-order Runge-Kutta method is used to solve the nonlinear dynamic equations in equations (35) to (40) of a gear-rotor-bearing system with strong nonlinear oil film force, nonlinear rub-impact force, and gear meshing force. The time step in the iterative solution procedure was assigned a value of π/300 and the termination criterion was specified as an error tolerance of less than 0.0001. The time series data corresponding to the first 800 revolutions of the two gears were deliberately excluded from the dynamic analysis to ensure that the analyzed data related to the steady-state conditions. The numerical results solved by the fourth-order Runge-Kutta method were used to generate the dynamic trajectories, Poincaré maps and bifurcation diagrams of the system in order to obtain a basic understanding of its dynamic behavior. The maximum Lyapunov exponent and the fractal dimension measure were then used to identify the onset of chaotic motion. For convenience, only the data of the displacements in the vertical direction were used to generate the diagrams. In contrast to the Fourier transform-based techniques and bifurcation diagrams, which provide only a general indication of the change from periodic motion to chaotic behavior, dimensional measures allow chaotic signals to be differentiated from random signals. Although many dimensional measures have been proposed, the most commonly applied measure is the correlation dimension d_G defined by Grassberger and Procaccia due to its computational speed and the consistency of its results. However, before the correlation dimension of a dynamic system flow can be evaluated, it is first necessary to generate a time series of one of the system variables using a time-delayed pseudo-phase-plane method. Assume an original time series of x_i = {x(iτ); i = 1,2,3, … N}, where τ is the time delay (or sampling time). If the system is acted upon by an excitation force with a frequency of ω, the sampling time, τ, is generally chosen such that it is much smaller than the driving period. The delay coordinates are then used to construct an n-dimensional vector X = (x(jτ), x[(j + 1) τ], x[(j + 2) τ], … ., x[(j + n-1) τ]), where j = 1,2,3, … (N-n + 1). The resulting vector comprises a total of (N−n + 1) vectors, which are then plotted in an n-dimensional embedding space. Importantly, the system flow in the reconstructed n-dimensional phase space retains the dynamic characteristics of the system in the original phase space. Simply stated, if the system flow has the form of a closed orbit in the original phase plane, it also forms a closed path in the n-dimensional embedding space. Similarly, if the system exhibits a chaotic behavior in the original phase plane, its path in the embedding space will also be chaotic. The characteristics of the attractor in the n-dimensional embedding space are generally tested using the function $\sum_{i, j = 1}^{N} H (r - | x_{i} - x_{j} |)$ to determine the number of pairs (i, j) lying within a distance |x_i – x_j| < r in ${x_{i}}_{i = 1}^{N}$ , where H denotes the Heaviside step function, N represents the number of data points, and r is the radius of an n-dimensional hyper-sphere. For many attractors, this function exhibits a power law dependence on r as r →0, i.e. $c (r) \propto r^{d_{G}}$ . Therefore, the correlation dimension, d_G, can be determined from the slope of a plot of [log c(r)] versus [log r]. Grassberger and Proccacia¹⁷ showed that the correlation dimension represents the lower bound to the capacity or fractal dimension d_c, and approaches its value asymptotically when the attracting set is distributed more uniformly in the embedding phase space. A set of points in the embedding space is said to be fractal if its dimension has a finite non-integer value. Otherwise, the attractor is referred to as a ‘strange attractor’. To establish the nature of the attractor, the embedding dimension is progressively increased, causing the slope of the characteristic curve to approach a steady state value. This value is then used to determine as to whether the system has a fractal structure or a strange attractor structure. If the dimension of the system flow is found to be fractal (i.e. to have a non-integer value), the system is judged to be chaotic.

The dimensionless rotating speed ratio s is one of the most important parameters, and is commonly used as a bifurcation control parameter when analyzing the dynamic characteristics of the system. Accordingly, the dynamic behaviors of the current gear-rotor-bearing system were examined using the dimensionless rotating speed ratio s as a bifurcation control parameter. The bifurcation diagrams for the vertical displacement of gear center and rotor center against the dimensionless rotating speed ratio are provided first. Figures 4 and 5 are bifurcation diagrams of a gear-rotor-bearing system with and without HSFD, respectively. According to the simulation results of the bifurcation diagrams, as shown in Figures 4 and 5, they reveal that the dynamic behaviors will perform more regular dynamic behaviors, and the amplitudes are also smaller for the case with HSFD in such strongly nonlinear dynamic systems. The use of HSFD can impressively suppress non-periodic dynamic behaviors and reduce amplitude even at the higher rotational speeds. A bifurcation diagram with a HSFD case (Figure 4) shows that the geometric center of the gear and rotor in the vertical direction performs a synchronous 1 T-periodic motion at low values of the rotating speed ratio, i.e. s < 0.80, and then the non-periodic motion, or the so-called chaotic motion, can be found as the dimensionless rotating speed ratio is increased over s = 0.80. It can be seen that the vibration amplitude is larger for higher rotating speed than low rotating speed. With the increase of the rotating speed, dynamic characteristics are disorderly and behave with non-periodic or even chaotic responses. To illustrate the non-periodic motions more clearly, phase diagrams, power spectra, Poincaré maps, Lyapunov exponents, and the fractal dimension are applied to identify the onset of chaotic motions. We highlight a case from Figure 5 to show the dynamic motion of the geometrical center of the gear and geometrical center of rotor without HSFD, and show in Figure 6. It represents the dynamic behavior shown with phase diagrams, power spectra, Poincaré maps, Lyapunov exponents, and the fractal dimensions of pinion center at s = 1.80. The motions shown in the phase diagrams performed a disordered dynamic behavior; the power spectra reveals numerous excitation frequencies; the returning points in the Poincaré maps is the so-called geometrically fractal structures; the maximum Lyapunov exponent is positive, and the fractal dimensions are also found to be 1.33 for $y_{p}$ at s = 1.80. Therefore, the chaotic motion is found at s = 1.80.

Figure 4.

Bifurcation diagrams of gear-rotor-bearing system with HSFD using dimensionless rotating speed ratio, s, as bifurcation parameter.

Figure 5.

Bifurcation diagrams of gear-rotor-bearing system without HSFD using dimensionless rotating speed ratio, s, as bifurcation parameter.

Figure 6.

Simulation results obtained for gear-rotor-bearing system without HSFD and s = 1.8 $(y_{g})$ . Phase diagram

The damping coefficient is also an important control parameter for observing the dynamic responses of the vibration system. From a physical point of view, the presence of damping coefficients in the vibration system will cause them to become a dissipative system. Figure 7 shows the bifurcation diagrams of gear geometric center and rotor geometric center with HSFD using a dimensionless damping coefficient, $ξ$ , as bifurcation parameter at s = 2.00. When compared to the above section, we know that the dynamic response is non-periodic or even chaotic at s = 2.00. The dynamic characteristics persist disorder or even chaotic under all the damping coefficients at this rotating speed, i.e. s = 2.00. As the damping coefficient increases, the amplitude of the vibration becomes smaller. One can find that the usage of damping can no longer suppress the aperiodic motion or the disordered motion at the higher rotating speeds. Therefore, the damping effect may not perform so well on such a strongly nonlinear system, and then the control tools will be important for those situations. The phase diagram is as shown in Figure 8, and it behaves in full disorder. In point of fact, numerous excitation frequencies can also be found in the power spectra, and some closed curves are also shown in the Poincaré maps. Meanwhile, the Lyapunov exponent performs negatively, and the fractal dimensions are found to be 0.41 for gear geometric center, $y_{g}$ at $ξ$ = 0.025. According to the identifying tools, quasi-periodic motion has been found in this case.

Figure 7.

Bifurcation diagrams of gear-rotor-bearing system with HSFD using dimensionless damping coefficient, $ξ$ , as bifurcation parameter (s = 2.00).

Figure 8.

Simulation results obtained for gear-rotor-bearing system with HSFD and $ξ$ = 0.025 ( $y_{g}$ ).

The dimensionless unbalance coefficient also plays an important role to study the dynamic response of the many rotating machineries. Therefore, it is valuable to study the various dynamic characteristics of a system for diagnosing some dynamic characteristics and faults by using an unbalance coefficient as the control parameters to plot the bifurcation diagrams. Figure 9 shows the bifurcation diagrams of the geometric center of the gear and rotor with HSFD using the dimensionless unbalance coefficient as a bifurcation parameter. It can be seen that the system response exhibits periodic motion at the interval of $0.001 \leq β \leq 0.243$ , but displays the strongly non-periodic motions when $β > 0.243$ . The result is unlike the usual chaotic mode (1 T => 2 T => 4 T => 8 T => 16 T => 32 T … => chaos or period => quasi-period => chaos), as it behaves with the periodic motions at a lower value of β, and suddenly develops into an aperiodic motion without any transition. In Figure 10, the irregular dynamic orbit has been shown in phase diagrams, and it also reveals that periodic motions don’t exist in those situations. Nevertheless, the irregular orbits cannot afford one with enough information about what types of dynamic response it will be. Therefore, the authors then use other distinguishing methods to identify the dynamic behaviors of the systems under individual situation, i.e. power spectra, Poincaré section, and the maximum Lyapunov exponent and fractal dimension. The power spectra show the numerous excitation frequencies and posit they are composed of many different frequencies. The Poincaré maps display two closed curves. The maximum Lyapunov exponent is negative and the fractal dimensions are found to be 1.39 for $y_{g}$ at $β = 0.52$ . The comparisons of the bifurcation diagrams of the gear geometric center using dimensionless rotating speed ratio, s, as the bifurcation parameter with different values of β and ξ are as shown in Figure 11(a) to (e). According to those bifurcation diagrams, dimensionless unbalance coefficient $β$ and dimensionless damping parameter $ξ$ are really important factors affecting the dynamic stability of the system. Comparison of Figure 11(a) to (c) shows more periodic motions when the dimensionless unbalance coefficient is low. Comparing Figure 11(b) and (e), greater value of the dimensionless damping parameter will suppress the non-periodic motions or even chaotic motions at low or medium rotating speeds, but it is possible to fail or suppress the non-periodic vibration at higher speeds.

Figure 9.

Bifurcation diagrams of gear-rotor-bearing system with HSFD using dimensionless unbalance coefficient, $β$ , as bifurcation parameter (s = 2.00).

Figure 10.

Simulation results obtained for gear-rotor-bearing system with HSFD and $β$ = 0.52 ( $y_{g}$ ).

Figure 11.

Comparison of bifurcation diagrams of gear geometric center using dimensionless rotating speed ratio, s, as bifurcation parameter with different values of $β$ and $ξ$ . (a) $β$ = 0.35 and $ξ$ = 0.005; (b) $β$ = 0.28 and $ξ$ = 0.005; (c) $β$ = 0.12 and $ξ$ = 0.005; (d) $β$ = 0.12 and $ξ$ = 0.0025; (e) $β$ = 0.28 and $ξ$ = 0.1.

With hydraulic active control

The initial state of the retaining spring is given to be $(X_{0}, Y_{0}) = (0, - 0.6)$ , and then a bifurcation diagram of the rotor trajectory in the horizontal direction with Kp = 0.01 and Kd = 0.01 is as shown in Figure 12. It can be found that the trajectory of the center of the rotor is in a disordered state at s > 5.4. The rotor trajectory at s = 6.0, the Poincaré plot and the time response are shown in Figure 13(a) to (d). This indicates that the trajectory of the center of the rotor is disordered. The structure of the Poincaré diagram in Figure 13(b) is confusing, and the time response of the rotor trajectory in the horizontal direction is irregular. In order to more profoundly understand the disordered dynamic characteristics of the rotor trajectory, the maximum Lyapunov exponent is applied to identify the chaotic motion in this state. A positive maximum Lyapunov exponent is obtained in Figure 14, which indicates that the rotor trajectory is in a chaotic motion state. Figure 15 shows the time response of control input $\bar{u} (1)$ and $\bar{u} (2)$ at s = 6.0. This shows that the dynamic of control inputs are also irregular at this time. To avoid the system operating in a chaotic motion at s = 6.0, an increased proportional gain $k_{p} = 0.1$ is applied from non-dimensional time $φ$ =1570. Figure 16 shows that the abrupt change in $k_{p}$ at $φ = 1570$ causes the rotor center trajectory to be regulated to a synchronous periodic motion after a transient response. The response of the control inputs is shown in Figure 17, as it can be seen that the time response of the control input $\bar{u} (1)$ and $\bar{u} (2)$ will become larger after an abrupt change in $k_{p}$ . The control inputs $\bar{u} (1)$ and $\bar{u} (2)$ are both steered to a periodic motion after a transient response.

Figure 12.

The bifurcation diagram of dynamic trajectory with Kp = 0.01 and Kd = 0.01.

Figure 13.

Aperiodic motion of rotor center at s = 6.0 with Kp = 0.01 and Kd = 0.01; (a) Rotor trajectory; (b) Poincare’ map in X(nT)-Y(nT) plane; (c), (d) Time response of rotor trajectory.

Figure 14.

The maximum Lyapunov exponent of rotor trajectory plotted as a function of the number of drive cycles at s = 6.0 with Kp = 0.01 and Kd = 0.01.

Figure 15.

The time response of $\bar{u} (1)$ and $\bar{u} (2)$ at s = 6.0 with Kp = 0.01.

Figure 16.

The time responses of rotor trajectories at s = 6.0 with Kp = 0.01 changes to Kp = 0.1 from non-dimensional time $φ$ = 1570.

Figure 17.

The time responses of $\bar{u} (1)$ and $\bar{u} (2)$ at s = 6.0 with Kp = 0.01 changes to Kp = 0.1 from non-dimensional time $φ$ =1570.

Conclusions

The numerical analysis of the nonlinear dynamic response of a gear-rotor-bearing system with nonlinear oil film force, nonlinear rubbing force, and gear meshing force is proposed in this study. The dynamic responses of the system are analyzed with reference to its dynamic trajectory, power spectrum, Poincaré maps, bifurcation diagram, maximum Lyapunov exponent, and fractal dimension. The bifurcation parameters of the dimensionless rotational speed ratio, the dimensionless unbalance parameter, and the dimensionless damping coefficient are used to study the nonlinear dynamic characteristics of the gear-rotor-bearing system. It can be observed that very rich aperiodic and chaotic motions can be found in the system. In this study, no general way of entering chaos was found (1 T => 2 T => 4 T => 8 T => 16 T => 32 T … => chaos or period => quasi-period => chaos). The authors also demonstrate that HSFD can suppress irregular dynamic behaviors and reduce the amplitude of rotating machines. In summary, the numerical results presented in this study provide a detailed understanding of the nonlinear dynamic responses of a gear-rotor-bearing system with some nonlinear effects. The numerical results will provide some useful viewpoints in designing those types of systems so that the chaotic behavior can be avoided, and thereby reduce the amplitude of vibration within the system, or even extend the life of the system. The authors have also proved that HSFD can be used to suppress the irregular or undesired vibrations found in the rotating machines.

Footnotes

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.

Appendix

References

Morgan

Cameron

Mechanism of lubrication in porous metal bearings. In:Proceedings of the conference on lubrication and wear, London, UK, May 1957, pp.151–157.

Savoulides

Breuer

KS.

Low-order models for very short hybrid gas Bearings. J Tribol 2001; 123: 368–375.

Adiletta

Bifurcating behaviour of a rotor on two-lobe wave squeeze film damper. Tribol Int 2015; 92: 72–83.

Faris

Khalid

Albagul

, et al. Investigation of the performance of a rotor-bearing system containing composite and non-composite squeeze film dampers. Mater Des 2012; 34: 340–345.

Zhang

Han

, et al. Multiple-objective design optimization of squirrel cage for squeeze film damper by using cell mapping method and experimental validation. Mech Mach Theory 2019; 132: 66–79.

Hsu

Wang

Shiau

, et al. Nonlinear dynamic study on effects of flywheel eccentricity in a turbine generator with a squeeze film damper. Proc Eng 2014; 79: 397–406.

Arghir

Frene

Finite-volume solution of the compressible Reynolds equation: linear and non-linear analysis of gas bearings. Prog Inst Mech Eng Part J: J Eng Tribol 2006; 220: 617–699.

San Andrés

Ryu

Flexure pivot tilting pad hybrid gas bearings: operation with worn clearances and two load-pad configurations. J Eng Gas Turb 2008; 130: 042506.

San Andrés

Ryu

Hybrid gas bearings with controlled supply pressure to eliminate rotor vibrations while crossing system critical speeds. J Eng Gas Turb 2008; 130: 062505.

10.

El-Shafei

Experimental and analytical investigation of hybrid squeeze film dampers. J Eng Gas Turb 1993; 115: 353–359.

11.

Chen

CL.

Research on chaos dynamic characteristics of hybrid squeeze film damper with electro-hydraulic active control. In: IEEE International Conference on Cybernetics & Intelligent Systems (CIS) and Robotics, Automation & Mechatronics, Chengdu, China, 3-6 June, 2008.

12.

Chen

Yau

YH.

Subharmonic and chaotic motions of a hybrid squeeze-film damper-mounted rigid rotor with active control. ASME J Vibr Acoust 2002; 124: 198–208.

13.

Amina

Ahmed

Performance analysis of four-pad hydrostatic squeeze film dampers loaded between pads under laminar and turbulent flow conditions. Tribol Mat Surf Int 2018; 12: 59–70.

14.

Chang-Jian

CW.

Gear dynamics analysis with turbulent journal bearings mounted hybrid squeeze film damper – chaos and active control analysis. J Comput Non Dynam 2015; 10: 1–11.

15.

Chang-Jian

CW.

Dynamic analysis of a hybrid squeeze film damper mounted rub-impact rotor-stator system. J Appl Math 2012; 2012: 1–21.

16.

Chang-Jian

Yau

HT.

Non-linear dynamic analysis of hybrid squeeze-film damper-mounted gear-bearing system and hydraulic active control. Proc Inst Mech Eng Part K J Multi-Body Dynam 2010; 224: 249–259.

17.

Grassberger P and Proccacia I. Characterization of Strange Attractors. Physical Review Letter 1983; 50: 346–351.