Bifurcation and Nonlinear Behavior Analysis of Dual-Directional Coupled Aerodynamic Bearing Systems

: Dual-directional coupled aerodynamic bearing (DCAB) systems have received considerable attention over the past few years. These systems are primarily used to solve air lubrication problems in high-precision mechanisms and equipment that run at a high rotational speed and require high rigidity and precision. DCABs have the advantages of axial and radial thrust and provide high rigidity, dual-directional support, and high load-carrying capacity. In DCAB systems, the nonlinearity of the air ﬁlm pressure and dynamic problems, such as critical speed, unbalanced air supply, or poor design, can cause the instability of the rotor-bearing system and phenomena such as nonperiodic or chaotic motion under certain parameters or conditions. Therefore, to investigate what conditions lead to nonperiodic phenomena and to avoid irregular vibration, the properties and performance of the DCAB system were explored in detail by using three numerical methods for verifying the accuracy of the numerical results. The rotor behavior was also studied by analyzing the spectral response, the bifurcation phenomenon, Poincar é maps, and the maximum Lyapunov exponent. The numerical results indicate that chaos occurs in the DCAB system for speciﬁc ranges of the rotor mass and bearing number. For example, when the rotor mass ( m r ) is 5.7 kg, chaotic regions where the maximum Lyapunov exponents are greater than 0 occur at bearing number ranges of 3.96–3.98 and 4.63–5.02. The coupling e ﬀ ect of the rotor mass and bearing number was also determined. This e ﬀ ect can provide an important guideline for avoiding an unstable state. the bearing number is to 5.02, the system returns to regular behavior. The aforementioned results indicate that distinctive changes occur in the behavior of the system as the bearing increases marginally from 3.96 to 3.98. The rotor motion is nonperiodic when Λ = 3.96 and periodic when Λ = 3.98. Thus, the DCAB system is very sensitive to the bearing number.


Introduction
Air lubrication theory was introduced in 1886 by Reynolds, who derived the well-known Reynolds equation that indicates the relationship among air pressure distribution, density, and velocity [1]. In 1958, Whipple [2] analyzed the characteristics of herringbone-grooved air bearings. In 1959, Whitley [3] published many designs for herringbone grooves on the basis of Whipple's research. In 1963, Vohr [4] conducted a study on spiral-grooved air bearings and obtained a modified Reynolds equation for air-lubricated pressure by assuming an infinite number of grooves. The theory proposed by Vohr came to be known as narrow groove theory. Studies have indicated that the type of grooves on the surface of air bearings considerably influences the bearing system. Absi [5] analyzed the properties of herringbone-grooved bearings, including their stability limit, their gas pressure distribution, and the influence of certain important parameters on them, by using the finite-element method. In 2005, 2 of 23 Coleman [6] developed an air-bearing system with high load-carrying efficiency and successfully applied it in laser-scanning-related technologies. Because of the increased demand for accuracy in the laser printing market, the technical demand for air-bearing scanners has increased each year. In 2009, Liu et al. [7] proposed a numerical method for solving the isothermal compressible Reynolds equation. They linearized the aforementioned equation by using the Newton method and obtained the pressure distribution of gas film lubrication through iterative continuous relaxation. The aforementioned research considers the theoretical basis of the air pocket externally pressurized bearing system, couples the effects of aerodynamics and aerostatic forces, and adds the rotation effect of the rotor. The results of Liu et al. indicated that their proposed numerical method was suitable for the calculation and analysis of the gas film pressure distribution and bearing capacity.
In 2014, Simek [8] constructed an air circulation machine with a rolling bearing rotor support and used it as a test platform to verify the performance of an air-bearing system. Simek's test results proved the feasibility of the examined air-bearing system. The bearing was a combination of an elastically supported tilt pad journal bearing and a thrust bearing with spiral grooves. The bearing had a maximum rotation speed of 65,000 rpm. At this speed, the root mean square value of the maximum amplitude generated by the rotor rotation was 1.5 µm, which was considerably lower than the limit value of 4.5 µm. The small rotor vibration value of the aforementioned air-bearing system proved that it had a very high operating efficiency and could maintain good rotor balance.
In 2016, Bao and Mao [9] proposed an innovative dynamic modeling method for air-bearing systems. In this method, the dynamic characteristics, direction of motion, and support are considered. The aforementioned dynamic modeling method was used to investigate an ultra-precision positioning dual-platform system, which was used in integrated circuit manufacturing equipment by using two sets of air bearings. The theoretical dynamic behavior of the ultra-precision positioning dual platform was also studied. Bao and Mao compared the theoretical dynamic behavior with the experimental results to verify the effectiveness, robustness, and correctness of their proposed method.
In 2018, Guo and Gao [10] designed an active noncontact journal bearing system with a two-way drive support function. They analyzed the basic design methods and principles of this system and evaluated the performance of the developed prototype. The aforementioned system includes a shaft positioner, noncontact journal bearing, and rotating motor that uses a vibration coupling mode. Noncontact rotation controls the pressure distribution of the air film through the coupled resonance mode. Moreover, the near-field sound force is used to generate a stable air film pressure so that the rotor can rotate and levitate stably in the bearing. This paper presents a dual-directional coupled air-bearing (DCAB) system for the application of supporting forces in the radial and thrust directions. The performance of the DCAB system was analyzed in terms of the vibration and dynamic performance of the bearing rotors.
Regarding the rotor motion behavior in an air-bearing system, Bently [11] found that the rotor in an oil film bearing system exhibits the second and third harmonic vibration behavior during the rotation process. Child et al. [12] proved that subharmonic vibration behavior occurs when the rotor rotates under the support of a bearing. Wang et al. [13][14][15] have analyzed the gas film pressure distribution in different air-bearing systems and have proposed effective solutions for investigating rotor dynamics. They proposed a numerical method with fast convergence for analyzing the Reynolds equation. The aforementioned authors also obtained the rotor behavior as well as the dynamic behavior of the rotor and journal center by analyzing the rotor trajectory, vibration spectrum, and bifurcation characteristics under different parameters. The results of the studies of Wang et al. indicated that the rotor and journal center exhibited periodic, quasi-periodic, and subharmonic motion under different operating conditions. Moreover, chaotic motion behavior was observed in some special air-bearing systems.
In the past 10 years, many researchers have successfully published machine learning schemes for the identification and classification of fault conditions generated by bearing systems. The proposed schemes include artificial neural networks, support vector machines, and principal component analysis. In particular, the application of deep learning methods has attracted considerable interest from the industry and academia [16,17]. For bearing fault detection, Ni et al. [18] proposed a method based on random matrix theory to evaluate the degradation of the rolling bearing system and increase the production safety. Their results proved that their method can be applied to rolling bearing systems for identifying different degradation and safety states. Moreover, the aforementioned method can be used to determine accurately the degradation rate of the bearing performance.
An air-film-bearing system, which is nonlinear, is considerably different from an oil-bearing system. Traditionally, the steady states of a deterministic nonlinear dynamic system are categorized as equilibrium points, periodic motion, and quasi-periodic motion. These three major states are commonly called attractors because a steady-state system is attracted to one of these states after transience. However, random chaotic motion can also develop in nonlinear dynamic systems. Once chaotic behavior begins, the outcome is unpredictable and may involve damage or even destruction. Therefore, chaos must be avoided, particularly in precision bearing systems. In this regard, this study mainly focused on the analysis of the properties of DCAB systems. The range of dynamic behaviors exhibited by this system under different operating conditions was also studied. In addition, the bifurcation properties of the nonlinear behavior produced by the system rotor were examined. These properties can be used to judge if the system displays chaos and predicts dynamic system trajectories accurately.

Design and Analysis of a DCAB System
This paper mainly focuses on the DCAB system used in the ultra-precision machinery industry. The aforementioned system is especially suitable for an environment that requires an extra-high rotation speed and a high load-carrying capability. As displayed in Figure 1, the key feature of a DCAB system is the provision of dual-directional support capability obtained through the combination of thrust and radial bearings. The industry is currently facing the problem of short bearing life caused by contact friction. Collisions between the bearing surfaces often occur when the rotor starts, stops, or is running. Collisions damage the bearing surfaces and can lead to a loss of support if the air cyclone is impeded. Instability becomes a serious problem when frequent collisions occur between the bearing surfaces.
Symmetry 2020, 12, x FOR PEER REVIEW 3 of 23 component analysis. In particular, the application of deep learning methods has attracted considerable interest from the industry and academia [16,17]. For bearing fault detection, Ni et al. [18] proposed a method based on random matrix theory to evaluate the degradation of the rolling bearing system and increase the production safety. Their results proved that their method can be applied to rolling bearing systems for identifying different degradation and safety states. Moreover, the aforementioned method can be used to determine accurately the degradation rate of the bearing performance.
An air-film-bearing system, which is nonlinear, is considerably different from an oil-bearing system. Traditionally, the steady states of a deterministic nonlinear dynamic system are categorized as equilibrium points, periodic motion, and quasi-periodic motion. These three major states are commonly called attractors because a steady-state system is attracted to one of these states after transience. However, random chaotic motion can also develop in nonlinear dynamic systems. Once chaotic behavior begins, the outcome is unpredictable and may involve damage or even destruction. Therefore, chaos must be avoided, particularly in precision bearing systems. In this regard, this study mainly focused on the analysis of the properties of DCAB systems. The range of dynamic behaviors exhibited by this system under different operating conditions was also studied. In addition, the bifurcation properties of the nonlinear behavior produced by the system rotor were examined. These properties can be used to judge if the system displays chaos and predicts dynamic system trajectories accurately.

Design and Analysis of a DCAB System
This paper mainly focuses on the DCAB system used in the ultra-precision machinery industry. The aforementioned system is especially suitable for an environment that requires an extra-high rotation speed and a high load-carrying capability. As displayed in Figure 1, the key feature of a DCAB system is the provision of dual-directional support capability obtained through the combination of thrust and radial bearings. The industry is currently facing the problem of short bearing life caused by contact friction. Collisions between the bearing surfaces often occur when the rotor starts, stops, or is running. Collisions damage the bearing surfaces and can lead to a loss of support if the air cyclone is impeded. Instability becomes a serious problem when frequent collisions occur between the bearing surfaces. A coupled air-bearing system with five degrees of freedom (x, y, z, θx, and θy,) was designed. Moreover, the performance of single-thrust and radial bearings was compared. This comparison was conducted to determine the optimal design criteria for a coupled air-bearing system and to obtain a complete picture of the different motion behaviors of this system. The Reynolds equations of a coupled air-bearing system, including the lubrication equations of the thrust and radial bearings A coupled air-bearing system with five degrees of freedom (x, y, z, θ x , and θ y, ) was designed. Moreover, the performance of single-thrust and radial bearings was compared. This comparison was conducted to determine the optimal design criteria for a coupled air-bearing system and to obtain a complete picture of the different motion behaviors of this system. The Reynolds equations of a coupled air-bearing system, including the lubrication equations of the thrust and radial bearings (presented in Equations (1) and (2), respectively), can be derived on the basis of narrow groove theory, [1,2]: where R represents the journal radius, h is the air film thickness, P is the air film pressure, µ is the viscosity coefficient, and rθz is the coordinate system. A perturbed system was analyzed using the perturbation method. In addition, the first-order perturbation equations for air film thickness and pressure were derived. The physical parameters of the perturbation (ξ) are x, y, z, θ x , θ y , .
x, rθ plane: The parameter z o in Equations (5) and (6) represents the coordinates of the center of mass. The pressure functions can be obtained through numerical calculation. The rigidity (K), damping (Ĉ), load bearing, and flow can then be obtained using the integrals of the pressure functions as follows [19]: The performance parameters of a bearing system in operation can be determined using Equations (7)- (10). The discriminant for determining if a system is stable can then be obtained using the Routh-Hurwitz method. After the optimized performance coefficients are determined, variables such as the optimal bearing dimensions, rotational speed parameter, number of grooves, and thrust disc diameter can also be obtained.
The analysis and design of the coupled air-bearing system are described in the aforementioned text. An optimal bearing system that can be operated under different conditions was also developed. The motion behaviors of the rotor in the coupled bearing system, including the dynamic trajectory, spectral response, bifurcation behavior, Poincaré map, and Lyapunov exponents, were also analyzed. These factors can be used as bases for developing effective control strategies so that the nonperiodic motion can be suppressed.

Equations of Rotor Dynamics
This study considered a flexible rotor of mass m r supported by DCABs mounted on rigid pedestals. The rotor rotates with two degrees of oscillated displacement in the transverse plane. The equations of rotor dynamics for the transient state are presented in Cartesian coordinates in Equations (11) and (12). Moreover, the equilibrium equation of forces at the journal center is presented in Equation (13). In Equation (13), (x 2 , y 2 ) and (x 3 , y 3 ) are the displacements of the rotor and journal centers, respectively; ω r is the rotor rotational speed; ρ is the rotor eccentricity; and K r is the rotor stiffness coefficient.
The following nondimensional groups are presented in Table 1: Table 1. Nomenclature of the variables used in this study.

Variable Meaning
x, y, z, θ x and θ y degrees of freedom for five directions R journal radius L bearing length h air film thickness P air film pressure P a atmospheric pressure µ viscosity coefficient r,θ,z coordinate system t time ξ physical quantity of perturbation K J ,K T rigidity of journal bearing and thrust bearing, respectivelŷ C J ,Ĉ T damping of journal bearing and thrust bearing, respectively m r rotor mass x 2 , y 2 displacements of rotor center x 3 , y 3 displacements of journal center ω r rotational speed of rotor ρ eccentricity of rotor K r stiffness coefficient of rotor f ax , f ay equilibrium equation of forces acting to the journal center

Mathematical Formulation of the Numerical Simulation
To determine the gas pressure distribution with the Reynolds equation, three numerical methods, namely the perturbation method, the finite-difference method (FDM), and a hybrid method, were used and their solutions were verified [14,15]. Equations (3) and (4) are used in the perturbation method to solve the Reynolds equations. In the FDM, Equations (1) and (2) are discretized using the central-difference method in the r-direction, θ-direction, and z-direction with a uniform mesh size. Moreover, the implicit backward-difference method is applied in the time domain. Equations (1) and (2) are then be transformed into Equations (17) and (18), respectively. Subsequently, the pressure distribution for the DCAB system at each time step is obtained through an iterative calculation process as follows: rθ plane: θz plane: A hybrid scheme was used with the differential transformation method (DTM) and FDM to determine the gas pressure distribution. The DTM is applied to the time domain of Equations (1) and (2). Moreover, this method is widely used for solving nonlinear differential equations due to its high convergence rate and calculation accuracy. Therefore, the DTM is first used to discretize the time domain of the Reynolds equations into intervals. The central-difference scheme is then applied to the r-, θ-, and z-coordinates of the locations. The principle of the DTM is presented in [14,15]. The DTM is used to transform Equations (1) and (2) with respect to τ as follows: rθ plane: θz plane: where Subsequently, the FDM is used to discretize Equations (19) and (20) with respect to the corresponding coordinate directions. Moreover, Equations (21) and (22) are substituted into Equations (19) and (20), respectively, to obtain the following equations: Symmetry 2020, 12, 1521 8 of 23 rθ plane: θz plane: where H is the time step.
To calculate the dynamic behavior of the rotor center effectively, the acceleration of the rotor center is first calculated through the balance equation. Subsequently, the velocity and displacement are calculated. The aforementioned iterative calculation process begins from the initial static equilibrium state. The initial displacement of the rotor is determined by the static equilibrium position. Moreover, the gap thickness between the rotor and the bearing is determined. Assume that the initial velocity of the rotor is 0. The calculation process of the overall iteration ( Figure 2) is described in the following text.
First step: After the time increment ∆τ has passed, recalculate the acceleration, velocity, and displacement of the rotor. Moreover, define the new gap thickness value h at different positions between the rotor and the bearing.
Second step: Substitute the new h value into Equations (23) and (24). Then, recalculate the new pressure distribution state between the rotor and the bearing.
Third step: Integrate the pressure distribution calculated in step 2 to obtain the bearing capacity. Fourth step: The new acceleration, velocity, and displacement values calculated in step 1; the pressure distribution calculated in step 2; and the bearing capacity obtained in step 3 are used as the new initial conditions for the next calculation. By using this new set of parameter conditions, return to the first step to recalculate the changes in the bearing system. By using the aforementioned steps, the pressure variation and distribution inside the bearing system can be determined. The aforementioned steps can also be used to determine the relationship between the air film thickness and the position change at each time point, the change in the air film bearing capacity, and the dynamic behavior of the rotor. First step: After the time increment Δτ has passed, recalculate the acceleration, velocity, and displacement of the rotor. Moreover, define the new gap thickness value h at different positions between the rotor and the bearing.
Second step: Substitute the new h value into Equations (23) and (24). Then, recalculate the new pressure distribution state between the rotor and the bearing.
Third step: Integrate the pressure distribution calculated in step 2 to obtain the bearing capacity. Fourth step: The new acceleration, velocity, and displacement values calculated in step 1; the pressure distribution calculated in step 2; and the bearing capacity obtained in step 3 are used as the new initial conditions for the next calculation. By using this new set of parameter conditions, return to the first step to recalculate the changes in the bearing system. By using the aforementioned steps, the pressure variation and distribution inside the bearing system can be determined. The aforementioned steps can also be used to determine the relationship between the air film thickness and the position change at each time point, the change in the air film bearing capacity, and the dynamic behavior of the rotor.
In this study, to understand in detail the dynamic behavior of the DCAB system, the number of bearings (related to the rotor speed) and the rotor mass were used as bifurcation parameters and a bifurcation diagram was drawn for stability analysis. The dimensionless bearing number (Λ) is expressed as follows [15]: Different numerical methods were used in this study to ensure the applicability and accuracy of the proposed hybrid method.

Comparison of the Results Obtained with Different Numerical Methods
In this study, three numerical methods were used to solve the Reynolds equation of the DCAB system. The simulation results indicate that the hybrid scheme combining the DTM and FDM produced nearly identical results to the perturbation method and FDM. However, the rotor center displacements obtained with the perturbation method exhibited unstable numerical phenomena In this study, to understand in detail the dynamic behavior of the DCAB system, the number of bearings (related to the rotor speed) and the rotor mass were used as bifurcation parameters and a bifurcation diagram was drawn for stability analysis. The dimensionless bearing number (Λ) is expressed as follows [15]: Different numerical methods were used in this study to ensure the applicability and accuracy of the proposed hybrid method.

Comparison of the Results Obtained with Different Numerical Methods
In this study, three numerical methods were used to solve the Reynolds equation of the DCAB system. The simulation results indicate that the hybrid scheme combining the DTM and FDM produced nearly identical results to the perturbation method and FDM. However, the rotor center displacements obtained with the perturbation method exhibited unstable numerical phenomena under specific operating conditions. The hybrid method exhibited a higher accuracy than the other two methods. Table 2 presents a comparison of the rotor center displacements (X2, Y2) obtained with the three adopted numerical methods for different m r and Λ values.
To determine and verify the stability of the numerical results, the influence of different time intervals on the results were studied. The rotor center displacements with periodic motions under specific conditions were also compared according to the values of the Poincaré sections ( Table 3). As indicated by the numerical results, the hybrid method exhibited outstanding convergence and applicability for the DCAB system. The rotor trajectories for different operation parameters and time-step values ( H) can be calculated to the fifth decimal place precisely when using the hybrid scheme. The rotor trajectories obtained with the hybrid scheme were superior to those obtained with the perturbation method and FDM ( Table 2). The values of the Poincaré sections of the rotor center in the DCAB system for different time intervals are presented in Table 3. The results obtained with the hybrid method had an accuracy of calculation to the fifth right digit under different time increments. Thus, a small time step did not have to be used in this study to save calculation time and obtain high-precision calculation results for the DCAB system. Therefore, π/300 was used as the time interval in the subsequent dynamic analysis. The rotor mass affects the strength of the air flotation effect between the rotor and the bearing, which is one of the key factors that influence the stability of the entire air-bearing system. Therefore, the influence of the rotor mass on the DCAB system was analyzed. In this section, the rotor mass (m r ) is considered as the bifurcation parameter and the bearing number (Λ) of the system is considered to be 2.5. Figure 3 indicate that the rotor center (X 2 , Y 2 ) has a regular trajectory for a low mass (m r = 11.2 kg). When the mass is increased to 20.92 kg, the regular motion changes to nonperiodic motion. When the rotor mass is maintained below 20.92 kg, the system is stable. The nonperiodic phenomenon becomes less obvious for rotor masses of 23.13, 24.31, 24.4, 24.52, 24.58, and 24.60 kg. When m r is 24.75 kg, the system suddenly becomes nonperiodic.

Dynamic Trajectory and Phase Plane Analysis
is considered as the bifurcation parameter and the bearing number (Λ) of the system is considered to be 2.5. Figure 3 indicate that the rotor center (X2, Y2) has a regular trajectory for a low mass (mr = 11.2 kg). When the mass is increased to 20.92 kg, the regular motion changes to nonperiodic motion. When the rotor mass is maintained below 20.92 kg, the system is stable. The nonperiodic phenomenon becomes less obvious for rotor masses of 23.13, 24.31, 24.4, 24.52, 24.58, and 24.60 kg. When mr is 24.75 kg, the system suddenly becomes nonperiodic.  Figure 4 illustrates the dynamic spectral responses of the rotor along the horizontal and vertical directions of the coordinate system. The results displayed in the aforementioned figures reveal that the rotor center exhibits T-periodic motion when the rotor mass is 11.2 kg. The spectral responses  Figure 4 illustrates the dynamic spectral responses of the rotor along the horizontal and vertical directions of the coordinate system. The results displayed in the aforementioned figures reveal that the rotor center exhibits T-periodic motion when the rotor mass is 11.2 kg. The spectral responses depicted in Figure 4b indicate that motion of the rotor center becomes quasi-periodic along the horizontal and vertical directions when m r is 20.92 kg. When the rotor mass is increased to 23.13, 24.4, or 24.58 kg, the system motion becomes T-periodic. At m r values of 24.31, 24.52, and 24.60 kg, the system exhibits subharmonic 3T-periodic motion. Moreover, when the rotor mass is increased to 24.75 kg, the system displays chaotic motion.

Spectral Analysis
(i)  Figure 4 illustrates the dynamic spectral responses of the rotor along the horizontal and vertical directions of the coordinate system. The results displayed in the aforementioned figures reveal that the rotor center exhibits T-periodic motion when the rotor mass is 11.2 kg. The spectral responses depicted in Figure 4b indicate that motion of the rotor center becomes quasi-periodic along the horizontal and vertical directions when mr is 20.92 kg. When the rotor mass is increased to 23.13, 24.4, or 24.58 kg, the system motion becomes T-periodic. At mr values of 24.31, 24.52, and 24.60 kg, the system exhibits subharmonic 3T-periodic motion. Moreover, when the rotor mass is increased to 24.75 kg, the system displays chaotic motion.   Figure 5 displays the influence of the rotor mass mr on the DCAB system. The rotor masses are set between 0.1 and 25.0 kg in the operation of the DCAB system. Figure 4a,c indicates that the rotor in the system exhibits T-periodic motion in the horizontal and vertical directions when mr is less than 20.92 kg. This phenomenon is confirmed by analyzing the Poincaré section [ Figure 6a]. When the rotor mass is increased to 20.92 kg, the T-periodic motion is replaced by quasi-periodic motion, as displayed in Figure 5b,d. The aforementioned state of motion is illustrated by the closed curve formed by the discrete points in the Poincaré section [ Figure 6b]. The aforementioned motion is observed for an mr range of 20.92-23.13 kg. When the rotor mass is decreased to 20.13 kg, the system again exhibits T-periodic motion [ Figure 6c]. This periodic motion is observed in the mr range of 20.13-24.31 kg. When the rotor mass is 24.31 kg, the original T-periodic motion bifurcates into 3T-periodic motion, as depicted in Figure  6d, in which three discrete points can be observed on the map. The 3T-periodic motion is observed within an mr range of 24.31-24.4 kg. The 3T-periodic motion changes to T-periodic motion when the rotor mass is increased to 24.4 kg, as displayed in Figure 6e. T-periodic motion is observed within an mr range of 24.4-24.52 kg. The aforementioned description indicates that the T-and 3T-periodic motions occur alternately over an mr range of 24.52-24.75 kg, as depicted in Figure 6f,h.  Figure 5 displays the influence of the rotor mass m r on the DCAB system. The rotor masses are set between 0.1 and 25.0 kg in the operation of the DCAB system. Figure 4a,c indicates that the rotor in the system exhibits T-periodic motion in the horizontal and vertical directions when m r is less than 20.92 kg. This phenomenon is confirmed by analyzing the Poincaré section [ Figure 6a]. When the rotor mass is increased to 20.92 kg, the T-periodic motion is replaced by quasi-periodic motion, as displayed in Figure 5b,d. The aforementioned state of motion is illustrated by the closed curve formed by the discrete points in the Poincaré section [ Figure 6b]. The aforementioned motion is observed for an m r range of 20.92-23.13 kg. When the rotor mass is decreased to 20.13 kg, the system again exhibits T-periodic motion [ Figure 6c]. This periodic motion is observed in the m r range of 20.13-24.31 kg. When the rotor mass is 24.31 kg, the original T-periodic motion bifurcates into 3T-periodic motion, as depicted in Figure 6d, in which three discrete points can be observed on the map. The 3T-periodic motion is observed within an m r range of 24.31-24.4 kg. The 3T-periodic motion changes to T-periodic motion when the rotor mass is increased to 24.4 kg, as displayed in Figure 6e. T-periodic motion is observed within an m r range of 24.4-24.52 kg. The aforementioned description indicates that the T-and 3T-periodic motions occur alternately over an m r range of 24.52-24.75 kg, as depicted in Figure 6f    When the rotor mass is increased to 24.75 kg, the original stable motion becomes nonperiodic, as illustrated in Figure 6i. Moreover, a number of irregular discrete points occur on the map. Nonperiodic motion occurs in the mr range of 24.75-25.0 kg. The types of rotor center behaviors for different rotor masses are presented in Table 4. In this study, maximum Lyapunov exponents (MLEs) were used to identify and verify chaotic behavior. Figure 7a-h indicates that when mr is 11.2, 20.92, 23.13, 24.31, 24.4, 24.52, 24.58, or 24.60 kg, the MLEs approach or are smaller than 0. Thus, the system exhibits nonchaotic behavior for the aforementioned mr values. When mr is 24.75 kg, as depicted in Figure 7i, the values of the MLEs are greater than 0; thus, the DCAB system exhibits chaotic motion. The MLE distribution of the DCAB system is displayed in Figure 8. Figure 8 indicates that the blue curve, which indicates chaotic motion (the MLEs are greater than 0), occurs between 24.75 and 25.0 kg. The red curve, which represents the When the rotor mass is increased to 24.75 kg, the original stable motion becomes nonperiodic, as illustrated in Figure 6i. Moreover, a number of irregular discrete points occur on the map. Nonperiodic motion occurs in the m r range of 24.75-25.0 kg. The types of rotor center behaviors for different rotor masses are presented in Table 4. In this study, maximum Lyapunov exponents (MLEs) were used to identify and verify chaotic behavior. Figure 7a-h indicates that when m r is 11.2, 20.92, 23.13, 24.31, 24.4, 24.52, 24.58, or 24.60 kg, the MLEs approach or are smaller than 0. Thus, the system exhibits nonchaotic behavior for the aforementioned m r values. When m r is 24.75 kg, as depicted in Figure 7i, the values of the MLEs are greater than 0; thus, the DCAB system exhibits chaotic motion. The MLE distribution of the DCAB system is displayed in Figure 8. Figure 8 indicates that the blue curve, which indicates chaotic motion (the MLEs are greater than 0), occurs between 24.75 and 25.0 kg. The red curve, which represents the stable situation, matches the aforementioned results and indicates that the rotor exhibits complicated types of motions.

Dynamic Behavior of a DCAB System: Using the Bearing Number as the Bifurcation Parameter
The bearing number (Λ) directly influences not only the pressure distribution inside the bearing but also the rotational speed and stability of the DCAB system. In this section, the bearing number Λ is considered as the bifurcation parameter and the rotor mass mr is considered as 5.7 kg. The influence of Λ on the rotor dynamic behavior of the DCAB system is described in the following text. Figure 9 indicates that the rotor displays regular behavior at low bearing numbers (Λ = 2.1 or 3.21). When Λ is increased to 3.96, the rotor center begins to exhibit irregular and asymmetric motion. When the bearing number is increased to 3.98, the dynamic motion changes to a more regular periodic motion. Moreover, when Λ = 4.63, nonperiodic motion reappears. When the bearing number is increased to 5.02, the system returns to regular behavior. The aforementioned results indicate that distinctive changes occur in the behavior of the system as the bearing number increases marginally from 3.96 to 3.98. The rotor motion is nonperiodic when Λ = 3.96 and periodic when Λ = 3.98. Thus, the DCAB system is very sensitive to the bearing number.

Dynamic Behavior of a DCAB System: Using the Bearing Number as the Bifurcation Parameter
The bearing number (Λ) directly influences not only the pressure distribution inside the bearing but also the rotational speed and stability of the DCAB system. In this section, the bearing number Λ is considered as the bifurcation parameter and the rotor mass m r is considered as 5.7 kg. The influence of Λ on the rotor dynamic behavior of the DCAB system is described in the following text. Figure 9 indicates that the rotor displays regular behavior at low bearing numbers (Λ = 2.1 or 3.21). When Λ is increased to 3.96, the rotor center begins to exhibit irregular and asymmetric motion. When the bearing number is increased to 3.98, the dynamic motion changes to a more regular periodic motion. Moreover, when Λ = 4.63, nonperiodic motion reappears. When the bearing number is increased to 5.02, the system returns to regular behavior. The aforementioned results indicate that distinctive changes occur in the behavior of the system as the bearing number increases marginally from 3.96 to 3.98. The rotor motion is nonperiodic when Λ = 3.96 and periodic when Λ = 3.98. Thus, the DCAB system is very sensitive to the bearing number.

Dynamic Behavior of a DCAB System: Using the Bearing Number as the Bifurcation Parameter
The bearing number (Λ) directly influences not only the pressure distribution inside the bearing but also the rotational speed and stability of the DCAB system. In this section, the bearing number Λ is considered as the bifurcation parameter and the rotor mass mr is considered as 5.7 kg. The influence of Λ on the rotor dynamic behavior of the DCAB system is described in the following text. Figure 9 indicates that the rotor displays regular behavior at low bearing numbers (Λ = 2.1 or 3.21). When Λ is increased to 3.96, the rotor center begins to exhibit irregular and asymmetric motion. When the bearing number is increased to 3.98, the dynamic motion changes to a more regular periodic motion. Moreover, when Λ = 4.63, nonperiodic motion reappears. When the bearing number is increased to 5.02, the system returns to regular behavior. The aforementioned results indicate that distinctive changes occur in the behavior of the system as the bearing number increases marginally from 3.96 to 3.98. The rotor motion is nonperiodic when Λ = 3.96 and periodic when Λ = 3.98. Thus, the DCAB system is very sensitive to the bearing number.  Figure 10 illustrates the rotor center dynamic responses along the horizontal and vertical directions of the coordinate system for different bearing numbers. The rotor center displays Tperiodic motion when the bearing number is 2.1 or 5.02. However, the spectral response indicates that when Λ = 3.96 or 4.63, the rotor center displays nonperiodic motion. Moreover, at bearing numbers of 3.21 and 3.98, the system motion becomes 2T-subperiodic.

Bifurcation Analysis
The bearing number Λ is a key parameter for exploring the behavior of the DCAB system ( Figure  11). The bearing number is set between 1.0 and 6.0, and this range complies with real operating conditions. Initially, when Λ = 2.1, the rotor center performs T-periodic motion in the horizontal and

Bifurcation Analysis
The bearing number Λ is a key parameter for exploring the behavior of the DCAB system ( Figure  11). The bearing number is set between 1.0 and 6.0, and this range complies with real operating conditions. Initially, when Λ = 2.1, the rotor center performs T-periodic motion in the horizontal and

Bifurcation Analysis
The bearing number Λ is a key parameter for exploring the behavior of the DCAB system ( Figure 11). The bearing number is set between 1.0 and 6.0, and this range complies with real operating conditions. Initially, when Λ = 2.1, the rotor center performs T-periodic motion in the horizontal and vertical directions. This phenomenon is revealed and verified by the Poincaré section depicted in Figure 12a, which only has a single point. The system exhibits T-periodic motion when the bearing number is between 1.0 and 3.21. When the bearing number is increased to 3.21, the system motion bifurcates into 2T-periodic motion. Figure 12b indicates that the system generates two discrete points, and the 2T-periodic motion continues in the Λ range of 3.21-3.96. When the bearing number is between 3.96 and 3.98, the system exhibits chaotic motion (see Figure 12c), in which irregular discrete points are observed). When Λ is equal to 3.98, the chaotic motion changes to 2T-periodic motion (Figure 12d). The system exhibits 2T-periodic motion when 3.98 ≤ Λ < 4.63. It displays chaotic motion again when 4.63 ≤ Λ < 5.02, as displayed in Figures 11 and 12e This unstable motion changes to stable T-periodic motion when Λ = 5.02. This stable motion is maintained until Λ = 6.0, as displayed in Figure 12f. The relationship between the rotor center behavior and Λ is complicated. All the system motions that occur when 1.0 ≤ Λ ≤ 6.0 are presented in Table 5.
Symmetry 2020, 12, x FOR PEER REVIEW 19 of 23 vertical directions. This phenomenon is revealed and verified by the Poincaré section depicted in Figure 12a, which only has a single point. The system exhibits T-periodic motion when the bearing number is between 1.0 and 3.21. When the bearing number is increased to 3.21, the system motion bifurcates into 2T-periodic motion. Figure 12b indicates that the system generates two discrete points, and the 2T-periodic motion continues in the Λ range of 3.21-3.96. When the bearing number is between 3.96 and 3.98, the system exhibits chaotic motion (see Figure 12c), in which irregular discrete points are observed). When Λ is equal to 3.98, the chaotic motion changes to 2T-periodic motion (Figure 12d). The system exhibits 2T-periodic motion when 3.98 ≤ Λ < 4.63. It displays chaotic motion again when 4.63 ≤ Λ < 5.02, as displayed in Figures 11 and 12e This unstable motion changes to stable T-periodic motion when Λ = 5.02. This stable motion is maintained until Λ = 6.0, as displayed in Figure 12f. The relationship between the rotor center behavior and Λ is complicated. All the system motions that occur when 1.0 ≤ Λ ≤ 6.0 are presented in Table 5.  vertical directions. This phenomenon is revealed and verified by the Poincaré section depicted in Figure 12a, which only has a single point. The system exhibits T-periodic motion when the bearing number is between 1.0 and 3.21. When the bearing number is increased to 3.21, the system motion bifurcates into 2T-periodic motion. Figure 12b indicates that the system generates two discrete points, and the 2T-periodic motion continues in the Λ range of 3.21-3.96. When the bearing number is between 3.96 and 3.98, the system exhibits chaotic motion (see Figure 12c), in which irregular discrete points are observed). When Λ is equal to 3.98, the chaotic motion changes to 2T-periodic motion (Figure 12d). The system exhibits 2T-periodic motion when 3.98 ≤ Λ < 4.63. It displays chaotic motion again when 4.63 ≤ Λ < 5.02, as displayed in Figures 11 and 12e This unstable motion changes to stable T-periodic motion when Λ = 5.02. This stable motion is maintained until Λ = 6.0, as displayed in Figure 12f. The relationship between the rotor center behavior and Λ is complicated. All the system motions that occur when 1.0 ≤ Λ ≤ 6.0 are presented in Table 5.  The MLEs were also used to determine and verify if chaotic behavior is affected by changes in the bearing number. Figure 13a,b,d,f indicates that the MLEs are smaller than 0 when Λ = 2.1, 3.21, 3.98, and 5.02; thus, the system exhibits nonchaotic motion for the aforementioned Λ values. Figure  13c,e indicates that the MLEs are greater than 0 when Λ = 3.96 and 4.63; thus, the system exhibits chaotic motion for the aforementioned Λ values. The variations in the MLEs with the bearing number are illustrated in Figure 14. The results in Figure 14 indicate that the chaotic regions (where the MLEs are greater than 0) occur when 3.96 ≤ Λ < 3.98 and 4.63 ≤ Λ < 5.02.   The MLEs were also used to determine and verify if chaotic behavior is affected by changes in the bearing number. Figure 13a,b,d,f indicates that the MLEs are smaller than 0 when Λ = 2.1, 3.21, 3.98, and 5.02; thus, the system exhibits nonchaotic motion for the aforementioned Λ values. Figure 13c,e indicates that the MLEs are greater than 0 when Λ = 3.96 and 4.63; thus, the system exhibits chaotic motion for the aforementioned Λ values. The variations in the MLEs with the bearing number are illustrated in Figure 14. The results in Figure 14 indicate that the chaotic regions (where the MLEs are greater than 0) occur when 3.96 ≤ Λ < 3.98 and 4.63 ≤ Λ < 5.02.  The MLEs were also used to determine and verify if chaotic behavior is affected by changes in the bearing number. Figure 13a,b,d,f indicates that the MLEs are smaller than 0 when Λ = 2.1, 3.21, 3.98, and 5.02; thus, the system exhibits nonchaotic motion for the aforementioned Λ values. Figure  13c,e indicates that the MLEs are greater than 0 when Λ = 3.96 and 4.63; thus, the system exhibits chaotic motion for the aforementioned Λ values. The variations in the MLEs with the bearing number are illustrated in Figure 14. The results in Figure 14 indicate that the chaotic regions (where the MLEs are greater than 0) occur when 3.96 ≤ Λ < 3.98 and 4.63 ≤ Λ < 5.02.

Establishment of the Steady and Nonsteady Operation Regions of the DCAB System
To determine if the behavior of the DCAB system is stable under various operating parameters, the rotor mass and bearing number were used as bifurcation parameters in the analysis conducted in this study. These two parameters influence the stability of the DCAB system. Moreover, their combined effect induces chaotic rotor behavior that cannot be predicted. Therefore, MLEs were used to confirm the stable and unstable operation regions for the DCAB system for different rotor masses and bearing numbers ( Figure 15). In Figure 15, the light pink to deep red areas represent unstable regions, for which the MLEs are greater than 1.0. The light blue to deep blue areas represent stable regions, for which the MLEs are less than 0. The results in Figure 15 indicate that the unstable (chaotic) regions are concentrated in areas with a high rotor mass. The aforementioned results can be effectively used to identify the stable and unstable regions for avoiding chaotic motion in a DCAB system. These results also provide guidance for the design of bearing systems with high stability and performance.

Conclusions
The aims of this study were to analyze the various dynamic behaviors of the flexible rotor in the DCAB system under different operating parameters and to identify the stable and chaotic operation

Establishment of the Steady and Nonsteady Operation Regions of the DCAB System
To determine if the behavior of the DCAB system is stable under various operating parameters, the rotor mass and bearing number were used as bifurcation parameters in the analysis conducted in this study. These two parameters influence the stability of the DCAB system. Moreover, their combined effect induces chaotic rotor behavior that cannot be predicted. Therefore, MLEs were used to confirm the stable and unstable operation regions for the DCAB system for different rotor masses and bearing numbers ( Figure 15). In Figure 15, the light pink to deep red areas represent unstable regions, for which the MLEs are greater than 1.0. The light blue to deep blue areas represent stable regions, for which the MLEs are less than 0. The results in Figure 15 indicate that the unstable (chaotic) regions are concentrated in areas with a high rotor mass. The aforementioned results can be effectively used to identify the stable and unstable regions for avoiding chaotic motion in a DCAB system. These results also provide guidance for the design of bearing systems with high stability and performance.

Establishment of the Steady and Nonsteady Operation Regions of the DCAB System
To determine if the behavior of the DCAB system is stable under various operating parameters, the rotor mass and bearing number were used as bifurcation parameters in the analysis conducted in this study. These two parameters influence the stability of the DCAB system. Moreover, their combined effect induces chaotic rotor behavior that cannot be predicted. Therefore, MLEs were used to confirm the stable and unstable operation regions for the DCAB system for different rotor masses and bearing numbers ( Figure 15). In Figure 15, the light pink to deep red areas represent unstable regions, for which the MLEs are greater than 1.0. The light blue to deep blue areas represent stable regions, for which the MLEs are less than 0. The results in Figure 15 indicate that the unstable (chaotic) regions are concentrated in areas with a high rotor mass. The aforementioned results can be effectively used to identify the stable and unstable regions for avoiding chaotic motion in a DCAB system. These results also provide guidance for the design of bearing systems with high stability and performance.

Conclusions
The aims of this study were to analyze the various dynamic behaviors of the flexible rotor in the DCAB system under different operating parameters and to identify the stable and chaotic operation

Conclusions
The aims of this study were to analyze the various dynamic behaviors of the flexible rotor in the DCAB system under different operating parameters and to identify the stable and chaotic operation areas of the bearing system. The FDM, perturbation method, and hybrid method were used in the numerical analysis to determine the gas film pressure distribution and gas film thickness function of the Reynolds equation. The flexible rotor dynamic balance equation was then used to solve the trajectory displacement of the rotor center. Vibration signal analysis was performed by analyzing the spectral response generated by the displacement of the rotor center, Poincaré section diagram, bifurcation diagram, and MLEs. The numerical analysis results indicate that the hybrid method had a higher accuracy and precision than the FDM and perturbation method. Furthermore, a fairly accurate numerical solution could be obtained without using a small time step. The vibration analysis results indicate that the trajectory of the rotor center changes significantly with changes in the rotor mass and the bearing number. Changes in the aforementioned parameters produce various vibration behaviors in the horizontal and vertical directions, including periodic, subharmonic, quasi-periodic and chaotic motions. In particular, when the rotor mass is high, the bearing system exhibits unstable nonlinear chaotic behavior. Moreover, changes in the bearing number also affect the trajectory of the rotor. Chaotic behavior occurs when 3.96 ≤ Λ < 3.98 and 4.63 ≤ Λ < 5.02 for a rotor mass of 4.2 kg. When the bearing number does not fall within the aforementioned ranges, the rotor exhibits relatively stable motion. The major contribution of this study is its finding that the rotor mass and bearing number influence the stability of DCAB system. Furthermore, this study confirmed that the rotor exhibits chaotic behavior for specific ranges of the rotor mass and bearing number. The MLE is a key index for defining stable and unstable regions (nonchaotic and chaotic regions, respectively) for the DCAB system over various values of the rotor mass and bearing number. The research results can be used as an analysis reference and important basis for avoiding unstable nonlinear behaviors in DCAB system design and application.