Flow Mechanisms and Lubrication Performance of Water-Lubricated Thrust Bearings with Herringbone Grooves

: Due to their excellent stability and zero leakage capability, thrust bearings with herringbone spiral grooves are frequently used in transmission mechanisms. However, the lubrication mechanism of thrust bearings has not been clearly understood and explained, preventing the optimization of the bearing performance. Thus, this paper is devoted to solving this problem by building a three-dimensional ﬁnite element ﬂow model. In this model, the change in viscosity temperature is considered using Roelands equation, and the turbulence and cavitation are taken into consideration. Using the established model, the inﬂuence of parameters such as spiral angle, groove width ratio, and rotational speed on the cavitation area of the thrust bearing are analyzed. The pressure contour and speed distribution are obtained inside the clearance, as well as the volume fraction of the gas phase at the end face. Finally, according to the analysis results, the optimum structural parameter for the herringbone spiral groove structure is proposed, which enables higher bearing stability and provides a reference for engineering practice.


Introduction
Herringbone groove thrust bearings (HGTBs) have broad application prospects in hydraulic machinery, machine tool spindles, and other industrial fields [1][2][3], due to the advantages of good stability, high stiffness, and zero leakage. In the bearing, oil and air are commonly used as lubricative mediums. However, the gas-lubricated HGTB has low stiffness, and the oil-lubricated HGTB has a rapid temperature rise. On the opposite, the water-lubricated [4,5] HGTB has excellent dynamic stability and low frictional power consumption. Therefore, water-lubricated HGTBs are more suitable as spindle supports for ultra-precision machine tools under high speeds [6][7][8]. Cavitation, however, commonly occurs in water lubrication, reducing bearing performance and service life. It is valuable for establishing a three-dimensional water-lubricated flow field model considering the cavitation effect, which is used to explore the influence of structural changes on the bearing performance.
In recent years, more and more scholars [9,10] have achieved bearing resistance reduction, and the bearing stability was improved. In some research, it was shown that changes in the lubricating properties of bearings were closely related to the design of the surface texture [11][12][13] (distribution, shape, and size). Guenat [14,15] used the experiment to study the effect of the herringbone groove geometry on lubrication performance deeply. Wang [16] developed a thermo-hydrodynamic model of a combined bearing, and the effect of the groove structure on the bearing characteristics was investigated. If the bearing worked at high speed, the cavitation inception number of water is much higher. Thus, cavitation effects cannot be ignored. Some researchers treat cavitation lubricants as gasliquid mixed fluids. According to the different calculation methods of the gas volume fraction, the cavitation models in the mixed fluid model can be divided into three types: (1) Model based on the R-P equation [17,18]; (2) Equation model based on gas solubility and surface tension of bubbles [19,20]; (3) Model of the Transport Equation based on Gas Volume Fraction [21,22]. Moreover, some researchers used CFD to analyze the cavitation. In this theory, the journal bearings were investigated numerically in three dimensions by using FLUENT software [23], and the viscosity changes have a significant impact on the pressure distribution, cavity volume, and load-bearing capacity. Furthermore, based on the multi-phase flow cavitation model and the dynamic mesh technology, the cavitation characteristics of the liquid film in the dynamic pressure mechanical seal gap were investigated [24]. Meanwhile, Xu et al. [25] determined the optimal mass transfer coefficient in the Lee phase transition mass transfer equation, and the influence of groove structure parameters on the vaporization characteristics was studied.
In the above literature summary, the herringbone groove structure has been widely used and studied in thrust bearing design. However, the effects of temperature change and cavitation on lubrication states have rarely been studied. In this paper, CFD is used for modeling and analysis. The cavitation effects of the herringbone groove under changes in viscosity and temperature are studied, and its sealing performance and load-bearing performance are analyzed. The lubrication mechanism of the herringbone groove seal is explored, and the influence of the groove structure on the cavitation effect is summarized.

Conservation Equations
Fluid domains are analyzed using CFD methods, and the relative equations can be found as following.
The mass conservation equation is expressed by Equation (1): where, ρ represents the fluid density, ρ → v is the gravitational body force, and t is the time. The momentum equation is presented as Equation (2): where, → v is the fluid velocity vector, p represents the static pressure, and τ is the stress. The energy equation is given in Equation (3): where, E represents the total energy of the fluid, h j is the enthalpy of component j, k eff is the effective heat transfer coefficient, k eff = k + k t , k t denotes the turbulent heat transfer coefficient, J is the diffusion flux of component j, and S h represents the volumetric heat source term.

Cavitation Equations
The cavitation phenomenon easily appears in the high-speed rotation of the water film. When the temperature changes drastically, the water rises to the boiling point, the originally small bubbles rapidly expand and escape from the surface of the liquid. In this study, the phase change process is mainly controlled by temperature. The volume of fluid (VOF) model is chosen for the tracing of the gas-liquid interface.
The interphase interface is tracked by solving the continuity equation. This equation is solved by using the volume fraction of one (or more) phases. For the q th phase, this equation has the following form: In Equation (4), α v is the phase q volume fraction, ρ q is the phase q density, . m qp is the mass transfer from phase q to phase p, and . m qp is the mass transfer from phase p to phase q. The source term on the right-hand side of Equation (4) (Sα p ) is zero.
Among various cavitation models, the Lee model is relatively mature and matches the VOF model well. Therefore, this study adopted the Lee model as the phase transition model. By simplifying Equation (4), the transport equation of the bubble phase is obtained.
In Equation (5), R e and R c represent the evaporation and condensation terms in the phase transition process, respectively. They are respectively determined by using the gas and liquid phase medium temperatures, the saturation vaporization temperature is considered, and the relevant equations are as follows: where, T represents the local saturation vaporization temperature; coeff is the evaporative condensation coefficient derived from the interfacial concentration under vesicular flow.
where, d is the vapor bubble diameter, and η represents called the coefficient of adaptation and characterizes the adsorption of vapor molecules on the surface of the liquid, approximated by 1 under equilibrium conditions. L is the latent heat of vaporization, M is the molar mass, and R is the gas constant.

Geometric Model
The simplified structure of the herringbone thrust bearing is shown in Figures 1 and 2. Compared with the general single-sided spiral groove structure, the lubricating medium at the inner and outer diameters is pumped into the groove when the bearing is in operation. This structure has excellent lubrication and sealing performance. Its structural expressions in the r-θ coordinate system are given by Equations (8) and (9):

CFD Model
The calculation model is simplified on the premise of assuring the accuracy of the calculation results and given the following assumptions: the effects of body force and inertial force are ignored, the speed slip is neglected, the effect of surface roughness is ignored.
A CFD model is built by considering the above assumptions. Solidworks is used for geometric modeling, and ANSYS is used for the meshing. Finally, the case is solved in Fluent. The film thickness direction is enlarged 100 times for facilitating observation. The grid and boundary conditions for the clearance water film are presented in Figure 3. The three-dimensional geometry of the clearance was discretized with Hexahedral elements, the quality of the grid was assessed by checking the elements' quality (0.98), aspect ratio (1.1), and skewness (0.04). Then, it is indicating that the grid was of high quality. Due to the existence of fluid viscosity, there will be a boundary layer when the fluid moves in the near-wall region. As shown in Figure 3b, grid refinements are made on both sides of the wall. The inlet and outlet of the three-dimensional model are set as pressure boundaries, the walls in contact with the rotating ring are set as rotating walls, and the others are set as static walls.

CFD Model
The calculation model is simplified on the premise of assuring the accuracy of the calculation results and given the following assumptions: the effects of body force and inertial force are ignored, the speed slip is neglected, the effect of surface roughness is ignored.
A CFD model is built by considering the above assumptions. Solidworks is used for geometric modeling, and ANSYS is used for the meshing. Finally, the case is solved in Fluent. The film thickness direction is enlarged 100 times for facilitating observation. The grid and boundary conditions for the clearance water film are presented in Figure 3. The three-dimensional geometry of the clearance was discretized with Hexahedral elements, the quality of the grid was assessed by checking the elements' quality (0.98), aspect ratio (1.1), and skewness (0.04). Then, it is indicating that the grid was of high quality. Due to the existence of fluid viscosity, there will be a boundary layer when the fluid moves in the near-wall region. As shown in Figure 3b, grid refinements are made on both sides of the wall. The inlet and outlet of the three-dimensional model are set as pressure boundaries, the walls in contact with the rotating ring are set as rotating walls, and the others are set as static walls. In Equations (8) and (9), R 1 is the inner diameter, R 2 represents the outer diameter, β is the spiral angle, θ denotes the angle, R in is the inner spiral, R out is the outer spiral, B g is the groove width, B r is the ridge width, B is the groove width ratio, B = B g /B r , H 2 is the water film depth, and H 1 is the bearing clearance.

CFD Model
The calculation model is simplified on the premise of assuring the accuracy of the calculation results and given the following assumptions: the effects of body force and inertial force are ignored, the speed slip is neglected, the effect of surface roughness is ignored.
A CFD model is built by considering the above assumptions. Solidworks is used for geometric modeling, and ANSYS is used for the meshing. Finally, the case is solved in Fluent. The film thickness direction is enlarged 100 times for facilitating observation. The grid and boundary conditions for the clearance water film are presented in Figure 3. The three-dimensional geometry of the clearance was discretized with Hexahedral elements, the quality of the grid was assessed by checking the elements' quality (0.98), aspect ratio (1.1), and skewness (0.04). Then, it is indicating that the grid was of high quality. Due to the existence of fluid viscosity, there will be a boundary layer when the fluid moves in the near-wall region. As shown in Figure 3b, grid refinements are made on both sides of the wall. The inlet and outlet of the three-dimensional model are set as pressure boundaries, the walls in contact with the rotating ring are set as rotating walls, and the others are set as static walls.

Viscosity-Temperature Relationship
For the high-temperature sealing medium, the change in medium viscosity with temperature must be considered. Roelands equation [26] is as follows: where μ is the dynamic viscosity of water, T is the water temperature, and μ0 is the viscosity at the reference temperature T0.

Solver Settings
A mixture model was selected and analyzed. Then, the implicit method was adopted. In this study, liquid water was set as the primary phase, and water vapor was set as the secondary phase. The two phases were embedded in the mass transfer mechanism. "Dispersed" was selected as the interface type. The Zwart-Gerber-Belamri model was used in the simulation.
A pressure-based solver was chosen for the numerical simulation. The SIMPLE algorithm was used, which is an algorithm based on the pressure-velocity coupling principle. Then, the default state was used in other settings. A tolerance of 10 −3 was set for the residual terms.

Verification of the Model
The work in Ref. 27 is similar to that of this paper, and the research method has been proved by experimental tests. Therefore, the results in Ref. 27 can be used to confirm the model established in this paper, and the result is given in Figure 4. It is found that the cavitation position obtained from the model established in this paper is consistent with the result in the literature [27].
The leakage was obtained experimentally and validated as a quantitative indicator. As shown in Figure 5, the results of current work calculation are analogous with theoretical values and tests in the literature [27]. The relative error between this study and Ref. 27 is less than 5%, which proves the reliability of the proposed model.

Viscosity-Temperature Relationship
For the high-temperature sealing medium, the change in medium viscosity with temperature must be considered. Roelands equation [26] is as follows: where µ is the dynamic viscosity of water, T is the water temperature, and µ 0 is the viscosity at the reference temperature T 0 .

Solver Settings
A mixture model was selected and analyzed. Then, the implicit method was adopted. In this study, liquid water was set as the primary phase, and water vapor was set as the secondary phase. The two phases were embedded in the mass transfer mechanism. "Dispersed" was selected as the interface type. The Zwart-Gerber-Belamri model was used in the simulation.
A pressure-based solver was chosen for the numerical simulation. The SIMPLE algorithm was used, which is an algorithm based on the pressure-velocity coupling principle. Then, the default state was used in other settings. A tolerance of 10 −3 was set for the residual terms.

Verification of the Model
The work in Ref. [27] is similar to that of this paper, and the research method has been proved by experimental tests. Therefore, the results in Ref. [27] can be used to confirm the model established in this paper, and the result is given in Figure 4. It is found that the cavitation position obtained from the model established in this paper is consistent with the result in the literature [27].

Results and Discussion
The subsequent simulation solution takes a single cycle model caused by the flow's periodicity. Then, computer resources are saved, and the accuracy of solution is improved. The geometric and working parameters of the face seal are shown in Table 1.  The leakage was obtained experimentally and validated as a quantitative indicator. As shown in Figure 5, the results of current work calculation are analogous with theoretical values and tests in the literature [27]. The relative error between this study and Ref. [27] is less than 5%, which proves the reliability of the proposed model.

Results and Discussion
The subsequent simulation solution takes a single cycle model caused by the flow's periodicity. Then, computer resources are saved, and the accuracy of solution is improved. The geometric and working parameters of the face seal are shown in Table 1.

Results and Discussion
The subsequent simulation solution takes a single cycle model caused by the flow's periodicity. Then, computer resources are saved, and the accuracy of solution is improved. The geometric and working parameters of the face seal are shown in Table 1.

Liquid Film Flow Mechanism
The flow, which is inside a micro gap, is difficult to capture. Thus, its flow characteristics cannot be accurately understood. In this paper, a cross-sectional approach was chosen to demonstrate micro gap flow. The method is as follows: the actual film thickness is used in the simulations, but it is scaled up in the direction of the film thickness when the detailed flow is charactered.
When HGTB is working, the lubricating medium for groove structure is pumped into the top. Under the obstruction of the ridge, a part of the fluid slows down and increases the pressure, and the other part of the fluid climbs over the ridge. Thus, the most complex area of flow is at the top of the groove. Figures 6 and 7 show the pressure distribution and velocity field distribution on the radius section at the top of the groove, respectively. As shown in Figure 6, the pressure in this section is concentrated at the top of the groove. Because of the extremely small gap size, there is no significant pressure change in the thickness direction of the film. From the results in Figure 7, eddy currents exist in the gap. The fluid is divided into two parts: one part creates eddy currents in the groove, and the other part goes to lubricate the bearing clearance.

Liquid Film Flow Mechanism
The flow, which is inside a micro gap, is difficult to capture. Thus, its flow characteristics cannot be accurately understood. In this paper, a cross-sectional approach was chosen to demonstrate micro gap flow. The method is as follows: the actual film thickness is used in the simulations, but it is scaled up in the direction of the film thickness when the detailed flow is charactered.
When HGTB is working, the lubricating medium for groove structure is pumped into the top. Under the obstruction of the ridge, a part of the fluid slows down and increases the pressure, and the other part of the fluid climbs over the ridge. Thus, the most complex area of flow is at the top of the groove. Figures 6 and 7 show the pressure distribution and velocity field distribution on the radius section at the top of the groove, respectively. As shown in Figure 6, the pressure in this section is concentrated at the top of the groove. Because of the extremely small gap size, there is no significant pressure change in the thickness direction of the film. From the results in Figure 7, eddy currents exist in the gap. The fluid is divided into two parts: one part creates eddy currents in the groove, and the other part goes to lubricate the bearing clearance.  In Figure 8, when the groove surface moves, the fluid in the groove is pushed from one side to the other. Fluid flow characteristics in the micro gap are abundant. The fluid is squeezed at the front face, and a pressure peak is generated. The fluid after the collision is divided into two parts: one forms a vortex in the groove, and the other flows to the bearing gap. The fluid in the clearance acts as a lubricating medium, and the bearing surfaces are pushed apart. Thus, the load-carrying capacity is increasing, and the friction coefficient is reduced [28].

Liquid Film Flow Mechanism
The flow, which is inside a micro gap, is difficult to capture. Thus, its flow characteristics cannot be accurately understood. In this paper, a cross-sectional approach was chosen to demonstrate micro gap flow. The method is as follows: the actual film thickness is used in the simulations, but it is scaled up in the direction of the film thickness when the detailed flow is charactered.
When HGTB is working, the lubricating medium for groove structure is pumped into the top. Under the obstruction of the ridge, a part of the fluid slows down and increases the pressure, and the other part of the fluid climbs over the ridge. Thus, the most complex area of flow is at the top of the groove. Figures 6 and 7 show the pressure distribution and velocity field distribution on the radius section at the top of the groove, respectively. As shown in Figure 6, the pressure in this section is concentrated at the top of the groove. Because of the extremely small gap size, there is no significant pressure change in the thickness direction of the film. From the results in Figure 7, eddy currents exist in the gap. The fluid is divided into two parts: one part creates eddy currents in the groove, and the other part goes to lubricate the bearing clearance.  In Figure 8, when the groove surface moves, the fluid in the groove is pushed from one side to the other. Fluid flow characteristics in the micro gap are abundant. The fluid is squeezed at the front face, and a pressure peak is generated. The fluid after the collision is divided into two parts: one forms a vortex in the groove, and the other flows to the bearing gap. The fluid in the clearance acts as a lubricating medium, and the bearing surfaces are pushed apart. Thus, the load-carrying capacity is increasing, and the friction coefficient is reduced [28]. In Figure 8, when the groove surface moves, the fluid in the groove is pushed from one side to the other. Fluid flow characteristics in the micro gap are abundant. The fluid is squeezed at the front face, and a pressure peak is generated. The fluid after the collision is divided into two parts: one forms a vortex in the groove, and the other flows to the bearing gap. The fluid in the clearance acts as a lubricating medium, and the bearing surfaces are pushed apart. Thus, the load-carrying capacity is increasing, and the friction coefficient is reduced [28].

Influence of Rotational Speed on the Flow Field
When the speed is changing, and other parameters remain constant. The fluid film pressure distribution can be seen in Figure 9. A high-pressure area appears at the tip of the herringbone groove. Low-pressure areas appear at the groove root of the outer diameter side of the herringbone groove, the groove root, and side wall of the inner diameter side groove. As shown in Figure 10, cavitation takes place at the root of the outer diameter side groove and the side wall of the inner diameter.
As the rotational speed increases, the maximum pressure of the fluid film magnifies, and the pressure distribution law of the fluid film is unchanged. The cavitation area gradually expands with the increasement of the rotational speed, which indicates that the cavitation area is extensible.

Influence of Rotational Speed on the Flow Field
When the speed is changing, and other parameters remain constant. The fluid film pressure distribution can be seen in Figure 9. A high-pressure area appears at the tip of the herringbone groove. Low-pressure areas appear at the groove root of the outer diameter side of the herringbone groove, the groove root, and side wall of the inner diameter side groove. As shown in Figure 10, cavitation takes place at the root of the outer diameter side groove and the side wall of the inner diameter.

Influence of Rotational Speed on the Flow Field
When the speed is changing, and other parameters remain constant. The fluid film pressure distribution can be seen in Figure 9. A high-pressure area appears at the tip of the herringbone groove. Low-pressure areas appear at the groove root of the outer diameter side of the herringbone groove, the groove root, and side wall of the inner diameter side groove. As shown in Figure 10, cavitation takes place at the root of the outer diameter side groove and the side wall of the inner diameter.
As the rotational speed increases, the maximum pressure of the fluid film magnifies, and the pressure distribution law of the fluid film is unchanged. The cavitation area gradually expands with the increasement of the rotational speed, which indicates that the cavitation area is extensible.

Spiral Angle
The spiral angle β directly affects the boosting effect of the thrust bearing. In the Figures 11 and 12, the effect of the spiral angle on the bearing's load-carrying capacity and leakage rate was investigated. The load capacity appears to increase with spiral angle from 50° to 70°. Then, it decreases as the spiral angle increases. When the spiral angle is about 68°, the bearing capacity reaches the maximum value. When the spiral angle is 0° or 90°, the dynamic pressure effect and step effect are greatly reduced, and the bearing capacity is almost lost. Thus, the leakage exhibits the downtrend.  As the rotational speed increases, the maximum pressure of the fluid film magnifies, and the pressure distribution law of the fluid film is unchanged. The cavitation area gradually expands with the increasement of the rotational speed, which indicates that the cavitation area is extensible.

Spiral Angle
The spiral angle β directly affects the boosting effect of the thrust bearing. In the Figures 11 and 12, the effect of the spiral angle on the bearing's load-carrying capacity and leakage rate was investigated. The load capacity appears to increase with spiral angle from 50 • to 70 • . Then, it decreases as the spiral angle increases. When the spiral angle is about 68 • , the bearing capacity reaches the maximum value. When the spiral angle is 0 • or 90 • , the dynamic pressure effect and step effect are greatly reduced, and the bearing capacity is almost lost. Thus, the leakage exhibits the downtrend.

Spiral Angle
The spiral angle β directly affects the boosting effect of the thrust bearing. In the Figures 11 and 12, the effect of the spiral angle on the bearing's load-carrying capacity and leakage rate was investigated. The load capacity appears to increase with spiral angle from 50° to 70°. Then, it decreases as the spiral angle increases. When the spiral angle is about 68°, the bearing capacity reaches the maximum value. When the spiral angle is 0° or 90°, the dynamic pressure effect and step effect are greatly reduced, and the bearing capacity is almost lost. Thus, the leakage exhibits the downtrend. Figure 11. Effect of the spiral angle on the load-carrying capacity. Figure 11. Effect of the spiral angle on the load-carrying capacity.

Groove Width Ratio
From Figure 13, as the groove width ratio increases, the axial load capacity of the thrust bearing gradually increases. When the groove width ratio is about 1, the bearing capacity reaches the maximum value, and then gradually decreases. If the groove width ratio is equal to 0 or infinite, the bearing end face becomes flat, and the dynamic pressure effect and step effect cannot be found. Therefore, there is a suitable groove width ratio to maximize the bearing capacity. As seen in Figure 14, when the groove width ratio increases, the leakage also increases and tends toward a stable value. The magnified groove width ratio corresponds indirectly increases the clearance of the bearing, and the leakage is also increasing.

Groove Width Ratio
From Figure 13, as the groove width ratio increases, the axial load capacity of the thrust bearing gradually increases. When the groove width ratio is about 1, the bearing capacity reaches the maximum value, and then gradually decreases. If the groove width ratio is equal to 0 or infinite, the bearing end face becomes flat, and the dynamic pressure effect and step effect cannot be found. Therefore, there is a suitable groove width ratio to maximize the bearing capacity. As seen in Figure 14, when the groove width ratio increases, the leakage also increases and tends toward a stable value. The magnified groove width ratio corresponds indirectly increases the clearance of the bearing, and the leakage is also increasing.

Groove Width Ratio
From Figure 13, as the groove width ratio increases, the axial load capacity of the thrust bearing gradually increases. When the groove width ratio is about 1, the bearing capacity reaches the maximum value, and then gradually decreases. If the groove width ratio is equal to 0 or infinite, the bearing end face becomes flat, and the dynamic pressure effect and step effect cannot be found. Therefore, there is a suitable groove width ratio to maximize the bearing capacity. As seen in Figure 14, when the groove width ratio increases, the leakage also increases and tends toward a stable value. The magnified groove width ratio corresponds indirectly increases the clearance of the bearing, and the leakage is also increasing. Figure 13. Effect of the groove width ratio on the load-carrying capacity. Figure 13. Effect of the groove width ratio on the load-carrying capacity.

Groove Width Ratio
From Figure 13, as the groove width ratio increases, the axial load capacity of the thrust bearing gradually increases. When the groove width ratio is about 1, the bearing capacity reaches the maximum value, and then gradually decreases. If the groove width ratio is equal to 0 or infinite, the bearing end face becomes flat, and the dynamic pressure effect and step effect cannot be found. Therefore, there is a suitable groove width ratio to maximize the bearing capacity. As seen in Figure 14, when the groove width ratio increases, the leakage also increases and tends toward a stable value. The magnified groove width ratio corresponds indirectly increases the clearance of the bearing, and the leakage is also increasing.

Conclusions
In this study, a finite element flow field model is established to understand and explain the flow mechanisms of HGTBs. The influence of structural parameters on the bearing performance was investigated. In this model, viscous temperature changes, turbulence, and cavitation are considered. The main conclusions are summarized as follows: (1) Water is pumped into the spiral groove from both ends by vicious and differential pressure forces. The lubricating medium is converged at the top of the herringbone groove. Then, one part of the fluid goes over the ridge, and the other part forms vortices in the groove. The lubrication performance of the thrust bearing is affected by the dynamic pressure and the step.
(2) The cavitation is concentrated at the root of the outer diameter side of the herringbone groove and the inner diameter wall. The cavitated area is extended by the increasing rotational speed. There is an optimum value for the groove width ratio B and the helix angle β, which can improve the bearing performance.
(3) Optimum values for the geometrical parameters of the thrust bearing have been determined using the load-carrying capacity and leakage. When the bearing clearance is 10 µm, the groove width ratio is taken as 1, and the spiral angle β can be selected from 65 • to 70 • .