Abstract
In recent years, spiral-grooved air bearing systems have attracted much attention and are especially useful in precision instruments and machines with spindles that rotate at high speed. Load support can be multidirectional and this type of bearing can also be very rigid. Studies show that some of the design problems encountered are dynamic and include critical speed, nonlinearity, gas film pressure, unbalanced rotors, and even poor design, all of which can result in the generation of chaotic aperiodic motion and instability under certain conditions. Such irregular motion on a large scale can cause severe damage to a machine or instrument. Therefore, understanding the conditions under which aperiodic behaviour and vibration arise is crucial for prevention. In this study, numerical analysis, including the Finite Difference and Differential Transformation Methods, is used to study these effects in detail in a front opposed-hemispherical spiral-grooved air bearing system. It was found that different rotor masses and bearing number could cause undesirable behaviour including periodic, subperiodic, quasi-periodic, and chaotic motion. The results obtained in this study can be used as a basis for future bearing system design and the prevention of instability.
1. Literature Review
The mathematical theory of gas lubrication was first introduced by Reynolds [1] in 1886 and he derived a partial differential equation related to pressure, density, relative motion, and speed, the well-known Reynolds equation. This established the foundation of fluid lubrication theory. In 1958, Whipple [2] was the first to study the properties of an air bearing with herringbone grooves. In 1959, Whitley and Williams [3] studied the model proposed by Whipple and also published a great number of methods for the design of the herringbone groove. In 1963, Vohr and Pan [4] studied air bearings with spiral grooves and derived a differential equation for the gas film pressure based on an assumption of the number of spiral grooves close to infinity and this became the “Narrow Groove Theory.” Our survey of the literature revealed that the dynamic characteristics, particularly the gas film pressure, gas film rigidity, and the damping effect of the gas film can have a significant impact on the entire bearing supporting system. In 1994, Bonneau and Absi [5] used the finite element method to study characteristics related to the herringbone groove, including stability analysis, pressure analysis, and the impact of several important parameters on the system. Recently, Hirayama et al. [6] discussed the optimization of groove dimensions in herringbone grooved journal bearings for the design of precision spindles with improved run-out characteristics. Schiffmann and Favrat [7] studied the comparison of the influences on the real gas and the perfect gas for herringbone grooved gas journal bearing system and suggested bearing designers to further consider the real gas effects. In 2011, Chen et al. [8] analyzed the groove-ridge discontinuity effect for herringbone grooved gas journal bearings and obtained the pressure distribution of the fluid film on the herringbone grooves. In 2013, Schiffmann [9] optimized the groove geometry of the gas bearing system and set up the design guidelines for enhanced herringbone grooved gas journal bearings. For the transient analysis for the gas bearing system, Hassini and Arghir [10, 11] proposed a novel simplified method for simulating large nonlinear displacements and stability analysis in gas-lubricated bearings. However, none of the references provided an analysis on the dynamic characteristics of a bearing. Furthermore, no series of complete and effective analytic methods for the motion behaviour of spiral groove air bearings and rotors could be found.
In 2013, Grigor’ev and Smirnov [12] used finite element method for calculating the characteristics of gas-lubricated spiral-grooved bearings based on integrating the Reynolds equation and the thrust bearing is also considered as a specific example. The results on the dependence of the carrying capacity on the compressibility number and the number of grooves are presented and compared with the asymptotic Narrow Groove Theory. With regard to studies of the dynamic behaviour of rotors, Bently [13] discovered in 1974 that a rotary oil film bearing system could exhibit two- or even three-stage subharmonic vibration and Child et al. [14, 15] used an analytic method to prove the existence of subharmonic vibration in a rotor-bearing system. However, these investigations were all directed towards oil film bearings, and there have been very few studies of rotor behaviour in an air bearing system. From 2007 to 2011, Wang [16–18] used the Finite Difference Method to investigate gas film pressure in a herringbone air bearing system, and he also used orbit graphs, spectrograms, and bifurcation diagrams to analyze the dynamic behaviour of the centres of the rotor and the journal. He learned that, under different operational conditions, the centres of the rotor and the journal demonstrated periodic, quasi-periodic, and subharmonic motion, and he also discovered that the system could generate nonlinear chaotic motion. These results proved useful for this present study. In 2015, Guangwei et al. [19, 20] studied the spiral-grooved opposed-hemisphere gas bearing system with a rigid rotor and focused particular attention on its whirl motion. They used finite element method combined with the Finite Difference Method to solve the time-dependent Reynolds equation and analyze the complicated dynamic behaviour of the rotor-bearing system by phase portraits, power spectra, Poincare maps, and bifurcation diagrams obtained from the numerical procedure. The results revealed that the conical whirl instability appears earlier than the cylindrical whirl instability with increasing rotational speed for the rotor-bearing system with no unbalance mass.
In addition, a spiral-grooved air bearing system is typically nonlinear. Conventionally, the steady-state condition of a deterministic nonlinear dynamic system can be regarded as being in one of three main states: equilibrium, periodic motion, and quasi-periodic motion. These may collectively be described as the “attractor” because, after a transient-state response, a system in a steady state can be attracted towards another state. However, recent research shows that a nonlinear dynamic system can exhibit not only regular behaviour but also random or chaotic motion. In a situation where a tiny uncertainty exists in the initial state of motion, the associated dynamic behaviour cannot easily be predicted. The geometric shape formed in the phase space may also differ significantly from that of a normal qualitative system. When chaotic behaviour arises it is very difficult to predict the outcome which may be totally unexpected and even include damage to the system. There is a definite need for a means to prevent the occurrence of such phenomena and the objective of this study is the analysis of the characteristics of a spiral-grooved air bearing and an evaluation of its dynamic behaviour in response to different operational criteria. The bifurcation characteristics associated with nonlinear behaviour generated by the rotor in the system are also discussed. Based on such analysis and study, determinations on whether the chaotic phenomena occur in the system and accurate prediction on the dynamic orbits of the system can be achieved.
2. Theoretical Analysis and Model Construction
2.1. System Governing Equation
Front opposed-hemispherical spiral-grooved air bearing (FOSAB) systems are used extensively in precision instruments, such as gyroscopes, and in other aerospace applications and are suitable for multidirectional load bearing. They can handle heavy loads making them extremely suitable for many practical engineering applications. In this study two different numerical analytic methods have been used to study a FOSAB system and the performance and dynamic characteristics of the bearing were investigated in some detail. The diagrams in Figures 1 and 2 show the structure of a FOSAB in which the centre and the centre of the sphere coincide. The -axis is vertical and downwards and the corresponding spherical coordinate is (, , ).
(a)
(b)
(c)
As shown in the Figures 1 and 2, the spherical radius of the bearing is , the average clearance is , and the rotational speed is . Parameters of the spiral groove include angle of groove , depth of groove , width of groove , and separation between grooves . The compressible Reynolds equation of the bearing of such type uses the spherical coordinate system in principle, which is modified as follows:wherein
is the dimensionless gas film pressure, is the dimensionless gas film thickness, is the gas film density, and is the gas film viscosity. To solve the bearing gas pressure distribution the traditional Finite Difference Method and a hybrid method [21, 22] were used. The hybrid method is a combination of the Differential Transformation Method (DTM) and the Finite Difference Method (FDM). For the nonlinear Reynolds equation time domain, the Differential Transformation Method was first used to discretize the time into intervals and the central difference scheme of the Finite Difference Method was then applied to the coordinates of the locations. The Differential Transformation Method was used to solve the nonlinear Reynolds equation, and the basic principle is as follows.
The differential transformation of the function of is defined as
Wherein refers to the transformation function in the transformation domain, refers to the original function in the time domain, refers to the time separation, and refers to the conversion parameter. The reverse equation of is
Calibration is needed to ensure precision during the calculation; therefore, when the second time interval is used, (14) and (15) are used to perform calibration.
If then the above equation becomes
If the calculated values on each side of the equals’ signs in (14) and (15) vary greatly, precision is insufficient. Therefore, the time interval needs to be reduced and recalculation is necessary.
The basic arithmetic operation for the differential transformation is as shown in Table 1.
The differential transformation of the governing equation and the boundary condition are obtained and the pressure function theorem domain is divided into number of subintervals; the pressure solution for each subinternal is expressed as an inverse transformation equation. If the theorem domains of the subintervals are , then the expression for the function in the first interval, based on (13), will be
If the initial value of is known, then, from the above equation, can be obtained. From , other discrete values of , in this interval can be obtained. In the first subinterval, the destination value of the function of refers to the initial value of the function in the second subinterval. The function in the second subinterval is
From , other , can also be obtained. In such order, the value of 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:
2.2. Numerical Analysis
The traditional Finite Difference Method is used first on the Reynolds equation and to apply the central difference scheme to the spatial coordinates. For the time domain, the implicit backwards difference method is used. Equation (1) is transformed into a differential equation, and the Finite Difference Method of Successive Over Relation (SOR) is used to reduce the five unknown numbers to three, following which iteration is used to obtain the correct pressure value at each time interval. Iteration at each time interval is performed in the direction only, and the three points of , , are unknown numbers. The two points of and that are originally unknown are replaced by the previous iteration values. Consequently, all the equations can be formed into a triple-diagonal matrix, which can be then solved by Gaussian Elimination.
After obtaining the pressure distribution between the bearing and the rotor, the bearing force of the gas film can be obtained by integration. The numerical calculations of the hybrid method can be done using (12)–(18), as shown in the following:
Next, the Differential Transformation Method is applied to the time domain in (19) to produce discrete time intervals, and the central difference scheme of the Finite Difference Method is then applied to the coordinates of the locations. The following can then be derived:
2.3. Dynamic Equations of the Rotor
In this system, the ends of a flexible rotor of mass are supported on two identical and symmetrical air bearings secured on a rigid base. Under such ideal conditions the motion of the rotor can be treated as two having two directions of movement. refers to the geometric centre of the rotor, to the geometric centre of the journal, to the eccentricity of the rotor, to the rigidity of the rotating axle, and to the rotational speed of the rotor.
In a steady state, the equation of motion of the rotor centre is as follows:
The force equilibrium equation at the journal centre is as follows:
Upon dimensionless transformation:
In addition, the dimensionless parameters are introduced as follows:
By substituting (25) into (21)~(24) and introducing the dimensionless parameters (26), the following equations can be obtained:
In this study iteration was used to calculate the motion behaviour of the rotor. First, the acceleration, speed, and displacement were calculated, step by step, with respect to time. The initial state of the system was taken as static, and the initial displacement (, ) of the rotor refers to the location of the rotor centre when the system was in static equilibrium; such equilibrium position also defines the initial gas film thickness value . In addition, the initial speed was set to zero.
The steps of calculation are as follows.
Step 1. At time , the initial condition is used; after a period of time , (31)~(36) can be used to calculate the new acceleration, speed, and displacement of the rotor:
Step 2. The displacement of the rotor centre obtained from Step 1 can be used to calculate the new displacement of the journal centre. The rotation of the journal causes changes in the gas film thickness ; consequently, after the new displacement value yields a new , this is substituted into (2) to calculate a new pressure value, and the gas film bearing force can be derived by integration.
Step 3. From the displacement and velocity values obtained in Step 1 and the pressure distribution from Step 2, a new set of initial conditions can be derived from the gas film force obtained. When the time is increased by , (i.e., ), the aforementioned new initial conditions are used again to further calculate all the variations in the system. 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, the dynamic orbits of the rotor, and so on 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 and mass of the rotor. 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 accuracy of the numerical analysis, two different numerical methods were used (as previously mentioned) to ensure that the hybrid method proposed in this paper was applicable and to increase the calculation precision.
3. Results and Discussion
3.1. Comparison between Results Generated by Different Numerical Methods
The results obtained by the hybrid method that combines both the Differential Transformation Method (DTM) and the Finite Difference Method (FDM) agree with the results of the traditional Finite Difference Method and Perturbation Method. This was confirmed by use of the two different methods as previously mentioned. The values of the Finite Difference Method and Perturbation Method obtained under certain parameters would cause unstable phenomena; in addition, the result also shows that the hybrid method (DTM and FDM) yields greater accuracy. Table 2 shows a comparison of the values of the trace of the rotor centre (, ), wherein refers to the time step for calculation. With regard to the analysis on the stability of the values, the impact of different time step on the values is also completed, showing that the stability of the values obtained from the hybrid method is indeed greater than those obtained from the traditional Finite Difference Method. Furthermore, precision of the displacement calculation can be accurate to five decimal places by hybrid method.
The results show that the hybrid method presented here demonstrates excellent convergence and applicability. The method can be used to effectively calculate the orbits of the system with different rotor mass and increased bearing number (rotational speed). Examination of the test results in Table 3 shows that the time step selected does not need to be highly precise to obtain sufficiently accurate numerical results and this further reduces the bifurcation characteristics of subsequent calculations. Therefore, is used as the time step for dynamic analysis calculation.
3.2. Analysis of the Dynamic Behaviour of a FOSAB System: Using Rotor Mass as a Bifurcation Parameter
In this section, analysis of the impact of rotor mass on the air bearing system is conducted; here the bearing number of the FOSAB system is set to , and the rotor mass is used as the bifurcation parameter.
3.2.1. Dynamic Orbits Analysis
An examination of Figures 3(a1, b1), 3(a2, b2), …, and 3(a7, b7) show that the orbits of the centres of the rotor (, ) and the journal (, ) will follow regular paths when the mass is small ( = 5.1 kg). However, when the mass is increased to = 7.06 kg, aperiodic motion occurs. This shows that, with a rotor mass of less than 7.06 kg, the system will remain relatively stable. When the rotor mass reaches = 8.06 kg, aperiodic motion is reduced but comes back when = 8.21 kg. However, the nonperiodicity stabilizes when = 8.52 kg. This stable behaviour does not last long, and the rotor centre becomes unstable again when its mass reaches 9.30 kg. Finally, at = 10.76 kg, the system becomes relatively stable.
3.2.2. Spectrum Analysis
Figures 3(c1, d1), 3(c2, d2), …, and 3(c7, d7) show the dynamic response of the rotor centre in the horizontal and vertical directions of rotation. The results show that when the rotor mass is = 5.1 and 10.76 kg, motion of the rotor centre is periodic; when = 7.06 and 9.30 kg, the spectrum response graphs (Figures 3(c2, d2), 3(c6, d6)) show that the motion of the rotor centre becomes chaotic. When the rotor mass is 8.06 and 8.52 kg, the system generates 2 subharmonic motion. However, motion of the system is quasi-periodic when the rotor mass is 8.21 kg.
3.2.3. Bifurcation Analysis
The effect of different rotor mass on the system was studied using rotor mass as an analytic parameter, as shown in Figure 4. The mass was set within a range of 1.0 to 12.0 kg according to actual operating conditions. Figure 5 shows that when < 7.06 kg, the rotor centre undergoes periodic motion in both the horizontal and vertical directions, and this is illustrated by the Poincaré map shown in Figure 5(a). Nevertheless, when the mass was increased to = 7.06 kg, the periodic motion of the rotor centre was replaced by chaotic motion, as shown in Figure 5(b). When the mass reached = 8.06 kg, the system returned to 2 subperiodic motion (see Figure 5(c)). Although this motion continued over the range 8.06 ≤ < 8.21 kg it did not last to the end. When increased to 8.21 kg, quasi-periodic motion appeared (see Figure 5(d)) and continued over the range 8.21 ≤ < 8.52 kg. This is illustrated by the closed curve formed by multiple discrete points shown on the Poincaré map in Figure 5(d). As the mass increased to = 8.52 kg, the quasi-periodic motion of the original system bifurcated into 2 subharmonic motion, as shown in Figure 5(e) where the map shows two discrete points and this type of motion is present throughout the range of 8.52 ≤ < 9.30 kg. Next, as the mass increased further to = 9.30 kg, the subharmonic motion of the original system became chaotic motion, as shown in Figure 5(f). Here the map generates two attraction subregions, and this is present throughout the range of 9.30 ≤ < 10.76 kg. Finally, when the rotor mass is increased to = 10.76 kg, the chaotic motion of the system has become a periodic motion that is relatively stable, as shown in Figure 5(g).
(a)
(b)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
The maximum Lyapunov exponents have been very useful for the analysis and verification of the chaotic behaviour of the system. From Figures 6(a), 6(c), 6(d), 6(e), and 6(g), it can be seen that when = 5.1, 8.06, 8.21, 8.52 and 10.76 kg, the maximum Lyapunov exponents are close to zero, and the system is nonchaotic. However, when = 7.06 and 9.30, as shown in Figures 6(b) and 6(f), the exponents are greater than zero and the system state of motion is indeed chaotic. A full view of the maximum Lyapunov exponents of the system is shown in Figure 7, and the result indicates that the areas of chaotic behaviour (exponent greater than zero) occur at two intervals, 7.06 ≤ < 8.06 kg and 9.30 ≤ < 10.76 kg. This is in clear agreement with the results. However, these results show that the motion of the rotor can be a relatively complicated model. A summary of the different kinds of system motion, with rotors of different mass, is given in Table 4.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
3.3. Analysis on Dynamic Behaviour of FOSAB System: Using the Bearing Number as the Bifurcation Parameter
For an air bearing system, the rotational speed referred to that bearing number of the rotor-bearing has a direct effect on pressure distribution inside the bearing, which in turn has a direct effect on the relevant performance and stability of the entire system. Accordingly, in this section, the bearing number is used as the bifurcation parameter, and the rotor mass is set to = 2.3 kg to analyze and study the dynamic behaviour of the system with respect to bearing number.
3.3.1. Dynamic Orbits Analysis
An examination of Figures 8(a1, b1), 8(a2, b2), …, and 8(a6, b6) shows that when the bearing number is small, = 2.3, 3.12, the behaviour of the journal and the rotor will be relatively regular. However, when is increased to 3.84, behaviour of the journal and the rotor centre will become irregular and asymmetrical. When the bearing number is further increased to 3.86, the dynamic orbits return to relatively regular periodic motion. However, irregular motion reappears when and only becomes regular again when the bearing number is increased to 4.89.
3.3.2. Spectrum Analysis
Figures 8(c1, d1), 8(c2, d2), …, and 8(c6, d6) show the dynamic responses of the rotor centre in the horizontal and vertical directions of rotation with different bearing numbers. The study shows that at bearing numbers of = 2.3, 3.86, and 4.89 the rotor centre demonstrates periodic motions. When = 3.84, the spectrum response diagram shows aperiodic motion of the rotor centre, but when the bearing number is = 3.12, the motion becomes 2 subharmonic motion. However, when = 4.55, the rotor centre motion becomes aperiodic.
3.3.3. Bifurcation Analysis
As shown in Figure 9, the bearing number can be used as a primary analytic parameter for studying the behaviour of an air bearing system and has been set within a range of 1.0 to 5.5 according to actual operational conditions. When the bearing number is small, such as , the rotor centre exhibits periodic motion in both the horizontal and rotational directions and this can be verified by the appearance of a single dot on the Poincaré map as shown in Figure 10(a). When the bearing number is increased to , the system will bifurcate and generate 2 periodic motion, clearly shown by the presence of two discrete dots on the map as shown in Figure 10(b). When the bearing number increases to , the system motion becomes chaotic, clearly shown by the irregular discrete points generated on the map in Figure 10(c). When , the chaotic motion becomes periodic motion, as shown in Figure 10(d), and the same occurs in the interval of ; see Figure 9. In addition, for the interval of , the system generates quasi-periodic motion, as shown in Figure 10(e). When , the system assumes -periodic motion that is relatively stable and this continues up to ; see Figure 10(f).
(a)
(b)
(a)
(b)
(c)
(d)
(e)
(f)
With regard to the question about whether a change of bearing number will cause the generation of chaotic behaviour, in this paper, similarly, the maximum Lyapunov exponents are also used for determination and further verification. From Figures 11(a), 11(b), 11(d), 11(e), and 11(f), it can be seen that when = 2.3, 3.12, 3.86, 4.55, and 4.89, the exponents obtained can all be smaller than or equal to zero and the system motion is nonchaotic. In view of the complexity of the relationship between the motion of the rotor centre and bearing number, all motions and behaviours generated with respect to the interval of are summarized in Table 5 and a full view of the maximum Lyapunov exponents for the system is shown in Figure 12. The result indicates that chaotic regions (exponents greater than zero) only occur when , which agrees with the results shown in Table 5.
(a)
(b)
(c)
(d)
(e)
(f)
4. Conclusion
The objective of this study was an analysis of the dynamic behaviour of a flexible rotor supported by a front opposed-hemispherical spiral-grooved air bearing system. The Finite Difference Method, Perturbation Method, and a hybrid method were used to solve the pressure distribution at the highest nonlinearity in the system, after which dynamic equations of the flexible rotor centre were used to obtain the orbits of the rotor centre. 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 centre change along with the mass and the bearing number and synchronously generate complicated motion in the horizontal and vertical directions and include normal and subharmonic vibration, as well as periodic, quasi-periodic, and chaotic motion. The hybrid method applied in this paper yields better accuracy and precision for the verification of values than the Finite Difference Method 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 behaviour of the rotor revealed that the rotor orbits demonstrate unstable nonlinear chaotic behaviour 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 behaviour in the interval 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 bearing system design and the prevention of instability.
Competing Interests
The author declares that they have no competing interests.
Acknowledgments
The author gratefully acknowledges the financial support provided to this study by the Ministry of Science and Technology of Taiwan under Grant nos. MOST-105-2221-E-167-010 and MOST-105-2622-E-167-007-CC3.