Multiscale simulation of flow in gas-lubricated journal bearings: A comparative study between the Reynolds equation and lattice Boltzmann methods

Abstract

On account of varying tiny cross-sectional clearance and fairly high bearing speed, the feature of lubricant gas of gas-lubricated journal bearings (GLJBs) is a complex multiscale problem. In this research, both macroscopic Reynolds equation and mesoscopic lattice Boltzmann methods (LBM) are adopted to successfully establish models for restoring flow characteristics of GLJBs from continuum, slip to transition flow regimes. The influences of eccentricity ratios, bearing speeds and Knudsen numbers (Kn) on pressure, thickness, velocity distributions and flow fields of lubricant flow are investigated. The distinctions between the results through two methods, which can be complementary, are compared and analyzed with maximum error ≯9.6%. It indicates that keeping the eccentricity ratio or bearing speed at appropriate high levels is good for improving its working performance. When Kn is enhanced, the maximum and minimum pressure peaks will decrease and increase, respectively, whereas, the impact degree of rarefaction effect on the pressure distribution will reduce after Kn exceeding a certain level. Besides, an interesting backflow phenomenon in the lubricant flow field is demonstrated. This investigation enriches numerical models for GLJBs at different flow regimes, and provides new fruitful insights into the flow characteristics at multiscale.

KEYWORDS:

• Gas-lubricated journal bearing (GLJB)
• multiscale flow
• Reynolds equation
• lattice Boltzmann method (LBM)
• boundary-condtion

1. Introduction

As an important type of rotating machinery support element, the gas-lubricated journal bearing (GLJB) possesses the advantages of strong adaptability, outstanding stability, perfect reliability, high temperature resistance, good starting/stopping properties and long working life (Samanta et al., 2019). It plays a crucial role in the stability and operation applications for rotating machinery such as turbo-expanders (Lai et al., 2018; Li et al., 2017), turbochargers (Andrés & Rodríguez, 2020; Guo et al., 2017), small micro-turbines (Sim et al., 2013; Zhang et al., 2018a) and engines for turboshaft propulsion (Lehn et al., 2017). Thus, the research on GLJBs has attracted much attention (Arghir & Benchekroun, 2019; Gu et al., 2020a, 2020b; Rom & Müller, 2019). Generally speaking, static characteristics, dynamic characteristics and stability are the dominating content of air bearing investigation (Chen et al., 2019; Baum et al., 2021; Guo et al., 2019). However, several difficult problems are still poorly understood, involving the multiscale features of lubrication flow in the variable cross-section clearance, the detailed mechanisms behind the macroscopic phenomena, the characteristics of bearing gas flow under high Knudsen number (Kn), etc. Therefore, it is of great significance to explore the flow characteristics of GLJBs at the multiscale level.

At the macroscale level, the theoretical foundation for studying the characteristics of lubrication films is the Reynolds equation, which possesses the advantages of convenient application, quick calculation and clear physical meaning (Bryant, 2013). The static characteristics of foil bearings were obtained using the Reynolds equations by Heshmat et al. (1983, 2000). Rao (2009) discovered that, if interface slip occurs, the bearing capacity will increase. It was pointed out that a compromise must be made in terms of convenience of use and simple programming when solving bearing problems by Larsen et al. (2016). Nielsen and Santos (2017) investigated the modelling problem taking into consideration the large droop effect in nonlinear transients. Previous studies of GLJBs mainly focus on the conventional Reynolds equation and its modifications. Nevertheless, it is worth noting that the clearance of the journal and bearing house ranges from microns to hundreds of microns, in which the flowing air film is a typical multiscale process. Exactly as in the findings from the real case studies of applications of computational fluid dynamics (CFD) models by Salih et al. (2019) and Ghalandari et al. (2019), the effects of the microchannel should not be neglected. However, owing to the continuum assumption, the Reynolds equation is no longer appropriate for high Kn operating conditions with rarefied gas effects. In addition, tiny variations in flow fields, such as the backflow phenomenon of lubricating air flow and the velocity streamline distribution, are scarcely captured or exhibited by the macroscopic method, either.

On the contrary, at the mesoscale level, the lattice Boltzmann method (LBM) has plentiful potential advantages, such as no restriction of the continuity hypothesis, complex boundaries adaptability, simple programming and easy parallelism, and specializes in describing multiphase flows, fluid mechanics, interfacial dynamics and interstitial fluids (Chen et al., 2016; Guiet et al., 2011; Jiao et al., 2021; Liu et al., 2016; Sergi et al., 2015). Many scholars have adopted the LBM to study the micro-flow in narrow and small spaces. For instance, Nie et al. (2002) certified that the LBM is capable of exhibiting the velocity slip and the nonlinear pressure drop along a microchannel. An LB model for rarefied gas flow at high Kn in a microchannel was derived by Gokaltun and Dulikravich (2007). Different LBM models to simulate the rotating slipstream in a coaxial cylinder were proposed by Watari (2011). Yang et al. (2019) successfully applied the LBM under an improved slip boundary. Inspired by the features and benefits of LBM, the present authors intend to introduce the LBM to further investigate the details of flow characteristics in the narrow clearance of GLJBs, and the flow in high Kn conditions, as a complement to the Reynolds equation. However, to the authors’ best knowledge, few LBM studies have been reported in the GLJB domain, especially for the lubricant air flowing in the GLJB clearance with the characteristics of strong shear action, annular variable cross-section and high velocity gradient. An available simulation approach for flow in GLJBs at the multiscale level still needs to be proposed.

In this study, both the macroscopic Reynolds equation method and the mesoscopic LBM are introduced to successfully establish numerical models of GLJBs at the multiscale level that are applicable in different flow regimes. We numerically predict the influences of eccentricity ratios, bearing speeds and Kn on the lubricating air film pressure distributions, thickness distributions, velocity distributions and flow field variations at the multiscale level, and a comparison study between the two methods is conducted. Additionally, we find and discuss the interesting backflow phenomenon in the lubricating air flow field and discuss its formation processes, which is novel to the best of our knowledge. Moreover, on account of clear physical meaning, the macroscopic model can be regarded as proof for the mesoscopic model; in contrast, owing to its ability to capture details at a lower scale and its applicability at higher Kn, the mesoscopic method can be seen as a complement to the macroscopic method. As a result, this investigation can be regarded as an improvement on current numerical models of GLJBs, providing new insights into the flow characteristics and making a contribution to engineering research on, and development of, GLJBs.

2. Models and methods

The structure of a GLJB is shown in Figure 1. The lubricant air in the journal and bearing house clearance is marked.

Figure 1. Schematic diagram of GLJB structure: (a) radial; (b) axial.

When the journal starts to rotate in a clockwise direction, the lubricating gas continuously accumulates in the convergent wedge and is compressed, leading to an increase of pressure in the convergent wedge. After the journal speed reaches a certain value, a lubricating air film is formed. The position of the journal is shifted to the left, as demonstrated in Figure 2, where the resultant pressure force generated by the lubricating air film is balanced by the load W.

Figure 2. Schematic diagram of axis balance position.

The equilibrium position of the axis O_j is described by the offset angle φ, which is defined as the angle between the line O_jO_b, and the eccentricity ratio ϵ, where $\epsilon =e\text{/}C$, $C=R_b-R_j$. The relevant geometric, operating and environmental parameters studied in the paper are summarized in Table 1.

Table 1. Geometric, operational and ambient bearing parameters.

Parameter name	Value
Bearing length, L (m)	25.0 × 10^−3
Shaft radius, R (m)	10.95 × 10^−3
Radial clearance, C0 (m)	0.05 × 10^−3
Ambient pressure, pa (Pa)	1.013 × 10^5
Dynamic viscosity of gas, μ (Pa·s)	2.16 × 10^−5
Eccentricity ratio, ϵ	0.3∼0.9
Bearing speed, ω (rpm)	5 × 10^4∼1.5 × 10^5

2.1. Macroscopic numerical models and boundary conditions

2.1.1. Macroscopic numerical models

The lubrication state of the GLJB is the aerodynamic lubrication. The basic equation of the theoretical basis for it is the Navier–Stokes (N-S) equation, and the most obvious feature of the lubrication problem is that the fluid film thickness is much smaller than the geometric sizes of the other two dimensions, which makes it possible to simplify the N-S equation as follows: (1) $\{\begin{array}{c}{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }x}}=\mu {\displaystyle \frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{y}^2}}\\ {\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }y}}=0\\ {\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }z}}=\mu {\displaystyle \frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{y}^2}}\end{array}$(1)

As we all know, the equation of continuity is (2) $\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\frac{\mathrm{\partial }(\rho u)}{\mathrm{\partial }x}+\frac{\mathrm{\partial }(\rho v)}{\mathrm{\partial }y}+\frac{\mathrm{\partial }(\rho w)}{\mathrm{\partial }z}=0$(2)

Combining Equation (1) with Equation (2), the general form of the Reynolds equation, which is the basic equation of gas lubrication, can be derived: (3) $\begin{array}{rl}& \frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(p{h}^3\frac{\mathrm{\partial }p}{\mathrm{\partial }x}\right)+\frac{\mathrm{\partial }}{\mathrm{\partial }z}\left(p{h}^3\frac{\mathrm{\partial }p}{\mathrm{\partial }z}\right)\\ & =6\mu U\frac{\mathrm{\partial }(ph)}{\mathrm{\partial }x}+12\mu \frac{\mathrm{\partial }(ph)}{\mathrm{\partial }t}\end{array}$(3) where h is the film thickness. Transferring the Reynolds equation from Cartesian coordinates into cylindrical coordinates with the relevant dimensionless parameters introduced, the Reynolds equation can be expressed as follows: (4) $\begin{array}{rl}& \frac{\mathrm{\partial }}{\mathrm{\partial }\theta }\left(P{H}^3\frac{\mathrm{\partial }P}{\mathrm{\partial }\theta }\right)+{\left(\frac{2R}{L}\right)}^2\frac{\mathrm{\partial }}{\mathrm{\partial }Z}\left(P{H}^3\frac{\mathrm{\partial }P}{\mathrm{\partial }Z}\right)\\ & =\mathrm{\Lambda }\frac{\mathrm{\partial }(PH)}{\mathrm{\partial }\theta }+2\mathrm{\Lambda }\gamma \frac{\mathrm{\partial }}{\mathrm{\partial }{\tau }_s}(PH)\end{array}$(4) where $\theta$ (rad) is the angular coordinate: $\begin{array}{rl}\theta & =\frac{x}{R},Z=\frac{z}{0.5L},P=\frac{p}{p_a},H=\frac{h}{C},\mathrm{\Lambda }=\frac{6\mu \omega }{p_a}{\left(\frac{R}{C}\right)}^2,\\ \gamma & =\frac{{\omega }_s}{\omega },\overline{t}={\omega }_{s}t\end{array}$

Z is the dimensionless bearing length, P is the dimensionless pressure, H is the dimensionless film thickness, $\mathrm{\Lambda }$ is the bearing number, $\mu$ (Pa/s) is the dynamic viscosity, $\omega$ (rpm) is the angular velocity of the bearing, $\gamma$ is the whirl frequency ratio, $\overline{t}$ is the dimensionless time and ${\omega }_s$ (rad/s) is angular frequency of oscillation.

As can be seen from Figure 2, the angle $\mathrm{\angle }{\text{O}}_j\text{P}{\text{O}}_b$ is quite small in $\mathrm{\Delta }{\text{O}}_j{\text{O}}_b\text{P}$ so that $\text{cos(}\mathrm{\angle }{\text{O}}_j\text{P}{\text{O}}_b\text{)}\approx \text{1}$. Thereby, the air film thickness equation is obtained: (5) $\begin{array}{r}h=(R_b-R_j)+e\cos \theta =C+e\cos \theta =C(1+\epsilon \cos \theta )\end{array}$(5) As the gas-lubricated bearing works stably, the clearance is quite narrow (tens of microns or less) and appears to have different characteristics from under no-slip conditions. This effect is related to Kn, a measure of the deviation of the gas from a continuum medium, which is defined as follows: (6) $\text{Kn}=\frac{\lambda }{h_0}$(6) where h₀ is the fluid’s characteristic size and $\lambda$ is the molecular free path, which is determined by the gas properties and state (Wang et al., 2016). (7) $\lambda =\frac{1}{\sqrt{2}\pi d^{2}n}=\frac{kT}{\sqrt{2}\pi d^{2}p}$(7) where d is the diameter of molecular, n is the number density of molecules, p is the gas pressure, T is the gas temperature and k is Boltzmann’s constant.

2.1.2. Macroscopic boundary conditions

The finite difference method (FDM) is used in solving the Reynolds equation. The fluid domain is meshed as shown in Figure 3. The z- and x-directions are divided into equal parts, as nz and nx, respectively. Therefore, there are nz + 1 nodes along the length of the bearing in the z-direction, and nx + 1 nodes listed along the length of the circumference in the x-direction. The solution domain is defined as (0: nx, 0: nz).

Figure 3. Schematic diagram of mesh of macroscopic method.

The boundary conditions are defined as follows: (8) $P\left(x,0\right)=1;P\left(x,nz\right)=1;$(8) In the circumferential direction, positions of x = 0 and x = nx coincide, so P(x, 0) = P(x, nz).

2.1.3. Macroscopic grids independence verification

Because the simulation results are affected by mesh generation modes, in order to verify the independence from macroscopic grid partitions, we gradually intensify the gridding of the model from 3 × 10³ to 1.9 × 10⁴ and the variation of maximum dimensionless pressure with the number of grids is shown in Figure 4, of which the relevant calculation parameters are shown in Table 2.

Figure 4. Macroscopic grids independence verification.

Table 2. Calculation parameters for macroscopic grids independence verification.

Bearing length,	Shaft radius,	Radial clearance,	Ambient pressure,	Dynamic viscosity of	Bearing load,	Bearing speed,
L (m)	R (m)	C (m)	pa (Pa)	gas, μ (Pa·s)	W (N)	ω (rpm)
25 × 10^−3	10.95 × 10^−3	0.05 × 10^−3	1.013 × 10^5	1.87 × 10^−5	25	6 × 10^4

It can be seen that the pressure keeps basically constant after the number of grids exceeds 7 × 10³ and the mean deviation of the numerical result is around 0.9%. Comprehensively considering the computing time and precision factor, we choose this as the grid independent solution.

2.2. Mesoscopic numerical models and boundary conditions

2.2.1. Mesoscopic numerical models

• Multi-relaxation-time LBM (MRT-LBM)

In the LBM, the two most used collision operators are the Bhatnagar–Gross–Ktook (BGK) (Qian & Orszag, 1993) and the MRT (Humières, 1994). Owing to the fact that the BGK collision operator has some disadvantages, such as that the viscosity depends on the physical results and the restriction of the density ratio (Zhang et al., 2018b), we choose the MRT in this simulation, of which the evolution equation is as follows (Lallemand & Luo, 2000): (9) $\begin{array}{rl}& f_i\left(\mathit{x}+{\mathit{e}}_i{\delta }_t,t+{\delta }_t\right)-f_i\left(\mathit{x},t\right)\\ & =-\sum _j{\left({\mathbf{\text{M}}}^{-1}\mathbf{\Lambda }\mathbf{\text{M}}\right)}_{ij}\left[f_j\left(\mathit{x},t\right)-f_j^{\left(eq\right)}\left(\mathit{x},t\right)\right]\end{array}$(9) where $f_i(\mathit{x},t)$ is the density distribution function for node x at time, ${\mathit{e}}_i$ is the discrete velocity, ${\delta }_t$ is the time step $(i,j=0,1,\dots ,b)$, b is the total number of discrete velocities, $\mathbf{\text{M}}$ is the orthogonal transformation matrix and $f_j^{(\text{eq})}$ is the equilibrium distribution function that relates the lubricant density $\rho$, velocity $\mathit{u}$ and temperature $T$, which can be explained as follows: (10) $f_i^{(\text{eq})}(\mathit{x},t)={\omega }_i\rho \left[1+\frac{({\mathit{e}}_i\cdot \mathit{u})}{{\mathit{c}}_s^2}+\frac{{({\mathit{e}}_i\cdot \mathit{u})}^2}{2{\mathit{c}}_s^4}-\frac{{\mathit{u}}^2}{{\mathit{c}}_s^2}\right]$(10) where e is the unit length lattice velocity, $\mathit{e}={\delta }_x/{\delta }_t$, ${\delta }_x$ is the lattice space step size, ${\delta }_t$ is the time step, ${\mathit{c}}_{\mathit{s}}$ is the lattice speed of sound, ${\mathit{c}}_{\mathit{s}}=\sqrt{RT}$ (generally set as $c/\sqrt{3}$ for isothermal flow), $R$ is the gas constant and ${\omega }_i$ is the weight coefficient: (11) ${\omega }_i=\{\begin{array}{c}4/9i=0\\ 1/9i=1,2,3,4\\ 1/36i=5,6,7,8\end{array}$(11) In this research, we choose the two-dimensional nine-velocity (D2Q9) model, in which ${\mathit{e}}_i$ is given by (Lallemand & Luo, 2000): (12) $\begin{array}{rl}& [{\mathit{e}}_0,{\mathit{e}}_1,{\mathit{e}}_2,{\mathit{e}}_3,{\mathit{e}}_4,{\mathit{e}}_5,{\mathit{e}}_6,{\mathit{e}}_7,{\mathit{e}}_8]\\ & =c\cdot \left[\begin{array}{ccccccccc}0& 1& 0& -1& 0& 1& -1& -1& 1\\ 0& 0& 1& 0& -1& 1& 1& -1& -1\end{array}\right]\end{array}$(12)

In Equation (9), the matrix M and diagonal matrix Λ are as follows (Lallemand & Luo, 2000; Li et al., 2018): (13) $\begin{array}{rl}\mathbf{\text{M}}& =\left[\begin{array}{ccccccccc}1& 1& 1& 1& 1& 1& 1& 1& 1\\ -4& -1& -1& -1& -1& 2& 2& 2& 2\\ 4& -2& -2& -2& -2& 1& 1& 1& 1\\ 0& 1& 0& -1& 0& 1& -1& -1& 1\\ 0& -2& 0& 2& 0& 1& -1& -1& 1\\ 0& 0& 1& 0& -1& 1& -1& -1& 1\\ 0& 0& -2& 0& 2& 1& 1& -1& -1\\ 0& 1& -1& 1& -1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& -1& 1& -1\end{array}\right]\end{array}$(13) (14) $\begin{array}{rl}\mathbf{\Lambda }& =\text{diag}({\tau }_{\rho }^{-1},{\tau }_e^{-1},{\tau }_{\varsigma }^{-1},{\tau }_j^{-1},{\tau }_q^{-1},{\tau }_j^{-1},{\tau }_q^{-1},{\tau }_v^{-1},{\tau }_v^{-1})\end{array}$(14) where${\tau }_v$ defines the gas dynamic viscosity with $v=c_s^2({\tau }_v-1/2)$. Based on the simulation tests and the literature, the relaxation times in the diagonal matrix are chosen as ${\tau }_{\rho }^{-1}={\tau }_j^{-1}=1.0$, ${\tau }_e^{-1}={\tau }_{\varsigma }^{-1}={\tau }_{\varsigma }^{-1}=0.8$, ${\tau }_q^{-1}=1.1$. The relaxation time ${\tau }_v^{-1}$ is determined by the shear viscosity.

The macroscopic variables, for instance lubricant density $\rho$ and velocity $\mathit{u}$, in dimensionless form are derived from the distribution function $f_i$ and can be defined as (15) $\rho =\sum _{i=0}^8{f}_i,\mathit{\rho }\mathit{u}=\sum _{i=0}^8{\mathit{e}}_i{\mathit{f}}_i$(15) The corresponding moment function is as follows: (16) $\mathit{m}=\mathbf{\text{M}}\mathit{f}=(\rho ,e,\epsilon ,j_x,{\mathit{q}}_{\mathit{x}},j_y,{\mathit{q}}_{\mathit{y}},{\mathit{p}}_{\mathit{x}\mathit{x}},{\mathit{p}}_{\mathit{x}\mathit{y}}{)}^{\text{T}}$(16)

The viscosity coefficient $\nu$ has a linear relationship with the $\lambda$ in Equation (6): (17) $v=\frac{1}{2}\overline{c}{\lambda }^{\ast }$(17) where $\overline{c}=\sqrt{8RT/\pi }$. According to Equation (6), $\lambda$ can be obtained by $\lambda =KnH$, $H$ is given by $H=N_H\mathrm{\Delta }x$ ($N_H$ is the number of grids, $\mathrm{\Delta }x$ is the grid spacing), and $v=c_s^2({\tau }_v-1/2)$ can be transformed into (18) $v=\frac{1}{2}\sqrt{\frac{8RT}{\pi }}\text{Kn}{N}_H\mathrm{\Delta }x$(18)

Here, the relationship between ${\tau }_v$ and Kn can be acquired as (Guo et al., 2006): (19) ${\tau }_v=\frac{1}{2}+\sqrt{\frac{\pi }{6}}Kn{N}_H$(19) To associate the relaxation time with lubricant density$\rho$, we set (20) ${\tau }_{\nu }^{\mathrm{\prime }}=\frac{1}{2}+\left({\tau }_v-\frac{1}{2}\right)/\rho$(20) Now, the LBM model suitable for continuous flow is set up, and the relationship between the Kn and the relaxation time ${\tau }_v$ is put forward as well. The problem on disposing of the slip boundary will be expounded later.

• Generalized form of interpolation supplemented LBM (GILBM)

In this section, it is our purpose to transform the MRT model governing equation (Equation 9) from uniform Cartesian coordinates into non-uniform curvilinear coordinates. That is, from the physical plane, expressed as $\mathit{x}=(x_1,x_2)=(x,y)$, into computational planes, described as $\mathit{\xi }=({\xi }_1,{\xi }_2)=(\xi ,\eta )$. The process of solving the lattice Boltzmann equation (LBE) can be divided into two steps: collision and advection.

By substituting $\mathit{x}$ in Equation (9) with $\mathit{\xi }$, we can obtain the collision term described in generalized coordinates (Imamura et al., 2005) as (21) $\begin{array}{rl}& f_i^{\ast }(\mathit{\xi },t)\\ & =f_i(\mathit{\xi },t)-\sum _j({\mathbf{\text{M}}}^{-1}\mathrm{\Lambda }\mathbf{\text{M}}{)}_{ij}[f_j(\mathit{\xi },t)-f_j^{(\text{eq})}(\mathit{\xi },t)]\end{array}$(21) After the transformation of the collision step, the advection term calculation is carried out and more discussion will be put forward. In Figure 5, the distributions of the vectors ${\mathit{e}}_1=\mathit{e}(1,0)$ on the physical and computational planes, respectively, are shown.

Figure 5. Discrete velocity vectors on the physical plane (a) and the computational plane (b).

It is obvious that ${\mathit{e}}_1$ is kept the same in the physical plane, whereas it changes with variations of the grid nodes in the computational domain. Then, the contravariant velocities ${\tilde{\mathit{e}}}_{\mathit{i},a}$ are established: (22) ${\tilde{\mathit{e}}}_{\mathit{i},\alpha }={\mathit{e}}_{\mathit{i},\beta }\frac{\mathrm{\partial }{\xi }_{\alpha }}{\mathrm{\partial }{x}_{\beta }}$(22) As a result, the advection term in Equation (9) is transformed into Equation (23): (23) $f_i\left(\mathit{\xi },t+\delta t\right)=f_i^{\ast }\left(\mathit{\xi }-\mathrm{\Delta }{\mathit{\xi }}_{up,i},t\right)$(23) where $\mathrm{\Delta }{\xi }_{\text{up},i}$ is the distance of particle advection and $\delta t$ is the time step. (24) $\mathrm{\Delta }{\xi }_{\text{up},i}=\underset{0}{\overset{\mathrm{\Delta }t}{\int }}\text{d}{\xi }_i=\underset{0}{\overset{\mathrm{\Delta }t}{\int }}{\tilde{\mathit{e}}}_i\text{d}t$(24) In Equation (23), because the right-hand side is not always on the lattice nodes on the computational plane, the value cannot be exactly captured. Thus it is essential to introduce the interpolation method to ensure calculation accuracy. For the inner nodes, the second-order upwind interpolation method is applied. In contrast, for the boundary-adjacent nodes, on account of the limitation quantity of nodes, the second-order central difference interpolation method is adopted. The stencils used for the interpolation function are displayed in Figure 6, in which the $g_{i,lm}$ are the stencils depending on the positions of $\xi -\mathrm{\Delta }{\xi }_{\text{up},i}$ (Imamura et al., 2005).

Figure 6. Stencils used for the interpolation function.

In order to keep the stability of numerical iteration, the time step $\delta t$ ought to satisfy the Courante–Friedrichs–Lewy(CFL) condition (Imamura et al., 2005). Inspired by the traditional CFD methods, the calculation formula can be written as (25) $\delta t=\lambda \underset{i,\alpha }{min}\left|\frac{1}{{\tilde{\mathit{e}}}_{i,\alpha }|\xi |}\right|$(25) where $\lambda$ is the CFL number ($0<\lambda \le 1$); herein, $\lambda =1$. After that, during the advection, the particle will not exceed the domains of adjacent nodes and the stability of interpolation will be ensured.

2.2.2. Mesoscopic boundary conditions

Boundary conditions have a significant effect on the mesoscale LBM for the gas flow (Chen & Tian, 2009; Dongari et al., 2007; Myong, 2004; Wang et al., 2018). When dealing with the slip boundary of microflows, the Maxwell slip model (Dongari et al., 2007) is one of the most commonly seen types of model. Nevertheless, it has some disadvantages, of which the tangential momentum accommodation coefficient (TMAC) is difficult to evaluate, for it depends on the fluid types, surface material properties and geometry shapes. Here, we decide to adopt the Langmuir slip model (Myong, 2004). Without taking the surface diffusion effect into consideration, for the first-order slip velocity $u_{\text{slip}}^{\alpha }$ can be written as (26) $u_{\text{slip}}^{\alpha }=(1-{\kappa }_1)u_f+{\kappa }_1{u}_w$(26) where ${\kappa }_1$ is the first-order slip coefficient, $u_f$ is the local velocity adjacent to the wall and $u_w$ is the local velocity at the wall. Combining Equation (25) with the Maxwell first-order slip model, ${\kappa }_1$ can be derived as (27) ${\kappa }_1=\frac{1}{1-C_{1}KnN}$(27) where $C_1$ is the first adjustment coefficient and $N=1/\delta x$ is the grid number.

As previously mentioned, the first-order slip model for the slip flow regime is set up. However, the accuracy of the first-order slip model does not satisfy the transition regime. We choose to introduce the second-order slip velocity $u_{\text{slip}}^{\beta }$ by the Langmuir model (Wang et al., 2018): (28) $\begin{array}{r}{u}_{\text{slip}}^{\beta }=\frac{1}{1-{\kappa }_2}[(1-{\kappa }_1-2{\kappa }_2)u_f+{\kappa }_1{u}_w+{\kappa }_2{u}_{ff}]\end{array}$(28) where ${\kappa }_2$ is the second-order slip coefficient and $u_{ff}$ is the local velocity adjacent to node $f$. Similarly to the procedure of Equation (27), we can derive the ${\kappa }_2$ as follows: (29) ${\kappa }_2=-C_2{\left(KnN\right)}^2{\kappa }_1$(29) where $C_2$ is the second adjustment coefficient.

According to the non-equilibrium extrapolation (NEE) scheme (Chen & Tian, 2009), the distribution function can be divided into two sections, the equilibrium and the non-equilibrium sections: (30) $f_i({\mathit{x}}_w,t)=f_i^{\text{eq}}({\mathit{x}}_w,t)+f_i^{\text{neq}}({\mathit{x}}_w,t)$(30) We can transform the Langmuir slip model into the NEE scheme. The first-order slip boundary and the second-order slip boundary can be expressed as Equations (31) and (32). Then, we have the slip boundary conditions for the slip flow and the transition flow regimes: (31) $\begin{array}{rl}{f}_i({\mathit{x}}_w,t)& ={\kappa }_1[f_i^{\text{eq}}(\rho ({\mathit{x}}_f,t),\mathit{u}({\mathit{x}}_w,t)-f_i^{\text{eq}}({\mathit{x}}_f,t)]\\ & +f_i({\mathit{x}}_f,t)\end{array}$(31) (32) $\begin{array}{rl}{f}_i({\mathit{x}}_w,t)& =\frac{1}{1-{\kappa }_2}\{{\kappa }_1[f_i^{\text{eq}}(\rho ({\mathit{x}}_f,t),\mathit{u}({\mathit{x}}_w,t))-f_i^{\text{eq}}({\mathit{x}}_f,t)]\\ & +(1-2{\kappa }_2)f_i({\mathit{x}}_f,t)+{\kappa }_2{f}_i({\mathit{x}}_{ff},t)\}\end{array}$(32)

2.2.3. Mesoscopic grids independence verification

Owing to the fact that the method of mesh generation has a great influence on the simulation results, we need to verify the independence of mesoscopic grid partition. The bearing clearance is meshed as shown in Figure 7. As we gradually intensify the grids from 1.38 × 10⁵ to 2.2 × 10⁶, the variation of maximum dimensionless pressure with the number of grids is shown in Figure 8, of which the relevant calculation parameters are shown in Table 3.

Figure 7. Schematic diagram of mesh of mesoscopic method.

Table 3. Calculation parameters for mesoscopic grids independence verification.

Shaft	Radial	Ambient	Bearing
radius, R (m)	clearance, C (m)	pressure, pa (Pa)	speed, ω (rpm)
10.95 × 10^−3	0.05 × 10^−3	1.013 × 10^5	6 × 10^4

Figure 8. Mesoscopic grids independence verification.

It can be seen that the pressure keeps basically constant after the number of grids exceeds 5.5 × 10⁵ and the mean deviation of the numerical result is less than 0.9%. Comprehensively considering the computing time and precision factor, we choose this as the grid independent solution.

3. Model validation

3.1. Macroscopic model validation

The curves are shown in Figure 9 of the lubricant air film pressure and film thickness distributions calculated by the macroscopic method compared with those in the literature (Peng & Khonsari, 2006; Xu, 2012), the data of which are highly consistent with the existing experimental results (Strom, 1988). The relative comparison case parameters are shown in Table 4. The maximum error values occur at the maximum pressure and thickness peaks, which are less than 3%. It can be concluded that the macroscopic method is valid and has high numerical accuracy.

Figure 9. Comparison of macroscopic method simulation results and experimental data in the literature: (a) dimensionless pressure distribution with angular coordinate; (b) dimensionless film thickness distribution with angular coordinate.

Table 4. Parameters of comparison examples in the literature (Peng & Khonsari, 2006; Xu, 2012).

Bearing length,	Shaft radius,	Radial clearance,	Ambient pressure,	Dynamic viscosity of	Eccentricity ratio,	Bearing speed,
L (m)	R (m)	C (m)	pa (Pa)	gas, μ (Pa·s)	ϵ	ω (rpm)	Kn
75 × 10^−3	50 × 10^−3	0.1 × 10^−3	1.013 × 10^5	1.932 × 10^−5	0.84	3 × 10^4	0.0033

3.2. Mesoscopic model validation

The radial velocity profile results of the microcylindrical Couette flow are displayed in Figure 10 compared with those obtained by the mesoscopic LBM and analytic method in the literature (Sun et al., 2005), respectively. The mesh division of LBM is shown in Figure 11.

Figure 10. Comparison of mesoscopic method simulation results and analytical results.

Figure 11. Grids of the microcylindrical Couette flow.

The inner journal with radius R₁ rotates clockwise with angular velocity $\omega$ and the outer cylinder house with radius R₂ keeps stationary, where R₁:R₂=3:5. The analytic model in the literature (Sun et al., 2005) is a combination of the N-S equation and the Maxwell’s diffusion model, expressed as follows: (33) $\frac{u_{\theta }}{\omega R_1}=\frac{1}{(A-B)R_1}\left(Ar-\frac{1}{R}\right)$(33) where (34) $\begin{array}{rl}A& =\frac{1}{R_2^2}\left(1-\frac{(2-{\sigma }_2)}{{\sigma }_2}\frac{2\lambda }{R_2}\right),\\ B& =\frac{1}{R_1^2}\left(1+\frac{(2-{\sigma }_1)}{{\sigma }_1}\frac{2\lambda }{R_1}\right)\end{array}$(34)

Then, Equation (33) can be rewritten as (35) $\begin{array}{rl}{u}_{\theta }& =\omega R\left(1-\frac{(2-{\sigma }_2)}{{\sigma }_2}\frac{2\lambda }{R_2}-\frac{R_2^2}{R}\right)\\ & \times [\left(1-\frac{(2-{\sigma }_2)}{{\sigma }_2}\frac{2\lambda }{R_2}\right)\\ & -{\frac{R_2^2}{R_1^2}\left(1+\frac{(2-{\sigma }_1)}{{\sigma }_1}\frac{2\lambda }{R_1}\right)]}^{-1}\end{array}$(35) where $u_{\theta }$ is the tangential velocity, ${\sigma }_1$ and ${\sigma }_2$ are the tangential momentum accommodation coefficients (TMACs) of the inner and outer cylindrical surfaces (set separately as one) and $R$ is the radial parameter in the cylindrical-polar coordinate $(R,\theta )$ reference frame. The error value of the maximum pressure is less than 3%, which indicates that the introduced mesoscopic LBM obtains reasonable feasibility and accuracy.

4. Results and discussion

4.1. Air film pressure, thickness, velocity distribution and the backflow phenomenon under different eccentricity ratios

It is shown in Figure 12 that when the eccentricity ratios vary from ϵ = 0.3 to ϵ = 0.9, with ω = 5 × 10⁴ rpm, the maximum and minimum lubricating air film pressure peaks will increase and decrease, respectively. Simultaneously, the pressure differences between the two peaks increase as the eccentricity ratios increase. It can be seen from Figure 13 that the variation trends of the maximum, minimum and pressure difference are stable before ϵ = 0.5, the slopes incline more obviously from ϵ = 0.5 to ϵ = 0.7, and the inclinations change rapidly after ϵ = 0.7. Additionally, the tendency of curves simulated by both methods are similar. Nevertheless, some differences exist: first, the width of the peaks simulated by LBM is narrower than that of the Reynolds equation (shown in Figure 12), which represents the fact that the area of the high pressure zone simulated by the former is less than that of the latter; also, the pressure variation is more concentrated. The peaks simulated by LBM appear slightly earlier than those of the Reynolds equation. Third, the minimum pressure peak obtained by LBM is distinctly lower than that of the Reynolds equation at ϵ = 0.9 (shown in Figure 13); however, the pressure differences of the maximum/minimum peaks of results from both methods are close. The errors in the maximum pressure by the two methods are 0.4%, 1.3%, 3.9%, 2.5% at ϵ = 0.3, ϵ = 0.5, ϵ = 0.7 and ϵ = 0.9, respectively. Meanwhile, those in the minimum pressure are 1.4%, 2.5%, 3.1%, 9% at ϵ = 0.3, ϵ = 0.5, ϵ = 0.7 and ϵ = 0.9, respectively.

Figure 12. Air film pressure distributions under different eccentricity ratios.

Figure 13. Variations of pressure under different eccentricity ratios.

These distinctions could be caused by, first, the transformation of the wedge-shape gas flow channel will be intensified by increasing the eccentricity ratio. The shearing force of the flow near the journal surface around the pressure peaks will change more extensively, the friction and viscous dissipation may increase more severely leading to the reduction of the pressure peak areas. Secondly, the gas bearing has a certain negative pressure zone with pressure lower than the ambient pressure, mainly around the minimum pressure peak position of the bearing, where lubricating gas can be sucked in to supply lubricant for bearing. With the increase of tbe eccentricity ratio, not enough lubrication air is sucked into the clearance in the negative pressure zone to fill up the loss caused by the high pressure output, resulting in a decrease of the minimum pressure peak. Thirdly, the slip boundary accelerates the decrease of the maximum air film pressure, making the maximum pressure peak overall narrower than that of the macroscopic result. Additionally, the continuum assumption could be the reason for a higher maximum air film pressure in the macroscopic methods, which is mentioned by Zhang et al. (2016).

To better understand the influence of eccentricity ratios on the bearing lubricant flow characteristics, the air flow streamlines in the maximum clearance are simulated by the mesoscopic LBM, which is correspondingly displayed in Figure 14. On account of the shear stress driving function on the lubricant air film caused by the rotating journal shaft, the velocity of the air film close to the journal surface is higher than that adjacent to the bearing house inner surface. With an increase in the eccentricity ratio, the area of the relative low velocity domain expands; in contrast, the area of the relatively high velocity domain shrinks. Then, at the eccentricity ratio ϵ = 0.3, the streamlines in the clearance are evenly distributed in the same direction. At ϵ = 0.5, vorticity streamlines appear around the bearing house inner surface, indicating backflow occurring on the side close to the bearing sleeve. As it reaches ϵ = 0.7, with the space of maximum clearance becoming larger, and the relative low-velocity air flow increasing, collisions between the relatively high-velocity and low-velocity air flows are aggravated, leading to the air congestion phenomenon in the wedge-shape channel becoming more severe and the vortices appearing more obviously. When the eccentricity ratio reaches ϵ = 0.9, both the vorticity streamlines become more intense and the area of the vortex becomes even larger, and a more intensified backflow phenomenon occurs. This can also explain the drop of the minimum pressure peak. Furthermore, we can explore the reasons for the form of Figure 12 further by means of Figure 15.

Figure 14. Change of air film backflow under different eccentricity ratios.

Figure 15. Flow velocity distribution under different eccentricity ratios.

Figure 16. Three-dimensional pictures of air film pressure distribution under different eccentricity ratios.

As shown in Figure 15, the dimensionless velocity varies with the dimensionless radial distance under different eccentricity ratios. At the low eccentricity ratio ϵ = 0.3, with the dimensionless radial distance shifting from the journal to the bearing house, the velocity of lubricant air flow declines approximately linearly, which is similar to the velocity distribution of a typical traditional Couette flow. When it reaches ϵ = 0.5, the air flow velocity at first goes down then lifts again around a relative position of 0.68, indicating that an inverse flow occurs. When it goes from ϵ = 0.7 to ϵ = 0.9, the reversal trends of the velocity flow appear more obvious, with the velocity difference turning larger. Because, the larger the eccentricity ratio is, the higher the maximum air film thickness will be, combined with the greater effect of the inverse pressure gradient, the shear force becomes insufficient, and more violent backflow occurs. Furthermore, the squeeze function of the increasing relatively-low velocity air flow to that of the decreasing relatively-high velocity air flow strengthens the phenomenon, which could be one of the primary causes of the friction loss in bearing operation reported by Serrano et al. (2020).

Displayed in Figures 16 and 17 are the corresponding three-dimensional cloud pictures of air film pressure and thickness distributions changing under different eccentricity ratios, from ϵ = 0.3 to ϵ = 0.9 with ω = 5 × 10⁴ rpm, by the macroscopic method. It can be seen that, as the angular coordinate θ gradually increases from 0° to 360°, the air film pressure distribution exhibits an increase–decrease–increase trend. In contrast, the air film thickness distribution presents a decrease–increase trend. With the eccentricity ratios increasing from ϵ = 0.3 to ϵ = 0.9, the maximum air film pressure peaks elevate and the minimum air film pressure peaks reduce. While the minimum air film pressure peaks reduce, the maximum air film pressure peaks elevate. Those appearances could be caused for the following reasons. First, for the eccentricity ratio representing the bearing load, the higher the eccentricity ratio is, the heavier the load on the journal will be, and the stronger the squeezing action on the lubricant air film will be. On account of the compressibility of gas, the minimum thickness of the air film shrinks, leading to an intensified cross-section effect variation. As a consequence, the maximum pressure peak is enhanced near the maximum air film pressure peak. Secondly, at high eccentricity ratios, there will be much more space near the maximum air film thickness position, which makes the minimum air film pressure peak reduce. The difference in pressure between the high-pressure and low-pressure zones expands, leading to more pressure gap existing between the pressure peaks. Moreover, the sharp jump in the maximum air film pressure peak at ϵ = 0.9 is caused by the severe compression of the air film in the clearance, and the interval is too narrow to flow through, causing the air pressure to increase abruptly, which is in agreement with the results reported in the literature (Shalash & Schiffmann, 2020). The study indicates that, by keeping the eccentricity ratio at an appropriately high level, the working performance of the GLJB will be improved.

Figure 17. Three-dimensional pictures of air film thickness distribution under different eccentricity ratios.

4.2. Air film pressure, velocity distribution under different bearing speeds

The change of the air film pressure distribution withdifferent bearing speeds is shown in Figure 18, which is simulated by both the macroscopic and mesoscopic methods at bearing speeds ranging from ω = 5 × 10⁴ rpm to ω = 1.5 × 10⁵ rpm and ϵ = 0.5. With increasing bearing speed, the peaks of the maximum and the minimum air film pressure will increase and decrease, respectively. The pressure difference between the two peaks of the air clearance enlarges, as displayed in Figure 19. By comparison, the pressure distribution obtained by the two methods are generally consistent with each other. However, the errors of the results change with the speeds—the errors in the maximum pressure are 1.3%, 1.2% and 0.81% at ω = 5 × 10⁴ rpm, ω = 1 × 10⁵ rpm and ω = 1.5 × 10⁵ rpm, respectively. Meanwhile, those in the minimum pressure are 0.08%, 2.8% and 9.6%, respectively.

Figure 18. Air film pressure distribution under different bearing speeds.

Figure 19. Variations of pressure under different bearing speeds.

We need to explore the reasons explaining this phenomenon. As shown in Figure 20, with the enhancement of journal speeds, the velocity of gas flow is accelerated and the radial gradient of the lubricating film airflow velocity is increased, thereby giving rise to a growth in the tangential shear force, and causing an increase in the viscous friction force. To some extent, the higher the speed, the greater the radial gradient of the lubricating film air flow velocity will be, and the larger the tangential shear force that will form. As a result of the shear force function, the gas flow velocity difference in the boundary layer near the journal surface is enlarged. As is shown in Figure 21, the bearing capacity of the GLJB will be enhanced. It can be concluded that an increase of bearing speed within a certain range is beneficial for the GLJB to undertake a higher bearing load at relative lower eccentricity ratio and larger attitude angle. However, it is worth noting that the bearing speed is not simply ‘the higher the better’. When the bearing speed is high enough, the slip function of lubrication flow near the journal surface is strengthened, and gas congestion can form around the pressure peaks, as mentioned in the literature (Zhang et al., 2016). Excessive bearing speed can lead to amounts of friction loss that will negatively influence the operating stability and working life of the GLJB. As a consequence, the working performance of the GLJB will not always be improved with a rise in the bearing speed.

Figure 20. Flow velocity distribution under different bearing speeds.

Figure 21. Variation of eccentricity ratio with bearing loads under different bearing speeds.

Figure 22. Variation of attitude angle with bearing loads under different bearing speeds.

Figure 21 illustrates the relationship between the eccentricity ratio and the bearing load at different bearing speeds. When the bearing load rises, the eccentricity ratios will first lift rapidly to a relative high level, then the increase in the slopes of the eccentricity ratio become gentle and linear. It is because, at the beginning of the loading effect, the air flow will be greatly compressed to small volumes, resulting the rapid growth of the eccentricity ratio. After that, due to the limited compressibility of air and the clearance space, gas compression becomes difficult and the distance between the journal and the bearing quite narrow, then the growth of the eccentricity ratio is slight. Moreover, to tolerate the same load on the bearing, a higher speed rotating bearing needs a lower eccentricity ratio than that of lower speed rotating bearing, which is a function of kinetic energy, which can counteract the effect of bearing load.

As shown in Figure 22, the attitude angle changes with the variation in bearing load at different bearing speeds. It can be observed that the larger the bearing speed, the less obvious the variation of the attitude angle will be. When the bearing load is higher than 40 N, all the attitude angels are less than 90°, which indicates that the effect of the difference in bearing load on the attitude angle is shrinking. According to fluid lubrication theory, the smaller the attitude angle, the more stable the bearing will be, and the less the vertical component of the force on the lubricant film will be, which is usually through the centrelines. Thus, to some extent, with an increase in bearing speed, the stability of the journal bearing will become better.

Figure 23. Air film pressure distribution under different values of Kn.

4.3. Air film pressure distribution under different Kn

Variations in the lubricating air film pressure distributions with Kn obtained by the LBM are displayed in Figure 23. As an important parameter to characterize channel size and fluid characteristics, Kn can be used to classify the gas flow situations into different flow regimes (Sun et al., 2018). Within the tiny journal-bearing house clearance, the Kn of lubricating gas flow could rise to a high value. Owing to the continuum assumption, the conventional macroscopic method is not available at high Kn conditions, and the advantage of the mesoscopic method is demonstrated. In this research, Kn varies from 0.0605 to 0.303 in Figure 23(a) and from 1 to 4 in Figure 23(b), involving the slip flow and the transition flow regimes, with the aim of exploring the effect of Kn on the bearing pressure distribution. It is shown in Figure 23 that, as Kn increases, the maximum and minimum pressure peaks will decline and rise, and the pressure gap between the two peaks will diminish. Additionally, the variation speed of the maximum pressure peaks will change from high to low with the enhancement of Kn. Comparing Figures 23(a) and 23(b), it can be seen that, at relative large Kn, the pressure transformation is different from that at relatively low Kn, although both of them belong to the transition regime. That is, if Kn is high enough, the influence of Kn on the pressure will be limited. This phenomenon could be caused by the increased Kn leading to a more pronounced rarefaction effect, making the wall slip effect of the bearing to strengthen, and the bearing's shear force on the gas to be decreased, resulting in decreased maximum pressure peaks. Moreover, if Kn is too high, the residual space for the rarefaction effect on the lubricant pressure will be too narrow, which will give rise to a very small difference in pressure peaks. To sum up, Kn will have an effect on the pressure distributions; however, in the tiny clearance of a GLJB, the degree of influence of the rarefaction effect on the pressure distribution will reduce with an increase in Kn. Furthermore, it indicates that the LBM can be introduced as an ideal method for research in high Kn conditions, which is of inherent compatibility with the slip boundary.

5. Conclusion

In this research, we propose a multiscale numerical investigation method combining the macroscopic Reynolds equation and the mesoscopic LBM to simulate the air film flow in the narrow clearance of a GLJB. At different bearing speeds, eccentricity ratios and slip conditions (Kn), the changing rules of the air film pressure, thickness and velocity distributions are investigated. The differences between the results of the two methods are compared and analysed. In addition, an interesting backflow phenomenon in the clearance is observed by the LBM and the reason for its formation is discussed. To sum up, conclusions are given as follows.

1. For the Reynold equation and the LBM, the former has a clear physical meaning and can be used as a proof for the latter in the simulation of GLJBs. In contrast, the latter can be applied to capture the velocity distributions and backflow phenomenon at the multiscale level as a complement to the former.

2. The multiscale model can be adopted to simulate the air flow of a GLJB at high Kn conditions in both the slip or transition flow regimes, which overcomes the disadvantage of the continuum hypothesis. The LBM can be an ideal method for research on the rarefaction effect.

3. When the eccentricity ratios increase from ϵ = 0.3 to ϵ = 0.9, the errors in the maximum pressure peaks by the two methods are 0.4%, 1.3%, 3.9% and 2.5%, respectively. In contrast, those in the minimum pressure peaks are 1.4%, 2.5%, 3.1% and 9%, respectively. When bearing speeds accelerate from ω = 5 × 10⁴ rpm to ω = 1.5 × 10⁵ rpm, the errors in the maximum pressure peaks are 1.3%, 1.2%, 0.81%, respectively. In contrast, those in the minimum pressure peaks are 0.08%, 2.8%, 9.6%, respectively.

4. This study indicates that, by keeping the eccentricity ratio and bearing speed at appropriately high levels, the working performance of GLJBs will be improved. However, the performance will not always be promoted with a rise in bearing speed. When Kn enhances, the maximum and minimum pressure peaks will decrease and increase, respectively; however, the extent of influence of the rarefaction effect on the pressure distribution will reduce after Kn exceeds a certain degree.

Multiscale numerical simulation methods can provide fruitful new insights in investigation, and be applied in theoretical support for improving lubricating characteristics, optimizing operating performance and promoting the engineering application of GLJBs. Although, at the multiscale level, the characteristics of gas flow in the micro clearance of GLJBs have been obtained successfully, experimental results are required to evaluate the utilized models more accurately. Future work should investigate various shapes and geometries for the bearing house. Furthermore, researching the influence of surface roughness could be helpful in acquiring insight into its impact so as to achieve better stability and carrying capacity.

Nomenclature

Symbols

C	=	Radial clearance (m)
Co	=	Radial average clearance (m)
c	=	Lattice velocity
cs	=	Lattice speed of sound
$\overline{c}$	=	Molecular average motion speed
D	=	Diameter of journal (m)
d	=	Effective molecular diameter (m)
e	=	Eccentricity
ei	=	Discrete velocity
Fi	=	External forcing term
$f_i$	=	Particle density distribution function
H	=	Dimensionless film thickness
h	=	Film thickness, gap height (m)
h0	=	Fluid characteristic size (m)
Kn	=	Knudsen number
k	=	Boltzmann's constant, 1.38 × 10−23 (J/K)
L	=	Bearing length (m)
M	=	Orthogonal transformation matrix
M−1	=	Inverse orthogonal transformation matrix
NH	=	Number of lattices
n	=	Number density of gas molecules
P	=	Dimensionless pressure
pa	=	Ambient pressure (Pa)
p	=	Pressure of the gas (Pa)
R	=	Radial pressure in the cylindrical-polar coordinate (m)
Rb	=	Radius of bearing (m)
Rj	=	Radius of journal (m)
S	=	Diagonal relaxation matrix
T	=	Temperature of the gas (K)
t	=	time (s)
$\overline{t}$	=	Dimensionless time
U	=	Tangential velocity of bearing (m/s)
u	=	Velocity vector
uf	=	Local velocity adjacent to the wall
uslip	=	Slip velocity
uw	=	Local velocity at wall
W	=	Bearing load (N)
x	=	Circumference direction
y	=	Radial direction
Z	=	Dimensionless bearing length
z	=	Axial direction

Greek

$\gamma$	=	whirl frequency ratio
Δx	=	Grid spacing
δt	=	Time step
ϵ	=	Eccentricity ratio
κ1, κ2	=	Slip coefficients
λ	=	Molecular free path (m)/CFL number
μ	=	Dynamic viscosity (Pa/s)
θ	=	Angular coordinate (rad)
Λ	=	Diagonal matrix
Λ	=	Bearing number
τ	=	Relaxation time that relies on the fluid kinematic viscosity
τs	=	Relaxation time
τv	=	Gas dynamic viscosity
ξ, η	=	Coordinates on the computational plane
φ	=	Attitude angle (rad)
ω	=	Bearing speed (rpm)
ωi	=	Weight coefficient
${\omega }_s$	=	Angular frequency of oscillation

Superscripts

*	=	Probability
eq	=	Equilibrium
1	=	First-order
2	=	Second-order

Subscripts

b	=	Bearing
j	=	Journal
i	=	Interacting atoms or component
o	=	Original/preload
s	=	Sound

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China [grant number No. 52006031].