Research on the nonlinear dynamic characteristics of opposed high-speed gas bearing systems

Abstract

In recent years, the opposed high-speed gas bearing system has been gradually valued and used in the field of precision machinery, especially for precision instruments and mechanisms requiring high speed, high precision, and high rigidity. Although the bearing capacity is not as good as the oil film bearings, it can provide a working environment where the rotor can generate high speed and low heat without deformation of the shaft, and the gas pressure distribution of clearance in bearing also has better stability. Due to the strong nonlinearity of the gas film pressure function of gas bearings and the fact that the actual shaft system possesses dynamic problems including critical speed, spindle imbalance or improper bearing design, it will cause the rotation process of the shaft to produce a nonperiodic motion and instability, and even chaotic motion under certain parameters. And these irregular movements can even cause machine damage or process delays when serious, so in order to understand the process of working under the conditions where the system will have a nonperiodic phenomenon and to avoid the occurrence of irregular vibration especially chaos. In this paper, the opposed high-speed gas bearing system feature will be discussed in detail with three different numerical analysis methods, i.e. the finite difference method, perturbation method, and mixing method. The relevant theories include dynamic trajectories, spectrum analysis, bifurcation diagram, Poincare map, and the maximum Lyapunov exponents. From the results of nonlinear dynamic behavior of the rotor center, periodic and nonperiodic motions occur at different rotor masses and bearing parameters, respectively. Especially, for the chaos of shaft exists at specific intervals and can be distinguished efficiently. Meanwhile, it is found to ensure that the bearing system can suppress the phenomena of chaos actively by adjusting the bearing parameters, and reduce the system loss caused by irregular vibration. It is expected to be an important basis for designing a precision shaft or mechanism and to enhance the stability and performance of bearing system.

Keywords

Opposed high-speed gas bearing nonperiodic chaos

Introduction

Nowadays, ultra-precision machining and measuring technologies for precision instruments and equipment components have been further developed to the nanometer level. In addition to the substantial increase in measurement accuracy requirements, the need has arisen for nanometer-scale displacement and positioning, such as precision processing machines and semiconductor wafer cutting machines, high-speed hard drives, etc. Besides, a system in need of clean, ultra-high speed, accurate positioning also has a large number of applications.^1,2 In order to meet the relevant requirements, such as requiring the rotation accuracy of the gauges to reach the nanometer level, and requiring the maintenance of a high degree of stability, the application of the opposed high-speed gas bearing system among the key components is gaining more and more attention. For example, the high accuracy supporting techniques for gyro motor becomes increasingly using the opposed high-speed gas bearing to replace the mechanical bearing, and performs higher stiffness, longer life and higher precision. Gyroscopes applied opposed high-speed lubrication bearing of gyro motor is currently the highest precision gyroscope. Especially for the wedge-shaped type of opposed high-speed gas bearing, it could work in high speed, high structure stability, strong anti-overload ability, good and centroid stability, etc.³ In the research on gas bearings, how to accurately calculate the gas film pressure distribution has always been an important issue. Nguyen et al.^4–6 has proposed to solve the gas bearing pressure distribution problem of complex boundaries, using the finite element method (FEM) after 1991. In contrast, the opposed bearings studied in this paper have grooves on their surface, and the related research in this area is quite rare. In 1958, Whipple⁷ was the first to study the characteristics of gas bearings for herringbone grooves. Vohr⁸ studied gas bearings with helical grooves in 1963 and deduced the differential equation of the gas film pressure on the assumption that the number of grooves is close to infinity. This theory is called “Narrow Groove Theory”. Until 1991, Absi and Bonneau⁹ used the FEM to solve the pressure distribution function of herringbone bearings and studied their related characteristics, including stability analysis, pressure distribution, and the influence of key parameters on the system. It can be seen from the above literature that the gas film dynamic characteristics of gas bearings have a great influence on the entire bearing support system, especially the pressure of the gas film, the rigidity of the gas film, and the damping effect. In the literature of the rotor dynamic behavior, in 1974 it was found by Bently¹⁰ through experimentation that a rotor bearing system has subharmonic vibration of order 2 and order 3. Choi and Noah¹¹ then used the analytic method to prove the existence of subharmonic vibration in rotor bearing systems.

Recently, many electromechanical and opto-electric devices are downsizing and printed-circuit boards used in these devices have also been made smaller and smaller. So, the diameters of the small holes for fixing electric parts to the printed-circuit board had to be reduced. Opposed conical bearings with spiral grooves were designed for a high-speed spindle to drill these small holes in printed-circuit boards. Compared with oil bearings, opposed conical bearing system has much lower power consumption at high speeds and without causing environmental problems. Moreover, opposed conical gas bearing has a larger load capacity than aerostatic bearings, as a higher supply pressure can be used to suppress the rotor instability. About the research of this gas bearing, in 1972, Kazimierski¹² studied the conical (wedge-shaped) gas bearing system and considered the rotor dynamic behaviors mounted in two hybrid conical gas lubricated bearings. It was found an optimum angle of conical gas bearings gives the maximum stability of the motion of the rotor with respect to the self-excited, cylindrical vibrations is important and can be obtained. Meanwhile, the existence of an optimum stiffness of elastic elements of support with respect to the creation of self-excited whirl vibrations of the rotor and foundation is proved. In 2002, Yoshimoto et al.¹³ applied the conical bearings with spiral grooves for high-speed spindles to investigate the performance of bearing system theoretically and experimentally. They first presented these bearing types for pressurized water fed to the inside of the rotating shaft and then introduced into spiral grooves through holes located at one end of each spiral groove. The stability of the proposed bearing was also theoretically predicted by the perturbation method, and the numerical results were compared with experimental results. In 1993, Agrawal¹⁴ studied the power-shaped multipad hydrodynamic conical gas bearings working in a compressible fluid that is suitable for high speed turbomachines. The results revealed that the bearing has higher load capacity both in the radial and axial directions than traditional aerostatic bearings. Especially, the radial load is three times that of an axial groove bearing under high rotational speed. This conical gas bearing is more stable and the friction factor is about one-third that of an axial groove bearing. It is concluded that the radial load is approximately a function of compressibility number and eccentricity only, and the axial load is approximately a function of compressibility number only. In 2007, Pan et al.¹⁵ presented the five degrees-of-freedom characterization of the stability of a conical gas bearing with the spiral-groove design for bearing numbers up to 500. The stability thresholds of bearing system with respect to each mode including axial, cylindrical, and conical modes are analyzed as functions of the bearing number and rotor mass, respectively.

Wang^16,17 used a mixing method to solve the gas film pressure of high-speed rotor gas bearings and grooved gas shaft systems. At the same time, the dynamic action of the rotor and journal centers was analyzed by the trajectory, spectrogram, and bifurcation diagrams. The study shows that the rotor center and journal center have periodic, quasi-periodic, or sub-harmonic motion under different operating conditions. Subsequent applicants successively analyzed other gas shaft systems and found that the system may even occur nonlinear chaotic motion under certain operating ranges.^18,19

In addition, in the related experiments, Young²⁰ used ultrasonic technology to detect abnormal rotation precision and spindle vibration instability produced by special operating conditions in 1992, but it was still in the experimental stage and lacked theoretical analysis and verification. There are still many scholars in the follow-up research, but most of them have focused on the discovery and experimental analysis of the abnormal vibration behavior of the oil film bearing system. There have only been a few developments on the experimental planning and testing platform technology for gas bearing systems. In 2002, Nikolaou²¹ designed a test platform and performed spectrum analysis on the vibration and noise signals of the spindles supported by the ball bearings, and attempted to compare them with the theoretical values, but the effect was limited. In the Taguchi Experimental Planning Section, Wu and Kung²² used the Taguchi method in 2005 for an experimental hybrid oil film bearing design test factor and level, containing the oil film viscosity, the amount of oil and orifice size, etc. to analyze normal and abnormal behavior of the rotor bearing operation, but it could not confirm whether this state is chaos phenomenon. In 2008, Al-Raheem²³ carried out anomaly diagnosis analysis of ball bearings through theoretical analysis, but it was still at the theoretical stage. However, the related research and development of the above-mentioned literature for opposed high-speed gas bearing systems is quite scarce and lacking in theoretical integration and verification. Therefore, integrating theoretical analysis technology will become an important direction for academic research and industrial production in the future, but it is still needed to do sustained research and development on the related theories. Today, with emphasis on speed and efficiency in production, to discuss shafting of high-precision and high stability, which is easy to detect has become an important research direction of contemporary industry and academia, and is also the main focus of this paper.

Research methods

This paper researches on the opposed high-speed gas bearing (OHGB) systems required for the precision machinery industry, especially suitable for ultra-high-speed areas and environments in need of a high load-carrying capacity. At present, the opposed-type bearing is a traditional H-type bearing structure, as shown in Figure 1. This type of bearing is a combination of plain thrust bearings and cylindrical journal bearings, sometimes referred to as “closed linear bearings”, which have matured to date and which have grooves in the bearing surfaces that produce a larger bearing ability. However, if the application has a tapered bearing surface or a spherical surface shaft, it is not suitable. For this purpose the present paper mainly studies two shafts, which are different from conventional bearings, comprising an opposed wedge-shaped gas bearing and an opposed dome gas bearing, in order to summarize the optimization design criteria for the opposed-type high-speed gas bearings.

Figure 1.

H-type bearing structure: (a) spiral groove type; (b) herringbone spiral groove type.

Governing equations of opposed wedge-shaped gas bearings

Opposed wedge-shaped gas bearings have the ability to support radial and axial loads, and can be used in the support of contact surfaces with special needs. There are two main forms of application: opposed obverse wedge-shaped and reverse wedge-shaped bearings, with the structure shown in Figure 2.

Figure 2.

Opposed wedge-shaped bearing structure: (a)obverse wedge shape; (b) reverse wedge shape.

Since the bearing system is mainly gas-fluid-supported, the following basic assumptions have to be fulfilled:

The flow of the gas film is in an isothermal, laminar flow state.

The gas viscosity coefficient does not change with the pressure and temperature, as a constant.

Therefore, the compressible Reynolds equation of the gas film pressure is as follows

\frac{\partial}{\partial θ} (\bar{P} {\bar{h}}^{3} \frac{\partial \bar{P}}{\partial θ}) + \frac{\partial}{\partial η} (\bar{P} {\bar{h}}^{3} \frac{\partial \bar{P}}{\partial η}) = Λ \frac{\partial (\bar{P} \bar{h})}{\partial θ} + 2 Λ \frac{\partial}{\partial τ} (\bar{P} \bar{h})

(1)

where

\bar{P}

is the dimensionless gas film pressure,

\bar{h}

is the dimensionless gas film thickness,

θ

and

η

are the coordinate systems,

Λ

is the bearing number relevant to speed, and ω is the rotational speed. Figure 2 shows such a bearing groove width, location, number, groove depth and slot angle are all different. The spiral groove parameters are groove angle α, groove depth δ, groove width b₁, table width b₂.

{\bar{G}}_{1}

and

{\bar{R}}_{g}

are the geometric parameters as shown in the following equation

{\bar{G}}_{1} = \frac{b_{1}}{b_{1} + b_{2}}, {\bar{R}}_{g} = \frac{R_{o} - R_{g}}{R_{o} - R_{j}}

(2)

This type of bearing in the wedge-shaped body parameters includes a half cone angle of $γ$ , the inner and outer radii R_j, R_o of the tank body, the axial length of the cone face and the axial length of the tank region L₁, L₂, and the radius R_g of the end of the spiral groove.

Opposed wedge-shaped bearings analyzed in this study are often used industrially. For realistic operating conditions, thereby simulating and analyzing the dynamic performance of the gas film pressure and the rotor, the groove angle is set to 161.7º–162.7º.

The thickness of the gas film can be expressed in terms of $θ$ in the groove and the ridge regions, respectively.

Groove region: h = C_{r} + C_{g} + ε_{x} \cos θ + ε_{y} \sin θ

(3)

Ridge region: h = C_{r} + ε_{x} \cos θ + ε_{y} \sin θ

(4)

Solution of instantaneous gas film supporting force

In Figure 1, the pressure distribution, $\bar{P} (θ, η)$ , of the gas film in the gap between the shaft and the bushing enables the estimation of the components of the force acting between the shaft and the bushing to be given as the integrals

f_{gx} = P_{a} R^{2} \int_{0}^{L / R} \int_{0}^{2 π} \bar{P} (θ, η) \cos θ d η d θ

(5)

f_{gy} = P_{a} R^{2} \int_{0}^{L / R} \int_{0}^{2 π} \bar{P} (θ, η) \sin θ d η d θ

(6)

Rotor dynamic equations and computation procedure

The analysis presented in this paper considers a rotor-bearing system comprising a perfectly balanced, flexible rotor of mass m_r supported symmetrically on opposed wedge-shaped bearings, which are mounted on rigid pedestals. With these idealizations, it is possible to consider one bearing in isolation. Hence, the system under investigation consists of a rotating rotor of mass m_r with two degrees of translatory oscillation in the transverse plane. The idealized rotor-bearing model is shown in Figure 3, in which O₁ is the geometric bearing center, O₂ is the geometric rotor center, O₃ is the geometric journal center, $ρ$ is the mass eccentricity of the rotor, K_p is the shaft stiffness, and $ω$ is the rotational speed of the shaft.

Figure 3.

Model of flexible rotor supported by two gas bearings.

In the transient state, the equations of motion of O₂ (i.e. the rotor center) can be written in Cartesian coordinates as

m_{r} \frac{d^{2} x_{2}}{d t^{2}} + K_{p} (x_{2} - x_{3}) = m_{r} ρ ω^{2} \cos ω t

(7)

m_{r} \frac{d^{2} y_{2}}{d t^{2}} + K_{p} (y_{2} - y_{3}) = m_{r} ρ ω^{2} \sin ω t

(8)

Regarding the journal center, O₃, the resulting forces in the horizontal and vertical directions are as shown in Figure 3. From the equilibrium of forces, the forces applied to the journal center, O₃, may be expressed as

f_{gx} = K_{p} (x_{2} - x_{3}) / 2

(9)

f_{gy} = K_{p} (y_{2} - y_{3}) / 2

(10)

The following transformation is introduced

X_{2} = x_{2} / C_{r}; X_{3} = x_{3} / C_{r}; Y_{2} = y_{2} / C_{r}; Y_{3} = y_{3} / C_{r}; τ = ω t; \frac{d}{dt} = ω \frac{d}{d τ}

(11)

Additionally, the following nondimensional groups are defined

F_{gx} = \frac{f_{gx}}{P_{a} R L} F_{gy} = \frac{f_{gy}}{P_{a} R L} M_{r} = \frac{m_{r} ω^{2}}{K_{p}} ξ = \frac{m_{r} ρ ω^{2}}{K_{p} C_{r}}

Substituting these nondimensional groups into equations (7) to (10) yields

F_{gx} = \frac{K_{p} C_{r}}{2 P_{a} R L} (X_{2} - X_{3})

(12)

F_{gy} = \frac{K_{p} C_{r}}{2 P_{a} R L} (Y_{2} - Y_{3})

(13)

M_{r} \frac{d^{2} X_{2}}{d τ^{2}} + (X_{2} - X_{3}) = ξ \cos τ

(14)

M_{r} \frac{d^{2} Y_{2}}{d τ^{2}} + (Y_{2} - Y_{3}) = ξ \sin τ

(15)

Mathematical formulation of numerical simulation

In the model of the bearing, the Reynolds’ equation is a nonlinear partial differential equation, so a direct numerical method is used for its solution. To solve the bearing gas pressure distribution the perturbation method,¹² traditional finite difference method (FDM),¹⁵ and a mixing method^16,17 were used. For traditional FDM, let us first discretize equation (1) using the central-difference scheme in the $θ$ and η directions, and the implicity-backward-difference scheme in time $τ$ . For simplicity, a uniform mesh size is used. Equation (1) can be transformed into the following form

\begin{array}{l} 3 {({\bar{h}}_{i, j}^{n + 1})}^{2} (\frac{{\bar{h}}_{i + 1, j}^{n + 1} - {\bar{h}}_{i - 1, j}^{n + 1}}{2 Δ θ}) \cdot (\frac{S_{i + 1, j}^{n + 1} - S_{i - 1, j}^{n + 1}}{2 Δ θ}) + {({\bar{h}}_{i, j}^{n + 1})}^{3} (\frac{S_{i + 1, j}^{n + 1} - 2 S_{i, j}^{n + 1} + S_{i - 1, j}^{n + 1}}{{(Δ θ)}^{2}}) \\ + 3 {({\bar{h}}_{i, j}^{n + 1})}^{2} (\frac{{\bar{h}}_{i, j + 1}^{n + 1} - {\bar{h}}_{i, j - 1}^{n + 1}}{2 Δ η}) \cdot (\frac{S_{i, j + 1}^{n + 1} - S_{i, j - 1}^{n + 1}}{2 Δ η}) + {({\bar{h}}_{i, j}^{n + 1})}^{3} (\frac{S_{i, j + 1}^{n + 1} - 2 S_{i, j}^{n + 1} + S_{i, j - 1}^{n + 1}}{{(Δ η)}^{2}}) \\ = 2 Λ \cdot {\bar{h}}_{i, j}^{n + 1} (\frac{{\bar{P}}_{i + 1, j}^{n + 1} - {\bar{P}}_{i - 1, j}^{n + 1}}{2 Δ θ}) + 2 Λ {\bar{P}}_{i, j}^{n + 1} (\frac{{\bar{h}}_{i + 1, j}^{n + 1} - {\bar{h}}_{i - 1, j}^{n + 1}}{2 Δ θ}) \\ + 4 Λ \cdot {\bar{h}}_{i, j}^{n + 1} (\frac{{\bar{P}}_{i, j}^{n + 1} - {\bar{P}}_{i, j}^{n}}{Δ τ}) + 4 Λ \cdot {\bar{P}}_{i, j}^{n + 1} (\frac{{\bar{h}}_{i, j}^{n + 1} - {\bar{h}}_{i, j}^{n}}{Δ τ}) \end{array}

(16)

where S=P².

Then, the pressure distribution at each time step is obtained using an iterative calculation process from the above traditional numerical method.

In this paper, the mixing method is used and it is a combination of the differential transformation method (DTM) and the FDM. Differential transformation¹⁹ is applied to deal with equation (1) and this method is also one of the most widely used techniques for solving differential equations due to its rapid convergence rate and minimal calculation error. A further advantage of this method over the integral transformation approach is its ability to solve nonlinear differential equations. The theoretical analysis is introduced in the following.

For the nonlinear Reynolds equation time domain, the DTM was first used to discretize the time into intervals and the central difference scheme of the FDM was then applied to the coordinates of the locations. The basic principle of the DTM is as follows.

The differential transformation of the function of f(t) is defined as

\bar{f} (k) = \frac{T^{k}}{k!} {[\frac{\partial^{k} f (t)}{\partial t^{k}}]}_{t = 0}

(17)

wherein

\bar{f} (k)

refers to the transformation function in the transformation domain,

f (t)

refers to the original function in the time domain, T refers to the time separation, and k refers to the conversion parameter. The reverse equation of

\bar{f} (k)

f (t) = \sum_{k = 0}^{\infty} \bar{f} (k) {(\frac{t}{T})}^{k}

(18)

Calibration is needed to ensure precision during the calculation, therefore, when the second time interval is used, equations (19) and (20) are used to perform calibration

f_{i - 1} (0) = f_{i} (0) + {\sum_{k = 1}^{\infty} {\bar{f}}_{i} (k) (\frac{- T_{i - 1}}{T_{i}})}^{k}

(19)

If $T_{0} = T_{1} = T_{2} = \dots$ then the above equation becomes

f_{i - 1} (0) = f_{i} (0) + \sum_{k = 1}^{\infty} {\bar{f}}_{i} (k) {(- 1)}^{k}

(20)

If the calculated values on each side of the equals’ signs in equations (19) and (20) vary greatly, precision is insufficient. Therefore, the time interval T needs to be reduced and re-calculation is necessary.

The basic arithmetic operation for the differential transformation is as shown in Table 1.

Table 1.

Basic arithmetic formula for the differential transformation.

Before transformation	After transformation
$af (t) \pm bg (t)$	$a \bar{f} (k) \pm b \bar{g} (k)$
$f (t) g (t)$	$\sum_{l = 0}^{k} \bar{f} (k - l) \bar{g} (l) = \sum_{l = 0}^{k} \bar{f} (l) \bar{g} (k - l)$
$f (t) / g (t)$	$\frac{\bar{f} (k) - \sum_{l = 0}^{k - 1} (\bar{f} (l) / \bar{g} (l)) \bar{g} (k - l)}{\bar{g} (0)} = \frac{\bar{f} (k) - \sum_{l = 0}^{k - 1} (\bar{f} (k - l) / \bar{g} (k - l)) \bar{g} (l)}{\bar{g} (0)}$
$\cos (at + b)$ , $a \neq 0$	$\frac{{(aT)}^{k}}{k!} \cos (\frac{π k}{2} + b)$
$\frac{d^{n} f (t)}{d t^{n}}$	$\frac{(k + n)!}{k! T^{n}} \bar{f} (k + n)$
$f^{m} (t)$	$\sum_{l = 0}^{k} \bar{f} (k - l) \sum_{n = 0}^{l} \bar{f} m - 2 (n) \bar{f} (l - n)$

The differential transformation of the governing equation and the boundary condition are obtained and the pressure function theorem domain are divided into n number of sub-intervals, the pressure solution pf for each sub-internal is expressed as an inverse transformation equation. If the theorem domains of the sub-intervals are T₀, T₁, T₂, T₃, …, then the expression for the function in the first interval, based on equation (18), will be

P_{o} (t) = \sum_{k = 0}^{\infty} {\bar{P}}_{o} (k) {(\frac{t}{T_{o}})}^{k}, 0 \leq t \leq T_{o}

(21)

If the initial value of $P_{o} (0)$ is known, then from the above equation, $P_{o} (0) = {\bar{P}}_{o} (0)$ can be obtained. From ${\bar{P}}_{o} (0)$ , other discrete values of ${\bar{P}}_{o} (k), k = 1, 2, \dots$ in this interval can be obtained. In the first sub-interval, the destination value $P_{o} (t = T_{o})$ of the function of $P_{o} (t)$ refers to the initial value $P_{1} (t = 0)$ of the function $P_{1} (t)$ in the second sub-interval. The function $P_{1} (t)$ in the second sub-interval is

P_{1} (t) = \sum_{k = 0}^{\infty} {\bar{P}}_{1} (k) {(\frac{t}{T_{1}})}^{k}, 0 \leq t \leq T_{1}

(22)

From ${\bar{P}}_{1} (0)$ , other ${\bar{P}}_{1} (k), k = 1, 2, \dots$ can also be obtained. In such order, the value of ${\bar{P}}_{i} (k)$ can be obtained. Then, the inverse differential transformation can be further used to obtain the solution at each time domain, and the equation for the series solution is as follows

P_{i} (t) = \sum_{k = 0}^{m} {\bar{P}}_{i} (k) {(\frac{t}{T_{i}})}^{k}, 0 \leq t \leq T_{i}

(23)

For solving the Reynolds equation for the current gas bearing system, the differential transformation method is used for taking transformation with respect to the time domain τ, and hence equation (1) becomes

\begin{array}{l} 3 I \otimes \frac{\partial \bar{h}}{\partial θ} \otimes \frac{\partial ξ}{\partial θ} + J \otimes \frac{\partial^{2} ξ}{\partial θ^{2}} + 3 I \otimes \frac{\partial \bar{h}}{\partial η} \otimes \frac{\partial ξ}{\partial η} + J \otimes \frac{\partial^{2} ξ}{\partial η^{2}} \\ = 2 Λ \bar{h} \otimes \frac{\partial \bar{P}}{\partial θ} + 2 Λ \bar{P} \otimes \frac{\partial \bar{h}}{\partial θ} + 4 Λ \bar{h} \otimes \frac{\partial \bar{P}}{\partial τ} + 4 Λ \bar{P} \otimes \frac{\partial \bar{h}}{\partial τ} \end{array}

(24)

where

ξ (k) = {\bar{P}}^{2} = \bar{P} \otimes \bar{P} = \sum_{l = 0}^{k} {\bar{P}}_{i, j} (k - l) {\bar{P}}_{i, j} (l)

(25)

I (k) = {\bar{h}}^{2} = \bar{h} \otimes \bar{h} = \sum_{l = 0}^{k} {\bar{h}}_{i, j} (k - l) {\bar{h}}_{i, j} (l)

(26)

J (k) = {\bar{h}}^{3} = \bar{h} \otimes \bar{h} \otimes \bar{h} = \sum_{l = 0}^{k} {\bar{h}}_{i, j} (k - l) \sum_{m = 0}^{l} {\bar{h}}_{i, j} (l - m) {\bar{h}}_{i, j} (m)

(27)

The FDM is then used to discretize equation (24) with respect to the θ and $η$ directions. Note that equation (24) is discretized using the second-order accurate central difference scheme for both the first and the second derivatives.

Substituting equations (25) to (27) into equation (24) yields

\begin{array}{l} 3 \sum_{l = 0}^{k} I_{i, j} (k - l) \sum_{m = 0}^{l} (\frac{{\bar{h}}_{i + 1, j} (l - m) - {\bar{h}}_{i - 1, j} (l - m)}{2 Δ θ}) \cdot (\frac{ξ_{i + 1, j} (m) - ξ_{i - 1, j} (m)}{2 Δ θ}) \\ + \sum_{l = 0}^{k} J_{i, j} (k - l) (\frac{ξ_{i + 1, j} (l) - 2 ξ_{i, j} (l) + ξ_{i - 1, j} (l)}{{(Δ θ)}^{2}}) \\ 3 \sum_{l = 0}^{k} I_{i, j} (k - l) \sum_{m = 0}^{l} (\frac{{\bar{h}}_{i, j + 1} (l - m) - {\bar{h}}_{i, j - 1} (l - m)}{2 Δ η}) \cdot (\frac{ξ_{i, j + 1} (m) - ξ_{i, j - 1} (m)}{2 Δ η}) \\ + \sum_{l = 0}^{k} J_{i, j} (k - l) (\frac{ξ_{i, j + 1} (l) - 2 ξ_{i, j} (l) + ξ_{i, j - 1} (l)}{{(Δ η)}^{2}}) \\ = 2 Λ \cdot \sum_{l = 0}^{k} [{\bar{h}}_{i, j} (k - l) (\frac{{\bar{P}}_{i + 1, j} (l) - {\bar{P}}_{i - 1, j} (l)}{2 Δ θ}) \cdot] + 2 Λ \sum_{l = 0}^{k} [{\bar{P}}_{i, j} (k - l) (\frac{{\bar{h}}_{i + 1, j} (l) - {\bar{h}}_{i - 1, j} (l)}{2 Δ θ}) \cdot] \\ + 4 Λ \sum_{l = 0}^{k} [(\frac{l + 1}{\tilde{H}}) {\bar{P}}_{i, j} (k - l) \cdot {\bar{h}}_{i, j} (l + 1)] + 4 Λ \sum_{l = 0}^{k} [(\frac{l + 1}{\tilde{H}}) {\bar{h}}_{i, j} (k - l) \cdot {\bar{P}}_{i, j} (l + 1)] \end{array}

(28)

where

\tilde{H}

is the time step. From equation (28),

{\bar{P}}_{i, j} (k)

is obtained for each time interval, where i and j indicate the node position and k indicates the kth term.

Computing the motions of the rotor center is an iterative procedure, which first determines the acceleration, then the velocity, and finally the displacement, step-by-step over time. The computation procedure begins by specifying an initial static equilibrium state. The initial displacement of the rotor ( $X_{o}$ , $Y_{o}$ ) corresponds to the static equilibrium position and defines the gap ${\bar{h}}_{i, j} (k)$ between the shaft and the journal bearing. The initial velocity of the rotor is assumed to be zero.

The iterative computation procedure can be summarized as follows:

Step 1: Following a time increment $Δ τ$ , the new values of the rotor acceleration, velocity, and displacement are calculated to obtain.

Step 2: The displacements of the rotor center obtained from Step 1 can then be determined and the corresponding change in the value of the gap ( $\bar{h}$ ) can be calculated. Substituting the new value of $\bar{h}$ into equation (28) gives the new pressure distribution in the gap between the shaft and the journal.

Step 3: The pressure distribution obtained from Step 2 is integrated to estimate the internal force.

Step 4: The displacement and velocity values computed in Step 1, the pressure distribution calculated in Step 2, and the internal force obtained in Step 3 are taken as the new initial conditions. Using this new set of conditions, the calculation procedure returns to Step 1 to compute the changes in the gas bearing system during the time interval $Δ τ \to 2 Δ τ$ .

These three steps allow determination of the pressure distribution, the relationship between the gas film thickness and change of position, variation of the gas film bearing force, and the dynamic orbits of the rotor, etc. at each time point. In addition, the maximum Lyapunov exponents can be used as a primary determination factor to predict periodic, quasi-periodic, aperiodic, and chaotic motion caused by such factors as the bearing number, mass of the rotor, and so on. By using the bearing number and the rotor mass of the air bearing as bifurcation parameters, a full view bifurcation diagram can be constructed. An analysis of stability can be conducted and the relationship between the maximum Lyapunov exponents and each of the previously mentioned factors can be interpreted. The parameters of the stable and the nonstable region of the system can be determined.

To further verify the accuracy of the numerical analysis, two different numerical methods were used (as previously mentioned) to ensure that the mixing method proposed in this paper was applicable and to increase the calculation precision.

Rotor dynamic analysis

In this study, the time series data of the first 1000 revolutions are excluded from the dynamic behavior investigation so as to ensure that the data used represent steady-state conditions. The resulting data include the orbital paths and velocity of the rotor center O₂ and the orbits of the journal center O₃. These data are used to generate the power spectra, Poincaré maps, and bifurcation diagrams.

Poincaré map

A Poincaré section (map) is a hypersurface in the state space that is transverse to the flow of a given dynamic system. In nonautonomous systems, a point on the Poincaré section is referred to as the return point of the time series at a constant interval of T, where T is the driving period of the exciting force. The projection of a Poincaré section on the X(nT)-Y(nT) plane is referred to as the Poincaré map of the dynamic system.

Power spectrum

The power spectra for the vertical and horizontal displacement of the journal center O₃, and the rotor center O₂, are obtained by using fast Fourier transformation to show the spectrum components of the motion. In this paper, the frequency axis of the power spectrum plot is normalized by using the rotor mass and rotational speed.

Bifurcation diagram

To generate a bifurcation diagram, the rotor mass (or speed) is varied with a constant step, and the state variables at the end of the integration process are used as the initial values for the next mass (or speed). This means that there is a tendency for the integration to follow a single response curve. Variation of the X(nT) and Y(nT) coordinates of the return in the Poincaré map with rotor mass(or speed) is then plotted to form the bifurcation diagram.

Results and discussion

Dynamic analysis of the bearing system with the results of different numerical methods

For the OHGB system, three different numerical methods including the FDM, perturbation method, and proposed method are used to compare the numerical simulation results for opposed wedge-shaped bearings (opposed obverse wedges and opposed reverse wedges). The result shows the proposed mixing method binding DTM and the FDM, the traditional FDM and perturbation method¹² have a better fit when the rotor has a large mass. This also shows that the mixing method raised in this paper reveals higher accuracy and better convergence than the other methods with respect to different time steps and different opposed wedge-shaped bearings. The comparison of rotor center trajectory values on Poincaré map (X2(nT) vs. Y2(nT)) is as shown in Table 2 and the behaviors of rotor center with specific operational parameters are all periodic motions. Then the simulation results in Table 2 can be compared objectively and the differences were observed by three numerical methods. For numerical stability analysis, the impact of the different time steps is studied and the displacement changes are compared under different operating conditions by the numerical value of the Poincaré section, as shown in Table 3. From Tables 2 and 3, DTM and FDM (mixing method) is applied to obtain the dynamic orbits under different operating parameters for the OHGB system. For different values of the time step $\tilde{H}$ , it is found that the system orbit performs five decimal places by DTM and FDM better than the other two methods with different types of OHGB including opposed obverse wedge shape and opposed reverse wedge shape bearings shown in Table 2. It can also be seen that the orbits of Poincaré maps for the opposed obverse wedge-shaped bearing are in agreement to approximately five decimal places for the different time increments obtained by DTM and FDM at different rotor masses as shown in Table 3.

Table 2.

Comparison of opposed wedge-shaped bearings (the opposite obverse wedge shape and opposite reverse wedge shape) system for rotor center displacement on Poincaré map with various numerical calculations.

DisplacementNumerical methods and operating conditions			X2(nT)		Y2(nT)
DisplacementNumerical methods and operating conditions			$\tilde{H} = 0 .001$	$\tilde{H} = 0 .01$	$\tilde{H} = 0 .001$	$\tilde{H} = 0 .01$
Opposed obverse wedge shape	FDM¹⁵	m_r =5.3 kgΛ=4.2	–0.1087137902	–0.1088537204	0.1489311312	0.1488991226
	DTM and FDM		–0.1088447339	–0.1088435108	0.1488009818	0.1488037843
	Perturbation method		–0.1089171212	–0.1088125852	0.1487818632	0.1487199569
Opposed reverse wedge shape	FDM	m_r =5.3 kgΛ=4.2	–0.1433787451	–0.1431442008	0.1544627618	0.1542364053
	DTM and FDM		–0.1432811874	–0.1432898964	0.1543399815	0.1543381615
	Perturbation method		–0.1434782833	–0.1433421190	0.1543261862	0.1545368265

DTM: differential transformation method; FDM: finite difference method.

Table 3.

Poincaré map numerical comparison for different time steps $\tilde{H}$ under different rotor mass (the results calculated by the DTM and FDM mixing method and the opposed obverse wedge-shaped bearing is taken as an example.).

m_r = 6.2 kg			m_r =8.14 kg
$\tilde{H}$	X2(nT)	Y2(nT)	$\tilde{H}$	X2(nT)	Y2(nT)
π/300	–0.1330085926	0.1459065217	π/300	–0.0949159566	0.2038387436
π/600	–0.1330049029	0.1459035148	π/600	–0.0949136681	0.2038321584

For the grid number analysis, the traditional FDM, perturbation, and mixing methods are used to solve equation (1) with uniform mesh size. Meanwhile, the three methods are calculated and compared with different values of the mesh size in space. It can be seen that the rotor center performs four decimal places by mixing method (DTM and FDM) better than the traditional FDM and perturbation method with different amount of grids shown in Table 4 (in the revised paper). It can also be seen that the rotor orbits are in agreement to approximately four decimal places for the different mesh sizes obtained by DTM and FDM at the values of bearing number (Λ = 4.2 and 7.1) as shown in Table 4.

Table 4.

Comparison of Poincaré maps of rotor center with different values of Λ and mesh size by FDM, mixing method, and perturbation method (time step: 0.01, m_r =5.3 kg).

Mesh size θ× $η$ Numerical methods		Λ =4.2		Λ= 7.1
		X2(nT)	Y2(nT)	X2(nT)	Y2(nT)
FDM	300 × 400	−0.1052676758	0.1461092741	0.6462281104	0.3293100827
DTM and FDM		−0.1088435108	0.1488037843	0.6426031498	0.3212320136
Perturbation method		−0.1039827871	0.1434840745	0.6463938891	0.3213155516
FDM	600 × 800	−0.1074930892	0.1422010995	0.6490348585	0.3228400422
DTM and FDM		−0.1088704818	0.1488144025	0.6426931214	0.3212108381
Perturbation method		−0.1039346372	0.1426370104	0.6473851651	0.3288238202

The opposed obverse wedge-shaped bearing is taken as an example. DTM: differential transformation method; FDM: finite difference method.

From the above numerical results, it can be seen that the mixing method of this study has the best accuracy than FDM and perturbation method. Meanwhile, it has better applicability for this OHGB system. Whether for different opposed-type bearings, the system trajectory can be calculated effectively and accurately by the mixing method. At the same time, in order to shorten the bifurcation characteristics of the subsequent calculation, the test results in Tables 2, 3, and 4 show that the simulations with sufficient accuracy can be obtained without too small time steps and too many grid numbers. Therefore, in the subsequent dynamic analysis section, π/300 is taken as a time step and 300 × 400 is taken as a grid number to calculate.

Rotor dynamic behavior analysis of the gas bearing system

Rotor impact on the system

Take the opposed obverse wedge-shaped bearing as an example, the bearing number in the system is (rotational speed) Λ = 4.5, and the mass of the rotor m_r is taken as a bifurcation parameter.

Dynamic trajectories map

From Figure 4.1(a), (b), 4.2(a), (b),…, 4.6(a), (b) it can be observed that under the case of the lower mass of the rotor center and of the journal center (m_r =5.3 kg), the trajectory shows the regular pattern, and when the mass raises to m_r =6.98 kg, the nonperiodic phenomenon replaces the regular motion. At the same time when the rotor mass has reached m_r =7.71 and 8.11 kg, the dynamic trajectory continues the cycle and the sub-cycle movement, but this cycle is in response to m_r =8.25 kg, and then turns into nonperiodic motion. And when m_r =9.39 kg, the bearing system remains unstable with nonperiodic behavior.

Figure 4.

The dynamic trajectory (Figure 4.1(a) to 4.6(a)) of the rotor center when the rotor mass is m_r = 5.3, 6.98, 7.71, 8.11, 8.25, 9.39, phase plan (Figure 4.1(b) to 4.6(b)), spectral response of the center of the rotor in the horizontal direction (Figure 4.1(c) to 4.6(c)), and the vertical direction (Figure 4.1(d) to 4.6(d)), the system bearing number is Λ = 4.5.

b. Spectrum analysis chart

Figure 4.1(c), (d), 4.2(c), (d),…, 4.6(c), (d) shows the dynamic response of the rotor center in the horizontal and vertical directions. The studies show when the rotor mass is m_r = 5.3 kg, the rotor center shows T-periodic motion, and when m_r = 6.98kg, the spectral response map (Figure 4.2(c) and (d)) displays the rotor center occurring as a quasi-periodic motion in both the horizontal and perpendicular direction. As the rotor mass changes to m_r = 7.71 and 8.11 kg respectively, the motion turns into T-periodic motion and the 3T cycle motion of subharmonic vibrations. With the mass increased to m_r = 8.25 and 9.39 kg, then the dynamic trajectory mutates into chaotic motion.

c. Bifurcation diagram

Shown in Figure 5, the rotor mass m_r is taken as an analysis parameter to explore the effect of different rotor masses on the system, while setting the mass range from 0.01 to 9.4 kg for the actual operating conditions. From Figures 5 to 7, we can see that system instability will occur earlier when the number of bearings (rotation speed) is larger. Therefore, the paper analyzes the influence of the number of bearings (rotation speed) later. At the same time, as it can be seen from Figure 5 when m_r < 6.97 kg, the rotor center of the system shows periodic motion in both the horizontal and vertical directions. This phenomenon of motion can be further proved by the Poincaré section (Figure 8(a)). However, as the mass of this cycle motion increased to m_r = 6.98 kg, the rotor center will be replaced by quasi-periodic movement, as shown in Figure 8(b), and this state of motion is maintained within the 6.97≤ m_r <7.71 kg range. But as the mass increases to m_r =7.71 kg, quasi-periodic motion quickly changes to T-periodic motion, and this state can be further confirmed by Figure 8(c).

Figure 5.

Bifurcation map of the rotor center to different rotor masses: (a) X2(nT); (b) Y2(nT), the number of system bearings Λ = 4.5.

Figure 6.

Bifurcation map of the rotor center to different rotor masses: (a) X2(nT); (b) Y2(nT), the number of the system bearing Λ = 6. 0.

Figure 7.

Bifurcation map of the rotor center to different rotor masses: (a) X2(nT); (b) Y2(nT), the number of the system bearing Λ = 2.0.

Figure 8.

Poincaré sectional view of the rotor center at different rotor masses: (a) m_r = 5.3; (b) 6.98; (c) 7.71; (d) 8.11; (e) 8.25; (f) 9.39 kg (system bearing number Λ = 4.5).

When the mass increases to m_r =8.11 kg, the system turns to the sub-cycle motion, as shown in Figure 8(d), while the movement is only maintained at 8.11≤ m_r <8.14kg range. As the mass continues to increase to m_r =8.18 kg, the system’s motion behavior is T and 3T intermittent movement. This phenomenon occurs in the 7.71≤ m_r <8.25kg range. When m_r =8.25kg, the system turns from 3T cycle to chaotic movement. The movement is only maintained at 8.25≤ m_r <8.93kg. Afterward the system presents the movement of Chaos and 3T vibration forms. The movement is only maintained at 8.25≤ m_r ≤9.4 kg range. When m_r =8.25 and 9.39 kg, the system is in a state of chaotic motion, which can be further proved by Figures 8(e) and (f). From the above results showing that the rotor movement patterns also showed a considerable degree of cyclical order, the following are the movement patterns presented by the system for rotors with different masses and makes a whole summary as shown in Table 5.

Table 5.

Rotor center dynamic behavior under the rotor mass change (system bearing number Λ = 4.5).

Rotor mass	[0.01, 6.97)	[6.97, 7.71)	[7.71, 8.11)	[8.11, 8.14)	[8.14, 8.18)
Dynamic behavior	T	Quasi	T	3T	T
Rotor mass	[8.18, 8.25)	[8.25, 8.93)	[8.93, 8.97)	[8.97, 8.98)	[8.98, 9.0)
Dynamic behavior	3T	Chaos	3T	Chaos	3T
Rotor mass	[9.0, 9.02)	[9.02, 9.07)	[9.07, 9.08)	[9.08, 9.1)	[9.1, 9.13)
Dynamic behavior	Chaos	3T	Chaos	3T	Chaos
Rotor mass	[9.13, 9.15)	[9.15, 9.16)	[9.16, 9.19)	[9.19, 9.2)	[9.2, 9.39)
Dynamic behavior	3T	Chaos	3T	Chaos	3T
Rotor mass	[9.39, 9.4]
Dynamic behavior	Chaos

Shown by Figure 8(a), (c) and (d), while m_r =5.3, 7.71, 8.11 kg, the motion of the bearing presented in a Poincaré sectional view is a single point or three discrete points, the system is in a single cycle movement or 3T periodic motion. If the motion state shown in Figure 8(b) is a single closed curve formed by multiple discrete points, it shows quasi-periodic motion, while Figure 8(e) and Figure 8(f) are irregular discrete points, which means the system is in chaotic movement.

For the confirmation of chaos behavior, in the execution of this paper, interpret and further verify the maximum Lyapunov exponents. As it can be seen from Figure 9(a) to (d), as m_r = 5.3, 6.98, 7.71, 8.11 kg, and the Lyapunov index tends to zero, then the system is nonchaotic behavior. And when m_r =8.25 and 9.39, as shown in Figure 9(e) and (f), the index is greater than zero, so the system’s motion state is indeed chaotic behavior. At the same time, the overall Lyapunov exponent for the system is also shown in Figure 10, which shows that the region of chaos (index greater than zero) occurs in several intervals, especially within the 8.25≤ m_r ≤9.4 kg range, which is consistent with the results summarized in Table 5.

Figure 9.

Lyapunov index of the rotor center at different rotor masses: (a) m_r = 5.3; (b) 6.98; (c) 7.71; (d) 8.11; (e) 8.25; (f) 9.39 kg (system bearing number Λ = 4.5).

Figure 10.

Maximum Lyapunov exponent with different rotor masses: (a) overall; (b) 8.0 ≤ m_r ≤ 9.4 kg (system bearings number Λ = 4.5).

Effect of the bearing number Λ on the system

Take an obverse wedge shape opposing bearing as an example. The rotor speed will affect the related performance of the entire system, and therefore the bearing number Λ is taken as a bifurcation parameter, the rotor mass is set as m_r =3.5 kg, and analyze the dynamic system behavior.

Dynamic trajectory map

From Figure 11.1(a), (b), 11.2(a), (b),…, 11.4(a), (b) it can be observed that, when there is a small number of bearings (Λ = 3.2), the behavior presented by the rotor is more regular, and when Λ increases to 5.07, the system still shows regular motion. When the bearing number reaches 6.51, the dynamic trajectory changes to non-regular motion, and when Λ = 6.97, the behavior of the rotor center is transformed into a regular phenomenon.

Figure 11.

Rotor center dynamic trajectory (Figure 11.1(a) to 11.4(a)) when the bearing number is Λ = 3.2, 5.07, 6.51, 6.97, journal center dynamic trajectory map (Figure 11.1(b) to 11.4(b)), spectral response of the center of the rotor in the horizontal direction (Figure 11.1(c) to 11.4(c)) and the vertical direction (Fig 11.1(d) to 11.4(d)) (rotor mass as m_r=3.5 kg).

b. Spectrum analysis map

Figure 11.1(c), (d), 11.2(c), (d),…, 11.4(c), (d) display the dynamic response of the rotor center in horizontal and vertical directions with a different number of bearings. Research shows when the number of bearings is Λ = 3.2, the rotor center presents T-periodic motion, and when Λ = 5.07, the spectral response graph (Figure 11.2(c), (d)) displays the rotor center presenting 2T subharmonic motion both in the horizontal and vertical directions. As the number of bearings increases to Λ = 6.51, the mode of motion transforms into nonlinear chaotic motion. At Λ = 6.97, the system turns into T-periodic motions, and its shape is the same both in the horizontal and vertical directions.

c. Bifurcation diagram

As shown in Figure 12, the number of bearings Λ is taken as the main analysis parameter, to study the effect of the bearing number Λ on systems, while setting the range of the number of bearings from 1.0 to 7.5 for the actual operating conditions. When the number of low bearings is Λ = 3.2, regardless of the horizontal and vertical directions, the system rotor center presents T-period movement, this phenomenon may be verified by Figure 13(a). When increasing the number of bearings to Λ = 5.07, the system in turn generates 2T movement, and from Figure 13(b) it is clear at this time that the system generates two discrete points. This 2T movement stands at just 5.07≤Λ < 6.51. When increasing the number of bearings to Λ = 6.51, the system turns to produce chaotic motion, this status may be further verified by Figure 13(c) and this motion state is maintained within the range 6.51 ≤ Λ < 6.97. As Λ increases to 6.97, the bearing system tends to be stable and exhibits T-period motions, as shown in Figure 13(d).

Figure 12.

The bifurcation diagram of the rotor center for different bearing numbers Λ: (a) X2(nT); (b) Y2(nT) (rotor mass at m_r=3.5 kg).

Figure 13.

Poincaré cross-sectional view of the center of the rotor in different number of bearing: (a) Λ = 3.2; (b) 5.07; (c) 6.51; (d) 6.97 (rotor mass at m_r=3.5 kg).

The relationship between the rotor center motion and the number of bearings is quite complex, so the range of 1.0 ≤ Λ ≤ 7.5 is summarized in Table 6 as follows.

Table 6.

Dynamic behavior of the rotor center with different bearing numbers (1.0 ≤ Λ ≤ 7.5, the rotor mass as m_r=3.5 kg).

Bearing number Λ	[1.0, 5.07]	[5.07, 6.51)	[6.51, 6.97)	[6.97, 7.5]
Dynamic behavior	T	2T	Chaos	T

For the confirmation of chaos behavior, in the execution of this paper, we further verify and interpret the maximum Lyapunov exponents. From Figure 14(a), (b), and (d) it can be seen that when Λ = 3.2, 5.07, 6.97, and the Lyapunov exponents are less than zero, the system is in nonchaotic behavior. When Λ = 6.51, in Figure 14(c), the index is greater than zero, so the system motion state is indeed chaotic behavior. While the whole picture of the influence of the system’s Lyapunov exponents on the changing bearings number is also shown in Figure 15, the results show a chaotic region (index greater than zero) occurs only in the 6.51≤Λ < 6.97 range, and the remaining area is a stable operating range, which is consistent with the results summarized in Table 6.

Figure 14.

Lyapunov exponents of rotor center behavior at bearing number: (a) Λ = 3.2; (b) 5.07; (c) 6.51; (d) 6.97 (rotor mass at m_r=3.5 kg).

Figure 15.

Overall system Lyapunov exponent with different bearing numbers (rotor mass at m_r = 3.5 kg).

Stable and unstable areas for OHGB system considering various rotor masses and bearing numbers

In order to know whether the behavior of the bearing system under various parameters is stable or not, this paper uses different rotor masses and bearing numbers as the analysis parameters and constructs both stable and unstable regions of the shaft through Lyapunov exponents, as shown in Figure 16. The white, red, and pink blocks in the figure belong to the unstable region while the other color blocks belong to the stable region. Therefore, the unstable regions of the system are mostly concentrated in the upper right corner of the figure (the rotor mass and bearing number are larger). And in the low number of bearings, the system is stable and has periodic motion, the results can effectively identify and establish the stable and unstable areas of the opposed-type high-speed gas bearing system in operation, and can be used as an important reference basis for the follow-up bearing system design.

Figure 16.

Full-size image of Lyapunov exponents of stable and unstable areas of opposed high-speed gas bearing systems.

Conclusion

The objective of this study was an analysis of the dynamic behavior of a flexible rotor supported by an OHGB system. The FDM, perturbation method, and a mixing method were used to solve the pressure distribution at the highest nonlinearity in the system, after which dynamic equations of the flexible rotor center were used to obtain the orbits of the rotor center. Analysis was then conducted on the orbit data to generate spectrum diagrams, Poincaré maps, bifurcation diagrams, and maximum Lyapunov exponents. This showed that the orbits of the rotor center change along with the mass and the bearing number and synchronously generate complicated motion in the horizontal and vertical directions, include normal and sub-harmonic vibration, as well as periodic, quasi-periodic, and chaotic motion. The mixing method applied in this paper yields better accuracy and precision for the verification of values than the FDM and perturbation method. Accurate solutions can be obtained without the need for very precise position and time step determination. Findings with respect to the motion and behavior of the rotor revealed that the rotor orbits demonstrate unstable nonlinear chaotic behavior when the mass of the rotor is high but the system is relatively stable at the other intervals and ranges. The bearing number affects the orbits of the rotor, which show chaotic behavior in the interval 3.84 ≤ Λ < 3.86 while the orbits in the other intervals and ranges show relatively organized and stable motion. The results obtained in this study can be used as a basis for future OHGB system design and the prevention of instability.

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 study was funded by the Ministry of Science and Technology of Taiwan (Grant numbers MOST-106–2221-E-167–012-MY2 and MOST-106–2622-E-167–010-CC3).

References

Marui

Ema

Kato

Contact rigidity at tool shank of turning tools. ASME J Eng Ind 1987; 109: 169–172.

Rivin

Trends in tooling for CNC machine tools: Machine system stiffness. Manuf Rev (USA) 1991; 4: 257–263.

Chen

Yang

, et al. Modeling cone self-acting gas lubrication bearing dynamics. Procedia Eng 2015; 126: 416–420.

Nguyen

SH.

P-version finite element analysis of gas bearings of finite width. ASME J Tribol 1991; 113: 417–420.

Malik

Bert

CW.

Differential quadrature solution for steady state incompressible and compressible lubrication problems. ASME J Tribol 1994; 116: 296–302.

Wang

Schellekens

PH.

Tri-conical gap-shaped externally pressurized gas bearing pads. Tribol Trans 1994; 37: 850–856.

Whipple

RT.

The inclined groove bearing. Harwell: Atomic Energy Research Establishment, 1958.

Vohr

Pan

CHT.

On the spiral-grooved, self-acting gas bearing, MTI Technical Report MTI63TR52, Prepared under Ofﬁce of Naval Research Contract Nonr-3780(00), Task NR061-131, 1963.

Absi J and Bonneau D. Caracteristiques Statiques et Dynamiques desPaliers a Chevrons. Proc. of 10 erne Congres Francois de Mecanique, Paris, 1991, Vol. 1, pp. 253--256.

10.

Bently

Forced subrotative speed dynamic action of rotating machinery. In: Petroleum mechanical engineering conference, Dallas. TX, USA, 1974.

11.

Choi

Noah

ST.

Nonlinear steady-state response of a rotor-support system. ASME J Vib Acoust Stress Reliab Des 1987; 109: 255–261.

12.

Kazimierski

Optimization of hybrid gas lubricated conical bearings. J Mech Eng 1972; 14: 312–318.

13.

Agrawal

GL.

Analysis of power-shaped multipad hydrodynamic conical gas bearings. Tribol Trans 1993; 36: 240–248.

14.

Yoshimoto

Oshima

Danbara

, et al. Stability of water-lubricated, hydrostatic, conical bearings with spiral grooves for high-speed spindles. ASME J Tribol 2002; 124: 398–405.

15.

Pan

CHT

Kim

Stability characteristics of a rigid rotor supported by a gas-lubricated spiral-groove conical bearing. ASME J Tribol 2007; 129: 375–383.

16.

Wang

Yau

HT.

Theoretical analysis of high speed spindle air bearings by a hybrid numerical method. Appl Math Comput 2010; 217: 2084–2096.

17.

Wang

CC.

Application of a hybrid method to the bifurcation analysis of a relative short gas journal bearing system with herringbone grooves. Ind Lubr Tribol 2011; 63: 307–319.

18.

Wang

CC.

Bifurcation and nonlinear dynamic analysis of united gas-lubricated bearing system. Comput Math Appl 2012; 64: 729–738.

19.

Wang

CC.

Bifurcation and nonlinear dynamic analysis of noncircular aerodynamic journal bearing system. Nonlinear Dyn 2013; 72: 477–489.

20.

Young

A high precision ultrasonic system for vibration measurements. Rev Sci Instrum 1992; 63: 5435–5441.

21.

Nikolaou

, et al. Demodulation of vibration signals generated by defects in rolling element bearings using complex shifted Morlet wavelets. Mech Syst Signal Process 2002; 16: 677–694.

22.

Kung

YT.

A parametric study on oil/air lubrication of a high-speed spindle. Precis Eng 2005; 29: 162–167.

23.

Al-Raheem

Roy

Ramachandran

, et al. Application of the Laplace-wavelet combined with ANN for rolling bearing fault diagnosis. J Vib Acoust 2008; 130: 221–229.