# Feasibility Analysis of Adopting the Hydrogen Hydrostatic Thrust Bearing
## 2 Numerical Model
### Computational Model  
The complete Navier-Stokes equation of three-dimensional compressible flow is numerically simulated as follows [29].  
\[\rho\bigg{(}\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{ \partial u}{\partial y}+w\frac{\partial u}{\partial z}\bigg{)}=\rho f_{x}- \frac{\partial P}{\partial x}+\mu(\frac{\partial^{2}u}{\partial x^{2}}+\frac{ \partial^{2}u}{\partial y^{2}}+\frac{\partial^{2}u}{\partial z^{2}}) \tag{1a}\] \[\rho\bigg{(}\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{ \partial v}{\partial y}+w\frac{\partial v}{\partial z}\bigg{)}=\rho f_{y}- \frac{\partial P}{\partial y}+\mu(\frac{\partial^{2}v}{\partial x^{2}}+\frac{ \partial^{2}v}{\partial y^{2}}+\frac{\partial^{2}v}{\partial z^{2}})\] (1b) \[\rho\bigg{(}\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+v \frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z}\bigg{)}=\rho f_{z} -\frac{\partial P}{\partial z}+\mu(\frac{\partial^{2}w}{\partial x^{2}}+\frac {\partial^{2}w}{\partial y^{2}}+\frac{\partial^{2}w}{\partial z^{2}}) \tag{1c}\]  
where, \(x\), \(y\) and \(z\) are directional coordinates, \(u\), \(v\) and \(w\) are velocities in the three directions, respectively, \(t\) is time, \(\mu\) is dynamic viscosity, \(\rho\) is gas density, \(P\) is gas pressure and \(f\) is external force per unit volume of fluid.  
Consider the following equation for the continuity of the flow.  
\[\frac{\partial(\rho u)}{\partial x}+\frac{\partial(\rho v)}{\partial y}+\frac {\partial(\rho w)}{\partial z}+\frac{\partial\rho}{\partial t}=0 \tag{2}\]  
To determine the flow pattern, the Reynolds number (\(Re\)) is used. The flow regime of the gas film is regarded as turbulent if the Reynolds number is greater than 2000, and the clearance thickness is the characteristic length [47]. \(Re\) is computed as [48]:  
\[Re=\frac{\omega rc_{0}}{\
u} \tag{3}\]  
where \(\
u\) is kinematic viscosity, \(\omega\) is rotational speed, \(c_{0}\) is film clearance, \(r\) is the radius and \(r=D_{1}/2\). The highest Reynolds number in the gas film reached 3356 at 50,000 rpm and 0.5 MPa; the flow pattern was considered to be turbulent flow.  
To explore the characteristics of the HHTB, the CFD turbulence model was computed using FLUENT(r). Considering the vortex effect, a \(k\)-\(\varepsilon\) model was used, which could better reflect complex flows involving impingement. The flow field of an HHTB is symmetrical, and calculated according to the structure of a gas supply orifice with periodic boundaries on both sides, as displayed in Figure 3. The specific boundary conditions are defined as follows:  
Figure 2: Hydrogen hydrostatic thrust bearing structure and simulation model: (**a**) front view and (**b**) sectional view.  
When the thrust plate is exposed to a small harmonic disturbance, the gas film clearance changes as well, which is expressed as follows:  
\[\Delta c=a\cdot sin\ (2\pi ft+\varphi) \tag{4}\]  
where \(a\) represents amplitude, \(f\) represents disturbance frequency and \(\varphi\) represents the initial phase. Simultaneously, the increment of the gas film bearing capacity also has harmonic motion [49]. The dynamic grid model was adopted to realize the periodic motion of the dynamic gas film.  
To comprehensively evaluate the dynamic characteristics of the HHTB, the damping dissipation energy, equivalent damping and dynamic stiffness coefficient were analyzed. The relationship between dynamic load capacity and dynamic stiffness can be formulated as follows:  
\[[K]=\sum_{i=1}^{n}[K]^{i},\{f\}=\sum_{i}^{n}\{f\}^{e} \tag{5}\]  
\[[K]\{\Delta c\}=\{f\} \tag{6}\]  
Equivalent damping was introduced to represent the energy dissipation capacity, as shown in the following [50]:  
\[D_{eq}=\frac{2\Delta E_{dis}}{\Omega l^{2}} \tag{7}\]  
where \(\Delta E_{dis}\) is the energy dissipation (the closed area of the load capacity-film clearance curve), \(l\) is the maximum deformation of the gas film, \(\Omega\) is the disturbance frequency and \(\Omega=2\pi f\).