Optimizing energy dissipation in gas foil bearings to eliminate bifurcations of limit cycles in unbalanced rotor systems

High-speed rotor systems mounted on gas foil bearings present bifurcations which change the quality of stability, and may compromise the operability of rotating systems, or increase noise level when amplitude of specific harmonics drastically increases. The paper identifies the dissipating work in the gas film to be the source of self-excited motions driving the rotor whirling close to bearing’s surface. The energy flow among the components of a rotor gas foil bearing system with unbalance is evaluated for various design sets of bump foil properties, rotor stiffness and unbalance magnitude. The paper presents a methodology to retain the dissipating work at positive values during the periodic limit cycle motions caused by unbalance. An optimization technique is embedded in the pseudo-arc length continuation of limit cycles, those evaluated (when exist) utilizing an orthogonal collocation method. The optimization scheme considers the bump foil stiffness and damping as the variables for which bifurcations do not appear in a certain speed range. It is found that secondary Hopf (Neimark–Sacker) bifurcations, which trigger large limit cycle motions, do not exist in the unbalanced rotors when bump foil properties follow the optimization pattern. Period-doubling (flip) bifurcations are possible to occur, without driving the rotor in high response amplitude. Different design sets of rotor stiffness and unbalance magnitude are investigated for the efficiency of the method to eliminate bifurcations. The quality of the optimization pattern allows optimization in real time, and gas foil bearing properties shift values during operation, eliminating bifurcations and allowing operation at higher speed margins. Compliant bump foil is found to enhance the stability of the system.

Abstract High-speed rotor systems mounted on gas foil bearings present bifurcations which change the quality of stability, and may compromise the operability of rotating systems, or increase noise level when amplitude of specific harmonics drastically increases. The paper identifies the dissipating work in the gas film to be the source of self-excited motions driving the rotor whirling close to bearing's surface. The energy flow among the components of a rotor gas foil bearing system with unbalance is evaluated for various design sets of bump foil properties, rotor stiffness and unbalance magnitude. The paper presents a methodology to retain the dissipating work at positive values during the periodic limit cycle motions caused by unbalance. An optimization technique is embedded in the pseudo-arc length continuation of limit cycles, those evaluated (when exist) utilizing an orthogonal collocation method. The optimization scheme considers the bump foil stiffness and damping as the variables for which bifurcations do not appear in a certain speed range. It is found that secondary Hopf (Neimark-Sacker) bifurcations, which trigger large limit cycle motions, do not exist in the unbalanced rotors when bump foil properties follow the optimization pattern. Period-doubling (flip) bifurcations are possible to occur, without driving the rotor in high response amplitude. Different design sets of rotor stiffness and unbalance magnitude are investigated for the efficiency of the method to eliminate bifurcations. The quality of the optimization pattern allows optimization in real time, and gas foil bearing properties shift values during operation, eliminating bifurcations and allowing operation at higher speed margins. Compliant bump foil is found to enhance the stability of the system.
Keywords Gas foil bearing Á Nonlinear dynamics Á Bifurcation elimination Á Energy dissipation Á Design optimization Á Numerical continuation

Abbreviations
Hellenic letters a f Dimensionless compliance of the bump foil per area, Gas Foil Bearings (GFBs) are bearing elements which support rotating shafts in high-speed applications and they have been evolving for more than half a century. The advantages of high reliability, simplicity, and environmental friendly characteristics (oil-free technology) make a rapidly evolving technology. The operating principle lies on the elasto-aerodynamic flow (EAD lubrication) where the rotating element (rotor) is supported without contacting the stationary elements (bearing) by the self-developing pressure of the gas (ambient air in most applications). Low power loss and wear can be achieved as a major advantage together with the enhanced stability characteristics for the entire rotor-bearing system. However, the low load capacity (due to the low viscosity of the gas) is a limitation for the wide applicability of GFBs. Several types of GFBs have been introduced in the past, with the most common and efficient being the bump-type foil bearing. The major difference detected in comparison with the conventional oil bearings is the presence of a thin gas film as a lubricant, which results to building up an aerodynamic, load-carrying lubrication wedge and eliminates the need for external pressurization [1-3, 3, 4]. For more than half a century GFBs have been included in various applications, among which are the air cycle machines in commercial aircraft [6], oil-free turbochargers [5], oil-free rotorcraft propulsion engines [6], and micro-gas turbines [7][8][9]. In most applications, GFBs tend to substitute rolling element bearings, increasing the reliability of the rotating machine, the mean time between failures [5,6,6], and the respective sound emissions. The stiffness of the foil element is lower than the effective stiffness of the gas film and a GFB adapts to the operating conditions through the deformation of the foil. Further to the enhanced stability characteristics, GFBs provide large range between 2nd and 3rd critical speeds, and high rotating speeds can be achieved.
Nonlinear dynamics of rotor-GFB systems, and its study with tools like continuation methods, is relatively new object. Continuation methods have been applied on the nonlinear dynamics of rotor systems on oil bearings of several types [42][43][44][45][46][47][48], and recently in GFB rotor systems as well [49]. Bifurcations of Hopf type have been investigated at both bearing types (with oil and gas) [50][51][52][53]. The strongly nonlinear aerodynamic forces render a variety of motions and stability quality in the system, including periodic, quasiperiodic and chaotic motions. Further to that, the system has the potential of totally different motions even at the same operating speed, according to initial conditions and operating parameters, such unbalance magnitude. Dynamic systems present stable and unstable solution branches and respective bifurcation sets, and this is the case in rotor-GFB systems too. Nonlinear dynamics of rotor-GFB systems were mostly used to correlate the quality of response of the system with advances in top/bump foil structure simulation, or alternative models in the aerodynamic simulation. Some of the recent work existing in the literature is presented hereby. Most of these works include advances in the simulation of aerodynamic lubrication and in the simulation of bump foil structure. The response of the rotor system at each case is analysed with respect to its nonlinear features.
Significant work has been conducted in solving the compressible Reynolds equation numerically. The authors in [12,22,54] used the Finite Difference Method (FDM) to substitute the time derivatives of pressure _ p ¼ dp=dt and of film thickness _ h ¼ dh=dt included in Reynolds equation by backward-difference approximations. In [25,26,28,29] a state variable w ¼ p Á h was introduced in order to solve the Reynolds equation and simultaneously acquire results for other state variables, at the same time step. In parallel, FDM and Galerkin reduction methods were used in order to perform the spatial discretization. In [19], a slightly different approach by using the FDM method and Galerkin reduction method without inserting the state variable w was presented. Instead, a backward-differences approximation for the spatial differentiation of the pressure and fluid film variables was included and the Reynolds equation for the time derivative _ p was solved. Further to the aerodynamic modelling and, the bump-type foil modelling has been investigated during the last decades as its properties are directly correlated to the aerodynamic performance of the bearing (fluid-structure interaction). Heshmat et al. [1] introduced the so-called simple elastic foundation model, consisting of linear elastic uncoupled springs. No viscous damping was taken into account until Peng and Carpino [11] and Ku and Heshmat [55] took into consideration the top-to-bump foil and bump-tohousing Coulomb friction damping. Later on, Peng and Carpino [56] introduced a 2D foil model structure characterized by linear stiffness and damping coefficients, while considering the elastic behaviour of the bump foil. Alternative approaches have been published by Lez et al. [54] and other researchers [32,49], where the foil bumps and their interaction are modelled by multi-degree freedom systems. Recently, Baum et al. [19] introduced another foil structure, consisting of rigid, massless, beam-like elements with one finite dimension in axial direction and no coupling of the elements in the circumferential one, where each bump foil is modelled by a nonlinear spring and a linear damper.
An ample number of tests and simulations have been executed throughout the years on the dynamics of complete systems with flexible rotors and gas foil bearings. Nonlinear behaviour has been detected by Lez et al. [54], which can be attributed to the strongly nonlinear bearing forces. Baum et al. [19] performed a run-up and a coast-down simulation and obtained the response of the system, in addition to a waterfall diagram depicting the vibrations arising throughout the whole process. After a certain speed, qualitative and quantitative changes in the rotor's behaviour occurred. Bhore and Darpe [27] conducted an analysis of a flexible rotor supported on two identical GFBs, where Poincaré maps, bifurcation and Fast Fourier transform (FFT) plots were extracted in order to observe the influence of different operating parameters on the system behaviour.
The above-mentioned work in the literature aims to the enhanced modelling and precise prediction of the nonlinear dynamics of rotor-GFB systems. In this paper, a rather simplistic model for bump foil properties of linearized stiffness and damping coefficients is utilized directly from the literature [27], and two design parameters are introduced as foil compliance a f and foil loss factor g. The rotor model follows the Jeffcott rotor model and the shaft's lateral stiffness k s is introduced as the third design parameter in the modelling.
In difference to the aforementioned literature, this paper aims to the insight of the local instability mechanisms, which trigger bifurcations, and drive the response of the system far from its elastic response (self-excited vibrations), usually close to the bearing clearance in journal positions and the rotor-stator clearance along the rotor. In this paper, it is found that the dissipating energy in the gas film should be directly correlated to the first (at lower speed) bifurcation detected in such a system, this being a Neimark-Sacker bifurcation in the unbalanced system, at most cases, or a period-doubling (flip) bifurcation in some others. Specifically, this paper benefits from the wellknown lemma that self-excited motions are triggered when negative damping is included in the system; similar notification was made in [57] for rotors on oil bearings, under linear harmonic analysis. The energy flow among the structural components of the system is evaluated for various design sets and this is used to an optimization scheme to avoid bifurcations in a certain speed range. The paper is organized as follows: In Sect. 2, the dynamic system of consisting of an elastic Jeffcott rotor mounted on two GFBs is composed for the autonomous (balanced), and the nonautonomous (unbalanced) case. The composition of the system renders a set of ordinary differential equations of 1st order (ODE set). The aerodynamic lubrication model and the bump foil structural model follow the existing literature [19,27]. The authors program the pseudo-arc length continuation method [58][59][60] with embedded orthogonal collocation method to provide the potential periodic solutions of the ODE set with unbalance excitation (elastic response or self-excited motions).
In Sect. 3, some indicative quality of nonlinear gas forces is given first for the static problem evaluating the equilibrium locus of the system (fixed point solutions) for several operating conditions. Then, the quality of motions developed in such system is discussed evaluating the bifurcation sets, the timefrequency decomposition of the response, Poincaré maps, and Floquet multipliers, for some design scenarios and at a certain speed range. The energy flow is evaluated in the system during the respective periodic motions, with primary interest on the dissipative work of gas forces. An optimization scheme of successive polls for the two GFB design variables is implemented from the literature [62] and retains the dissipative work in the gas film positive in the respective speed range. Different rotor stiffness and unbalance magnitude are considered in the results.
The paper concludes that there are specific values of bump foil stiffness and damping, evaluated at discrete rotating speeds in the operating speed range, which retain the dissipative work in the gas film in positive values, and therefore self-excited motions (secondary Hopf bifurcations) are not triggered. The specific stiffness and damping of the bump foil is the outcome of an optimization scheme, checked on its efficiency with different objective functions, and for different values of unbalance magnitude and rotor stiffness (flexible rotor or rigid rotor). The bifurcationfree range of operating speed achieved through the optimization procedure is approximately two times larger than that of the reference design (without optimization to be embedded).

Physical and analytical model of a rotor on gas bearings
A representation of a symmetric flexible Jeffcott rotor carrying a disc mass m d at its centre, and two journal masses m j at its ends is depicted in Fig. 1. The simplified symmetric rotor system is selected in this study to avoid the influence of gyroscopic effect in the stability of the system. Figure 2a depicts a gas foil bearing consisting of three major parts: the rigid housing, the bump foil, and the top foil. The shaft is considered elastic, with lateral stiffness k s . The coordinates O d x d ; y d ð Þ and O j x j ; y j À Á represent the geometric centres of the disc and of the journals, respectively, and make up the four degrees of freedom (4 DoFs) of the rotor (both journals perform identical motions due to symmetry). The rotational axes of the journals are considered parallel to the axes of the bearings, an assumption necessary to neglect misalignment and gyroscopic affect. The geometrical centre of the bearing is denoted by O b , and its nominal radius is defined as R þ c r , see Fig. 2a. Notice the fixation point of the foil at v ¼ h 0 .
The physical model of the rotating system is depicted in Fig. 2b where the gas forces are evaluated in the next Section as a function of the bump foil stiffness k f and damping c f .

Analytical model of the gas bearing
The assumptions introduced in the elastoaerodynamic lubrication problem are: a) isothermal gas film, b) laminar flow of the gas, c) no-slip boundary conditions, d) continuum flow, e) negligible fluid inertia, f) ideal isothermal gas law p=q ¼ constant ð Þ , g) negligible entrance and exit effects, and h) negligible curvature R % R þ c r ð Þ . The Reynolds equation for compressible gas flow under these assumptions is given in Eq. 1 [19], and it is an implicit function of time and of journal and top foil kinematics.
Analytical solution for Eq. (1) cannot be defined; a common approach to evaluate the pressure distribution is the FDM. The pressure domain is converted into a grid of i ¼ 1; . . .; N X þ 1 and j ¼ 1; . . .; N Z þ 1 points, where i and j are the indexes in the circumferential and axial direction, respectively, see Fig. 3a.
The Reynolds equation is first rewritten defining the first time derivative of the pressure in Eq. (2a) and after some math in Eq. (2b). Then, the discrete Reynolds equation is defined in the grid points expressing the partial derivatives with finite differences. o os The dimensionless parameters of gas pressure p, gas film thickness h, spatial coordinates (circumferential and axial, respectively)h and z, time s, rotating speed X, and ratio j ¼ R=L b are included in the elastoaerodynamic lubrication problem of Eq. (1). The gas film thickness function is defined in Eq. (3) for the continuous and the discrete pressure domain Boundary and initial conditions of the problem are defined in continue. Ambient pressure is assumed at the starting and ending angle of the foil (periodic boundary condition) in Eq. (4), in the continuous and the discrete pressure domain, respectively.
Taking into account the symmetry of the lubrication problem, instead of assuming the gas pressure equal to the ambient p o at the axial ends, , the boundary condition can be written in Eq. (5) (for the continuous and the discrete pressure domain). In this way the lubrication problem is solved in the half domain, reducing the evaluation cost severely.
Last, the initial conditions for the dimensionless form of the problem are defined in Eq. (6) (in the continuous and the discrete pressure domain).
After evaluating gas pressure p(as p i;j ), the nonlinear gas forces are determined in Eq. (7) In this way the aerodynamic problem renders N X Á N Z ODEs of 1st order with respect to the time derivative of the point pressure in Eq. (8) The vectors x and q may be perceived as x ¼ x j y j x d y d f g T representing the journal motion (coupled to the disc motion through the rotor's equations of motion) and q ¼ (coupled to the journal motion through the Reynolds equation due to the gas film thickness function).
It is important to mention that it is quite common that sub-ambient pressure arises in GFBs. The subambient pressure can cause the top foil to separate from the bumps into a position in which the pressure on both sides of the pad is equalized. Heshmat et al. [1,2] introduced a set of boundary conditions accounting for this separation effect. More specifically, a simple Gümbel boundary condition is imposed, meaning that sub-ambient pressures are discarded when integrating the pressure in Eq. 7 to obtain the bearing force components F B;X ,F B;Y essentially leaving the sub-ambient regions ineffective. In terms of numerical calculations, the assumption made by Heshmat [1,2] can be simply explained as following: in case fluid pressure p is lower than the ambient p 0 , then the former should be considered equal to p 0 and the foil deformation at these points will be zero q i ¼ 0 for p i \1 ð Þ . The simplified model for the bump foil structure is depicted at Fig. 2b. The structure consists of N X À 2 linear massless elements of stiffness k f (compliance a f ¼ 1=k f ) and damping c f . The springs and dampers mount the corresponding N X À 1 top foil stripes of area Dx Á L b (or dimensionless area Dx Á 1), see Fig. 3b. The top foil of the bearing is not covering a complete cylinder; a single gap can be found at h ¼ h 0 , see Fig. 2a, where the top foil is clamped to the bearing housing. Therefore, the moving top foil stripes are N X À 1, see Fig. 3b. The top foil stripes are assumed to remain parallel to the bearing longitudinal axis during their lateral motion; therefore, no axial coordinate is required for the top foil motion. The geometry of the foil structure and its properties, shown in Figs. 2a and 3b, render the dimensionless compliance . The motion of each of the top foil stripe is excited by the mean gas pressure p m;i acting on the top of it, creating the gas force F B i ð Þ, see Figs. 2b and 3b. The mean gas pressure p m;i is defined in Eq. (9) (in the continuous and discrete pressure domain, respectively), for dimensional and dimensionless form.
The foil stiffness and damping coefficient are given as k f ¼ 1=a f and c f ¼ gk f for foil motion synchronous to the excitation. The N X À 1 ODEs that describe the radial displacement q i of the stripe i are defined in Eq. (10).
The ODEs in Eq. (10) may be written as in Eq. (11) to be used in continue.
2.2 Analytical model of the flexible rotor The equations of motion for the Jeffcott rotor shown in Fig. 1 are defined in Eq. (12) for the journal and the disc, in the two main directions, horizontal and vertical.
The ODEs in Eq. (12) may be written in Eq. (13), in the state space representation, to be used in continue. x In Eq. (12) k s is the dimensionless shaft stiffness coefficient, and n,r are dimensionless parameters defined in Eq. (14).
In addition, in Eq. (12), F U;X and F U;Y are the dimensionless unbalance forces defined in Eq. (15a) for constant rotating speed X, and in Eq. (15b) for linearly varying rotating speed X ¼ a Á s with constant acceleration a.
Dimensionless unbalance eccentricity e ¼ e u =c r follows in this paper the ISO unbalance grades (Ggrades) for low, medium, and high unbalance as G1, G2.5 and G6.3 correspondingly. The unbalance located in the disc is of magnitude u ¼ m d þ 2m j À Á e u at each case, and the corresponding eccentricity e u is given by Eq. (16), where the service speed of the system is selected at X r ¼ 2500 rad/s.

Composition and solution of the dynamic system
Equations (8,11, and 13) compose a coupled ODE set which is composed by the discretized Reynolds equation in the f B equations, the foil motion in f F equations, and the rotor motion in f R equations. The coupled nonlinear ODE set is defined in Eq. (17) expressing a non-autonomous dynamic system which will be studied with respect to the bifurcation parameter X. The ODE set is characterized as nonautonomous due to the explicit time presence in the equations of unbalance forces, see Eq. (15). The state vector s and the respective functions f are defined in Eq. (18).
The total number of equations in Eq. (17) Þþ8 with the first term coming from the pressure equations, the second term coming from the foil equations, and the third term from the rotor equations in state space.
The ODE set in Eq. (17) renders the time response of the physical system when time integration is applied [62]. The system is numerically stiff and special algorithms are applied in time integration [62]. Furthermore, the Reynolds equation can be reduced in size applying an order reduction method [19], improving the computational cost. The time integration can handle both cases of unbalance equations, for constant rotating speed or for run-up, see Eq. (15).
An orthogonal collocation method [59] is applied for the computation of limit cycle motions produced by the ODE set in Eq. (17) at a constant X; Eq. (15a) applies for unbalance forces at this case. Numerical continuation of limit cycles has been programmed by the authors according to pseudo-arc length continuation method [58,61,63,64] with embedded collocation scheme [59]. The formulation of the method is defined also in Appendix A1. The need for the implementation of collocation method for the evaluation of periodic limit cycles is based on the advantage of it to handle large state vectors (many degrees of freedom) which in the current system are sourced on the pressure states defined by the finite difference mesh in the solution of the aerodynamic lubrication problem. An alternative of e.g. the shooting method would not be efficient in this case as this requires perturbations of the initial state at a limit cycle, at all degrees of freedom, with respect to all state variables; every perturbation is followed by time integration for the time of one period, increasing the computational time dramatically; this would be not time efficient in the problem of this work. A finite difference method could have been used to solve the boundary value problem for the definition of periodic limit cycles in this problem. The authors chose to implement the collocation method as defined by Doedel in [59], as in the literature this is recommended in dynamic systems of high order.
As the collocation method cannot handle nonautonomous ODE systems, Eq. (17) has to be converted to autonomous. This is achieved by coupling the ODE set of Eq. (17) with a two DoF oscillator, see Eq. (19), whose unique solution is a harmonic motion of frequency X, see Eq. (20) [59].
The size of the final autonomous ODE set is N þ 2 and is defined in Eq. (21) with the unbalance forces to be defined at constant rotating speed, in Eq. (22).
3 Results and discussion

Static performance of the gas bearing
The nonlinear feature of gas forces is presented in this Section. The vertical displacement of each journal (equal at both journals due to symmetry) when load W is applied vertically at each journal, further to the gravity load, is depicted in Fig. 4 for low, moderate, and high bump foil compliance, and for different design and operating conditions of the GFBs. In Fig. 5, the corresponding equilibrium locus of the journal inside the radial clearance is depicted for the various cases mentioned. One may notice that the nonlinear feature of the gas forces becomes less or more intense depending on the design configuration of the bearing and the respective operating speed X. The gas foil bearings configured for the results in Figs. 4 and 5 are of length to diameter ratio equal to 1, while the ratio of the radial clearance to the journal radius differs among the cases R=c r ¼ 500; 750; 1000 to produce different Sommerfeld number S. In Fig. 5 one may notice that the gas film pressure at the axial middle of the bearing (pressure profile) may exceed two times the ambient pressure when bump foil is relatively stiff and specific load is high. The dynamic viscosity l ¼ 0:018mPa Á s of the gas equals this of the ambient air at 20 o C, and the ambient pressure is set at p 0 ¼ 100kPa % 1atm. Radial clearance ratio is R=c r ¼ 500 in the results produced at all next Sections while all other parameters are retained as defined hereby.

Quality of motion of the dynamic system
The dynamic system defined in Eq. (21) in autonomous and in Eq. (18) in non-autonomous version is investigated on its potential to develop a variety of bifurcation sets with respect to the key design parameters, namely rotor stiffness k s , foil compliance a f , foil loss factor g, and unbalance magnitude u. In Fig. 4 Vertical displacement of the journal as a function of vertical specific load for a a f ¼ 0:01, b a f ¼ 0:1, and c a f ¼ 1. Different rotating speed is considered to allow for influence of compressibility this paper, the key design parameters are defined within specific intervals, composing the case studies which are presented in continue. The design parameters follow a variation of ''low'', ''reference'', and ''high''. This is interpreted to the rotor stiffness values k s ¼ 0:3; 1; 3(flexible to rigid rotor), foil compliance values a f ¼ 0:01; 0:1; 1(stiff to flexible foil), foil loss factor g ¼ 0:005; 0:05; 0:5(low to high foil damping), and unbalance magnitude u G1 ð Þ; u G2:5 ð Þ; u G6:3 ð Þ(low to high unbalance). A reference system is defined with the design parameters to be k s ¼ 1, a f ¼ 0:1, g ¼ 0:05, and u G6:3 ð Þ. The validation of the developed code is presented in Fig. 6 considering the fixed-point equilibrium state and the respective stability property, and the extension of a stable limit cycle at a selected rotating speed, comparing with the recent literature [53]. In Fig. 6a one may note that the fixed-point solution at the rotating speed of X ¼ 12000 RPM, and with the properties of Table 1, converges to this evaluated in [53] with the same system properties. The divergence in the transient response (path followed by the journal towards the asymptotically stable fixed point) is due to the different time integration solver implemented in the numerically stiff system.
A bearing load of F B;Y ¼ 10N is then applied according to [53] to evaluate the stable fixed-point solutions and limit cycles as presented in Fig. 6b at the rotating speed of 13,000 RPM, and in Fig. 6c-d for a range of rotating speed; the fixed-point solutions agree well between the two models. However, the extent of limit cycles is not identical and this is due to the fact that the shooting method is used in [53] to predict the limit cycle solutions while in current work, numerical continuation and collocation method are implemented. In this example, Hopf bifurcation occurs in 12,030 RPM according to the current model, while in [53] this is at 11,475 RPM, see Fig. 6c-d; this is a difference of less than 5% and can be explained by a different perturbation in the evaluation of Jacobian matrix. Note that fixed point stability is calculated with a linear stability analysis at both works.
The time transient response of the reference system is evaluated applying time integration in Eq. (21) defined for variable rotating speed (run-up) of constant rotating acceleration; the result is used only for comparison. The ODE system is subjected to numerical continuation algorithm which evaluates limit cycles at the different rotating speeds (bifurcation parameter). In this paper, limit cycles refer in periodic motions caused either by elastic response or selfexcitation.
In Fig. 7a the time history of the journal motion in the vertical plane (evaluated by time integration) is presented together with the maximum and minimum values of the limit cycle at each rotating speed. The rotating speed is retained constant at each limit cycle evaluation (several discrete speed values are used to cover the entire range of rotating speed), and the unbalance forces are implemented through Eq. (15a).  Fig. 7a with PD, SN, and NS bifurcations to be presented. The frequency content of the time history obtained during the run-up (by time integration) is depicted in Fig. 7b where time-frequency decomposition is applied (STFT). The journal trajectory at a selected rotating speed is depicted in Fig. 7c together with the respective limit cycles evaluated with the collocation method. The transient response and the Poincaré map depicted also in Fig. 7c consider the last 100 periods of the response evaluated for 500 periods applying time integration [62]. The Floquet multipliers in Fig. 7d provide information regarding the quality of bifurcations mentioned above, at the respective speeds.
In general, periodic, quasi-periodic and chaotic motions are expected to be generated by the system studied in this paper. A drawback of the current continuation scheme programmed by the authors is that only periodic motions can be evaluated by the collocation method. The quasi-periodic motions can only be evaluated by time integration in this paper and depicted as transient response. For the evaluation and continuation of quasi-periodic limit cycles, the reader may consider the recent work in [66].
Further bifurcation sets for the respective design sets are presented in continue in order to study the system motion for the various combinations of design  parameters and to reveal the types of motions and of bifurcations which may occur. These are defined considering one design parameter, changing each time, to the next higher and lower value of the respective design variable. Figure 8 considers systems of different rotor stiffness, and two cases are presented in Fig. 8a and c, for k þ s ¼ 3(rigid rotor) and k À s ¼ 0:3(flexible rotor). One may notice the different bifurcation sets compared to the reference case. Period-doubling bifurcation is not noticed in this case.
In Fig. 9, systems of different foil compliance are considered, and two cases are presented in Fig. 9a and c, for a þ f ¼ 1(flexible foil) and a À f ¼ 0:01(rigid foil). One may notice the different bifurcation sets compared to the reference case, and the previous case (Fig. 8). The type of bifurcations is same to this at the reference case, but the speed in which they appear is different.
In Fig. 10, systems of different foil damping are considered, and two cases are presented in Fig. 10a and c, for g þ ¼ 0:5(high foil damping) and g À ¼ 0:005(low foil damping). The influence of foil damping in the bifurcation set is not severe, compared to the reference design (see Fig. 7).
In Fig. 11, systems of different unbalance are considered, and three cases are presented in Fig. 11a, for u G0 ð Þ (balanced rotor-autonomous system), u G1 ð Þ(low unbalance), and u G2:5 ð Þ (medium unbalance). It is worth noticing that the autonomous system of u G0 ð Þ loses local stability of fixed-point equilibria through an Andronov-Hopf bifurcation, at similar speed where the unbalanced systems lose local stability through secondary Hopf (Neimark-Sacker bifurcations). Further to that, in the unbalanced systems, the higher the unbalance is, the lower the speed of NS bifurcations is. Stable limit cycles close to radial clearance occur with higher amplitude in the balanced system, than in the unbalanced systems. In the balanced system the limit cycles of amplitude Fig. 7 Reference system: system of k s ¼ 1, a f ¼ 0:1, g ¼ 0:05, and u G6:3 ð Þ. Transient response and continuation of limit cycles during run-up in a Vertical direction, and b Horizontal direction. c STFT of the response time history, d Floquet multipliers of the corresponding limit cycles close to radial clearance will lose stability through a NS bifurcation at high speeds.
In the figures above, Poincaré maps are constructed at selected speeds to depict the different quality of motions (periodicity), and the respective frequency content can be considered through the FFT and the STFT in each figure. At lower speeds, the system oscillates periodically with 1T period at all cases (T ¼ 2p=X is the driving period hereby). This renders one point in Poincaré map, and additional harmonics of integer multiple in the FFT, see Fig. 8c. This status remains till a periodic doubling occurs at some cases, e.g. Figure 10c, and the system oscillates with period 2T immediately after; this renders two points in Poincaré map. At the highest speeds investigated, only quasi-periodic or chaotic motions where evaluated, these through time integration. The respective quality and frequency content is depicted in Figs. 9c and 11c, where Poincaré map is constructed with several points without forming a shape. The frequency content at such speeds includes higher and lower harmonics of synchronous response without integer multiple.

Energy flow and optimization for bifurcation elimination
The energy flow among the structural components of the system is evaluated in this Section aiming to reveal any correlation of the non-conservative work with the existence of bifurcations. The reference design cases presented in Sect. 3.2 were used also in this Section. The scope of this Section is to apply an optimization In order to evaluate the work of the non-conservative forces, the energy balance in the system is studied first. The storage of potential energy occurs in the stiffness elements of the system, these being the flexible rotor of lateral stiffness k s , and the flexible bump foil consisting of N x elements of stiffness k f . Note that k f expresses the stiffness coefficient per top foil area. The gas film possesses an unknown stiffness property and therefore potential energy is considered to be stored in the gas through the work of the conservative part of the gas forces. Kinetic energy storage occurs in the three masses of the system, these being the disc mass m d and the two journal masses m j ; note that bump foil is considered massless and no inertia effects are considered in the gas flow inside the gas bearings. Dissipation of energy occurs in the gas film forces (viscous damping) and in the bump foil damping elements of damping coefficient per area c f . The system executes planar motions with the energy offered through the work of the unbalance forces. The unbalance force is a result of rotation, and a motor is supposed to retain the rotating speed constant. The lateral vibrations are uncoupled to the energy offered by the motor.
In a closed trajectory (limit cycle motion) of the three rotor masses and the bump foil elements, the collocation method offers all values of the state vector s while _ s can be easily found by Eq. (21a). Let the time intervals during a limit cycle motion be N t , defining N t þ 1 discrete time points (these are the time intervals which are further divided by collocation points; this is not of interest hereby as N t renders a short enough time Fig. 9 System of k s ¼ 1,a À f ; a þ f , g ¼ 0:05,u ¼ G6:3 ð Þ. a Continuation of limit cycles, b STFT of the response time history for a þ f , c Trajectory, Poincare map, and FFT at X ¼ 0:9, d Floquet multipliers of the corresponding limit cycles interval). The corresponding gas forces at each discrete time point are evaluated by the state vectors and its time derivative _ s using Eq. (7), these being F B;X i ð Þ and F B;Y i ð Þ, for i ¼ 1; 2; :::N t þ 1. Similarly, the bump foil forces are evaluated in the radial direction as F f ; j i ð Þ for the j th element of the N x in total (see Figs. 2b and 3b), q j to be contained in the already knowns and _ s. The unbalance forces are defined as F U;X i ð Þ and F U;Y i ð Þ from Eq. (15a). The work of the bearing forces is evaluated in Eq. (23a), the work of the bump foil forces is evaluated in Eq. (23b), and the work of unbalance forces in evaluated in Eq. (23c). Fig. 10 System of k s ¼ 1 a f ¼ 0:1g À ; g þ ,u G6:3 ð Þ. a Continuation of limit cycles, b STFT of the response time history for g þ , c Trajectory, Poincare map, and FFT at X ¼ 0:57, d Floquet multipliers of the corresponding limit cycles for g þ In Eq. (23a) the portions W cb and W kb cannot be evaluated separately, but only in total. However, it is expected that W kb % 0 as this is the work of the conservative portion of the bearing forces in the closed trajectory. In Eq. (23b) W cf is evaluated as W cf ¼ W f ; it is found that W kf ¼ 0 as this is the work of conservative forces in a closed trajectory. In a closed trajectory, the kinetic energy T 1 and T 2 of the system at the beginning and at the end of the trajectory, respectively, must equal the work of the non-conservative forces plus the work offered to the system. Therefore, Eq. (24) holds, with T 1 ¼ T 2 (periodic motion) [67].
Equation 24 renders W cb as W fu is known by Eq. (23c) and W cf is known by Eq. (23b). The value of W cb is also validated by Eq. (23a). The work of the non-conservative forces W cb and W cf together with the work of the unbalance force W fu is depicted in Fig. 12 for selected design sets from the previous Section.
In Fig. 12a-d, both stable and unstable limit cycles are considered with the respective notation. At all cases, it is found that Neimark-Sacker bifurcations are triggered simultaneously to the reverse (from positive values to negative) of the dissipating work in the gas film W cb , meaning that energy is not dissipated in the gas film (when W cb \0) and self-excitation takes place. The respective limit cycles for the cases in Fig. 12 can be found in Sect. 3.2. In Fig. 12a-d, the arrows depict the path that would be followed during the run-up of the system (evaluated by time integration of the motion equations).
For the understanding of the sensitivity of nonconservative work W cb to the bump foil properties, this is evaluated for several values of bump foil properties, Fig. 11 System of k s ¼ 1,a f ¼ 0:1g ¼ 0:05. a Continuation of limit cycles, b STFT of the response time history, c Trajectory, Poincare map, and FFT at X ¼ 0:97, d Floquet multipliers of the corresponding limit cycles Fig. 12 Evaluation of energy flow at the respective limit cycles for a) k s ¼ 3, ð Þ, and d) k s ¼ 1, a f ¼ 0:1, g ¼ 0:05, u G2:5 ð Þ Fig. 13 Dissipated energy in the gas film in one limit cycle for various values of foil compliance a f and foil loss factor g, at a X ¼ 0:2, b X ¼ 0:4, c X ¼ 0:6 and at three different speeds (low speeds), where local stability is guaranteed according to the analysis in Sect. 3.2.
The work W cb is plotted in Fig. 13 for the respective cases of bump foil design and for the reference system regarding rotor stiffness and unbalance magnitude. Figure 13 depicts the sensitivity of W cb in regards to bump foil properties. It is worth noticing that at the lowest rotating speed X ¼ 0:2 the dissipated energy has a maximum for specific value of foil compliance a f % 1 and for low values of loss factor g. This is not the case at the higher speed X ¼ 0:4, where W cb has a maximum still for a f % 1, but for high values of loss factor g, see Fig. 13b. Increasing the speed further at X ¼ 0:6, W cb has the highest values for several pairs of a f and g, the former receiving the lower values, and the latter at higher values, see Fig. 13c.
The conclusion made from Fig. 13 is that the dissipated energy in the gas film has different progress (related to the design of the foil) as the rotating speed changes. Therefore, retaining the dissipated energy in the gas film at highest values requires a continuous adaptation of foil properties in regards to the rotating speed.
Further rotor designs and unbalance magnitudes were checked with relatively similar conclusion on the sensitivity of dissipated energy to design characteristics. Figure 14 depicts that the dissipating energy in the bump foil dampers (through work W cf ) follows similar progress in regards to the bump foil design, at all three speeds. This is increased when foil becomes flexible and of low damping, as this combination allows higher displacement.
Specific values of bump foil properties which render the highest dissipation energy in the gas film are sought applying the optimization pattern of successive poles [65]. The optimization scheme is implemented in several discrete speeds. The optimization requires the minimization of an objective function OBJ, which is defined as the inverse of dissipated energy in the gas film, OBJ ¼ 1=W cb . This is the 1st optimization problem.
Starting from random input values for foil compliance a f and foil loss factor g, the optimization pattern renders after some iterations the values of a f and g that maximize W cb at every speed X. The related progress of W cb during optimization is depicted in Fig. 15 at selected speeds. The optimization pattern works in the following sequence, see also Appendix A2: (a) a low rotating speed X(in which local stability is guaranteed) is selected for the system and random foil properties are assigned to the model, as a (e) A prediction for the limit cycle at the next speed, say X þ DX, is computed by pseudo-arc length continuation, given the limit cycle e s n , and a n f ,g n , as initialization. The limit cycles n is plotted in Fig. 16 for different rotor stiffness scenarios (flexible-rigid), and in Fig. 17 for different unbalance scenarios (low-high). The respective values a In Figs. 16a and 17a one may notice the displacement of the limit cycles in vertical direction as the foil compliance receives higher and lower values. The Floquet multipliers in Fig. 16d and 17d are retained within the unit circle as the rotating speed increases and tend to cross the unit circle at the points À0:9; AE0:25 ð Þ . As the optimization pattern aims to retain stability, the Floquet multipliers appear to turn to a new direction as the rotating speed increases, and they continue moving along the path of the unity circle, still inside the circle, up to the points À0:2; AE0:9 ð Þ . At the highest rotating speeds checked, the optimization cannot retain the dissipative energy positive and this results in a bifurcation. The Floquet multipliers do not lie inside the unit circle anymore, resulting in Neimark-Sacker bifurcation. The respective values of the dissipative energy are presented in Figs. 16a and One may notice that the optimization pattern renders flexible foil (higher compliance) in order to allow motion of the foil and more energy dissipation in the system. However, this is not the case for the rigid rotor (k s ¼ 3) in Fig. 16b.
An alternative objective function was investigated in the 2nd optimization problem. At each limit cycles i , the highest magnitude of the Floquet multipliers was defined as objective function, OBJ ¼ max l j À Á , neglecting the one existing always at the unit circle, point þ1; 0 ð Þ. In this way, all Floquet multipliers are retained inside the unity circle, l j \1. The respective limit cycles are presented in Fig. 18 for different rotor stiffness and in Fig. 19 for different unbalance magnitude. The evaluation time for the limit cycles was lower at this case (compared to the optimization considering the dissipative work), as the computation of dissipative work is not required. It is interesting to note how the optimization algorithm retains the Floquet multipliers inside the unity circle as the rotating speed increases, see Figs. 18d and 19d. In Figs. 18b and 19b one may notice the very similar trend of bump foil compliance to obtain high values, at all design cases. Similar trend is noticed for the bump foil loss factor, when comparing Figs. 18c and 19c; the trend is to obtain lower values at all design cases. Utilizing the Floquet multipliers as objective function, the rotor does not experience discontinuities in the location of limit cycles, as occurred when dissipative work was used in the objective function (1 st optimization problem, see Figs. 16a and 17a). The 2 nd optimization problem reveals that a stable system requires rather compliant bump foil with low damping.

Conclusions
The bifurcation sets of a rotating shaft on gas foil bearings were initially presented in this paper for flexible and rigid rotors with high and low unbalance, and for various bump foil stiffness and damping properties, at a certain range of rotating speed. The periodic limit cycle motions were evaluated for these sets of designs applying a pseudo-arc length continuation method with embedded orthogonal collocation. At that stage the potential of the system to develop different motions (periodic, quasi-periodic, chaotic) and different types of stability was investigated.
The work of the non-conservative (and nonlinear) damping force of the gas film was evaluated at each limit cycle, even when unstable, as the collocation method allows for this possibility. As the rotating speed of the system increased, the dissipative work generated by the gas film forces was found to be directly correlated to the self-exciting mechanism which triggers secondary Hopf bifurcations (Neimark-Sacker bifurcations) leading the journal motion towards limit cycle motions of high extent (operability Fig. 17 Elimination of bifurcations at a speed range for the system of k s ¼ 1 and u G1 ð Þ, u G2:5 ð Þ, u G6:3 ð Þ: a Journal motion limit cycles, and corresponding values for b Compliance a f , c Loss factor g, d Floquet multipliers, and e Dissipating work in the gas film of the system may be compromised at this case). It is found that the loss of this local stability occurs simultaneously with the reversal in the energy flow in the gas film, meaning that the dissipative work changes sign (from positive to negative) when the NS bifurcation takes place.
Using this notation, an optimization pattern was applied at discrete speeds during a virtual run-up of the system, and aimed to maximize the dissipated work in the gas film (retain it positive); this was the 1st optimization problem. The input variables for the optimization were the values of bump foil stiffness and damping. In this manner, bifurcations of limit cycles were avoided at a wide range of rotating speed, when specific bump foil stiffness and damping were applied through the optimization pattern.
A generic outcome from the optimization patterns is that bifurcations are avoided when bump foil stiffness is reduced; bump foil damping appeared to influence the existence of bifurcations much less than stiffness. Low damping values should be preferred in the bump foil though, according to the 2nd optimization pattern. The operating speed range without bifurcations to take place was found almost double, compared to the other sets (without optimization). The optimization procedure was repeated for several design scenarios of rotor stiffness and unbalance magnitude, and similar efficiency was noticed regarding bifurcation elimination. Research on design solutions to implement the change of foil damping and stiffness in real system operation belongs to ongoing work.
For the theoretical completeness of the work, an optimization pattern utilizing the magnitude of Floquet multipliers in the objective function (retaining Floquet multipliers inside the unit circle) was checked, and similar results for the bump foil properties were Fig. 18 Elimination of bifurcations at a speed range for a system with unbalance u G2:5 ð Þ and stiffness k s ¼ 0:3,k s ¼ 1, and k s ¼ 3: a Journal motion limit cycles, and corresponding values for b Compliance a f , c loss factor g, d Floquet multipliers, and e Dissipating work in the gas film found for a bifurcation-free operating range; this was the 2nd optimization problem.
Funding Open access funding provided by HEAL-Link Greece.
Data availability All data generated during the study are available from the corresponding author by request.
Code availability Custom code was used for the composition of the dynamic system, the collocation method, and the pseudoarc length continuation of limit cycles and of fixed points. Time integration when needed was performed by [62]. Optimization was performed by [65].

Declarations
Conflict of interest The authors declare that they have no conflict of interest or competing interest with regard to the publication of this manuscript. All data generated during the study are available from the corresponding author by request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Fig. 19 Elimination of bifurcations at a speed range for the system of k s ¼ 1 and u G1 ð Þ, u G2:5 ð Þ, u G6:3 ð Þ: a Journal motion limit cycles, and corresponding values for b Compliance a f , c Loss factor g, d Floquet multipliers, and e Dissipating work in the gas film where x; _ x 0 h i denotes the scalar product and _ x 0 is the time derivative of the previous solution. If the pseudo arc length is used as a continuation parameter then X also becomes an unknown and an additional equation is required (pseudo-arc length condition), where ðÞ 0 denotes the derivative with respect to arc length d Á ð Þ=ds. Setting u ¼ x; T; X ð Þand writing Eqs.
(27 and 28) as F u ð Þ ¼ 0 f g the system to solve becomes and Newton's method for solving (31) is which is iterated until a suitable convergence criterion is satisfied. The arc length derivatives du=ds can be calculated either by backwards differences or by solving The final step is to discretize in time and calculate A. To this end the method of orthogonal collocation at Gauss points with piece-wise polynomials was used. An overview of the method applied to nonlinear BVPs with periodicity boundary conditions with unknown period is given below.
The time interval ½0; 1 is discretized into N subintervals. For the ith sub-interval the collocation equations must be assembled at the required time nodes where h i ¼ t iþ1 À t i the length of the time sub-interval i and q j are chosen as the zeroes of an mth order Legendre polynomial. At the above time nodes, an initial solution x ij must be provided along with the function evaluation T Á f x ij ; X À Á (abbreviated henceforth as f ij ), Jacobian of ox x ij ; X À Á , and parameter derivative (for the case of pseudo arc length continuation) of oX x ij ; X À Á . Equivalently, the values x ij can be extracted from the solutions at the global time nodes t i as l¼1 a jl f il ; ð36Þ where a j1 ; a j2 ; . . .; a jk are the quadrature weights. Then the quasi linearized two point BVP (equivalent to Newton's method) can be written as shown below where r ij ¼ Tf x ij ; X À Á À _ x ij . Applying parameter condensation to eliminate the local unknowns x ij at every time interval t i we can write the derivatives T for the local unknowns as a function of the global unknowns x i . Substituting (36) in (37) yields The above can be rewritten as where where a ij are the quadrature weights for the time intervals t iÀ1 t i ½ . The solution of the linear system can be achieved by various methods. Iterative methods are applied in this work. Floquet multipliers are evaluated as the eigenvalues of the matrix C 1 Á C 2 Á Á Á C N when the iterative solution of the n Â N system is achieved (right hand side less than a maximum). Calculating the Floquet multipliers in this way severely reduces the evaluation time compared to other methods (e.g. shooting method). The normal form coefficient is calculated with different method for the different type of bifurcation occurring. For more detailed information, the reader may refer to [63,64].

A2. Implementation of the optimization pattern
Patternsearch [65] finds a sequence of points x 0 ; x 1 ; x 2 ; :::; that approach an optimal point. The value of the objective function either decreases or remains the same from each point in the sequence to the next. The pattern search begins at the initial point x 0 provide as initialization by the user.

Iteration 1
At the first iteration, the mesh size is 1 and the Generalized Pattern Search (GPS) algorithm adds the pattern vectors to the initial point x 0 to compute the following mesh points: 1; 0 ½ þ x 0 , 0; 1 ½ þ x 0 , À1; 0 ½ þx 0 , 0; À1 ½ þx 0 . The algorithm computes the objective function at the mesh points in the order shown above. The algorithm polls the mesh points by computing their objective function values until it finds one whose value is smaller than the value at x 0 . If it finds that value then, the poll at iteration 1 is successful and the algorithm sets the next point in the sequence equal to x 1 .

Iteration 2
After a successful poll, the algorithm multiplies the current mesh size by 2. Because the initial mesh size is 1, at the second iteration the mesh size is 2. The mesh at iteration 2 contains the following points: 2 Á 1; 0 ½ þx 1 , 2Á 0; 1 ½ þx 1 , 2Á À1; 0 ½ þx 1 , 2 Á 0; À1 ½ þx 1 . The algorithm polls the mesh points until finds one whose value is smaller than the value at x 1 . If it finds that value then, the poll at iteration 2 is again successful. The algorithm sets the second point in the sequence equal to x 2 . As the poll is successful, the algorithm multiplies the current mesh size by 2 to get a mesh size of 4 at the third iteration.

An unsuccessful poll
By the fourth iteration, the current point is x 3 and the mesh size is 8, so the mesh consists of the points: 8 Á 1; 0 ½ þx 3 , 8Á 0; 1 ½ þx 3 , 8Á À1; 0 ½ þx 3 , 8 Á 0; À1 ½ þx 3 . At this iteration, none of the mesh points has a smaller objective function value than the value at 9 3, so the poll is unsuccessful. In this case, the algorithm does not change the current point at the next iteration. That is x 4 ¼ x 3 . At the next iteration, the algorithm multiplies the current mesh size by 0.5, so that the mesh size at the next iteration is 4. The algorithm then polls with a smaller mesh size. Patternsearch modifies poll points to be feasible at each iteration, meaning to satisfy all bounds and linear constraints. In this paper, the boundaries that are used to implement the optimization for a f and g are 0:005\a f \2 and 0:0001\g\10, respectively.

Stopping conditions for the pattern search
The algorithm stops when any of the following conditions occurs: 1. The mesh size is less than the MeshTolerance option which was set at 10 À6 . 2 The number of iterations performed by the algorithm reaches the value of the MaxIterations option which was set at 300. 3. The total number of objective function evaluations performed by the algorithm reaches the value of the MaxFunctionEvaluations option. 4. The time, in seconds, the algorithm runs until it reaches the value of the MaxTime option. In our optimization that was free. 5. After a successful poll, the distance between the point found in the previous two iterations and the mesh size are both less than the StepTolerance option. 6. After a successful poll, the change in the objective function in the previous two iterations is less than the FunctionTolerance option and the mesh size is less than the StepTolerance option.