Dynamic characteristics of water-lubricated journal bearings

The increasing ecological awareness and stringent requirements for environmental protection have led to the development of water lubricated journal bearings. For the investigation of water-lubricated journal bearings, a new structured mesh movement algorithm for the CFD model is developed and based on this method, the nonlinear transient hydrodynamic force model is established. Then, with consideration of velocity perturbation, a method to determine dynamic coefficients and linear hydrodynamic forces is promoted. After validation of static equilibrium position and stiffness coefficients, a comparative linear and nonlinear hydrodynamic force analysis of multiple grooves water-lubricated journal bearings (MGWJBs) is conducted. The calculation results indicate that the structured mesh movement algorithm is suitable for the dynamic characteristics investigation of water-lubricated journal bearings. And the comparative study shows that there is a considerable difference between the linear and nonlinear hydrodynamic forces of MGWJBs. Further investigation should be carried to evaluate the dynamic responses of rotor supported by MGWJBs under difference force models.


Introduction
The increasing ecological awareness and stringent requirements for environmental protection have led to the development of water lubricated journal bearings in many applications where oil was used as the lubricant [1]. The key applications of water-lubricated journal bearing include marine propulsion systems, water pumps and hydraulic turbines [2]. To facilitate the hydrodynamic process, the constant supply of water is commonly fed into the bearing clearance through longitudinal grooves to maintain the thin fluid film between the journal and bearing.
Commonly, the dynamic characteristics of a rotor supported by water lubricated journal bearings include unbalance response and whirling instability. Those characteristics can be investigated by two methods: the linearized method and nonlinear transient analysis. For the linearized method, a small perturbation is given to the journal around static equilibrium position and the stiffness and damping coefficients can be determined.
Then, the unbalance response and stability threshold speed are related to those coefficients by solving the equations of motion. Using the linearized method, Stability analysis of MGWJBs and hybrid water-lubricated journal bearings have been investigated by Majumdar et al. [3] and Ren et al. [4]. Nonlinear transient analysis model of the rotor-journal bearings system gives the orbital trajectory within the clearance circle by solving the equations of hydrodynamic lubrication and motions iteratively. This model can provide a better understanding of transient flow field and stability of rotor. For transient analysis model of MGWJBs, Pai et al. [5,6] solved the time-dependent Reynolds equation and equations of motion to predict the transient behavior of the rotor. However, a comparative linear and nonlinear dynamic analysis of cylindrical journal bearing has shown that a linear model is found to deliver acceptable results at a relatively small shaft unbalance [7]. Meanwhile, linear models can predict the imbalance response only when operating conditions are below the instability threshold speed at high eccentricities [8]. In addition, the difference of linear and nonlinear hydrodynamic forces is relative to bearing structures [9]. Therefore, understanding the difference between linear and nonlinear hydrodynamic forces of water-lubricated journal bearings is an important issue that can reveal the reason for the inherent stability and provide more precise and detailed prediction of bearing behaviors in the future.
Up to now, there are two possible modeling approaches for numerical solution of hydrodynamic flow issue of waterlubricated journal bearings: on the one hand, as listed above, there are models based on classical Reynolds equation [3][4][5][6][7][8][9]. However, low viscosity of water compared to oil increases Reynolds number drastically and makes these bearings prone to significant fluid inertia effects. Dousi et al. [10][11][12] highlighted that the inertia, neglected by classical Reynolds equation, has a considerable effect on the dynamic coefficients and stability. To eliminate this limitation, on the other hand, there are the models based on the extended Reynolds equation considering inertia effect or CFD programs solving full Navier-Stokes equations. Although, extending Reynolds equation has the advantage of short computation time, it may need more work to meet requirements when a complex flow geometry is used or when a more detailed analysis is required. With the rapid development of computer technology, more and more researchers use CFD programs to predict the performance of water-lubricated journal bearings.
The CFD method has already been proved to be a very useful tool in the lubrication analysis of water-lubricated journal bearings. Cabrera et al. [13] constructed a twodimensional CFD model and the calculation results showed similarity in experimental results. Gao et al. [14,15] established a 3D CFD model to analysis load-carrying capacity of water-lubricated journal bearing under hydrodynamic lubrication condition. Based on a same method, Zhang et al. [16] provide a method for determining the stiffness coefficients of hydrodynamic plain journal bearings lubricated by water. And then, the load carrying capacity of misaligned condition was analyzed [17]. Using CFD and FSI method, an elastohydrodynamic (EHD) model was established for water-lubricated journal bearings by Wang et al. [18]. Then, this method was used to analyze the performance of worn hydrodynamic waterlubricated plain journal bearings [19].
However, for the transient or dynamic characteristic investigation of journal bearings using CFD techniques, only few studies have been carried on. Guo et al. [20]. developed a CFD model for dynamic coefficients. Gertzos et al. [21] and Cheqamahi et al. [22] proposed a transient analysis model of journal bearing lubricated by Bingham lubricant and turbocharger's bearing using dynamic mesh method, respectively. Liu et al. [23] and Lin et al. [24] established the model for oil lubrication journal bearings using CFD-FSI methodology. However, due to the magnitude difference between the dimensions of clearance and journal bearing, using of unstructured grids and dynamic mesh models in CFD will stop the calculation process due to numerical failure or the negative volume of grids in the transient condition. Meanwhile, the viscosity of the water is about 30 ∼ 40 times less than mineral oils [1]. It contributes to a lower hydrodynamic load carrying capacity of a water-lubricated plain bearing. And a larger eccentricity ratio such as 0.6 or 0.7 is regarded as a proper value [15]. This also resulted in a great increase in the difficulty of nonlinear transient CFD analysis model.
Recent years, Li et al. promoted a structured mesh movement algorithm to calculate the 3D transient flow of classic circular journal bearings [25]. Then, this method is revised to study the transient flow field of tilting pad journal bearings [26] and the self-circulating oil bearing [27,28]. However, the water-lubricated bearings have not been covered.
In the present work, two weakness of the CFD model of water-lubricated journal bearings are overcome. Firstly, a new structured mesh movement algorithm for the CFD model of water-lubricated journal bearings is developed and based on this method, the nonlinear transient hydrodynamic force model is established. Secondly, with consideration of velocity perturbation, an efficient method to determine dynamic coefficients and linear hydrodynamic forces is promoted based on the calculation of 3D transient flow field. Figure 1 shows assembly diagram of multiple grooves water-lubricated journal bearing. Commonly, the bearing is submersed in water and water is fed from one end through longitudinal grooves. Then, hydrodynamic pressure is generated in the non-grooved region to maintain a thin water film between shaft and bearing. The 3D diagram of the fluid domain including grooves domain and clearance domain is shown in Figure 2. When steady states reached (Fig. 3), the journal is displaced from the bearing with a center distance (e), which is referred to the journal eccentricity. The parameters of the journal bearing used in the numerical analysis are shown in Table 1.

Physical and meshing models of a bearing
To obtain a structured grid distribution at clearance domain, the grooves is split by the outer circular surface keeping connection to entire fluid domain (Fig. 2). The clearance fluid domain could be meshed with a hexahedral type element. And the grooves domain is meshed using "cooper" method. Divisions used across the film thickness are 10 [16], and the interval size used in the circumferential and axial direction is 0.5. The grid distribution at one of the grooves and clearance domain are shown in Figure 4.

Structured mesh movement algorithm
Transient simulations depend on the grid quality. Using unstructured grids and the existing dynamic mesh models will stop the calculation process due to numerical failure or the negative volume of grids as mentioned before. A new mesh moving method based on structured grids stated above is proposed. The following process is used to update the generated mesh: -find the coordinate (x pt , y pt ) of each grid p at time step t and determine which fluid domain p belongs to: ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi > R j grooves domain: The nodes in the grooves domain need not to be updated because the grooves domain are treated as rigid and the nodes position would not change with the time step. And the journal is moving the under the action of hydrodynamic force, so the nodes in the clearance domain would be updated; calculate the displacement of journal center (Dx jt , Dy jt ) at time step t by solving equations of motion and store the coordinate in memory temporarily; judge which the radial reticulate layer N i each node p at clearance domain belongs to at time step t; calculate the position (x pt+1 , y pt+1 ) of each node p at clearance domain at next time step t + 1 according to which layer it belongs to: where N total is the total number of radial reticulate layers. Thus completing the motion updating.
To describe the capabilities of the method clearly, Figure 5 models a simplified MGWJBs with larger clearance than in reality. It can be concluded that, at a larger eccentricity ratio, the existing mesh movement algorithm will result poor mesh quality, even stop calculating. And the grid number is decreased dramatically, which will produce a considerable calculation error. The above-mentioned mesh movement algorithm is based on structured mesh distribution, shown as Figure 5e. The grids after updating is still keeping a uniform distribution and the grids number is unchanged. The updating process is stopped manually without any mesh distortion or numerical failure. Therefore, the structured mesh movement algorithm proposed in this paper is suitable for the transient flow in MGWJBs.

Assumption and governing equations
A rigid aligned bearing with the geometry of Figure 2 is considered. The flow is assumed as laminar and water properties are treated as constant. The multiphase flow of the lubricate with cavitation is described by "mixture model". The "mixture model" solves the continuity equation and momentum equation for the mixture, and the volume fraction equation for the secondary phases.
The continuity equation for the mixture is: where v m is the mass averaged velocity: and r m is the mixture density: where a k is the volume fraction of phase k. The momentum equation for the mixture can be expressed as: where n is the number of phases,F is a body force, and m m is the viscosity of the mixture: where v dr;k is the drift velocity for secondary phase k: Volume fraction equation for secondary phase k can be derived from the corresponding continuity equation: In a water-lubricated journal bearing, the clearance between journal and bearing is not uniform because of eccentricity of the shaft. Where the clearance is divergent, a subambient region may be formed. Once the pressure drops below vaporization pressure, a gaseous phase begins to fill the divergent region, then, cavitation occurs. In the current study, the mass transfer between the liquid and vapor lubricate is based on the "full cavitation model" [29], which treats the mass transfer as source term of equation (9). The transfer rate is calculated as follows: ð10Þ A numerical solution of equations (3)-(11) using a finite volume method (FVM) with a pressure-implicit with splitting of operators (PISO) scheme for pressure-velocity coupling equation satisfying the boundary conditions gives the pressure distribution. A time-dependent solution "unsteady" is choose to model the transient state of lubrication filed when the rotor is moving. Double precision calculations are employed to avoid the negative influence of the large aspect ratio. The PRESTO! scheme is used for pressure interpolation because of high-speed rotating flows in the journal bearing. In addition, other convection terms are discretized by the first-order upwind method to improve convergence rate.
The boundary conditions of the inlet and outlet are respectively "pressure inlet" and "pressure outlet" with 50 kPa and 42 kPa. Outer surface is defined as a stationary wall. Inner journal surface is defined as a moving wall with an absolute rotational speed. And the rotational speed is defined as follows: where N is rotational speed.

Linear and nonlinear transient hydrodynamic forces
For a small amplitude motion Dx; Dy; D _ x; D _ y ð Þ , the expansion of hydrodynamic force around the equilibrium position in a Taylor's series, with only the linear terms retained, can be linearized as the following form: where F x and F y are the hydrodynamic force components in the xand y-axis directions, respectively; F x 0 and F y 0 are the force components at the static equilibrium position.    Those dynamic coefficients have significant effects on the stability of rotor systems and to our knowledge, are still a challenge to experimental determined. Therefore, it is necessary to develop an efficient numerical method. In the reference [16], a method to the determination of stiffness coefficients of hydrodynamic water-lubricated plain journal bearing is proposed by neglecting velocity perturbation. And the hydrodynamic force is linearized as: However, it is difficult to evaluate the influence of velocity perturbation on the dynamic characteristic. Therefore, a linear velocity perturbation was used to determinate the full dynamic coefficients in this study. Firstly, the journal is moved from the static equilibrium position in x-direction and y-direction with the velocity perturbation of D _ x and D _ y, respectively. The corresponding hydrodynamic forces can be written as: ( F x2 ¼ F x0 À k xy Dy À c xy D _ y F y2 ¼ F y0 À k yy Dy À c yy D _ y y-direction: ð16Þ Then, equations (15) and (16) were fitted by linear relation: Finally, the dynamic coefficients can be calculated as: c xx c xy c yx c yy However, linear and nonlinear forces are in a good agreement for small journal amplitudes only. The transient hydrodynamic force obtained using linearized dynamic coefficients is not accurate enough at large journal amplitudes. Solving the equations of motion (Eqs. (21)-(23)) successively, the acceleration (€ x t , € y t ), velocity ( _ x t , _ y t ) and displacement (x t , y t ) can be calculated at each time step. And (x t , y t ) and ( _ x t , _ y t ) are set equal to zero at initial time step.
Then, the nonlinear transient hydrodynamics forces can be directly computed from: The flow chart of the numerical procedure is described in Figure 6. A CFD package FLUENT was used to analyze flow field and the structured mesh movement algorithm was achieved by UDFs. The hydrodynamic force was fitted with linear relationship in MATLAB.

Validation of CFD model and structured mesh movement algorithm
To facilitate the hydrodynamic process, the constant supply of water is commonly fed into the bearing clearance through longitudinal grooves to maintain the thin fluid film  between the journal and bearing. However, the dynamic mesh zone, which is called as clearance domain in Figure 2, has a similar working condition and flow state with plain water-lubricated journal bearing. Meanwhile, the plain water lubricated journal bearings is relatively full investigated by experimental measurements and numerical simulations. So, validation is carried for a plain water lubricated journal bearing. The detailed description of apparatus is shown in reference [15].
Calculation is carried out in the transient state with a constant load (e = 0). Under the action of nonlinear hydrodynamic force, the journal is moved from the concentric position of the bearing and finally stabilized at a certain position, which is called static equilibrium position and is shown Figure 7. With increasing of speed, the static equilibrium position is coming closer to concentric position of the bearing. The calculation results (expressed as "CFD") are compared with the experimental data (expressed as "Exp.") for a validation about nonlinear hydrodynamic force. They are all in good agreement as shown in Table 2.
The pressure and vapor phase contours at static equilibrium position are shown in Figures 8 and 9, respectively. It is obvious that, with the increasing of speed, the maximum pressure is decreasing and positive        area is increasing. The cavitation region is narrower and the gaseous phase is gathered in the rotational direction, which are in good agreement with references [30,31]. In summary, it can be concluded that the 3D transient flow calculation of water-lubricate bearings, based on structured mesh movement algorithm, is reliable.

Validation of dynamic coefficients determination method
As shown in Figure 10, the perturbation of velocity can lead to the increasing of hydrodynamic force by damping effect. The reference [16] neglected this damping effect. As the rotational speed is increasing, the narrower difference implies that the damping effect is weaken. However, as key application of water-lubricated journal bearing, the rotational speed of the marine propulsion systems is commonly around 1000 r/min. Therefore, the velocity perturbation will have a considerable influence of the dynamic characteristic of water-lubricated journal bearings. Although, a larger perturbation velocity can result a more obvious damping effect, the influence of solution initialization value is more obvious, as shown in Figures 11  and 12. This influence will lead to the linearization of hydrodynamic force more away from the given static equilibrium position. Because the dynamic coefficients are varied with static equilibrium position, a larger perturbation velocity will bring errors to the determination of those coefficients.
By neglecting the initial value, hydrodynamic forces are fitted by linear relationship. For a comparison, initial value at the displace perturbation 0.0-1.0 mm is neglected. And fitted by linear relation, the hydrodynamic force under different perturbation velocity is shown in Figures 13 and  14. The square of correlation coefficient is keeping larger than 0.9990, which means the fitting is relatively precise. Using the equations (19) and (20), the full dynamic coefficients can be determined in Table 3. Under different perturbation velocity, the difference of dynamic coefficients, which should be the same in theory, may be the results of initial value.    After a balancing consideration of damping effect and the influence of solution initial value, 2eÀ4 m/s is chosen as the perturbation velocity. The comparison with the reference [16], in terms of stiffness coefficients, is shown in Figure 15. For a constant eccentricity ratio, the stiffness coefficients are proportional to the rotational speed and show a good agreement with the reference value. The difference may be resulted by different pressure-velocity coupling scheme. SIMPLEC is used for the reference and PISO for present work. Because PISO is more stable for transient calculation. Another validation is carried between linear and nonlinear hydrodynamic forces and is shown in Figure 16. Linear and nonlinear hydrodynamic forces are calculated by equations (13) and (24), respectively. When the rotor is whiling in a circular with a radius of 1 mm, a good agreement is achieved between linear and nonlinear hydrodynamic force model.
However, linear and nonlinear forces are in a good agreement for small journal amplitudes. The transient hydrodynamic force obtained using linearized dynamic coefficients is not accurate enough at large journal amplitudes, as shown in Figures 17 and 18. Therefore, in the design phase, a designer should keep in mind that the inherent nonlinear effect of the hydrodynamic force of the journal at large journal amplitudes has a considerable influence for designing an efficient water-lubricated journal bearing.

Comparison of two hydrodynamic force models of MGWJBs
Traditionally, in the investigation of the dynamic characteristic of MGWJBs, the linear hydrodynamic forces which are calculated from dynamic coefficient are used [6,32]. However, the nonlinear hydrodynamic force gives a more detailed understanding of the dynamic behavior of MGWJBs. To evaluate the influence of nonlinear effect of hydrodynamic force for the rotors supported by MGWJBs, further comparisons of hydrodynamic forces were explored for four-axial grooves and six-axial grooves. The pressure contours of bearing surface at the static equilibrium positon are shown in Figure 19. For the comparative study of linear and nonlinear hydrodynamic forces of MGWJBs, the motion of the journal is set to circles with different radiuses to simulate whirling effects. The whirling speed is equal to the journal rotational speed.
During the whirling process, the variation of the hydrodynamic forces in time domain was shown in Figures 20 and 21. By conducting Fourier frequency transform (FFT), the spectrogram of the hydrodynamic forces can be seen in Figures 22 and 23: as shown in time domain (Figs. 20 and 21), there is a considerable amplitude difference between the linear and nonlinear hydrodynamic forces, even with a smaller whirling radius (R = 0.5 mm). However, there is no phase difference. Both hydrodynamic force models achieved the maximum and minimum point almost at the same time; as shown in frequency domain (Figs. 22 and 23), the frequency of the linear and nonlinear is equal to whirling frequency. However, the amplitude difference between nonlinear and linear model in y-direction is larger than that in x-direction. With the increasing of orbit size, amplitude of linear hydrodynamic force is much larger than nonlinear hydrodynamic force, especially for y-direction.

Conclusion
In the present work, two weakness of the CFD model of water-lubricated journal bearings are overcome. Firstly, a new structured mesh movement algorithm for the CFD model of multiple grooves water-lubricated journal bearings is developed and based on this method, the nonlinear transient hydrodynamic force model was established. Secondly, with consideration of velocity perturbation, an efficient method to determine dynamic coefficients and linear hydrodynamic forces is promoted based on the calculation of 3D transient flow field. After validation of static equilibrium position and stiffness coefficients, a  comparative linear and nonlinear hydrodynamic force analysis of multiple grooves water-lubricated journal bearings is conducted. The conclusions drawn from the study can be summarized as follows: a self-developed structured mesh movement algorithm is used to insure the high quality of the grid in the transient calculation process of water-lubricated journal bearings. The grids after updating is still keeping a uniform distribution and the grids number is unchanged. The updating process is stopped manually without any mesh distortion or numerical failure, showing the structured mesh movement algorithm proposed in this paper is suitable for the transient flow of water-lubricated journal bearings; perturbation of velocity can lead to the increasing of hydrodynamic force by damping effect, which was neglected by previous study. Based on the transient flow calculation, a method to determine the stiffness and damping coefficients is promoted using linear fitting. The stiffness coefficients show a good agreement with the reference value; there is a considerable difference between the linear and nonlinear hydrodynamic forces, even with a smaller whirling radius (R = 0.5 mm). And with the increasing in the orbit size, amplitude of linear hydrodynamic force is much larger than nonlinear hydrodynamic force, especially for y-direction.
Although the transient hydrodynamic force models are established and a comparative study is carried out, the dynamic responses of rotor supported by MGWJBs under difference force models still have not been investigated. The proposed method will be employed for further studies of the dynamic characteristic of grooves water-lubricated journal bearings-rotor system.