Investigation on the dynamic behaviors of a rod fastening rotor based on an analytical solution of the oil film force of the supporting bearing

Abstract

The dynamic characteristics of rod fastening rotor supported by a finite journal bearing are investigated in this study. To model the dynamic behaviors of the bearing-rotor system, the oil film force of bearing is calculated by approximately solving the Reynolds equation with the variables separation method and Sturm–Liouville theory, and then a motion equation is developed with consideration of the contact and gyro effects of the disks of the rotor. To solve the motion equation with small error and excellent stability, an improved Newmark method is proposed. On this basis, the dynamics characteristics of the rod fastening rotor are analyzed for different rotor speeds, disk eccentricities, shaft bearing stiffness, and contact stiffness. And the orbits of the rod fastening rotor and integral rotor are compared. The numerical results indicate that the analytical solution of the oil film force has higher computational efficiency than the finite difference method. The rod fastening rotor shows higher stability than the integral rotor, and exhibits rich dynamic behaviors, such as periodic, qusi-periodic, period-2, period-4, and period-6.

Keywords

Rod fastening rotor Newmark method Sturm–Liouville theory approximate solution stability

Introduction

Rotor-bearing system is an important component in the power equipment, and its stability affects greatly the safety and normal service performance of the equipment. However, the oil film force of the bearing is usually nonlinear, and consequently causes complex and unstable dynamic behaviors of the rotor when the rotor is working under the conditions of high speed and overload. Therefore, the dynamic stability of the rotor-bearing system caused by the nonlinear oil film force has investigated by many researchers.^1–6

Yang et al.¹ established the model of rotor with transverse crack by considering parametric uncertainties. And the equation of the system was solved by the Harmonic Balance method. Then, the influences of transverse crack and uncertain parameters on the responses of the rotor were studied. Wang et al.² studied the stability of a flexible liquid-filled rotor based on Fourier series. The analytical expression of the dynamic pressure of the liquid was derived, and the dynamic stability of the rotor was analyzed by the numerical method. Desavale et al.³ established a comprehensive empirical model with dimensionless parameters. Based on the dimensional analysis approach, the diagnosis of misalignment and bearing looseness were studied. The calculating results of the empirical model agreed well with the experimental results. Barbosa et al.⁴ predicted the rotor dynamics by using the Kriging surrogate models. The models can reduce the calculation time effectively. Zhang et al.⁵ proposed an adjustable elliptical journal bearing to suppress the vibration of amplitude of the rotor, and the results indicated that the bearing can effectively reduce the vibration. Santos et al.⁶ studied the dynamic characteristic of a rigid-rotor system experimentally. The damping and natural frequencies were obtained, and the maximum servovalve voltage and radial injection pressure were found for the rotor system.

In the studies above, the integral rotor was taken as the object. However, besides the integral rotor, the rod fastening rotor is also very commonly used in the heavy gas turbine and so on. For the rod fastening rotor, the rod bolts are used to press all disks tightly together. Therefore, the rod fastening rotor will show some different dynamic behaviors from the integral rotor. In this context, the dynamics characteristics of the rod fastening rotor has been studied by some researchers. Hei et al.⁷ established the model of rod fastening rotor supported by fixed-tilting pad journal bearing. The database method was employed to solve the oil film force of the fixed-tilting pad journal bearing. The dynamics behaviors of the rotor system were investigated by the orbit diagrams, the time series, the frequency spectrum diagrams. But, the gyro effects of the disks were not taken into account when modelling. Zhang⁸ studied the tangential contact stiffness of the disks of the rod fastening rotor based on the fractal contact theory. The correctness of the fractal contact theory was verified by experiment. At the same time, it was found that the fractal contact theory has its scope of application. The study only solved the tangential contact stiffness, but the dynamics behaviors of the rod fastening rotor were not analyzed. Li⁹ established the model of rod fastening rotor by considering the contact friction force and pre-tightening force of the disks. The effects of the contact friction force and pre-tightening force on the rotor dynamics were studied. The results shown that the contact stiffness and the natural frequency of the system increased with the increase of pre-tightening force. Hu et al¹⁰ derived the motion equation of the rod fastening rotor by Lagrange equation. The influence of the initial deflection caused by the unbalanced pre-tightening force on the dynamic behaviors of the rotor was studied. Their work showed that the initial deflection has a great effect on the rotor dynamics. Meanwhile, Hu et al.¹¹ also studied the dynamic characteristics of the rub-impact rod fastening rotor. The bifurcation diagram, vibration waveform, frequency spectrum, shaft orbit, and Poincaré map were used to analyze the dynamics responses of the rotor system. However, in the works of Hu et al.¹¹, the influence of the gyro effect on the dynamic behaviors of the rod fastening rotor was not considered, and a short bearing model was assumed.

For the rod fastening rotor, the hydrodynamic sliding bearing is usually taken as the supporting part, and the oil film force of the bearing is very important for the analysis of the dynamic behaviors of the rotor. So, the investigation of the oil film force was implemented by many researches. In many studies, the model of infinite long bearing and infinite short bearing^12–14 were adopted widely. Considering journal bearing parametric uncertainties, Ramos et al.¹² studied the dynamics of the rotor system, and the rotor is supported by the short fluid film bearings. Wei et al.¹³ studied the nonlinear dynamic behaviors of the multi-disk rotor-bearing-seal system which supported by short bearing. The response of the system was calculated by the fourth order Runge–Kutta method. Liu et al.¹⁴ studied the dynamics characteristics of the rod fastening rotor which supported by the infinite long bearing. The shooting method and path-following technique were adopted to calculate the dynamics responses of the system.

The infinite long bearing and infinite short bearing models are very simple and idealized models. In practice, the slide bearings are finite long bearings. So, the two idealized models cannot describe the actual bearing. In order to solve the oil film force of finite long bearing accurately, numerical methods are proposed. Smolík et al.¹⁵ studied the shape of the threshold curve for a rigid rotor supported by journal slide bearing. In order to solve the oil film force of journal bearing, the four methods were proposed: the infinitely short approximation, the infinitely short approximation, the finite differences method, and the finite elements method. Tuckmantel¹⁶ investigated the vibration signature of rotor-coupling-bearing system by solving the oil film force of the bearing with finite volume method. Lu et al.¹⁷ calculated the oil film forces and their Jacobis by the variational constraint approach. However, the numerical methods (such as the finite volume method, finite elements method, finite differences method, and variational constraint approach) have large calculation cost.

In this study, an approximation analytical solution of oil film force is proposed for the supporting bearing of the rod fastening rotor by using the method of separation of variables and Sturm–Liouville theory. This method not only ensures the calculation accuracy, but also improves the calculation efficiency. Meanwhile, a model of the rod fastening rotor is established by considering the Gyro effect and contact effect of the disks, and an improved Newmark method is proposed to solve the dynamic responses of the rotor system. On this basis, the comparison of obit is implemented between the integral rotor and rod fastening rotor. The influence of rotating speed, eccentricity, bending stiffness, and contact stiffness on the dynamic behaviors of the rotor system is also investigated.

Equation of the rod fastening rotor system

The geometry model of the rod fastening rotor system supported by oil film hydrodynamic journal bearings is given in Figure 1. In this study, the contact effect and gyro effect of the disks are considered, and the contact effect of disks is taken as a bending spring with nonlinear stiffness.¹⁸

Figure 1.

Schematic diagram of the bearing-rod fastening rotor system.

In Figure 1, O₁ and O₂ are the centers of the two disks, l is the length of the whole rotor, a and b are the length of the two shafts, a₁ is the distance between A and O₁, b₁ is the distance between B and O₁, a₂ is the distance between A and O₂, b₂ is the distance between B and O₂, m_O₁ and m_O₂ are the mass of disk O₁ and O₂, m_A and m_B are the mass of the two shafts (the length are a and b respectively), e_O₁ and e_O₂ are the eccentricities of the disk O₁ and O₂, f_x and f_y are the nonlinear oil film forces of bearing in negative x and y directions, respectively.

Figure 2.

Coordinate of finite length journal bearing.

The dynamic equation of the rod fastening rotor system can be written as follows

M \ddot{q} + G \dot{q} + Kq = f + Q + W + f_{R}

(1)

where M is the mass matrix, G is the gyroscopic matrix, K is the stiffness matrix, f is the vector of nonlinear oil film forces of bearing, Q is periodic exciting forces (i.e., unbalanced or steam excitation with the same phase as rotating speed) vector acting on rotor, W is the vector of gravity, f_R is the vector of nonlinear restoring force caused by the nonlinear stiffness, and q is the displacement vector of the rotor. The 12 degrees-of-freedom of the system are [x_A, y_A, x_O₁, y_O₁, θ_x₁, θ_y₁, x_O₂, y_O₂, θ_x₂, θ_y₂, x_B, y_B]. x_A and y_A are the displacements of the rotor at bearing station A in x and y directions, x_O₁ and y_O₁ are the displacements of the center of disk O₁ in x and y directions, θ_x₁ and θ_y₁ are the angles that disk O₁ rotates around the x and y axes, x_O₂ and y_O₂ are the displacements of the center of disk O₂ in x and y directions, θ_x₂ and θ_y₂ are the angles that disk O₂ rotates around the x and y axes, x_B and y_B are the displacements of the rotor at bearing station B in x and y directions.

\begin{array}{l} q = {[x_{A}, y_{A}, x_{O}_{1}, y_{O}_{1}, θ_{x}_{1}, θ_{y}_{1}, x_{O}_{2}, y_{O}_{2}, θ_{x}_{2}, θ_{y}_{2}, x_{B}, y_{B}]}^{T} \\ f = {[- f_{x A}, - f_{y A}, 0, 0, 0, 0, 0, 0, 0, 0, - f_{x B}, - f_{y B}]}^{T} \end{array}

\begin{array}{l} W = {[0, m_{A} g, 0, m_{O}_{1} g, 0, 0, 0, m_{O}_{2} g, 0, 0, 0, m_{B} g]}^{T} \\ Q = {[0, 0, m_{O}_{1} e_{O}_{1} ω^{2} \sin (ω t), m_{O}_{1} e_{O}_{1} ω^{2} \cos (ω t), 0, 0, m_{O}_{2} e_{O}_{2} ω^{2} \sin (ω t), m_{O}_{2} e_{O}_{2} ω^{2} \cos (ω t), 0, 0, 0, 0]}^{T} \\ f_{R} = {[0, 0, k_{R} {(x_{O}_{2} - x_{O}_{1})}^{3}, k_{R} {(y_{O}_{2} - y_{O}_{1})}^{3}, 0, 0, - k_{R} {(x_{O}_{2} - x_{O}_{1})}^{3}, - k_{R} {(y_{O}_{2} - y_{O}_{1})}^{3}, 0, 0, 0, 0]}^{T} \end{array}

M = [\begin{matrix} m_{A} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ m_{A} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ m_{O 1} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ m_{O 1} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ J_{d 1} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ J_{d 1} & 0 & 0 & 0 & 0 & 0 & 0 \\ m_{O 2} & 0 & 0 & 0 & 0 & 0 \\ m_{O 2} & 0 & 0 & 0 & 0 \\ J_{d 2} & 0 & 0 & 0 \\ J_{d 2} & 0 & 0 \\ m_{B} & 0 \\ m_{B} \end{matrix}]

G = [\begin{matrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & - J_{z 1} ω & 0 & 0 & 0 & 0 & 0 & 0 \\ J_{z 1} ω & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & - J_{z 2} ω & 0 & 0 \\ J_{z 2} ω & 0 & 0 & 0 \\ 0 & 0 \\ 0 \end{matrix}]

K = [\begin{matrix} {\begin{matrix} (k_{11} \frac{b_{1}^{2}}{l^{2}} + \frac{k_{44}}{l^{2}} \\ - \frac{2 k_{14} b_{1}}{l^{2}} \\ + k_{11}^{'} \frac{b_{2}^{2}}{l^{2}} + \frac{k_{44}^{'}}{l^{2}} \\ - \frac{2 k_{14}^{'} b_{2}}{l^{2}}) \end{matrix}} & 0 & (- \frac{k_{11} b_{1}}{l} + \frac{k_{14}}{l}) & 0 & 0 & (\frac{k_{44}}{l} - \frac{k_{14} b_{1}}{l}) \\ 0 & {\begin{matrix} (k_{22} \frac{b_{1}^{2}}{l^{2}} + \frac{k_{33}}{l^{2}} \\ - \frac{2 k_{23} b_{1}}{l^{2}} + k_{22}^{'} \frac{b_{2}^{2}}{l^{2}} \\ + \frac{k_{33}^{'}}{l^{2}} - \frac{2 k_{23}^{'} b_{2}}{l^{2}}) \end{matrix}} & 0 & (- k_{22} \frac{b_{1}}{l} + \frac{k_{23}}{l}) & (\frac{k_{33}}{l} - \frac{k_{23} b_{1}}{l}) & 0 \\ (\frac{k_{14}}{l} - \frac{k_{11} b_{1}}{l}) & 0 & (k_{11} + k_{b}) & 0 & 0 & k_{14} \\ 0 & (\frac{k_{23}}{l} - \frac{k_{22} b_{1}}{l}) & 0 & (k_{22} + k_{b}) & k_{23} & 0 \\ 0 & (\frac{k_{33}}{l} - \frac{k_{23} b_{1}}{l}) & 0 & k_{23} & k_{33} & 0 \\ (\frac{k_{44}}{l} - \frac{k_{14} b_{1}}{l}) & 0 & k_{14} & 0 & 0 & k_{44} \\ (- \frac{k_{11}^{'} b_{2}}{l} + \frac{k_{14}^{'}}{l}) & 0 & - k_{b} & 0 & 0 & 0 \\ 0 & (- \frac{k_{22}^{'} b_{2}}{l} + \frac{k_{23}^{'}}{l}) & 0 & - k_{b} & 0 & 0 \\ 0 & (\frac{k_{33}^{'}}{l} - \frac{k_{23}^{'} b_{2}}{l}) & 0 & 0 & 0 & 0 \\ (\frac{k_{44}^{'}}{l} - \frac{k_{14}^{'} b_{2}}{l}) & 0 & 0 & 0 & 0 & 0 \\ {\begin{matrix} (\frac{k_{11} a_{1} b_{1}}{l^{2}} - \frac{k_{44}}{l^{2}} \\ + \frac{k_{14} (b_{1} - a_{1})}{l^{2}} \\ + \frac{k_{11}^{'} a_{2} b_{2}}{l^{2}} - \frac{k_{44}^{'}}{l^{2}} \\ + \frac{k_{14}^{'} (b_{2} - a_{2})}{l^{2}}) \end{matrix}} & 0 & (- \frac{k_{11} a_{1}}{l} - \frac{k_{14}}{l}) & 0 & 0 & (- \frac{k_{44}}{l} - \frac{k_{14} a_{1}}{l}) \\ 0 & {\begin{matrix} (\frac{k_{22} a_{1} b_{1}}{l^{2}} - \frac{k_{33}}{l^{2}} \\ + \frac{k_{23} (b_{1} - a_{1})}{l^{2}} \\ + \frac{k_{22}^{'} a_{2} b_{2}}{l^{2}} - \frac{k_{33}^{'}}{l^{2}} \\ + \frac{k_{23}^{'} (b_{2} - a_{2})}{l^{2}}) \end{matrix}} & 0 & (- \frac{k_{22} a_{1}}{l} - \frac{k_{23}}{l}) & (- \frac{k_{33}}{l} - \frac{k_{23} a_{1}}{l}) & 0 \end{matrix}

\begin{matrix} (- \frac{k_{11}^{'} b_{2}}{l} + \frac{k_{14}^{'}}{l}) & 0 & 0 & (\frac{k_{44}^{'}}{l} - \frac{k_{14}^{'} b_{2}}{l}) & {\begin{matrix} (\frac{k_{11} a_{1} b_{1}}{l^{2}} - \frac{k_{44}}{l^{2}} \\ + \frac{k_{14} (b_{1} - a_{1})}{l^{2}} \\ + \frac{k_{11}^{'} a_{2} b_{2}}{l^{2}} - \frac{k_{44}^{'}}{l^{2}} \\ + \frac{k_{14}^{'} (b_{2} - a_{2})}{l^{2}}) \end{matrix}} & 0 \\ 0 & (- k_{22}^{'} \frac{b_{2}}{l} + \frac{k_{23}^{'}}{l}) & (\frac{k_{33}^{'}}{l} - \frac{k_{23}^{'} b_{2}}{l}) & 0 & 0 & {\begin{matrix} (\frac{k_{22} a_{1} b_{1}}{l^{2}} - \frac{k_{33}}{l^{2}} \\ + \frac{k_{23} (b_{1} - a_{1})}{l^{2}} \\ + \frac{k_{22}^{'} a_{2} b_{2}}{l^{2}} - \frac{k_{33}^{'}}{l^{2}} \\ + \frac{k_{23}^{'} (b_{2} - a_{2})}{l^{2}}) \end{matrix}} \\ - k_{b} & 0 & 0 & 0 & (- \frac{k_{11} a_{1}}{l} - \frac{k_{14}}{l}) & 0 \\ 0 & - k_{b} & 0 & 0 & 0 & (- \frac{k_{22} a_{1}}{l} - \frac{k_{23}}{l}) \\ 0 & 0 & 0 & 0 & 0 & (- \frac{k_{33}}{l} - \frac{k_{23} a_{1}}{l}) \\ 0 & 0 & 0 & 0 & (- \frac{k_{44}}{l} - \frac{k_{14} a_{1}}{l}) & 0 \\ (k_{11}^{'} + k_{b}) & 0 & 0 & k_{14}^{'} & (- \frac{k_{11}^{'} a_{2}}{l} - \frac{k_{14}^{'}}{l}) & 0 \\ 0 & (k_{22}^{'} + k_{b}) & k_{23}^{'} & 0 & 0 & (- \frac{k_{22}^{'} a_{2}}{l} - \frac{k_{23}^{'}}{l}) \\ 0 & k_{23}^{'} & k_{33}^{'} & 0 & 0 & (- \frac{k_{33}^{'}}{l} - \frac{k_{23}^{'} a_{2}}{l}) \\ k_{14}^{'} & 0 & 0 & k_{44}^{'} & (- \frac{k_{44}^{'}}{l} - \frac{k_{14}^{'} a_{2}}{l}) & 0 \\ (- \frac{k_{11}^{'} a_{2}}{l} - \frac{k_{14}^{'}}{l}) & 0 & 0 & (- \frac{k_{44}^{'}}{l} - \frac{k_{14}^{'} a_{2}}{l}) & {\begin{cases} (\frac{k_{11} a_{1}^{2}}{l^{2}} + \frac{k_{44}}{l^{2}} \\ + \frac{2 k_{14} a_{1}}{l^{2}} + \\ \frac{k_{11}^{'} a_{2}^{2}}{l^{2}} + \frac{k_{44}^{'}}{l^{2}} \\ + \frac{2 k_{14}^{'} a_{2}}{l^{2}}) \end{cases}} & 0 \\ 0 & (- \frac{k_{22}^{'} a_{2}}{l} - \frac{k_{23}^{'}}{l}) & (- \frac{k_{33}^{'}}{l} - \frac{k_{23}^{'} a_{2}}{l}) & 0 & 0 & {\begin{matrix} (\frac{k_{22} a_{1}^{2}}{l^{2}} + \frac{k_{33}}{l^{2}} \\ + \frac{2 k_{23} a_{1}}{l^{2}} \\ + \frac{k_{22}^{'} a_{2}^{2}}{l^{2}} + \frac{k_{33}^{'}}{l^{2}} \\ + \frac{2 k_{23}^{'} a_{2}}{l^{2}}) \end{matrix}} \end{matrix}}

where f_xA and f_yA are the nonlinear oil film forces of bearing at the station A in negative x and y directions, respectively. f_xB and f_yB are the nonlinear oil film forces of bearing at the station B in negative x and y directions, respectively. J_d₁ and J_d₂ are the polar moments of the inertia for disks O₁ and O₂. J_z₁ and J_z₂ are the equatorial moments of the inertia for disks O₁ and O₂. k_b is the bending stiffness of the rod. k_R is the nonlinear stiffness of contact interface. ω is the rotating speed of the rotor. k₁₁ and k₂₂ are the bending stiffness of shaft AO in the x and y directions. k₃₃ and k₄₄ are the stiffness caused by the swing angle which disk O₁ oscillated around the x and y axes, k₁₄ and k₂₃ are the cross stiffness of the shaft AO. kʹ₁₁ and kʹ₂₂ are the bending stiffness of shaft BO in the x and y directions. kʹ₃₃ and kʹ₄₄ are the stiffness caused by the swing angle which disk O₂ oscillated around the x and y axes. kʹ₁₄ and kʹ₂₃ are the cross stiffness of the shaft BO. The rod fastening rotor is a symmetric rotor in this paper, so the relationship of the stiffness are k₁₁ = k₂₂ = kʹ₁₁ = kʹ₂₂, k₃₃ = k₄₄ = kʹ₃₃ = kʹ₄₄, k₁₄ = k₂₃ = kʹ₁₄ = kʹ₂₃.¹⁹

For the convenience of calculating and derivation, the following dimensionless variables are introduced.

\begin{array}{l} X_{i} = x_{i} / c, Y_{i} = y_{i} / c, {\bar{m}}_{i} = m_{i} / m (i = A, O_{1}, O_{2}, B), \bar{l} = l / c, \bar{ω} = ω \sqrt{c / g}, τ = ω t, \\ ε_{j} = e_{j} / c (j = O_{1}, O_{2}), σ = μ ω B r^{3} / (m g c^{2}) \\ F = f / (σ m g), K_{11} = K_{22} = k_{11} c / (m g), K_{33} = K_{44} = k_{33} / (mgc), \\ K_{23} = K_{14} = k_{23} / (m g), \bar{J} = J / (m c^{2}) \end{array}

where c is radius clearance, r is the radius of the shaft journal, B is the width of the bearing, and σ is Sommerfeld number.

Nonlinear oil film force

Solution of nonlinear oil film force

The geometry of the journal bearing is shown in Figure 2, and the cross section of the bearing is shown in Figure 3. In the two figures, O_b is the center of the bearing, O_j is the center of the shaft journal, θ is the deviation angle, Φ is the angle goes from negative y-axis direction along the clockwise direction, φ is the angle which goes from extension of O_jO_b along the clockwise direction, f_r and f_t are oil film forces in the negative radial and the tangential directions, f_x and f_y are oil film forces in negative x and y directions, h is the thickness of oil film, R is the radius of the bearing, and W is the load of the bearing.

Figure 3.

The cross section of the finite length journal bearing.

By assuming that the lubricant is in-compressible, the Reynolds equation for the lubrication of the journal bearing can be written as

\frac{1}{R^{2}} \frac{\partial}{\partial φ} (\frac{h^{3}}{6} \frac{\partial p (φ, z)}{\partial φ}) + \frac{\partial}{\partial z} (\frac{h^{3}}{6} \frac{\partial p (φ, z)}{\partial z}) = ω \frac{\partial h}{\partial φ} + 2 \frac{\partial h}{\partial t}

(2)

where μ is the viscosity of the lubricating oil film. p is the distribution of pressure of oil film, h is the thickness of oil film (

h = c + e \cos φ

), R is the radius of bearing, and z is the axial coordinate.

For the convenience of derivation and calculation, the dimensionless variables are introduced as follows

\begin{array}{l} λ = z / (\frac{B}{2}), H = 1 + ε \cos φ, ε = e / c, τ = ω t, ε^{'} = d ε / d τ, θ^{'} = d θ / d τ, \\ ψ = c / R, P = p / p_{0} (p_{0} = 2 ω μ / ψ^{2}) \end{array}

(3)

where λ is the dimensionless axial coordinate, H is the dimensionless thickness of the oil film, ε is the dimensionless eccentricity, P is the dimensionless oil film pressure distribution, ψ is clearance ratio, and

ε^{'}

is the radial velocity.

Substituting the dimensionless variables into equation (2), the dimensionless Reynolds equation can be written as

\frac{\partial}{\partial φ} (H^{3} \frac{\partial P}{\partial φ}) + {(\frac{D}{B})}^{2} \frac{\partial}{\partial λ} (H^{3} \frac{\partial P}{\partial λ}) = - 3 ε \sin φ + 6 ε^{'} \cos φ + 6 ε θ^{'} \sin φ

(4)

where

ε θ^{'}

is the dimensionless tangential velocity.

To solve equation (4), the method of separation of variables and Sturm–Liouville theory are used. Based on the method of separation of variables, the dimensionless oil film pressure P can written as

P (φ, λ) = P_{p} (φ, λ) + P_{h} (φ, λ)

(5)

where P_p(φ, λ) is the special solution, P_h(φ, λ) is the general solution, and their expressions are given as

P_{p} (φ, λ) = P_{p 1} (φ) + P_{p 2} (λ)

(6)

P_{h} (φ, λ) = P_{h 1} (φ) P_{h 2} (λ)

(7)

In equation (5), it can be seen that P_p(φ, λ) and P_h(φ, λ) determine the oil film pressure distribution P. Therefore, P_p(φ, λ) and P_h(φ, λ) need to be solved respectively, and their solutions are introduced in the next section.

Calculation of special solution

Substituting equation (6) into equation (4), the dimensionless Reynolds equation can be written as

\frac{d^{2} P_{p 2} (λ)}{d λ^{2}} = {(\frac{B}{D})}^{2} (- \frac{d^{2} P_{p 1} (φ)}{d φ^{2}} + \frac{3 ε \sin φ}{H} \frac{d P_{p 1} (φ)}{d φ} - \frac{3 ε \sin φ}{H^{3}} + \frac{6 ε^{'} \cos φ}{H^{3}} + \frac{6 ε θ^{'} \sin φ}{H^{3}})

(8)

The expression for the left item of equation (8) is a function of λ, and the expression for the right item of equation (8) is a function of φ.

Let

\frac{d^{2} P_{p 2} (λ)}{d λ^{2}} = C

(9)

{(\frac{B}{D})}^{2} (\frac{d^{2} P_{p 1} (φ)}{d φ^{2}} - \frac{3 ε \sin φ}{H} \frac{d P_{p 1} (φ)}{d φ} + \frac{3 ε \sin φ}{H^{3}} - \frac{6 ε^{'} \cos φ}{H^{3}} - \frac{6 ε θ^{'} \sin φ}{H^{3}}) = C

(10)

P_p₂(λ) can be obtained by integrating equation (9) twice.

P_{p 2} (λ) = C λ^{2} / 2 + c_{1} λ + c_{2}

(11)

In equation (11), C is arbitrary constant, c₁ and c₂ are integral constants. Let C = 0, c₁ and c₂ can be obtained by solving the boundary conditions. The boundary conditions of the special solution in the axial direction are as follows

P_{p 2} |_{λ = \pm 1} = 0

(12)

Based on equation (12), the constant c₁ and c₂ are obtained, and c₁ = c₂ = 0.

Let

\frac{d P_{p 1} (φ)}{d φ} = {\bar{P}}_{p 1} (φ)

(13)

Substituting equation (13) into equation (10), ${\bar{P}}_{p 1} (φ)$ can be written as

{\bar{P}}_{p 1} (φ) = \frac{1}{{(1 + ε \cos φ)}^{3}} (c_{3} + 6 ε^{'} \sin φ + 3 ε \cos φ - 6 ε θ^{'} \cos φ)

(14)

where c₃ is undetermined constant. By integrating equation (13), P_p₁(φ) can be expressed as

P_{p 1} (φ) = \int {\bar{P}}_{p 1} (φ) d φ + c_{4}

(15)

where c₄ is undetermined constant.

In order to obtain the expression of P_p₁(φ), Sommerfeld transformation is introduced as follows

\begin{array}{l} 1 + ε \cos φ = \frac{1 - ε^{2}}{1 - ε \cos α}, d φ = \frac{{(1 - ε^{2})}^{1 / 2}}{1 - ε \cos α} d α, \sin α = \frac{\sin φ {(1 - ε^{2})}^{1 / 2}}{1 + ε \cos φ}, \cos α = \frac{ε + \cos φ}{1 + ε \cos φ}, \\ α = 2 \arctan (\sqrt{\frac{1 - ε}{1 + ε}} \tan \frac{φ}{2}) \end{array}

(16)

Based on Sommerfeld transformation, expression of P_p₁(φ) can be written as follows

\begin{array}{l} P_{p 1} (φ) = c_{3} \frac{α - 2 ε \sin α + 0.5 ε^{2} \sin α \cos α + 0.5 ε^{2} α}{{(1 - ε^{2})}^{5 / 2}} + \frac{3 ε^{'}}{ε {(1 + ε \cos φ)}^{2}} \\ + (\frac{3 ε}{{(1 - ε^{2})}^{5 / 2}} + \frac{- 6 ε θ^{'}}{{(1 - ε^{2})}^{5 / 2}}) (ε^{2} \sin α + \sin α - \frac{\sin α \cos α}{2} ε - \frac{3 ε α}{2}) + c_{4} \end{array}

(17)

The boundary conditions of the special solution in the circumferential direction are as follows

φ = φ_{a} = 0, P_{p 1} (φ_{a}) = 0; φ = φ_{c} = π, {\frac{\partial P_{p 1}}{\partial φ} |}_{φ = φ_{c}} = 0

(18)

where φ_a is the starting position of the oil film pressure and φ_c is the termination position of the oil film.

Based on the boundary conditions (i.e., equation (18)), c₃ and c₄ can be expressed as follows

c_{4} = - \frac{3 ε^{'}}{ε {(1 + ε)}^{2}}

(19)

\frac{((\frac{6 ε θ^{'}}{{(1 - ε^{2})}^{5 / 2}} - \frac{3 ε}{{(1 - ε^{2})}^{5 / 2}}) (ε^{2} \sin α_{c} + \sin α_{c} - \frac{ε \sin α_{c} \cos α_{c}}{2} - \frac{3 ε α_{c}}{2}) + \frac{3 ε^{'}}{ε {(1 + ε)}^{2}} - \frac{3 ε^{'}}{ε {(1 + ε \cos φ_{c})}^{2}}) {(1 - ε^{2})}^{5 / 2}}{α_{c} - 2 ε \sin α_{c} + 0.5 ε^{2} \sin α_{c} \cos α_{c} + 0.5 ε^{2} α_{c}} = c_{3}

(20)

The special solutions of the oil film force in the radial and tangential directions can be obtained by integrating the P_p(φ, λ), and they can be written as follows

{\begin{matrix} F_{r}^{p} = - \int_{- 1}^{1} \int_{φ_{a}}^{φ_{c}} P_{p 1} (φ) \cos φ d φ d λ, \\ F_{t}^{p} = - \int_{- 1}^{1} \int_{φ_{a}}^{φ_{c}} P_{p 1} (φ) \sin φ d φ d λ \end{matrix}

(21)

Calculation of general solution

By substituting P_h(φ, λ) into equation (4), then the equation (4) can be written as follows

H^{3} \frac{\partial^{2} P_{h} (φ, λ)}{\partial φ^{2}} - 3 H^{2} ε \sin φ \frac{\partial P_{h} (φ, λ)}{\partial φ} + {(\frac{D}{B})}^{2} H^{3} \frac{\partial^{2} P_{h} (φ, λ)}{\partial λ^{2}} = 0

(22)

Substituting equation (7) into equation (22), the following expression can be obtained

- \frac{P_{h 2}^{″} (λ)}{P_{h 2} (λ)} = {(\frac{B}{D})}^{2} \frac{P_{h 1}^{″} (φ)}{P_{h 1} (φ)} - {(\frac{B}{D})}^{2} \frac{3}{H} ε \sin φ \frac{P_{h 1}^{'} (φ)}{P_{h 1} (φ)}

(23)

In order to solve the equation (23), the following expressions should be satisfied

- \frac{P_{h 2}^{″} (λ)}{P_{h 2} (λ)} = κ

(24)

{(\frac{B}{D})}^{2} \frac{P_{h 1}^{″} (φ)}{P_{h 1} (φ)} - {(\frac{B}{D})}^{2} \frac{3}{H} ε \sin φ \frac{P_{h 1}^{'} (φ)}{P_{h 1} (φ)} = κ

(25)

For the journal bearing, the following boundary conditions should be satisfied

P_{p} (φ, \pm 1) = P_{p 1} (φ) + P_{p 2} (\pm 1) = - P_{h} (φ, \pm 1) = - P_{h 1} (φ) P_{h 2} (\pm 1)

(26)

Based on the boundary conditions above, the axial boundary conditions for the general solution is obtained

{\begin{matrix} {P_{h 2} |}_{λ = 1} = 1, \\ {P_{h 2} |}_{λ = - 1} = 1 \end{matrix}

(27)

In order to solve equation (24), three situations of κ need to be considered (i.e., κ > 0, κ = 0, and κ < 0).

when κ = 0

Equation (24) can be written as

{P^{″}}_{h 2} (λ) = 0

(28)

By integrating the equation (28) twice, the expression of P_h₂(λ) can be obtained as follows

P_{h}_{2} (λ) = c_{5} λ + c_{6}

(29)

where c₅ and c₆ are integral constants. c₅ and c₆ can be obtained by the axial boundary conditions (i.e., equation (27)), and c₅ = 0, c₆ = 1. Substituting c₅ and c₆ into equation (29), it is found that the function of P_h₂(λ) is not related to λ. So, the situation of κ = 0 is not suitable.

2. when κ > 0

Let κ = k², the solution of equation (24) is as follows

P_{h 2} (λ) = c_{7} \cos (k λ) + c_{8} \sin (k λ)

(30)

where c₇ and c₈ are constant. There is no solution to the equation (24) combining with the boundary conditions, so the equation (30) is also not suitable.

3. when κ < 0

Let κ = −k², the solution of equation (24) can be expressed as follows

P_{h 2} (λ) = c_{9} e^{k λ} + c_{10} e^{- k λ}

(31)

where c₉ and c₁₀ are constant, they can be obtained by using the axial boundary conditions of the general solution, and c₉ = c₁₀ = e^k/(1 + e² ^k ). The equation (31) is the solution of equation (24). The final expression of P_h₂(λ) can be obtained

P_{h 2} (λ) = \frac{e^{k}}{1 + e^{2 k}} (e^{k λ} + e^{- k λ})

(32)

To solve the equation (25), the following transformations are adopted

H_{h} = c_{h} φ, P_{h 1} (φ) = f (φ) / H_{h}, H_{h}^{2} = H^{3}

(33)

where H_h is an approximation function of the thickness of the oil film, c_h is the coefficient of the approximation function, and f(φ) is the transformation function of P_h₁(φ).

Equation (25) can be expressed as Strum–Liouville equation based on equation (33), and then the equation (25) can be written as follows

f^{″} (φ) - k {(\frac{D}{B})}^{2} f (φ) = 0

(34)

The boundary conditions of the equation (34) are as follows

f (φ_{a}) = 0, f (φ_{c}) = 0

(35)

Equations (34) and (35) are taken as Strum–Liouville equation, and the equation (34) can be expressed as follows

{(p (φ) f^{'} (φ))}^{'} - q (φ) f (φ) - k {(\frac{D}{B})}^{2} ρ (φ) f (φ) = 0

(36)

where p(φ) is the coefficient function of Strum-Liouville equation, q(φ) is the potential function of Strum-Liouville equation, ρ(φ) is the weight function of Strum–Liouville equation, and p(φ) = 1, q(φ) = 0, ρ(φ) = 1.

The eigenvalue and eigenfunction of equation (36) can be obtained directly as follows

k_{i} \frac{D}{B} = \frac{i π}{(φ_{c} - φ_{a})}, i = 1, 2 \dots

(37)

f_{i} (φ) = sin (k_{i} \frac{D}{B} = (φ - φ_{a})) = sin (\frac{i π (φ - φ_{a})}{(φ_{c} - φ_{a})})

(38)

where i is the number of eigenvalue, k_i is the ith eigenvalue, and f_i(φ) is the ith eigenfunction. Because of the orthogonality of eigenfunction, the following expression can be obtained

{\bar{f}}_{i} (φ) = η_{i} f_{i} (φ) (η_{i} = \frac{1}{\int_{φ_{a}}^{φ_{c}} ρ (φ) f_{i} {(φ)}^{2} d φ})

(39)

Then, general solution can be expressed as follows

P_{h} (φ, λ) = \sum_{i = 1}^{\infty} σ_{i} f_{i} (φ) P_{h 2 i} (λ)

(40)

where

σ_{i} = η_{i} \cdot \int_{φ_{a}}^{φ_{c}} - ρ (φ) \cdot p_{p 1} (φ) \cdot f_{i} (φ) d φ P_{h 2_{i}} (λ) = \frac{e^{k}}{1 + e^{2 k_{i}}} (e^{k_{i} λ} + e^{- k_{i} λ})

The general solutions of the oil film force in the radial and tangential directions can be obtained by integrating P_h(φ, λ), and they are expressed as

{\begin{matrix} F_{r}^{h} = - \int_{- 1}^{1} \int_{φ_{a}}^{φ_{c}} P_{h} (φ, λ) \cos φ d φ d λ, \\ F_{t}^{h} = - \int_{- 1}^{1} \int_{φ_{a}}^{φ_{c}} P_{h} (φ, λ) \sin φ d φ d λ \end{matrix}

(41)

Nonlinear oil film force

According to the equations (21) and (41), the nonlinear oil film forces in the radial and tangential directions can be expressed as follows

{\begin{matrix} F_{r} = F_{r}^{P} + F_{r}^{h} \\ F_{t} = F_{t}^{P} + F_{t}^{h} \end{matrix}

(42)

The expressions of the nonlinear oil film force in the x and y directions are as follows

{\begin{matrix} F_{x} = F_{t} \cos θ + F_{r} \sin θ \\ F_{y} = - F_{t} \sin θ + F_{r} cos θ \end{matrix}

(43)

Verification of the proposed method

To verify the proposed solution method of nonlinear oil film force of the bearing. The nonlinear oil film forces of bearing calculated by the proposed method, the finite difference method,²⁰ and the infinite long bearing model are compared. In the verification, the parameter of the bearing is listed as follows: B/d = 1, ε = 0.5, xʹ = 0.1, yʹ = 0.2, and i = 30. The comparison result is shown in Figure 4. From Figure 4, it can be seen that the oil film force calculated by the proposed method is in good agreement with the result calculated by the finite difference method. Meanwhile, the calculation time of the proposed method and the finite difference method are also compared, and shown in Figure 5. From Figure 5, it can be seen that the calculation time of the proposed method is smaller than that of the finite difference method. The results of Figures 4 and 5 indicate that the proposed method not only ensures the calculation accuracy, but also improves the calculation efficiency.

Figure 4.

Nonlinear oil film forces versus eccentric angle θ. (a) F_x versus θ; (b) F_y versus θ.

Figure 5.

Comparison of calculation efficiency of two method.

Solving method of nonlinear dynamics

Improved Newmark method

The Newmark method is usually based on an assumption of constant acceleration. It is also called constant average acceleration hypothesis. The solution process of the Newmark method is as follows

At time t + Δt, the equation (1) can be written as follows

M {\ddot{q}}_{t + Δ t} + G {\dot{q}}_{t + Δ t} + K q_{t + Δ t} = P_{t + Δ t}

(44)

where P = f + Q + W + f_R .

According to the constant average acceleration hypothesis, the following equations can be obtained

{\dot{q}}_{t + Δ t} = {\dot{q}}_{t} + a_{6} {\ddot{q}}_{t} + a_{7} {\ddot{q}}_{t + Δ t}

(45)

where a₆ and a₇ are integral parameters of the Newmark method, and

a_{6} = Δ t (1 - γ)

a_{7} = γ Δ t

. Δt is the step of time, γ is the control parameter of the Newmark method, and γ ≥ 0.5.

q_{t + Δ t} = q_{t} + {\dot{q}}_{t} Δ t + [(0.5 - α) {\ddot{q}}_{t} + α {\ddot{q}}_{t + Δ t}] Δ t^{2}

(46)

where α is the control parameter of the Newmark method, and α ≥ 0.25(0.5 + γ).²

Based on constant average acceleration hypothesis, ${\ddot{q}}_{t + Δ t}$ and ${\dot{q}}_{t + Δ t}$ can be expressed by $q_{t + Δ t}$ , ${\ddot{q}}_{t}$ , ${\dot{q}}_{t}$ , and $q_{t}$ in the following forms

{\ddot{q}}_{t + Δ t} = a_{0} (q_{t + Δ t} - q_{t}) - a_{2} {\dot{q}}_{t} - a_{3} {\ddot{q}}_{t}

(47)

{\dot{q}}_{t + Δ t} = a_{1} (q_{t + Δ t} - q_{t}) - a_{4} {\dot{q}}_{t} - a_{5} {\ddot{q}}_{t}

(48)

where a₀, a₁, a₂, a₃, a₄, and a₅ are integral parameters of the Newmark method, and

a_{0} = \frac{1}{a Δ t^{2}}, a_{1} = \frac{γ}{a Δ t}, a_{2} = \frac{1}{a Δ t}, a_{3} = \frac{1}{2 a} - 1, a_{4} = \frac{γ}{a} - 1, a_{5} = (\frac{γ}{2 a} - 1) Δ t

By substituting the equations (47) and (48) into equation (44), and then the equation (44) can be written as

\tilde{K} q_{t + Δ t} = {\tilde{P}}_{t + Δ t}

(49)

where the equivalent stiffness matrix (

\tilde{K} = K + a_{0} M + a_{1} G

) is irrelevant to time, and the equivalent load vector

{\tilde{P}}_{t + Δ t} = P_{t + Δ t} + M [a_{0} q_{t} + a_{2} {\dot{q}}_{t} + a_{3} {\ddot{q}}_{t}] + G [a_{1} q_{t} + a_{4} {\dot{q}}_{t} + a_{5} {\ddot{q}}_{t}]

is relevant to time.

q_{t + Δ t}

{\dot{q}}_{t + Δ t}

{\ddot{q}}_{t + Δ t}

can be obtained by the equations (47) to (49).

q_{t + Δ t}

{\dot{q}}_{t + Δ t}

{\ddot{q}}_{t + Δ t}

satisfy the constrained condition of motion, but cannot satisfy the dynamic equation simultaneously.

Substituting the values of $q_{t + Δ t}$ , ${\dot{q}}_{t + Δ t}$ into the equation (44), ${\ddot{q}}^{'}_{t + Δ t}$ can be obtained

{\ddot{q}}^{'}_{t + Δ t} = M^{- 1} (P_{t + Δ t} - G {\dot{q}}_{t + Δ t} - K q_{t + Δ t})

(50)

It should be noted that $q_{t + Δ t}$ , ${\dot{q}}_{t + Δ t}$ , ${\ddot{q}}^{'}_{t + Δ t}$ satisfy the dynamic equation, but they do not satisfy the motion constrained condition simultaneously.

In order to made the displacement, velocity, and acceleration satisfy the motion constrained condition and dynamic equation simultaneously, the Newmark method is improved as follows.

From the equations (47) and (50), it can be seen that Newmark method will produce an imbalance acceleration ( $δ {\ddot{q}}_{t + Δ t}$ ) at each calculation step, and $δ {\ddot{q}}_{t + Δ t}$ can be obtained by equations (47) and (50).

δ {\ddot{q}}_{t + Δ t} = {\ddot{q}}_{t + Δ t} - {\ddot{q}}^{'}_{t + Δ t}

(51)

By substituting the equation (51) into equation (50), the following equation can be obtained

M {\ddot{q}}_{t + Δ t} + G {\dot{q}}_{t + Δ t} + K q_{t + Δ t} = P_{t + Δ t} + δ P_{t + Δ t}^{*}

(52)

where

δ P_{t + Δ t}^{*}

is the increment of dynamic load, and

δ P_{t + Δ t}^{*} = M δ {\ddot{q}}_{t + Δ t}

In order to eliminate the increment of dynamic load in equation (52), an increment equation of dynamics can be obtained, and it is as follows

M δ {\ddot{q}}_{t + Δ t}^{*} + G δ {\dot{q}}_{t + Δ t}^{*} + K δ q_{t + Δ t}^{*} = - δ P_{t + Δ t}^{*}

(53)

where

δ {\ddot{q}}_{t + Δ t}^{*}

δ {\dot{q}}_{t + Δ t}^{*}

, and

δ q_{t + Δ t}^{*}

are the increments of displacement, velocity, and acceleration respectively, which are caused by

- δ P_{t + Δ t}^{*}

. They can be calculated by the following equations

δ q_{t + Δ t}^{*} = - {\tilde{K}}^{-1} δ P_{t + Δ t}^{*}

(54)

δ {\dot{q}}_{t + Δ t}^{*} = a_{1} δ q_{t + Δ t}^{*}

(55)

δ {\ddot{q}}_{t + Δ t}^{*} = a_{0} δ q_{t + Δ t}^{*}

(56)

δ {\ddot{q}}_{t + Δ t}^{*}

δ {\dot{q}}_{t + Δ t}^{*}

, and

δ q_{t + Δ t}^{*}

satisfy the motion constrained condition. Because the calculate disturbance

- δ P_{t + Δ t}^{*}

is very small,

δ {\ddot{q}}_{t + Δ t}^{*}

δ {\dot{q}}_{t + Δ t}^{*}

, and

δ q_{t + Δ t}^{*}

will also satisfy the equation (53). This means that the

δ {\ddot{q}}_{t + Δ t}^{*}

δ {\dot{q}}_{t + Δ t}^{*}

, and

δ q_{t + Δ t}^{*}

satisfy the motion constrained condition and the dynamic equation simultaneously.

By adding equation (52) to equation (53), and the increment of dynamic load is eliminated. And then the displacement, velocity and acceleration at the time (t + Δt) can be written as follows

q_{t + Δ t}^{*} = q_{t + Δ t} + δ q_{t + Δ t}^{*}

(57)

{\dot{q}}_{t + Δ t}^{*} = {\dot{q}}_{t + Δ t} + δ {\dot{q}}_{t + Δ t}^{*}

(58)

{\ddot{q}}_{t + Δ t}^{*} = {\ddot{q}}_{t + Δ t} + δ {\ddot{q}}_{t + Δ t}^{*}

(59)

q_{t + Δ t}^{*}

{\dot{q}}_{t + Δ t}^{*}

{\ddot{q}}_{t + Δ t}^{*}

are taken as the true solution to equation (44).

From the analysis above, it can be seen that the solution of the improved Newmark method not only satisfies the motion constrained condition, but also satisfies dynamics equation basically.

Verification of the improved Newmark method

In order to verify the improved Newmark method, a model of the rotor bearing system²¹ is adopted. The comparison of the calculating results among Newmark method, improved Newmark method and Runge–Kutta method is implemented.

In the verification, the parameters of the rotor-bearing system are as follows: the width of the bearing B = 0.16 m, the width to diameter ratio B/d = 1, the viscosity of lubricating oil μ = 0.02626 Pa s, the radius clearance c = 0.000288 m, the parameter of system s = 0.2662. The equivalent dimensionless mass of the rotor at the two bearing stations are ${\bar{m}}_{a}$ = 0.3825 and ${\bar{m}}_{b}$ = 0.3825, the dimensionless mass of the disk are ${\bar{m}}_{O}$ = 1.235, the dimensionless stiffness of two shafts are ${\bar{k}}_{1} = {\bar{k}}_{2}$ = 4.3, and the eccentricity of the two disks is $ε_{O}$ = 0.2.

Figure 6 shows the bifurcation diagrams which calculated by the proposed method (i.e., improved Newmark method) and Newmark method. From Figure 6, it can be seen that the dynamic behaviors of the rotor can be investigated for a wide range of speeds when the improved Newmark method is adopted. Figure 7(a) shows the orbit of the rotor calculated by the improved Newmark method, Newmark method and Runge–Kutta method when $\bar{ω} = 1$ . From Figure 7(a), it can be seen that the orbit calculated by the improved Newmark method is very closed to the calculating result of the Runge–Kutta method. Because the error of each step is compensated, the improved Newmark method has higher precision than the Newmark method. The comparison of calculation efficiency between the Runge-Kutta method and the improved Newmark method has been implemented in Figure 7(b). From Figure 7(b), it can be seen that the improved Newmark method has higher computational efficiency than the Runge–Kutta method.

Figure 6.

The bifurcation diagram calculated by improved Newmark method and Newmark method. (a) Bifurcation diagram calculated by the proposed method; (b) Bifurcation diagram calculated by Newmark.

Figure 7.

The comparison of obits and calculation efficiency between Runge–Kutta method and the improved Newmark method. (a) Obits calculated by the three method for $\bar{ω} = 1$ (b) Comparison on the calculation efficiency between the Runge-Kutta method and the improved Newmark method.

Numerical examples and results

The rod fastening rotor is adopted, and the oil film force is calculated by the proposed method. The improved Newmark method is adopted to calculated the dynamic response of the rod fastening rotor-bearing system. The parameter of the rotor system are as follows: r = 0.075 m, l = 3 m, a = 1.5 m, b = 1.5 m, a₁ = 1.3 m, b₁ = 1.7 m, a₂ = 1.7 m, b₂ = 1.3 m, m_A = m_B = 206.75 kg, h = 0.4 m (the width of disks), R = 0.16 m, e_O_1 =e_O₂ = 0.3c (i.e., ε = 0.3), m_o₁ = m_o₂ = 250.93 kg, J_z₁ = J_z₂ = 3.2118 kg·m², J_d₁ = J_d₂ = 1.6059 kg·m², 2 m = m_A + m_B + m_o₁ + m_o₂, m = 457.668 kg, B = 0.15 m, c = 0.000195 m, μ = 0.02 N·s/m². The stiffness of the shafts are as follows: k₁₁ = k₂₂ = k′₁₁ = k′₂₂ = 1.0116 × 10⁷ N/m, k₁₄ = k₂₃ = −3.7733 × 10⁶ N, k′₁₄ = k′₂₃ = 3.7733 × 10⁶ N, k₃₃ = k₄₄ =k′₃₃ = k′₄₄ = 2.0847 × 10⁷ N·m, k_b = k_R = 5k₁₁/3. The number of eigenvalue i = 30.

The comparison of the dynamic behaviors is implemented between the integral rotor and rod fastening rotor. Figure 8(a) and (b) show the orbits of the rotor at bearing station for $\bar{ω} = 1$ and $\bar{ω} = 1.1$ . From Figure 8(a), it can be seen that the orbits are all periodic orbits, and the orbits amplitude of the rod fastening rotor is smaller than that of integral rotor. In Figure 8(b), it can be seen that the motion of the rod fastening rotor is periodic motion, but the motion of the integral rotor is not periodic motion. The result in Figure 8 demonstrated that the rod fastening rotor is more stable.

Figure 8.

The results of the two models for $\bar{ω} = 1, 1.1$ . (a) The results of the two models for $\bar{ω} = 1$ ; (b) The results of the two models for $\bar{ω} = 1.1$ .

Taking the dimensionless speed as control parameter

Taking the dimensionless speed $\bar{ω}$ as control parameter, and the dimensionless speed range changes from 0.9 to 1.35. Figure 9 shows the bifurcation diagram of the rotor at bearing station and disk 1 station. From Figure 9, the rotor shows rich dynamic behaviors when the dimensionless speed of the rotor changes. Specially, the rotor is stable when the rotating speed is low, and the rotor becomes unstable when the dimensionless speed of the rotor increases.

Figure 9.

Bifurcation diagram of y versus $\bar{ω}$ at bearing station and disk stations. (a) Bifurcation diagram of y at bearing station; (b) Bifurcation diagram of y at disk 1.

When $\bar{ω} = 0.9$ , Figure 10 shows the orbit, time series and spectrum of the rotor at the bearing station. From Figure 10(a), it can be seen that the orbit of the rotor is a closed curve, which is similar to a circle. It means that the motion of the rotor is periodic and the motion of the rotor is stable.

Figure 10.

The orbit, time series and spectrum of the rotor at bearing station for $\bar{ω} = 0.9$ . (a) The orbit of the journal at bearing station A; (b) Time series; (c) The spectrum.

Figure 11 shows the swing angles of the two disks. In Figure 11, it can be seen that the swing angles of the disks also change periodically. With the increase of the rotating speed, the stable periodic response of the rotor system changes to period-2 response. Figure 12 shows the orbit, Poincaré map, time series and spectrum of the rotor at bearing station when $\bar{ω} = 1.1$ . From Figure 12(a), it can be seen that the motions of the rotor are typical half-frequency whirl motions, and there are two fixed point when the orbit pass through the Poincaré section. Therefore, the motion of the rotor is period-2 motion. The swing angles of the disks are shown in Figure 13, and it can be seen that the swing of the disks is consistent with the motion of the rotor.

Figure 11.

The swing angle of disks 1 and 2 around the x and y axis for $\bar{ω} = 0.9$ . (a) The swing angle of disks 1 around the x and y axis; (b) The swing angle of disks 2 around the x and y axis.

Figure 12.

The orbit, Poincaré map, time series and spectrum of the rotor at bearing station for $\bar{ω} = 1.1$ . (a) The orbit of the journal at bearing station A; (b) Poincaré map; (c) Time series; (d) The spectrum.

Figure 13.

The swing angle of disks 1 and 2 around the x and y axis for $\bar{ω} = 1.1$ . (a) The swing angle of disks 1 around the x and y axis; (b) The swing angle of disks 2 around the x and y axis.

With the increase of the rotating speed, the period-2 motion of the rotor turns to period-4 motion. Figure 14 shows the orbit, Poincaré map, time series and spectrum of the rotor at bearing station when $\bar{ω} = 1.18$ . It can be seen that there are four closed curves in Figure 14(a), and there are four fixed points in the Poincaré section in Figure 14(b). Figure 15 shows the swing angles of disks. In Figure 15, it is demonstrated that the swing of the disks is consistent with the motion of the rotor.

Figure 14.

The orbit of rotor, Poincaré map, time series and spectrum of the rotor at bearing station for $\bar{ω} = 1.18$ . (a) The orbit of the journal at bearing station A; (b) Poincaré map; (c) Time series; (d) The spectrum.

Figure 15.

The swing angle of disks 1 and 2 around the x and y axis for $\bar{ω} = 1.18$ . (a) The swing angle of disks 1 around the x and y axis; (b) The swing angle of disks 2 around the x and y axis.

For $\bar{ω} = 1.26$ , Figure 16 shows the orbit, Poincaré map, time series and spectrum of the rotor at bearing station. From Figure 16(a), it can be seen that the orbit of the rotor is disarray, and the fixed point of the Poincaré section is also disarray in Figure 6(b). The frequency spectrum is continuous spectrum in Figure 16(d). Figure 17 shows the swing angles of the disks. From Figure 17, it can be seen that the swing angles of the disks are irregular. From the describe above, it can be seen that the chaos phenomenon is exhibited (i.e., the periodic-4 motion of the rotor turns to chaos motion) with the rotating speed increase continuously.

Figure 16.

The orbit, Poincaré map, time series and spectrum of rotor at bearing station for $\bar{ω} = 1.26$ . (a) The orbit of the journal at bearing station A; (b) Poincaré map; (c) Time series (d) The spectrum.

Figure 17.

The swing angle of disks 1 and 2 around the x and y axis for $\bar{ω} = 1.26$ . (a) The swing angle of disks 1 around the x and y axis; (b) The swing angle of disks 2 around the x and y axis.

Taking the dimensionless bending stiffness of the shafts as control parameter

The dimensionless bending stiffness of the shafts contain K₁₁, K₂₂, Kʹ₁₁, Kʹ₂₂. In this section, K₁₁ = K₂₂ = Kʹ₁₁ = Kʹ₂₂ = K. Taking the dimensionless bending stiffness K as control parameter, the dynamics responses of the rod fastening rotor system are investigated when $\bar{ω} = 1.18$ .

When the dimensionless stiffness K changes from 3.0 to 6.0, Figure 18 shows the bifurcation diagram of the rotor at bearing station and disk 1station respectively. From Figure 18, it can be seen that the motion types of the rotor system are period-4 motion, period-2 motion, periodic motion, quasi-periodic motion, period-6 motion, quasi-periodic motion, and periodic motion when the dimensionless stiffness K are 3.0, 3.5, 3.8, 4.3, 4.5, 4.7, and 6.0.

Figure 18.

Bifurcation diagram of y versus K at bearing station and disk stations. (a) Bifurcation diagram of y at bearing station; (b) Bifurcation diagram of y at disk 1 station.

When the dimensionless bending stiffness of the shaft is low, the motion of the rotor is period-4 motion. When K = 3.0, the orbit, Poincaré map, time series and spectrum of period-4 motion of the rotor are shown in Figure 14. With the increase of the dimensionless bending stiffness of the shaft, period-4 motion of the rotor turns to a period-2 motion. When K = 3.5, Figure 19 shows the orbit and Poincaré map of the period-2 motion of the rotor.

Figure 19.

The orbit, Poincaré map of rotor at bearing station for K = 3.5. (a) The orbit of the journal at bearing station A; (b) Poincaré map.

With the dimensionless bending stiffness of the shaft increase continuously, the period-2 motion bifurcates to stable periodic motion. When K = 3.8, Figure 20 shows the orbit of periodic motion of the rotor. Then, the quasi-periodic bifurcation (i.e., the periodic motion turns to quasi-periodic motion) appears with the increase of the dimensionless bending stiffness of the shaft. When K = 4.3, Figure 21 shows the orbit, Poincaré map, time series and spectrum of the rotor. From Figure 21(a), it can be seen that the quasi-periodic orbit is in a bounded region, and the projection of the Poincaré map is a closed curve which is shown in Figure 21(b). From Figure 21(c) and (d), it can be seen that the amplitude of the orbit is irregular, and frequency spectrum contains multiple frequency components. Figure 22 shows the swing angle of the disks, and the swing angle of the disks are irregular.

Figure 20.

The orbit of rotor at bearing station for K = 3.8.

Figure 21.

The orbit, Poincaré map, time series and spectrum of the rotor at bearing station for K = 4.3. (a) The orbit of the journal at bearing station A; (b) Poincaré map; (c) Time series (d) The spectrum.

Figure 22.

The swing angle of disks 1 and 2 around the x and y axis for K = 4.3. (a) The swing angle of disks 1 around the x and y axis; (b) The swing angle of disks 2 around the x and y axis.

When K = 4.5, Figure 23 shows the orbit, Poincaré map, time series and spectrum of the rotor. From Figure 23(a), it can be seen that a period-6 motion is appeared for the rotor, and the orbit contains six closed curves. Moreover, there are six fixed points in the Poincaré section, as shown in Figure 23(b). From Figure 23(c) and (d), it can be seen that there are multiple amplitudes and frequency components. Figure 24 shows the swing angle of the disks. With the increase of dimensionless bending stiffness of the shaft, the motion of the rotor changed from period-6 motion to quasi-periodic motion, and then, the quasi-periodic motion turns to stable periodic motion. When K = 4.7, Figure 25 shows the orbit, Poincaré map of the rotor. In Figure 25, it can be seen that the orbit of the rotor is in a bounded area, and projection of the Poincaré map is a closed curve. Figure 26 shows the orbit of the rotor when K = 6, and the orbit of the rotor is a closed curve.

Figure 23.

The orbit, Poincaré map, time series and spectrum of the rotor at bearing station for K = 4.5. (a) The orbit of the journal at bearing station A; (b) Poincaré map; (c) Time series; (d) The spectrum.

Figure 24.

The swing angle of disks 1 and 2 around the x and y axis for K = 4.5. (a) The swing angle of disks 1 around the x and y axis; (b) The swing angle of disks 2 around the x and y axis.

Figure 25.

The orbit, Poincaré map of rotor at bearing station for K = 4.7. (a) The orbit of the journal at bearing station A; (b) Poincaré map.

Figure 26.

The orbit of rotor at bearing station for K = 6.0.

Taking the dimensionless eccentricity of the disks as control parameter

Taking the dimensionless eccentricity of the disks as control parameter, the dynamic behaviors of the rotor system are investigated. In the investigation, the dimensionless rotating speed $\bar{ω} = 1.18$ , and the range of dimensionless eccentricity of the disks changes from 0.01 to 0.47.

Figure 27 shows the bifurcation of the rotor at the bearing and disk 1 station. From Figure 27, it can be seen that the rotor will exhibits quasi-periodic motion, periodic motion, period-2 motion, period-4 motion, period-2 motion, and periodic motion with the increase of the dimensionless eccentricity.

Figure 27.

Bifurcation diagram of y versus ε at bearing station and disk stations. (a) Bifurcation diagram of y at bearing station; (b) Bifurcation diagram of y at disk 1 station.

When the dimensionless eccentricity of the disks is very low, the motion of the rotor is quasi-periodic motion. Figure 28 shows the orbit and Poincaré map of the rotor at the bearing station for ε = 0.01. Then, quasi-periodic motion will turn to periodic motion. When ε = 0.08, Figure 29 shows the periodic orbit of the rotor at the bearing station. With the increase of the dimensionless eccentricity of the disks, the periodic motion bifurcates to a period-2 motion. When ε = 0.15, Figure 30 shows the orbit and Poincaré map of the rotor at the bearing station. The bifurcation phenomenon occurs again when the dimensionless eccentricity increases continuously, and the period-2 motion will turn to period-4 motion. When ε = 0.3, Figure 14 shows the calculating results. The period-4 motion bifurcates to period-2 motion when the dimensionless eccentricity increase continuously. Figure 31 shows the orbit and Poincaré map of the rotor at the bearing station when ε = 0.42. Then, the motion of the rotor will turn to stable periodic motion with the increase of the dimensionless eccentricity of the disks. Figure 32 shows the orbit of the rotor at bearing station when ε = 0.46.

Figure 28.

The orbit, Poincaré map of rotor at bearing station for ε = 0.01. (a) The orbit of the journal at bearing station A; (b) Poincaré map.

Figure 29.

The orbit of rotor at bearing station for ε = 0.08.

Figure 30.

The orbit, Poincaré map of rotor at bearing station for ε = 0.15. (a) The orbit of the journal at bearing station A; (b) Poincaré map.

Figure 31.

The orbit, Poincaré map of rotor at bearing station for ε = 0.42. (a) The orbit of the journal at bearing station A; (b) Poincaré map.

Figure 32.

The orbit of rotor at bearing station for ε = 0.46.

Taking the dimensionless contact stiffness of the disks as control parameter

Letting K_b=K_R, the influence of the dimensionless contact stiffness K_b and K_R on the dynamics of the rod fastening rotor bearing system is studied. The parameters are the same as the section “Taking the dimensionless speed as control parameter.” For K_b=K_R=0.1, 0.7, 2, 20, Figure33(a) to (d) show the bifurcation diagram at bearing station. From Figure 33, it can be seen that the types of the motions of the rotor system are similar when $\bar{ω} < 1.195$ , and the motions of the rotor system are periodic motion, period-2 motion, and period-4 motion. The motion of the rotor system turns to chaos motion when $\bar{ω} \geq 1.195$ . From Figure 33(a), it can be seen that the chaos motion of the rotor system turns to period-4 motion, and then, the period-4 motion turns to period-2 motion with the increase of the rotating speed. Finally, the period-2 motion turns to chaos motion. The range of rotating speed in which the motions of the rotor system are period-4 motion and period-2 motion is from 1.25 to 1.29, and this range of rotating speed shrinks gradually with the increase of the contact stiffness which are shown in Figure 33(b) and (c).

Figure 33.

Bifurcation diagram of y versus $\bar{ω}$ at bearing station and disk stations for K_b=K_R=0.1,0.7,2,20. (a) Bifurcation diagram of y for K_b = K_R = 0.1; (b) Bifurcation diagram of y for K_b = K_R = 0.7; (c) Bifurcation diagram of y for K_b = K_R = 2; (d) Bifurcation diagram of y for K_b = K_R = 20.

At the same time, it is found that the rotor collides with the wall of the bearing in the speed range of 1.21–1.245 when K_b=K_R=0.1. The speed range shrinks gradually with the increase of the dimensionless contact stiffness, and the range disappears completely when K_b=K_R=20. For K_b=K_R=0.1 and $\bar{ω} = 1.24$ , Figure 34 shows the orbit and time series of the rotor at the bearing station and swing angle of the disk 1. From Figure 34(a) and (b), it can be seen that the orbit of the rotor at the bearing station is regular, and the orbit becomes disordered with the increase of the time. The corresponding time series is shown in Figure 34(c). Figure 34(d) shows the swing angle of the disk 1. From Figure 34(d), it can be seen that the swing angle becomes large abruptly. The instability of the rotor is due to the sudden increase of the swing angle of the disk 1. Therefore, it is very important to choose the appropriated contact stiffness.

Figure 34.

The orbit, time series of the rotor at bearing station and the swing angle of disks 1 around the x and y axis for K_b=K_R=0.1, $\bar{ω} = 1.24$ . (a) The orbit of the journal for 25–0 circle; (b) The orbit of the journal for 25-last; (c) Time series; (d) The swing angle of disks 1 around the x and y.

Conclusions

In this study, to investigate the dynamic behaviors of a rod fastening rotor. A motion model of rod fastening rotor is developed by considering the contact effect and gyro effect of the disks. An approximate analytical solution of the oil film force for the supporting bearing of the rotor is proposed based on the method of separation of variables and Sturm–Liouville theory, and the Newmark method is improved to solve the motion model of the rotor bearing system. The main results are as follows:

The proposed approximate analytical solution can solve the oil film force of the bearing fast, and the results calculated by the proposed approximate analytical solution agree well with the results calculated by the finite difference method. Moreover, compared with the Newmark method, the improved Newmark method has high precision by compensating the disturbance error, and the improved Newmark method has higher calculation efficiency than Runge–Kutta method.

Based on the approximate analytical solution of the oil film force and the improved Newmark method, the dynamics behaviors of the rod fastening rotor system and the integral rotor system are studied. The results indicate that the rod fastening rotor system is more stable than the integral rotor system.

Taking the dimensionless rotating speed, bending stiffness of the shaft and the eccentricity of the disks as control parameter, the rod fastening rotor system exhibits rich and complex nonlinear dynamic behaviors. Furthermore, the influence of the contact stiffness of the disks on the dynamic motion of the rotor system is analyzed. The result indicates that the contact stiffness of the disks has great effect on the dynamic behaviors of the rotor system when the rotor operates under non-periodic motion state. Meanwhile, a large swing angle will be observed for the disks when the contact stiffness is inappropriate, and it will consequently result in collision between the rotor and bearing.

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) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Natural Scientific Research Program of Shaanxi Province of China (No. 2019JQ-928), Research Program of Weinan of Shaanxi Province of China (Nos.2016KYJ-1–2, 2018-ZDYF-JCYJ-60), Research Program of Shaanxi Railway Institute of China (No. KY2018-53).

ORCID iD

Di Hei

References

Yongfeng

Qinyu

Yanlin

, et al. Dynamic characteristics of cracked uncertain hollow-shaft. Mech Syst Signal Process 2019; 124: 36–48.

Wang

Yuan

Dynamic stability analysis for a flexible rotor filled with liquid based on three-dimensional flow. ASME J Fluids Eng 2018; 20: 2253–2267.

Desavale

RG.

Dynamics characteristics and diagnosis of a rotor-bearing’s system through a dimensional analysis approach: an experimental study.

ASME J Comput Nonlinear Dyn 2019; 14: 014501–014501.

Barbosa

MPF

Alves

WM.

Kriging-based surrogate modeling for rotordynamics prediction in rotor-bearing system[C]. In: Cavalca

Weber

(eds). Proceedings of the 10th international conference on rotor dynamics-IFToMM. Cham: Springer, 2019. doi: 10.1007/978-3-319-99270-9_22.

Zhang

, et al. Vibration suppression mechanism research of adjustable elliptical journal bearing under synchronous unbalance load. Tribol Int 2018; 132: 1–286. doi: https://doi.org/10.1016/j.triboint.2018.10.043.

Santos

Svendsen

PK.

Noninvasive parameter identification in rotordynamics via fluid film bearings-linking active lubrication and operational modal analysis.

ASME J Eng Gas Turbines Power 2017; 139: 062507.

Hei

Zhang

, et al. Nonlinear dynamic behaviors of a rod fastening rotor supported by fixed-tilting pad journal bearings. Chaos Solitons Fractals 2014; 69: 129–150.

Zhang

Experimental study on tangential contact stiffness of plate joint surface of the rod fastening rotor. DTEEES 2018; 10: 12783. dio: 10.12783/dteees/epme2018/24631.

Research on the dynamic characteristics of the rod fastening rotor. In: Proceedings of the 2018 7th international conference on energy and environmental protection. Paris: Atlantis Press, 2018. dio: 10.2991/iceep-18.2018.173.

10.

Liu

Zhao

, et al. Nonlinear dynamic behaviors of circumferential rod fastening rotor under unbalanced pre-tightening force. Arch Appl Mech 2016; 86: 1621–1631.

11.

Liu

Zhao

, et al. Nonlinear dynamic response of a rub-impact rod fastening rotor considering nonlinear contact characteristic. Arch Appl Mech 2016; 86: 1869–1886.

12.

Ramos

Ferraz

Daniel

, et al. Dynamic analysis of rotating systems considering uncertainties in the bearings’ parameters. In: Proceedings of the 10th international conference on rotor dynamics-IFToMM. Cham: Springer, 2019. doi: 10.1007/978-3-319-99272-3_32.

13.

Wei

Chen

Jiao

, et al. Computational analysis of nonlinear dynamics of a multi-disk rotor-bearing-brush seal system. Int J Str Stab Dyn 2018; 18: 350–362.

14.

Liu

Wang

, et al. Nonlinear dynamic characteristics of a three-dimensional rod-fastening rotor bearing system. Proc IMechE Part C: J Mechanical Engineering Science 2015; 229: 882–894.

15.

Smolík

Rendl

Dyk

, et al. Threshold stability curves for a nonlinear rotor-bearing system. J Sound Vibrat 2019; 442: 698–713.

16.

Wenzel da Silva Tuckmantel

Cavalca

KL.

Vibration signatures of a rotor-coupling-bearing system under angular misalignment. Mech Machine Theor 2019; 133: 559–583.

17.

Zhang

, et al. Dynamic behaviours of the rotor non-linear system with fixed-tilting-pad journal bearings support. IMechE Part J: J Engineering Tribology. 2010; 224: 1037–1047.

18.

Qian

Cheng

Chen

, et al. Analysis on bistable response of a disk-rod-fastening rotor. J Aerospace Power 2011; 26: 1563–1568. (In Chinese).

19.

Zhong

, et al. Dynamics of the rotor. Beijing: Tsinghua University Press, 1987. (in Chinese).

20.

Feng

The dynamic characteristics of tilting-pad journal bearing. Xi’an: Xi’an University of Technology, 2011. (in Chinese).

21.

Hei

Zheng

Dynamics analysis of flexible rotor system based on Wilson-θ method. China Mech Eng 2011; 22: 2918–2922. (in Chinese).