A superlinear iteration method for calculation of finite length journal bearing's static equilibrium position

Solving the static equilibrium position is one of the most important parts of dynamic coefficients calculation and further coupled calculation of rotor system. The main contribution of this study is testing the superlinear iteration convergence method—twofold secant method, for the determination of the static equilibrium position of journal bearing with finite length. Essentially, the Reynolds equation for stable motion is solved by the finite difference method and the inner pressure is obtained by the successive over-relaxation iterative method reinforced by the compound Simpson quadrature formula. The accuracy and efficiency of the twofold secant method are higher in comparison with the secant method and dichotomy. The total number of iterative steps required for the twofold secant method are about one-third of the secant method and less than one-eighth of dichotomy for the same equilibrium position. The calculations for equilibrium position and pressure distribution for different bearing length, clearance and rotating speed were done. In the results, the eccentricity presents linear inverse proportional relationship to the attitude angle. The influence of the bearing length, clearance and bearing radius on the load-carrying capacity was also investigated. The results illustrate that larger bearing length, larger radius and smaller clearance are good for the load-carrying capacity of journal bearing. The application of the twofold secant method can greatly reduce the computational time for calculation of the dynamic coefficients and dynamic characteristics of rotor-bearing system with a journal bearing of finite length.


Introduction
A journal bearing is the main supporting component required for the stability and sustainable dynamic characteristics of the rotor system [1][2][3][4][5]. Many researchers are devoted to analysing the lubrication mechanism in a journal bearing. On the basis of fluids continuity equation and viscous fluid motion equation, Reynolds [6] proposed his equation, which laid the theoretical basis of fluid lubrication mechanism for journal bearing. The dynamic characteristics of journal bearing received considerable attention when Newkirk and Taylor found the unstable vibration phenomenon in journal bearing caused by the oil film [7,8]. They called this unstable phenomenon 'oil whip'. Concurrently, Stodola regarded the fluids oil film as a simple spring support, but the model could not have explained the observed finite amplitude of oscillation of a shaft operating at a critical speed [9]. Hagg and Sankey [10,11] represented the dynamic characteristics of a journal bearing by means of two positive stiffness and two positive damping coefficients. Subsequently, Lund and Sternlicht [12][13][14] proposed the eight dynamic coefficients model, which improved the calculating model of a journal bearing and widely adopted in the current calculations of the rotor's dynamic characteristics (table 1).
The introduction of the dynamic characteristic coefficients was greatly convenient for the coupled solution of a rotor-bearing system. Considering the difficulty of direct solution of the analytical Reynolds equation, the short theory was applied for the engineering calculations. Alnefaie [15] researched the startup and steady-state dynamics of a rotor supported by fluid film bearings. In this set-up, the bearing was viewed as short-plain cylindrical bearing, for which different damping ratios caused super-harmonic oscillations. Using the approximation for a short bearing, Chang-Jian and Chen [16,17] obtained the nonlinear bearing force by direct integration, which was used to facilitate the dynamic vibration of nonlinear suspension rotor. Li et al. [18] calculated the dynamical characteristic coefficients of a journal bearing by adopting the theory of a narrow bearing and Gümbel boundary conditions. In this research, the dynamic performance of multi-stage rotor system was also simulated for the varying dynamic coefficients and geometric parameters of journal bearings. Adiletta et al. [19] proposed the nonlinear dynamic model of short journal bearing under the hypothesis that the lubrication film was laminar and isothermal. The model is commonly used for a rotor system because of its good convergence and accuracy [4,[20][21][22][23].
In fact, the short bearing theory was just suitable for the bearing with a small length-diameter ratio. Muszynska and Bently [24,25] proposed a simplified nonlinear fluid dynamics model that considered the circulating velocity as the key factor affecting the dynamic characteristics of a fluid film. Although the model overcame the deficiency of a short bearing compared with other bearing models, it could be adopted only for continuous fluid film, thus it is widely used to describe the nonlinear sealing force in a rotor system [26][27][28]. The finite difference method (FDM) [29,30], partial derivative method (PDM) [31,32] and finite element method (FEM) [33,34] are the common solving methods for eight dynamic characteristics of a journal bearing with finite length. These dynamic coefficients are important for testing the dynamic model of the rotor system with a journal bearing of finite length.
The equilibrium position is the key for solving the dynamic coefficients of a journal bearing with finite length. In order to obtain the dynamic coefficients of this bearing, the equilibrium position has to be found first [35]. Unfortunately, there are various degrees of calculation efficiency problem in the current numerical methods that determine the static equilibrium position. Therefore, in order to improve the computational efficiency of the dynamic coefficients for journal bearings of finite length and the dynamic characteristics of rotor system, in the present study, effort was made to numerically determine the equilibrium position of a journal bearing by the twofold secant method. Moreover, comparative study is done among the results of convergent iteration steps solved by the twofold secant method, secant method and dichotomy to reveal the high efficiency of the proposed superlinear iteration method. The effects of journal length, clearance and rotating speed on equilibrium position and pressure distribution were researched by the twofold secant method. Finally, analysis of the load-carrying capacity of journal bearings with different geometric parameters is also studied.
F r ,F t dimensionless radial and tangential components of dimensionless lubricant film force h max maximum thickness of lubricant medium ( h min minimum thickness of lubricant medium (

2.
Forces in the fluids film of the journal bearing 2. 1. The dimensionless Reynolds equation Figure 1 shows the geometry of a journal bearing and the coordinate system. A journal and a bearing rotate at an angular velocity of Ω j and Ω b , respectively. The motion of the journal's centre point o j is determined by the external vertical load and the pressure in a lubricant film. The flow of lubricant medium between the journal and the bearing obeys the generalized Reynolds equation. In cylindrical coordinates, for cylindrical journal bearing under turbulent conditions, it reads as [36]: 1 The hypothesis of iso-viscous and incompressible Newtonian lubricating liquid in a journal bearing is assumed. Considering the bearing is always fixed on a base, that is can be simplified according to [34,37,38]: Introducing the following dimensionless parameters: As shown in figure 1, the first two items on the right-hand side of equation (2.4) represent the rotating effects of the angular velocity Ω j and dθ/dt; the last term represents the squeezing effect of dε/dt. For the stable motion, dθ/dt = dε/dt = 0, then equation (2.4) can be simplified into: The boundary condition (2.6a) and periodic characteristic (2.6b) for the cylindrical journal bearing of finite length can be concluded as:pz

Forces in the lubricant film
The FDM is applied to solve the dimensionless steady-state Reynolds equation. The discretization of the film domain is usually done in two-dimensional problems [39,40]. The difference grid is shown in figure 2, where the middle point of the adjacent points is the half-integer point.
The first-and second-order derivatives and equation (2.5) are both discretized by the second-order central difference scheme. Then, equation (2.5) takes the following difference form [36]:  n + 1 The resulting discretized form of dimensionless Reynolds equation takes the form: where (2.10) Analysis of equations suggest that the pressure at any grid point (i, j) is expressed in terms of pressure of four adjacent points. In order to obtain the pressure values in every discretized point quickly and accurately, equation (2.5) is solved using the successive over-relaxation iterative (SOR) method [29], supplied with the boundary conditions (2.6a) and periodic characteristic (2.6b).
When the lubricant film pressure distribution is calculated, components (F t ,F r ) of the dimensionless forces in the lubricant film, shown in figure 3, can be calculated using the compound   Simpson quadrature formula: (P 2 sin ϕ 2i ), (P 3 sin ϕ 2i+1 ),T 4 =P 4 sin ϕ m+1 , (2.12) (P 2 cos ϕ 2i ), Furthermore, the forces in the lubricant film and acting angle φ can be obtained [36]: TheF x andF y can be calculated by transformation of the coordinates according to equation (2.11):

Static equilibrium position
In the equilibrium state, the journal bearing of finite length must satisfy the following equilibrium condition: In other words, the equilibrium position can be determined by finding the attitude angle θ and eccentricity ε satisfying equation (3.1). According to the theory of fluid hydraulic lubrication, the only pressure distribution of the lubricant film can be obtained if the geometric parameters and eccentricity ε of a journal bearing used in equation (2.5) are given. The attitude angle θ adjusts the pressure distribution in the circumferential direction. Therefore, one must find θ and ε for the known external vertical load and rotating speed.
In the present work, the twofold secant method for solving nonlinear equation is applied to efficiently and accurately obtain the attitude angle θ and eccentricity ε. Compared to the secant method and bisection method, it has higher computational efficiency. If the initial values (θ 0 , θ 0 ) and (ε 0 , ε 0 ) are given, the iterative format of twofold secant method [41] for the static equilibrium position reads as follows: The convergence criteria of equation (3.2a) [42] and equation (3.2b) give: and There are two independent variables in the iterative process; therefore, it is impossible to acquire the attitude angle and eccentricity simultaneously. Thus, the iterative process of attitude angle is embedded in the iterative process for eccentricity. Moreover, the equilibrium position depends on the rotating speed [15], i.e. it changes for the different rotating speed. Therefore, the equilibrium position can be determined

Results and discussion
The static equilibrium position of a journal bearing was obtained by the twofold secant method. The used calculating parameters of the journal bearing are shown in table 2 ( figure 1 for reference). The convergence process, the inner pressure distribution and the load-carrying capacity were presented using the superlinear iteration method. The influence of the geometric parameters and the working conditions on the characteristics of journal bearing were also analysed.

Comparison of three convergence methods
In this section, in order to present the advantage of the twofold secant method for identification of the static equilibrium position, the convergence process and iterative steps of three methods are shown for different length L b , radial clearance c and rotating speed Ω j . Considering the accuracy and efficiency of the dichotomy and the scant method for nonlinear equations [42][43][44], these two methods were selected for calculation of the equilibrium position, and the results were compared with the results obtained by the twofold secant method. Figure 5 presents the process of identification of the journal's equilibrium position when journal length was set to 0.1 m. The other parameters of the journal bearing are listed in table 2. It can be seen that the convergence curve of the secant method is most twisty and its iterative range is also largest compared with those of twofold secant and dichotomy methods. There is only one fold point for the twofold secant method and the dichotomy method in the seeking process. Thus, these two methods effectively reduce the seeking scope at the beginning of the iteration process. However, for the dichotomy approach, the high density of the points shows poor convergence performance near the equilibrium position. Moreover, all three curves converged to the same equilibrium position point (ε = 0.489, θ = 59°). As this figure shows, the calculating efficiency of the twofold secant method is superior to the secant method and dichotomy. Figure 6 represents the convergence process of attitude angle corresponding to the different iterative number of eccentricity shown in figure 5. It can be seen from figure 6a-c that the iterative steps of attitude angle, i.e. the number of iterations in bottom loop, increase from the twofold secant method to dichotomy. The maximum numbers of iterative steps of the twofold secant method, secant method and dichotomy are 2, 3 10, respectively. Therefore, the convergence efficiency for attitude angle decreases from the twofold secant method to dichotomy. Also, compared with the iterative process of middle loop and bottom loop using the secant method, the numbers of iterations for attitude angle are much less than those for eccentricity. This is because the convergence efficiency for the secant method is mainly affected by the monotonicity of iterative results.
The convergence process for the equilibrium position when the clearance is 0.8 mm and the rotating speed is 2000 r min −1 is plotted in figures 7 and 8, respectively. From the two figures, it is clearly shown that the convergence trajectory for all three methods is nearly the same. The convergence curve of the secant method has more folds than the other two curves, and the convergence efficiency of the twofold secant method is also higher than that of the secant method and dichotomy. Besides, for the dichotomy, when the iterative values approach the point of convergence, the convergence speed becomes slow sharply. This accounts for the low convergence efficiency of the dichotomy for equilibrium position compared with the other two methods.

The effects of L b , c and Ω j on the equilibrium position
The numerical results calculated by the twofold secant method are compared with those of previous reference [45]. The comparative results are listed in table 5. It is evident through this table that the equilibrium position and lubrication film force calculated by the superlinear iteration method are consistent with those of previous reference. The maximum relative error and minimum relative error are 3.14% and 0.06%, respectively. The small difference of the numerical results implies that the twofold secant method proposed in the paper is accurate and feasible. In addition, in order to investigate the influence of geometric parameters and working condition on the equilibrium position, the bearing length L b , clearance c and rotating speed Ω j were selected as the independent variables to be used in further calculations by the twofold secant method.   The parameters used in the calculation of different equilibrium positions, shown in figure 9a, are listed in table 2 without the bearing length L b . It can be seen that the equilibrium position gradually approaches the geometric centre of the bearing when the length increases from 0.1 to 0.2 m. This reflects on the eccentricity and attitude angle-ε decreases from 0.489 to 0.124, while θ increases from 58.95°to 83. 43°. This phenomenon illustrates that higher bearing length facilitates the load-carrying capacity of the journal bearing. This fact is consistent with the results of Gengyuan et al. [46].
The corresponding fitting curve of the eccentricity and the attitude angle is plotted on figure 9b. The straight line, fitted by the least square method, indicates that the relationship of ε and θ is linear to the change of the bearing length. Figure 10 presents the distribution curves of the dimensionless pressure for corresponding equilibrium positions shown in figure 9 (z = 0). Figure 10 obviously shows that the peak pressure  The pressure distribution in the entire flow field clearly shows that the positive pressure emerges in the convergent wedge (ϕ = 0°→ 180°). Theoretically, once the lubrication fluid crosses the minimum clearance and enters the divergence wedge (ϕ = 180°→ 360°), the pressure quickly decreases and the values change from positive to negative on account of the expansion effect.
In fact, the actual lubrication fluid cannot bear the tensile stress and the liquid film will fracture in the divergence wedge, this is why half Sommerfeld boundary condition was widely applied in the   calculations of Reynolds equation [38,[47][48][49]. In addition, the pressure is distributed symmetrically on the left and right side of central plane (z = 0) in the axial direction.
The changing trajectory of equilibrium positions obtained by the twofold secant method and the corresponding fitting curve for different clearances and rotating speeds are shown in figures 12 and 13, respectively. Similar to figure 9a, the changing trajectory presents parabolic characteristics. The difference in the three trajectory curves could be explained by the fact that the equilibrium position is more sensitive for both bearing length and clearance, and less sensitive for rotating speed. This is shown in the changing range of equilibrium positions: in figures 9a and 12a it is larger than that shown in figure 13a. It also can be seen that with the increase of clearance, the eccentricity increases from 0.119 to 0.513. By contrast, the eccentricity decreases from 0.472 to 0.279 as the rotating speed changes from 2000 to 4500 r min −1 . The calculated results imply that smaller clearance and higher rotating speed is beneficial to load-carrying capacity. Figures 12b and 13b, respectively, present the relationship between the eccentricity and the attitude angle for different clearances and rotating speeds. Similarly, the negative slope of the two fitting curves demonstrates that eccentricity and attitude angle are in the inverse proportional relationship. Whereas compared with the slope of the fitting curve for bearing length shown in figure 9b, the slopes of clearance and rotating speed change from −68.17 to −60.94 and −61.57, respectively. These findings mean that the variation of attitude angle caused by the bearing length is greater than that of the clearance or the rotating speed for the same variation of the eccentricity.
The  these findings with the conclusion obtained from figure 10, it is clear that the larger eccentricity can cause higher pressure under the same geometric parameters and working condition. The angular position of the maximum pressure decreases with reduction of the clearance or rise of the rotating speed. Moreover, the increments for the maximum pressure and angle position increase as the bearing length decreases or the clearance increases, as shown in figures 10 and 14, respectively.

The effects of L b , c and R b on the forces in lubricant film
The geometric parameters, including L b , c and R b , play an important role in lubricant film force, directly affect the load-carrying capacity and the equilibrium position of the journal bearing and finally determine the calculated dynamic coefficients. In order to determine the equilibrium position using the proposed model, the lubricant film force needs to be transformed into dimensionless form [47,50]. In this section, the main geometric parameters L b , c and R b are taken into consideration for investigation of the lubricant film forces for different eccentricities. Figure 16a-c presents the changing trend for the radial film force F r , the tangential film force F t and the total lubricant film force F, respectively, for the different bearing lengths and eccentricity changing from 0.1 to 0. 9. The negative values of F r illustrate that the actual acting direction of the radial film force is opposite to the direction shown in figure 3. The rising curves given in figure 16a-c also demonstrate that when the eccentricity is fixed, the values of film force with larger L b are obviously greater than those with smaller L b ; the same results are also true for the fixed L b and changing eccentricity. That is, larger bearing length and eccentricity have a positive role in improving the load-carrying capacity of journal bearing. In addition, the curves of the film's forces can be roughly divided into nonlinear region and linear region at L b = 0.15 m on the basis of the slope change. The relationship of the acting angle and the bearing length with the different eccentricities is shown in figure 16d. Actually, the changing trend of the acting angle reflects the F t -F ratio according to the equation (2.16). As mentioned in figure 16d,    The aim of table 7 is to present the relationship between the bearing radius and the lubricant film force under different eccentricities. In this case, the forces grow with fixed eccentricity when the bearing radius increases from 0.05 to 0.15 m. This fact means that the journal bearing with larger radius (R b = 0.15 m) has better load-carrying capacity compared with those of smaller radius (R b = 0.05 m).

Conclusion
A novel superlinear iteration convergence method-twofold secant method, was applied for identification of the equilibrium position of journal bearing of finite length. The Reynolds equation for stable motion was solved numerically and effects of the geometric parameters and the operating condition on equilibrium position, pressure distribution and load-carrying capacity were studied. Main conclusions are as follows: (1) The number of iterative steps of the proposed method-twofold secant method-required for obtaining equilibrium position are obviously less than for the present methods: secant method and dichotomy. The efficiency of calculation of dynamic coefficients and dynamic characteristics of journal bearing with finite length using the twofold secant method is faster. (2) The trajectory of equilibrium positions is parabolic and the relationship curves between attitude angle θ and eccentricity ε fitted by the least square method indicates that the attitude angle linearly decreases when the eccentricity increases. In addition, a variation of the attitude angle is more sensitive to the bearing length than to the clearance or the rotating speed. (3) The inner pressure is distributed symmetrically in the axial direction on the left and right side relative to the central plane (z = 0). The large bearing length (L b = 0.2 m), higher rotating speed (Ω j = 3500 r min −1 ) and small clearance (c = 0.3 mm) shift the angle position of maximum pressure away from the position of minimum clearance between journal and bearing for the same external vertical loads. (4) The geometric parameters, including L b , c and R b , have a major impact on the load-carrying capacity. The load-carrying capacity can be enhanced by enlarging the bearing length, its radius or reducing the clearance. Besides, the curves of load-carrying capacity can be divided into nonlinear region and linear region according to the slopes as the eccentricity increases from 0.1 to 0. 9.