Bifurcation analysis of rotor/bearing system using third-order journal bearing stiffness and damping coefficients

Hydrodynamic journal bearings are used in many applications which involve high speeds and loads. However, they are susceptible to oil whirl instability, which may cause bearing failure. In this work, a flexible Jeffcott rotor supported by two identical journal bearings is used to investigate the stability and bifurcations of rotor bearing system. Since a closed form for the finite bearing forces is not exist, nonlinear bearing stiffness and damping coefficients are used to represent the bearing forces. The bearing forces are approximated to the third order using Taylor expansion, and infinitesimal perturbation method is used to evaluate the nonlinear bearing coefficients. The mesh sensitivity on the bearing coefficients is investigated. Then, the equations of motion based on bearing coefficients are used to investigate the dynamics and stability of the rotor-bearing system. The effect of rotor stiffness ratio and applied load on the Hopf bifurcation stability and limit cycle continuation of the system are investigated. The results of this work show that evaluating the bearing forces using Taylor’s expansion up to the third-order bearing coefficients can be used to profoundly investigate the rich dynamics of rotor-bearing systems.

order bearing coefficients can be used to profoundly investigate the rich dynamics of rotor-bearing systems.

List of symbols
The angular coordinate The rotational speed of the journal rad/s

Introduction
Journal bearing is one of the crucial elements used in industry. It has many applications in heavy-duty machinery whether moderate or high speed such as reciprocating engines [1,2], turbomachines [3,4], centrifugal pumps [5,6] and turbocharger [7][8][9]. Many of these machines are required to spin at very high speeds to improve their efficiency. Therefore, it is important to investigate and evaluate the threshold whirling speed, because above this speed oil whirl may be occurred [10,11]. However, some applications are reported to operate at speeds higher than this threshold value in the stable condition [12]. In addition, the experimental work of Muszynska [13], Deepak and Noah [14] showed that a stable whirling is occurred after threshold speed. Therefore, it is important to investigate the nonlinear dynamics behavior of journal bearings to improve the future design of rotating machinery. The dynamic modelling of journal bearings requires the evaluation of the bearing forces. Since the fluid film pressure distribution inside the journal bearing can be described by Reynolds equation, the bearing forces can be obtained from the integration of the bearing fluid film pressure distribution. In case of short and long bearings, Reynolds equation can be approximated and simplified so that the bearing forces can be obtained analytically, see for example [6,15,16]. Nishimura et al. [6] investigated the nonlinear dynamics and stability of a vertical flexible rotor supported by short journal bearing. Castro et al. [15] investigated the nonlinear dynamics of rotors supported by journal bearings based on short bearing approximation. They investigated two types of instability, which are 'oil whirl' and 'oil whip' during run-up and run-down. In case of finite length bearings, where the analytical solution for Reynolds equation is not available, the bearing forces can be obtained by direct integration of Reynolds equation.
In dynamic analysis, this integration is required to be done each time step, which is computationally expensive especially if very fine mesh is considered [17][18][19].
Investigating the Hopf bifurcation stability is important because it represents the transformation from spin motion to periodic whirl motion. Although the bearing forces based on linear first-order approximation coefficients are sufficient to obtain the threshold speed, they are not sufficient to judge the stability of periodic solution and whether it is subcritical or supercritical bifurcation. In the case of subcritical bifurcation, if the amplitude of perturbation is located outside the limit cycle, the stability of the rotor bearing system is lost and the perturbation will trigger oil whip; subsequently, the system becomes unstable. If the perturbation is located inside the limit cycle, the rotor will stabilize at a point and the system will remain stable. In the case of supercritical bifurcation and above the threshold speed, the limit cycles are stable and they attract the rotor centre if the perturbation was inside or outside the limit cycle, which physically means that the rotor will whirl in this limit cycle [14,37]. This motivates the researchers to use analytical solution for bearing forces where possible or to use higher-order bearing coefficients.
Dake et al. [38] investigated the nonlinear dynamics and stability of a rotor supported by hydrodynamic journal bearings. They investigated the stability based on both linear bearing coefficients and on bearing forces from the analytical solution. Smolik et al. [39] investigated the stability of rotor bearing system based on the threshold stability curves. They used different approaches to calculate the bearing forces, infinitely short, infinitely short approximation with correction factor for finite length bearing, finite difference, and finite element. They also investigated the run-up time and acceleration.
Miraskari et al. [40] investigated the nonlinear dynamics of flexible rotor supported by journal bearings. The second-order nonlinear bearing coefficients are obtained using the perturbation method. Then, they used a shooting method to evaluate the stability of the periodic solution. Chasalevris [41] investigated the Hopf bifurcation of flexible rotor supported by journal bearings. Chasalevris investigated the lemon bore and partial arc types of bearings. Chasalevris used Hopf bifurcation theory introduced by Hassard et al. [42] to investigate the stability of periodic solution.
The previous literature shows that many studies for the nonlinear dynamics of journal bearing are based on the analytical solution which is only available in case of short or long bearings. Little studies are done based on the second-order linear coefficient. Moreover, very little studies have been done on higher-order bearing coefficients. The main purpose of this study is to evaluate the third-order bearing coefficients using infinitesimal perturbation method. The obtained coefficients are used to investigate the stability of a flexible rotor/disc supported on two symmetric hydrodynamic journal bearings.
After this introduction, the numerical steps to obtain the journal bearing third-order bearing coefficients are presented on the mathematical model section. Then, the dynamic model of flexible rotor supported on two symmetric journal bearings is presented. The results of journal bearing coefficients and their verification and mesh sensitivity are presented on the results and discussion section. Moreover, the stability and continuation analysis of the dynamical model are discussed at the dynamical results subsection. Finally, the most important outcomes and conclusions are presented in the conclusion section.

Analytical model
In this section, the pressure distribution of the journal bearing is evaluated and then used in the calculations of the stiffness and damping coefficients of the oil film journal bearing. Then, these coefficients are used in the dynamic modeling of the elastic rotor-bearing system. The model is for full circular journal bearing o b is the center of the bearing, o js is the center of the journal or shaft at the equilibrium position and o j is the per-turbed journal center. R is the radius of the bearing and R j is the radius of the journal. The radial clearance is c = R − R j . The height of the oil film between the journal and the bearing at any angle is h. e 0 , θ 0 are the eccentricity and the attitude angle at the steady-state equilibrium. O js is the center of the journal or shaft at the equilibrium position in terms of dimensionless coordinates and O j is the perturbed journal center in terms of dimensionless coordinates , see Fig. 1c.
The pressure distribution inside the fluid film can be described by Reynolds equation as follows: where h and p are the oil film thickness and pressure, respectively. φ is the angular coordinate with respect to vertical axis, t is the time in sec., is the rotational speed of the shaft in rad/sec.. x and y are the sectional bearing coordinates as shown in Fig. 1a. z is the axial coordinate as shown in Fig. 1b. The fluid film viscosity and density are μ in N s/m 2 and ρ in kg/m 3 , respectively.
Reformulating the Reynolds equation Eq. (1) for the dimensionless representation, it can be written as: where H is the oil film dimensionless thickness, P = p 6μ c R 2 is the dimensionless pressure, Z = z R is dimensionless location in the axial direction. The dimensionless time is τ = t. To evaluate the bearing coefficients, a small perturbation is introduced from the steady-state point as shown in Fig. 1c. The height in dimensionless can be presented as in Eq. (3). After perturbation, the eccentricity ratio, attitude angle and dimensionless height can be obtained from Eqs. (4), (5) and (6), respectively. Since the analysis in the present article is based on the third-order approximation, all the differences with order higher than three are omitted.
for which ε is the eccentricity ratio and θ is the attitude angle, for a small perturbation For small δθ, substitute cos (δθ ) = 1− δθ 2 2 and sin (δθ ) =δθ − δθ 3 6 , then Eq (6) can be written as Equation (7) can be expanded as follows: for which Now it is important to transform the oil film height Eq. (8) from polar coordinates to Cartesian coordinates. From the geometry, Fig. 1c, the components δ X and δY can be written as: Expanding Eqs. (13)(14) and reducing higher-order terms, we get the following relation: Reordering Eq. (15), we can obtain the following : Substituting the transformation relations in Eq. (16) into Eqs. (9)(10)(11)(12) results in The dimensionless pressure distribution P(X, Y, X , Y ) in the fluid film can be expanded using Taylor expansion as follows: where is the steady-state equilibrium position, and ∂ X ∂Y ∂Y and so on.
To obtain the steady-state pressure P 0 the Reynolds equation in Eq. (2) is solved using finite difference relations. The derivation of the first-, second-and thirdorder pressure gradients of Eq. (20) is discussed in "Appendix A".

Nonlinear bearing forces
As the induced fluid film pressure, the dimensionless bearing forcesF X andF Y are approximated to the third order using Taylor series, for whichF X = F x W and F Y = F y W , where F x is the bearing force in the x direc-tion, F y is the bearing force in y direction and W is the bearing load.
where ζ = X or Y . The bearing forces in X, Y directions in Eq. (21) can be written using the bearing first-, second-and thirdorder coefficients as shown below: for which X = x c and Y = Y c are the dimensionless coordinates, c is the radial clearance.
are the dimensionless linear stiffness and damping coefficients, respectively, where (i, j = X, Y ). k i j and c i j are stiffness and damping dimensional coefficients, respectively.
are the secondorder dimensionless stiffness and damping coefficients, respectively, where i, j, k = X, Y . k i jk and c i jk are the second-order stiffness and damping dimensional coefficients, respectively.
are the thirdorder dimensionless stiffness and damping coefficients where i, j, k, l = X, Y . k i jkl and c i jkl are the thirdorder stiffness and damping dimensional coefficients.
The bearing forces can be calculated by integrating the pressure over the area as shown below: Comparing Eqs. (20), (22) to (24), the steady-state force can be obtained as: It is worth noting that through the present analysis the only steady-state load in y direction is the weight and no static load in x direction. Therefore, F x0 = 0 and F y0 = W which also means that Using this condition and by an iterative method, the steady-state attitude angle can be obtained. Also, the eight first-order coefficients can be obtained from the following equations.
where ζ = X or Y which means that by replacing ζ by either X or Y the eight equations for the eight coefficients can be obtained from Eqs. (28)(29)(30)(31). In addition, the fourteen second-order coefficients can be calculated from where each of ζ and η can be either X or Y . This means that Eqs. (32)(33)(34)(35) represent sixteen equations. These equations are reduced to fourteen equations because Similarly, the twenty third-order coefficients can be evaluated from the following equations: where each of ζ, η and γ can be either X or Y . This means that thirty-two equations can be obtained from Eqs. (36)(37)(38)(39), but these equations are reduced to only twenty because of similarity of several coefficients such as K XY X X = K X XY X = K X X XY .

Rotor-bearing model
In this section, the analytical model for elastic rotor supported on two symmetric journal bearings is investigated. The model is four degree of freedom in x and y directions. In this model, both the mass of rotor inside the journal m j and disk mass m d are considered, as shown in Fig. 2. A static load 2W is applied on the disc.
The equations of motion for the presented system can be presented as: For which m j is the mass of shaft inside and around the journal and m d is the sum of the disc and remaining part of the shaft mass, k s is the shaft stiffness, F x and F y are the bearing forces in x and y directions, respectively, W is the bearing load and x j , y j , x d , y d are journal mass geometrical center and disc geometrical center current positions with respect to the journal steady-state position. m = m d + 2m j is the total mass of the rotor.
Transforming the equations of motion Eq. (40) into dimensionless form, we get the following equations: W . M J and M D are the dimensionless masses of the journal and disc, respectively. Through the current paper, the following mass ratios are considered 2M J = 0.1M and M D = 0.9M, whereM is the dimensionless total mass of the rotor. The dimensionless shaft stiffness is K S . The rotational speed is , and the dimensionless time is τ . The dimensionless static equilibrium forces Substituting Eqs. (22) and (27) in Eq. (41), and converting the resulting equations into the state space form, the final resulting equations can be written as: where Finally, Eq. (42) can be written as: . This is occurred above what is called the threshold speed M ≥M th , where T is the dimensionless periodic time of solution and M th is the dimensionless threshold mass of the system. This will be discussed in detail in the results section.

Results and discussion
In this section, the first-, second-and third-order bearing coefficients are calculated using infinitesimal perturbation analysis. The bearing coefficients are presented and listed for full circular journal bearing. Perturbation analysis is used to investigate the validity range of these coefficients. Afterwards, the bearing coefficients are used to investigate the dynamics of flexible rotor supported by two symmetric journal bearings.

Perturbation analysis
In this section, a perturbation analysis is used to validate the obtained bearing coefficients. This analysis is executed using un-grooved full circular journal bearing with L/D ratio of 1. This is done by applying displacement and velocity perturbations to the equilibrium point and then evaluated the obtained forces using four different methods. These methods are based on direct integration of Reynolds equation (F Re. ), first-order coefficients (F 1 ), second-order coefficients (F 2 ) and third-order coefficients (F 3 ). The component of these forces in X and Y directions is evaluated and plotted in  Fig. 3a, b. However, by increasing the values of perturbation as in Fig. 3c, a significant deviation between the forces evaluated based on the nonlinear coefficients and that based on the direct Reynolds integration forces (F Re. ) is recognized. The last case in Fig. 3d shows that at very low and high Sommerfeld number the deviation between the forces is larger than when S ∈ [0.1 − 0.6]. Also, the x component of force approximation based on the third-order coefficients was closer to the forces based on direct integration of Reynolds equation at S < 0.55 and y component of force approximation was closer to the forces based on direct integration of Reynolds equation when S < 0.3.

Journal bearing coefficients and mesh sensitivity
In this section, the bearing linear and nonlinear coefficients based on first-, second-and third-order coefficients are evaluated. This is done using infinitesimal perturbation as discussed in the analysis section. The presented results are for un-grooved full circular journal bearing with L/D ratio of 1. The linear bearing coefficients are eight and plotted in Fig. 4 a. The nonlinear second-order bearing coefficients are fourteen and are plotted in Fig. 4b . The third-order coefficients are twenty and are plotted in Fig. 5. The effect of mesh size on the bearing coefficients is also investigated graphically for all the coefficients. Several mesh sizes are investigated, but only four cases are presented to avoid over density of the figures. The presented mesh sizes are selected to vary from coarse to fine mesh. These mesh sizes are 50×100, 100×200, 150×600 and 200×800. In both Figs. 4 and 5, the color is used to differentiate the bearing coefficients and the line type is used to differentiate the mesh size as shown in the subfigures legend and the figure general legend. The results of Fig. 4 show that for the plotted coefficients the results of mesh 100 × 400, 150 × 600 and 200 × 800 are identical over the full range of Sommerfeld number. The results of Fig. 5 show that the bearing coefficients obtained using 132 T. A. El-Sayed, H. Sayed mesh 100 × 400 converge with that 200 × 800 in most coefficients. However, for K X X X X , C X X X X , K Y X X X and C Y X X X mesh 100 × 400 is not sufficient. In addition, the figure shows that mesh 150×600 and 200×800 results are converged for all coefficients. Therefore, mesh 150 × 600 was considered as minimum mesh used in coefficients evaluation during the present paper. Numerical values of the obtained parameters using mesh size of 150 × 600 are presented in Tables 2 and 3.

Verification with previous results
Here, the first-, second-and third-order bearing coefficients are compared with the results of Miraskari et al. [27]. However, the reference coefficients are calculated based on Reynolds short bearing equation. Eighteen selected bearing coefficients are used for the comparison. These are divided as follows: six first-order bearing coefficients which are plotted in the first row of  Fig. 6 and six second-order bearing coefficients located in the second row of Fig. 6 and six third-order bearing coefficients plotted at the last row of Fig. 6. The finite bearing coefficients are obtained at L/D = 0.5 to approach the short bearing results, and these results are plotted using solid lines. The results of Miraskari et al. [27] are plotted using dashed lines. It is worth noting that the reference results are processed before plotting because the dimensionless parameters used in the reference paper are not the same with that considered in that paper. Figure 6 shows that the results of both methods are close to each other with a small deviation. This deviation can be explained by the fact that the short bearing Reynolds equation is a simplification to the general Reynolds equation.

Stability analysis
In this section, the calculated bearing coefficients are used to investigate the dynamic analysis of the elastic rotor system presented in Fig. 2. Initially, the stability of the autonomous dynamical system in Eq. (42) is investigated. This is done by evaluating the dimensionless threshold massM th = mc 2 th /W , which implicitly includes the threshold speed of the system. The dimensionless threshold massM th can be evaluated using the eigenvalues of the Jacobian matrix J x,M in Eq. (44) at the fixed points. The fixed points can be Finite Bearing L/D=0.5 Short Bearing L/D=0.5 [Miraskari] obtained by equating Eqs. (43) to zero. This results in x * = 0 0 0 0 0 0 2/K s 0 T .

J x,M
where Evaluating the Jacobian matrix at fixed points x = x * results in Evaluating the Jacobian matrix eigenvalues (λ 1 , λ 2 , · · · , λ 8 ) then reorder them in descending order such as Re (λ 1 ) ≥ Re (λ 2 ) ≥ · · · ≥ Re (λ n ). The value of theM where Re (λ 1 ) = 0 is the threshold valueM th . This analysis is repeated for the available range of Sommerfeld number.
These results are presented in Fig. 7a, b. In Fig. 7 a, the threshold mass is plotted versus the Sommerfeld  Figure 7b is considered as assembly of two-dimensional cross sections in Fig. 7a at four values of K s . To investigate the stability of the periodic solution when the rotor speed exceeds the threshold speed the Monodromy matrix is evaluated from the Floquet multipliers [40,43,44]. For the details of the Monodromy analysis, see Appendix D. If the eigenvalues of Monodromy matrix is less than one, the periodic solution is stable limit cycle, i.e., supercritical Hopf bifurcation, and if the eigenvalue of the Monodromy matrix is greater than one, the periodic solution is unstable, i.e., subcritical Hopf bifurcation. This is plotted in Fig. 7 using solid and dashed lines where the solid lines indi-cate supercritical Hopf bifurcation, and the dashed lines represent subcritical bifurcation solution. The results of the figure show that for the present dynamical model decreasing the stiffness of the rotor shaft increases the possibility of subcritical Hopf bifurcation at large value of Sommerfeld number(small load).

Time response comparative study
In this section, selected points from the stability diagrams I, II, III, IV, V and VI are investigated using orbit plots of the shaft geometrical center. The input data for these cases are shown in Table 1 and presented in Fig. 7b. The comparative analysis is done using four different analyses. The first analysis is based on the linear first-order approximation of the bearing forces (LA). This analysis can be done using Eq. (42) while eliminating the second-and the third-order bearing coefficients. The second analysis is based on the nonlinear second-order approximation of the bearing forces, and this is done while reducing the third-order coefficients to zero (NLA2). The third analysis is based on the nonlinear third-order force approximation (NLA3) as stated in Eq. (42). The last analysis is based on the nonlinear force evaluated directly from Reynolds equation (RNL). In this analysis, the Reynolds equation is solved every time step to evaluate the bearing forces.
The six selected operating points are listed in Table 1. All the cases are for L/D = 1. The first two cases I, II are at S = 0.148 and for K s = 1, the second two cases III, IV are at S = 0.507, K s = 1 and the last two cases V, VI are at S = 0.507 and K s = 20. The cases I, III and V are below the threshold speed, and the cases II, IV and VI are above the threshold speed.
For case I and case III, where both of them are below the threshold speed the four analyses give the same response as shown in the first and third rows of Fig. 8. Case II and case IV are above the threshold speed for supercritical and subcritical cases according to the Monodromy bifurcation analysis. For case II, the response based on the linear analysis shows instability as the rotor hits the journal bearing. Both analyses based on the second-order coefficients (NLA2) and third-order coefficients (NLA3) show a stable limit cycles. However, the analysis based on direct solution of Reynolds (RNL) shows bifurcation leads to unstable limit cycle, then jumps to a larger stable limit cycle. This behavior for journal bearings is previously reported in [45,46].    Fig. 9 Orbit plot for the journal geometrical center at two operating conditions V and VI as shown in Table 1 For case IV which is subcritical according to the Monodromy analysis. The first-order analysis shows unstable behavior above threshold speed, which is basically expected because analysis based on first-order (LA) is insensitive to bifurcation type [40]. In this case, the analyses based on NLA2 and NLA3 do not show the same response. The analysis based on second-order coefficients shows a stable limit cycle, while the analysis based on third-order coefficients shows unstable limit cycle. In this case, the nonlinear analysis RNL agrees with (NLA3) in showing unstable limit cycles. However, the RNL response shows a jump from unstable limit cycle to bigger stable limit cycle as shown in Fig. 8IV-d. The explanation for the difference between the third-order coefficient and RNL can be explained by the fact that the coefficients based on NLA3 are only valid in the neighborhood of the steady-state point and as the rotor receding from this point, these coefficients become not accurate as shown in the perturbation analysis.
Another two cases, V and VI, having higher rotor dimensionless stiffness are investigated, see Fig. 9. The dimensionless stiffness in this condition is K s = 20. The stability diagram shows that all bifurcations for this case are supercritical. The orbit plot for the journal geometrical center is plotted at S = 0.507. The first case V is below the threshold speed, and the second case VI is above the threshold speed. The results show that for case V below the threshold speed all the four analyses show the same response. For case VI which is higher than the threshold speed the analysis based on first-order coefficients shows growing of the response until the journal hits the bearing. The analyses NLA2, NLA3 and RNL show stable limit cycles. The size of the limit cycles is not the same for the three cases. This can be explained by the fact that the bearing coefficients are evaluated to be correct in the vicinity of the steady-state point.

Continuation analysis
In this section, numerical continuation is used to investigate the dynamics of the present model. The continuation analysis is based on the equations of motion based on the third-order nonlinear bearing coefficients. MatCont module is used for the present study which is a toolbox written in MATLAB language and it used for continuation analysis [47]. The continuation analysis is used to investigate the effect of changing the applied static load (Sommerfeld number) and the effect of changing the rotor stiffness as shown in Figs. 10 and 11, respectively. . The first row represents the limit cycle continuation, the second row represents a section in limit cycle continuation plot and the last row represents a top view of the limit cycle continuation plot To apply the continuation analysis, numerical integration is used to identify an equilibrium point below the threshold curve. From this point, continuation analysis is applied at variableM until reaching a Hopf bifurcation point which represent the threshold speed. From evaluating the Lyapunov exponent at this point, the type of Hopf bifurcation can be identified whether supercritical or subcritical. Then, limit cycle continuation is applied to investigate the rotor dynamics behind the threshold speed. MATCONT is used to identify the stability of limit cycles based on the modulus of Floquet multiplier [47]. If at least one multiplier has modulus greater than one, then the limit cycle is unstable. In Figs. 10 and 11, the stable limit cycles are presented by solid lines and the unstable limit cycles are presented by dotted lines.
In Fig. 10, the continuation analysis is applied at four values of eccentricity ratios (a) ε = 0.24 (column one), (b) ε = 0.45 (column two), (c) ε = 0.55 (column three), (d) ε = 0.65 (column four). These values are corresponding to four different static loads with Sommerfeld numbers 0.507, 0.216, 0.148 and 0.0983, respectively, see the four blue vertical dotted lines 7 b. The stiffness ratio of the rotor is K s = 1. Increasing the eccentricity ratio at the equilibrium is associated with the increase in the applied load. From the four selected cases, case a represents subcritical bifurcation atM th = 3.35. The limit cycle continuation in that case shows an unstable limit cycles with the decrease inM. . The first row represents the limit cycle continuation, the second row represents a section in limit cycle continuation plot and the last row represents a top view of the limit cycle continuation plot AtM=3.331, a limit point cycle is occurred (purple circle in Fig. 10) and the limit cycle becomes stable limit cycles. Further increase ofM = 3.333 , another limit point cycle is occurred, and unstable limit cycle continues with the decrease ofM. At cases b and c where ε = 0.45 and ε = 0.55, supercritical bifurcation is occurred followed by stable limit cycle continuation. Further increase ofM results in limit point cycle (purple circle in Fig. 10) atM = 3.664 andM = 4.934 in cases b and c, respectively. After this limit point cycle unstable limit cycles continuation is started. For case d which is supercritical bifurcation, a stable limit cycles are occurred. No limit point cycle is monitored within the investigated range.
In Fig. 11, a continuation analysis is done for four cases having ε = 0.5 and four different values of relative rotor stiffness (a) K s = 1, (b) K s = 5, (c) K s = 10 and (d) K s = 20, see the red dotted vertical dotted line in Fig. 7b. Results in Fig. 11 show the point of Hopf bifurcation for the four cases is atM th = 3.576, 8.673, 10.47 and 11.66, respectively. This indicates that the threshold speed increases with the increase of rotor stiffness ratio. In Fig. 11, although the type of Hopf bifurcation at the four different rotor stiffness ratio values is supercritical and the continuation analysis shows that at K s =1, 5 and 20, limit point cycles are monitored with the increase of the dimensionless massM. These limit point cycles are colored in purple. After these limit point cycles, the continuation analysis shows that the limit cycles change from stable to unstable limit cycles. However, at K s = 10, no limit point cycle is found within the investigated range. The continuation results presented in Fig. 11 show that increasing the rotor flexibility affects greatly the dynamics and stability of the rotor-bearing system.

Conclusion
In the present paper, the linear first and nonlinear second-and third-order bearing coefficients are obtained using infinitesimal perturbation method and finite difference method. The main conclusions from the present work are as follows: • A mesh of 150 × 600 is sufficient to obtain an accurate third-order bearing coefficients. • The perturbation analysis shows that the perturbation force based on the third-order coefficients are more accurate than that based on first-and secondorder coefficients. • The rotor stiffness ratio and applied load are a crucial parameters affecting the type of rotor-bearing stability whether subcritical or supercritical. Additionally, they affect the continuation of limit cycles behind the threshold speed. • The orbit plot results show that the analysis based on the third-order bearing coefficients is closer to the analysis based on the complete solution of Reynolds equation than the analyses based on the first-and second-order bearing coefficients.
At the end of our conclusions, it is recommended to experimentally validate the present analysis with experimental work and to investigate the integration of controller to avoid the rotor-bearing instability.

Funding Information
The authors declare that they do not receive any funds from any organization for this research.

Data availability statement
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
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/.

Appendix A: Derivative of dimensionless pressure gradients
To obtain the first-order pressure gradients the Reynolds Eq. (2) is differentiated with respect to X, Y, X and Y and rearranged as shown in Eq. (A.1). Similarly, the second-order pressure gradients can be obtained by differentiating the Reynolds first derivative equations in Eq. (A.1) with respect to X, Y, X and Y as shown in Eq. (A.2). Likewise, the third-order pressure gradients in Eq. (20) can be obtained by differentiating the second-order derivative of Eq. (A.2) with respect to X, Y, X and Y as shown in Eqs. (A.3-A.12). Equations (A.1-A.12) are discretized using finite difference over the bearing domain, Fig. 1b, and then solved for the first-, second-and third-order pressure gradients. It is worth noting that during the solution of these equations the first-, second-and third-order dimensionless height derivation are required which are listed in Appendix B.

(B.19)
Appendix C: Bearing first-, second-and third-order coefficients Table 2 List of the first-and second-order bearing coefficients for fully circular journal bearing at L/D  Table 3 List of the third-order bearing coefficients for fully circular journal bearing at L/D

Appendix D. Monodromy matrix derivation
In the current manuscript, the type of Hopf bifurcation is investigated by evaluating the Floquet multipliers which are the eigenvalues of the Monodromy matrix. Shooting method is used to evaluate the Monodromy matrix as in [43,44]. In this method the initial value problem is converted to boundary value problem to obtain a periodic solution. The initial conditions with minimal period are searched as shown below.
The periodic solution of Eq. (43) at τ = 0 is for which η is a vector of initial conditions; both the values of vector η and the period T are unknown yet. According to the shooting method [43,44], an initial guess value for η and T is selected such η = η 0 and T = T 0 ; then, the deviations from the actual values can be written as: Substituting Eqs. (D.2) into (D.1) and using firstorder Taylor expansion results in, see ref. [43]. The function ∂x ∂η (T 0 , η 0 ) can be determined with solving the ordinary differential Eqs. (D.7) with initial conditions, Eq. (D.8) coupled with initial value problem, Eq. (43) with the assumed initial conditions as follows: x (0) = η 0 .
(D.9) Equation (D.7) is a 8×8 matrix differential equation which means 64 differential equation. These equations are coupled with the eight differential equations in the dynamical system of Eq. (43). These 72 coupled differential equations should be integrated to τ = T 0 to evaluate ∂x ∂η (T 0 , η 0 ) and x (T 0 , η 0 ). Equation (D.3) has nine unknowns for the system of equation (43), eight for the vector η and the ninth unknown is T , so an additional equation is needed. Therefore, the orthogonality condition f T . η = 0 is imposed to enforce the perturbations η to be normal to the vector field. Therefore, the complete set of equations can be written as Initially, the values of η 0 and T 0 are assumed. Afterwards, Eq. (D.10) can be solved for η κ+1 and T κ+1 . Then, a new guess of η κ = η κ + η κ+1 and T κ = T κ + T κ+1 . For more details, see ref. [43].
The subscript κ in Eq. (D.10) is the iteration number. At each iteration κ , the coefficient of matrices is solved using ordinary differential Eqs. (43) and (D.7). The final solution is converged when the error ( ) is small. At this instant, the matrix ∂x ∂η η k , T k is the Monodromy matrix. Finally, to evaluate the stability of the periodic solution, the eigenvalues of the Monodromy matrix are calculated which are the Floquet multipliers. If the modulus of all eigenvalues is less than or equal to one, ρ j ≤ 1, the bifurcation type is supercritical, else, the system bifurcation is subcritical. The steps for evaluating the bifurcation analysis using Monodromy matrix are summarized in Fig. 12.