Dynamic Characteristics of High-speed Water-Lubricated Spiral Groove Thrust Bearing Based on Turbulent Cavitating Flow Lubrication Model

The water-lubricated bearings are usually the state of turbulent cavitating flow under high-speed conditions. And the distribution of cavitation bubbles and the interface effect between the two phases have not been included in previous studies on high-speed water-lubricated bearings. In order to study the influence of interface effect and cavitation bubble distribution on the dynamic characteristics of high-speed water-lubricated spiral groove thrust bearings (SGTB).A turbulent cavitating flow lubrication model based on two-phase fluid and population balance equation of bubbles was established in this paper. Stiffness and the damping coefficients of the SGTB were calculated using the perturbation pressure equations. An experimental apparatus was developed to verify the theoretical model. Simulating and experimental results show that the small-sized bubbles tend to generate in the turbulent cavitating flow when at a high rotary speed, and the bubbles mainly locate at the edges of the spiral groove. The simulating results also show that the direct stiffness coefficients are increased due to cavitation effect, and cross stiffness coefficients and damping coefficients are hardly affected by the cavitation effect. Turbulent effect on the dynamic characteristics of SGTB is much stronger than the cavitating effect


Introduction
Water-lubricated bearing has been widely applied due to its low viscosity, low temperature rise, no pollution and etc. Spiral groove thrust bearing (SGTB) has advantages in hydrodynamic effect and stability. Nowadays, the water-lubricated spiral groove thrust bearing has been used in high-speed spindle. However, since inception cavitation number of water is larger than that of oil, cavitation phenomenon in water-lubricated bearings is more serious than the oil-lubricated bearing at high speed, and the water-lubricated bearing tends to fall into turbulent fluid condition when rotary speed exceeds a certain value.Hence,both cavitating effect and turbulent effect should be considered to model the high-speed waterlubricated bearing.
The cavitating and turbulent effects of water or oil bearings with different configurations have been investigated in previous studies. In the early studies on bearing cavitation, the different boundary conditions for lubrication equation including the Swift-Stiebe condition [1], the JFO condition [2], the Elrod condition [3] have been adopted to describe the cavitation effect of the lubricating film. With those conditions, solution of a two-phase flow in cavitated region can be avoided. However, the cavitating flow is not involved in those models. To overcome this problem,the cavitation lubricant is treated as two-phase mixed fluid. The two-phase mixed fluid models have been developed, those models can be divided into three types based on calculation approaches of gas phase volume fraction: (1) the model based on the R-P equation [4,5] (2) the model based on gas solubility and surface tension of bubble [4,[6][7][8][9] (3) the model based on transport equation of the gas phase volume fraction [10][11][12][13][14][15].
In addition, with the rapid development of computational fluid dynamics (CFD) technology,some researchers have proposed performance analysis of bearings with the cavitation using the CFD software [7,11,[16][17][18][19]. However, the interfacial effect (the momentum and energy transfer effect between phases at the interface of two phases) of the two-phase flow is not included in the above lubrication models of cavitating flow. Actually, the momentum and energy exchange between the two phases produces at the interface, so the interface effect is a very important feature for two-phase flow,which can not be ignored when modeling the high-speed water-lubricated bearing.
In addition, the size distribution of the bubbles cannot be obtained using the present mixed fluid model or the CFD model. In our previous study, lubrication models considering the interface effect of two-phase flow have been proposed [20][21][22][23]. However, those models are established under heat isolation, and the bubble population equilibrium equation is solved directly. Since the bubble size probability density function is an unknown quantity, so it is difficult to determine the bubble size probability density function at the initial moment and at the interface. Furthermore, this kind of algorithm tends to unstable. This study aims to establish a lubrication model for the water-lubricated spiral groove thrust bearings including both two-phase interfacial and fluid turbulent and bearing heat conduction effects. A generalized Reynolds equation considering the cavitation interface and turbulence, and an energy equation considering bearing heat conduction are derived. A population balance equation considering the breakage and coalescence of bubbles is introduced to describe the evolution of the size distribution of bubbles. A discrete interval method is employed to solve the bubble population balance equation. The influence of cavitation and turbulence on the dynamic characteristics of the SGTB is analyzed using the established model. The size distribution of cavitation bubbles and stiffness coefficient of water film are measured using a self-developed experimental device. A comparison of the simulating results with experimental ones is given to verify the proposed model.

Generalized Reynolds Equation Based on the Two-Phase Flow Theory
This paper studies the water-lubricated SGTB with the turbulent cavitating flow film ( Figure 1). The classical Reynolds equation is not suitable for describing the bearing lubrication state, a Reynolds equation which contains turbulent and two-phase interfacial effects is needed, so the following basic assumptions are made for the modeling.
(1) The bubbles in the cavitation flow be regarded as spheres.
(2) The flow field is a small perturbation of the turbulent Couette flow [24]. (3) The turbulent shear stress is defined by the law of wall, Reichardt [25] empirical formula is used to define the eddy viscosity coefficient. (4) The average viscosity along the film thickness is adopted in generalized turbulent Reynolds equation.
The spiral area is an irregular area in the polar coordinate system, a spiral coordinates system (η, ζ ) are adopted to improve numerical calculation accuracy. The polar coordinates (θ , r) are converted into the spiral coordinates system (η, ζ ) as follows: The spiral patterns before and after the coordinates conversion are shown in Figure 2.
A generalized turbulent Reynolds equation with cavitation interface effects and inertial effects in the (η, ζ ) coordinates is derived based on two-phase flow theory [26]: Figure 1 Structure schematic of inward-pumping SGTB where The following expression can be given for the logarithmic spirals: The circumferential velocity distribution of the water film is written as where f eq (η, ς) is the equilibrium distribution function of the cavitation bubble diameter; f eq (η, ς)dς represents the number of bubbles with diameter between (ς , ς + dς) per unit volume of water at spatial coordinate η under the equilibrium state. The first term on the right-hand side of Eq. (2) is the hydrodynamic effect, the second term on the right side of Eq. (2) is the interface hydrodynamic effect, the third and fourth terms on the right side of Eq. (2) is the inertia hydrodynamic effect. The interface momentum transfer function includes the following three parts: where M 1 (η) is the interface momentum transfer term due to mass transfer, M 2 (η) is the interface momentum transfer term due to viscous drag, M 3 (η) is the interface momentum transfer term due to surface tension. Because the film thickness is thin, the liquid velocity in M 1 (η), M 2 (η) are approximated by the average velocity. M 1 (η), M 2 (η) and M 3 (η) are written as: where C Ds = 24 R e . Assume that the range of bubble diameter is divided into M smaller intervals, the diameter distribution of the bubbles is uniform in a smaller interval. So the equilibrium distribution function of bubble diameter in the kth smaller diameter interval (ς k , ς k+1 ) can be expressed as: where n k eq (η) is the number of bubbles in the kth diameter interval (ς k , ς k+1 ) . And the integral in Eq. (12) can be approximated as: The boundary conditions of Eq. (2) are as follows: The pressure at the groove-ridge boundary can be solved by using continuous condition of flow at the groove-ridge boundary. The finite volume at the grooveridge boundary is shown in Figure 3.
In a finite volume at the groove-ridge boundary, the continuous condition of the flow can be written as x k n k eq (η).

Turbulent Cavitating Flow Energy Equation
The temperature field of the cavitating flow water film can be obtained by solving the two-phase flow energy equation under high-speed conditions, and the following assumptions are adopted to simplify the energy equation.
(1) The circumferential convection term is only considered in the energy equation, because the circumferential velocity is much larger than the radial velocity. (2) Turbulent pulsation heat conduction term is ignored. (3) The viscous dissipation term along the film thickness is only considered in the energy equation.
According to the above assumption, the energy equation in the (η, ζ ) coordinates in the spiral coordinate system can be simplified to where The interface term E int in the energy equation is expressed as The following boundary conditions are adopted for the energy Eq. (21).
+2σ uπ M k=1 x k n k eq (η). The viscosity temperature relationship is written as The water film temperature field can be obtained solved by simultaneously solving the energy Eq. (21), the heat conduction equation of the stationary ring and the heat conduction equation of the thrust disk.The surface temperature T s (δ s ) of the stationary ring and surface temperature T d (η, 0) of the thrust plate are set as the boundary condition of the energy Eq. (21).

Heat Conduction Equation in Stationary Ring
The circumferential and radial heat conduction of the stationary ring are ignored, the heat conduction equation is simplified as follows: The boundary conditions of Eq. (29) are given as follows.
(1) At the surface between the stationary ring and the fluid (z s = δ s ): (2) At the surface between the stationary ring and the ambient (z s = 0): The expression of the temperature at the surface of the stationary ring can be obtained from Eq. (29) and boundary conditions (30), (31). The temperature at the surface of the stationary ring is written as

Heat Conduction Equation in Thrust Disk
For the rotary thrust disk, the heat transfer effects in circumferential and radial directions are also ignored, the heat conduction equation is simplified as follows: The boundary conditions of Eq. (33) are given as follows: (1) At the surface between the thrust disk and the fluid (z d = 0): (2) At the surface between the thrust disk and the ambient (z d = δ d ):

Force Balance Equation of Bubble
The force balance equation of bubbles is expressed as The bubble velocity u b in the cavitating flow is calculated by the force balance Eq. (36). Substitute F B , F G , F D into Eq. (36), the bubble velocity u b can be expressed as [27]

Population Balance Equation of Bubbles
Interface effect is an important phenomenon for cavitating flow, the momentum, mass and energy transfer occur through the interface between the gas-liquid. The interfacial area per unit volume liquid is positively correlated with the size distribution of the bubble, the bubble size distribution can be described by defining a probability density function f (r, ς , t) , and the internal coordinates ς are taken as bubble diameter in the model, and the f (r, ς, t)dς represents the number of bubbles with between (ς , ς + dς) at per volume liquid. Thus, integration of f over bubble diameter results in the total number of bubbles per volume liquid.The breakage, coalescence are two main factors affecting bubble size distribution in the cavitating flow, and bubble size distribution is predicted by the model of the breakage and coalescence processes of bubbles, this leads to the so-called population balance equation (PBE).
The PBE is a transport equation, the evolution of function f can be described by the PBE in the spiral coordinates (η, ζ ) , the PBE is written as follows: where, the first term on the right-hand side of Eq. (44) represents breakup sink term of bubbles with diameter ς per unit time, the second term on the right-hand side of Eq. (44) represents breakup source term of bubbles with diameter ς per unit time, the third term on the righthand side of Eq. (44) represents coalescence sink term of bubbles with diameter ς per unit time,the fourth term on the right-hand side of Eq. (44) represents coalescence x k n k eq (η). (43) source term of bubbles with diameter ς per unit time, the fifth term on the right-hand side of Eq. (44) represents bubbles with diameter ς source term.
The initial condition for Eq. (44) is given The breakage frequency b(ς ) is given as [28]: The daughter bubble size redistribution function h b (ξ , ς) is given as [28] The coalescence closure is given as [28]: The value of parameter β c is taken as 2.48. The expression of source term S c (ς ) is where The redistribution function of bubble diameter ς is given as The diameter distribution density function of the gas nucleus is given by experimental data [29]: where the value of parameter α is taken as − 10/3.

Calculation of Dynamic Characteristics of the SGTB
As shown in Figure 4, the origin of Cartesian coordinate is located at the center of the stationary ring. The possible independent motions of the thrust disk are axial movement along z axis and angular movements around x and y axes. The perturbed displacements and velocities of the thrust disk with respect to the steady equilibrium position are defined as (�z, �ϕ x , �ϕ y ) and (�ż, �φ x , �φ y ) . The assume initial position of the thrust disk is (h 0 , ϕ x0 , ϕ y0 ) , the film thickness in quasi-equilibrium can be expressed as The perturbations will have influence on the pressure distribution in the water film through Reynolds equation. The high order terms are ignored, the transient pressure of water film can be expressed as Substituting Eqs. (56) and (57) into Eq. (2), the zeroth and first-order perturbation generalized Reynolds equations in the (η, ζ ) coordinates system are obtained for the SGTB respectively, where the operator Rey() is defined as and q j j = z, ϕ x , ϕ y ,ż,φ x ,φ y can be found in Appendix.
The stiffness and damping coefficients in the (η, ζ ) coordinates system can be found by the following integrals where

Numerical Calculation Method of Theoretical Model
The generalized Reynolds Eq. (2), energy Eq. (21), and PBE Eq. (44) are the main three governing equations for the bearing. The pressure field, temperature field, and bubble diameter distribution are obtained by solving these three equations respectively, and pressure field; temperature field and bubble diameter distribution are coupled to each other, and an iterative algorithm needs to be employed to solve these distributions. Eq. (2) and Eq. (21) are discretized using finite difference, the pressure field is obtained by solving discretization equations of Eq.
(2) with the SOR iterative method, the temperature field is obtained by solving discretization equations of Eq. (21) with the stepping method.
(58) Rey() The shear stress parameter h * c of the Couette flow is included in turbulent velocity distribution Eq. (2), The governing equation for parameter h * c is written as [24] Eq. (62) is a nonlinear equation with respect to h * c , the parameter h * c can be solved using the Newton iteration method. The Newton iteration formula is expressed as h * c where

Numerical Solution of Bubble PBE
A discrete interval method [30] is employed to solve the bubble population balance Eq. (44). The procedure of the discrete interval method is as follows.
(1) Assume that the entire bubble diameter range is divided into M smaller intervals. The bubble diameter is approximated uniform distribution at the kth representative diameter range (ς k , ς k+1 ) , then the f (η, ς) can be expressed by the bubble number density n k (η) at the kth representative diamete range (ς k , ς k+1 ) . Thus, where (2) The integrating the continuous Eq. (44) over a smaller range (ς k , ς k+1 ) . Thus, Since coalescence or breakage of bubbles of these diameters results in the formation of new bubbles whose diameters do not match with any of the node diameters, the new bubble diameter will be assigned to the adjoining node by introducing weight factor γ , the weight factor γ can be determined by the total mass conservation of the bubble and the number of bubbles conservation. The expression is written as Substituting Eq. (65) into Eq. (67), and the integral term at the right end of Eq. (67) is reconstructed using the mean value theorem on breakage frequency and coalescence frequency, so Eq. (67) can be closed, a closed set of equations with regard to n k (η, t) can be written as where Coefficient α k,i is given as Eq. (69) is a set of first-order differential equations. Suppose that at a fixed time t, a set of nonlinear equations with regard to n t k,l can be obtained by discretizing the second term on the left side of Eq. (69). The discrete equation at node l is written as follows: where the operator Q (n k,l ) are defined as Eq. (74) can be solved using the Newton-SOR iterative method.In addition, since h b (x k , ς) is a given function, the coefficient α k,i can be calculated using the Gaussian integration method, and the value of the coefficient α k,i is stored in a matrix, it will be called directly in the (70) iterative calculation process, and the calculation time can be greatly saved by using the iterative process. The term of the time derivativein the transient Eq. (69) is specially processed with finite difference,and a weighted averaging method is adopted for the time-differencing. It can be written as [31] where, n t+1 k,l = n k,l (t + �t) , weight factor 0 ≤ ω t ≤ 1. Aset of nonlinear equations with regard to n t+1 k,l can be rewritten as The nonlinear Eq. (77) is also solved using the Newton-SOR iterative method. And downhill condition is adopted in numerical iterative procedure to improve convergence. The downhill condition is expressed as where m is iteration number. The overall algorithm flow diagram is given by Figure 5. Table 1 list the related parameters of the water-lubricated SGTB for the numerical calculation.In order to verify this theoretical model, the simulating results of bubble distribution and axial stiffness coefficient are compared with experimental ones from our previous test study in Refs. [22,23]. The experimental rotational speed are specified as 15000 r/min and 18000 r/min respectively, and the corresponding Reynolds numbers (1256 and 1507) are greater than the critical Reynolds number R e (R e = 1000) [32] of the bearing. And the cavitation number (σ = 0.14 to 0.2) of the water-lubricated bearing is smaller than the cavitation inception number (σ i ≈ 0.7) [33] of water under the above specified experimental rotational speed and standard atmospheric pressure, and it shows that cavitation occurs in experimental waterlubricated bearings. Figure 6 shows the photographs of bubble distributions taken during the experiments under the rotational speeds of 15000 r/min and 18000 r/min.

Bubble Distribution
F q n t+1 k,l = n t+1 k,l + �tω tQ n t+1 The experimental data of bubble diameter distribution can be obtained through the photo 6 using the experimental processing method provided in Ref. [22]. In order to improve the reliability of the experimental data of the small bubble diameter distribution, the multiple bubble photos are taken under a specified experimental rotation speed condition.The experimental data of bubble diameter can be obtained from multiple bubble photos using the image processing method of Ref. [22]. The average experimental value of the bubble diameter is calculated by statistically averaging the experimental data of the bubble diameter. A comparison of the theoretical results of the bubble size distribution with the experimental statistics results is shown in Figure 7. It can be seen that, for a specified rotational speed of 15000 r/min, the results predicted by the model are in good agreement with the experimental statistics results over the entire bubble size range; while, for a specified rotational speed of 18000 r/min, the results predicted by the model are in good agreement with the experimental statistics results in the large size bubble range, but the deviation between the Table 1 Parameters of the water-lubricated SGTB

Item Value
Inner radius r in (mm) 7.5 Outer radius r out (mm) 20 Basecircle radius r b (mm) is not easy to identify the edge of the small bubble outline in the experimental photos. This deviation can be attributed to the statistical error of the experimental data of the small bubble diameter due to the above reasons.

Axial Stiffness Coefficient
A comparison of the predicting axial stiffness coefficients with the experimental ones is shown in Figure 8. It can be seen that the predicting stiffness coefficients with axial load are generally in agreement with the experimental ones, but the experimental values are larger than the predicting ones in the entire range of applied load. This may be explained by the fact that the actual water flowrate of the bearing is less than the theoretical one, because the water is tended to be prevented from entering the clearance of inward-pumping SGTB due to the centrifugal effect. Figure 9 shows the influence of cavitation effect on the stiffness coefficients of the SGTB with different spiral angles. The results show the K * zz , K * ϕ x ϕ x , K * ϕ y ϕ y predicted by the cavitation flow model are larger than that of the noncavitation flow model when the spiral angle is less than 50°. This may be explained by the fact that the bubbles cause the interfacial hydrodynamic effect for the the cavitating flow. In addition, the relative deviation between the two models decreases with the increasing of spiral angle, so the interface hydrodynamic effect of the bubbles decreases with the spiral angle.It means that the influence of the cavitation effect on the water film stiffness is weakened with the increasing of spiral angle. The  K * ϕ x ϕ y , K * ϕ y ϕ x are of opposite sign, and K * ϕ x ϕ y and K * ϕ y ϕ x of cavitating flow model with the angle generally coincide with those of non-cavitating flow model. It shows that the influence of cavitation effect on K * ϕ x ϕ y , K * ϕ y ϕ x are negligi-ble. In addition, the K * ϕ x z , K * zϕ x , K * ϕ y z , K * zϕ y with different spiral angle sare approximately equal to zero for the cavitating flow or the non-cavitating flow. Figures 10 and 11 show the influence of cavitation effect on the water film stiffness coefficients for different groove depths and rotating speeds. The result also show the K * zz , K * ϕ x ϕ x , K * ϕ y ϕ y predicted by the cavitation flow model are larger than that of the non-cavitation flow model at a specific groove depth and rotating speed. However, the deviation of K * zz , K * ϕ x ϕ x , K * ϕ y ϕ y between the two models increase first and then decrease with the groove depth. Therefore, the interface hydrodynamic effect changes with the groove depth. A possible explanation is that the interface hydrodynamic effect is weakened and the cavitation effect is increased when the groove depth is increased in a certain range (20-35 μm), and the interface hydrodynamic effect is increased and the cavitation effect is weakened when the groove depth exceeds a given value. And the deviation of K * zz , K * ϕ x ϕ x , K * ϕ y ϕ y between the two models gradually increase as the rotating speed increases, This also can be explained by the fact that the interface hydrodynamic effect enhances and the cavitation effect is increased as the rotational speed increases.

Stiffness Coefficient
In addition, the variation of the K * ϕ x ϕ y , K * ϕ y ϕ x with the groove depth is identical with that with the spiral angle. Similarly, it also shows that the influence of cavitation effect on K * ϕ x ϕ y , K * ϕ y ϕ x are negligible for different groove depths and rotating speeds. The K * ϕ x z , K * zϕ x , K * ϕ y z , K * zϕ y predicted by the two models are approximately equal to zero for different groove depths and rotating speeds. Figures 12 and 13 show the influence of cavitation effect on the water film damping coefficients for different spiral angles and groove depths. It can be seen that the cavitation effect on damping coefficients is very weak. Because the damping coefficient is determined by the perturbed pressure, and the perturbed pressure with cavitating flow is only affected by the fraction C g of bubble volume, and the variation of bubble volume fraction is very small for the cavitating flow. In addition, the C * zz , C * ϕ x ϕ x , C * ϕ y ϕ y decreases with the increasing of the spiral angle and groove depth, and the C * zϕ x , C * ϕ x z , C * zϕ y , C * ϕ y z , C * ϕ x ϕ y , C * ϕ y ϕ x predicted by the two models are also approximately equal to zero for different spiral angles and groove depth. Figure 14 shows the influence of turbulence effect on the water film stiffness coefficients for different spiral angles. It can be seen that the K * zz ,K * lated by the turbulent cavitating flow model are larger than those of the laminar cavitating flow model at a specified spiral angle. This may be explained by the fact that additional Reynolds shear stress generates due to pulsating velocity in the turbulent flow, and equivalent viscosity of the turbulent water film increases, and the force of the turbulent water film needed to produce unit displacement increases, so the stiffness of the water film increases. In addition, it can be seen that the deviation of stiffness coefficient between turbulence model and laminar model gradually decreases with the increase of the spiral angle. This is due to the fact that the outward pumping hydrodynamic effect increases first and then decreases with the increase of the spiral angle. Therefore, when the spiral angle exceeds a certain value, the outward pumping hydrodynamic effect becomes weak and the additional Reynolds shear stress of the turbulent flow reduces. This indicates that the influence of turbulence effect on the water film stiffness coefficients is weakened when the spiral angle exceeds a certain value. And the K * zϕ x , K * ϕ x z , K * zϕ y , K * ϕ y z calculated by the two models are approximately equal to zero for different spiral angles. Figures 15 and 16 show the influence of turbulence effect on the water film stiffness coefficients for different groove depths and rotating speeds. It can also be seen that the K * zz , K * ϕ x ϕ x , K * ϕ y ϕ y , K * ϕ x ϕ y , K * ϕ y ϕ x calculated by the turbulent flow model are also larger than that of the laminar flow model for different groove depths and rotating speeds. And the deviation of the calculation results of the two models decreases with the increase of groove depth, and the deviation of the calculation results of the two models increase with the increase of rotating speed, and the K * zϕ x , K * ϕ x z , K * zϕ y , K * ϕ y z calculated by the two models are also approximately equal to zero for different groove depths and rotating speeds. Figures 17 and 18 show the influence of turbulence effect on the water film damping coefficients for different spiral angles and groove depths.It can also be seen that the C * zz , C * ϕ x ϕ x , C * ϕ y ϕ y calculated by the turbulent flow model are also larger than that of the laminar flow model for different spiral angles and groove depths. And the deviation between C * zz , C * ϕ x ϕ x , C * ϕ y ϕ y calculated by the two models is gradually decreased with the increase of spiral angle. This indicates that the turbulence effect is gradually weakened when the spiral angle increases. And the deviation between C * zz , C * ϕ x ϕ x , C * ϕ y ϕ y calculated by the two modes hardly changes with the increase of groove depth, and this indicates that the turbulence effect is not sensitive to changes in groove depth. In addition, the C * zϕ x , C * ϕ x z , C * zϕ y , C * ϕ y z , C * ϕ x ϕ y , C * ϕ y ϕ x calculated by the two models are also approximately equal to zero for different spiral angles and groove depth. Figure 19 shows the influence of turbulence effect on the water film damping coefficients for different rotating speed. It can be seen that the C * zz , C * ϕ x ϕ x , C * ϕ y ϕ y calculated by the two models increases gently with the increase of rotating speed, it means that C * zz , C * ϕ x ϕ x , C * ϕ y ϕ y is insensitive when rotating speed is increased. And the deviation between C * zz , C * ϕ x ϕ x , C * ϕ y ϕ y calculated by the two modes increase alse with the increase of rotating speed. The C * zϕ x , C * ϕ x z , C * zϕ y , C * ϕ y z , C * ϕ x ϕ y , C * ϕ y ϕ x calculated by the two models are also approximately equal to zero for different rotating speeds. Figure 20 shows the influence of turbulence effect on the bubble number for different radius. It can be seen that the bubbles mainly locate at the groove-ridge edges and the outer rim of bearing. This may be explained by the fact that the bubbles are formed due to the pressure drop near the groove-ridge edges. The result also show that the bubbles are approximately evenly distributed in dam region of bearing. In addition, the number of bubbles with laminar flow is more than the number of bubbles with turbulent flow at a specified circumferential position. Because the Reynolds shear stress due to the turbulent flow causes the equivalent viscosity of the turbulent liquid to increase, consequently, the turbulent flow effect on dynamic pressure is enhanced, which contributes to the generation of the bubbles.

Conclusions
In this paper, a theoretical model of water lubrication based on two-phase flow considering both turbulence and cavitation bubble effect is established. In the turbulent state, the evolution of the size distribution of the bubbles is described by the bubble transport equation considering the mechanism of breakage and coalescence. The interval discretization method is used to solve the bubble transport equation. The size and spatial equilibrium distribution of the bubbles under turbulent flow are predicted by solving the generalized Reynolds equation and the bubble transport equation simultaneously. The dynamic characteristics of water-lubricated SGTB under turbulent cavitating flow is analyzed. Based on discussion, the following conclusion can be drawn: (1) The bubbles are mainly concentrated near the spiral-groove edges of the SGTB. (2) The small-sized bubbles is much more than the large-sized bubbles for the water-lubricated SGTB under the cavitation flow condition. (3) The cavitation has a significant effect on the direct stiffness coefficients, and the cavitation hardly affect the cross-coupled stiffness coefficients and the damping coefficients of SGTB.
(4) Compared to the cavitation, the turbulence has a much greater effect on the dynamic characteristics of the high-speed water-lubricated SGTB.