Research on Static and Dynamic Characteristics of Gas-lubricated Microbearings for Microfluidic Devices

Aerodynamic journal microbearings with microtribological phenomena significantly influence the operating stability of microfluidic devices. The modified Reynolds equations including different rarefaction models are derived and are solved by the partial derivative method and relaxation iteration algorithm. The effects of Knudsen number and bearing parameters on the static and dynamic characteristics of microbearings are investigated in detail. The results show that the rarefaction effect plays a crucial role in the ultra-thin gas film lubrication. The maximum film pressure of Fukui–Kaneko (FK) model is lowest and the result in Boltzmann model is largest in small Knudsen number regions. As the Knudsen number increases further, the curve of FK model is coincident with that of the Boltzmann model in the transition regime. The direct stiffness coefficients of Boltzmann model increase with the increase of eccentricity ratio and aspect ratio, whereas the effect of Knudsen number on the damping coefficients yields a relatively complicated trend.


Introduction
Gas journal microbearing is of great importance in various microelectromechanical systems (MEMS) devices such as ultrasonic motors, micro-gas turbines, microaccelerometers, inertial navigation systems and magnetic hard disk drives (HDD) [1][2][3]. Due to the continuous development of rotating machinery for high precision, miniaturization, low energy consumption and long life span. Aerodynamic microbearings are widely used owing to their advantages of simple structure, low friction resistance, no contamination and minimal temperature limitation in comparison with the oil-lubricated bearings and electromagnetic bearings [4][5][6][7]. In these microfluidic applications, the lubricating films between the journals and bearing pads usually become very small and are found to be even smaller than the molecular mean free path [8,9], and the continuum flow of Reynolds equation is no longer valid for describing the rarefied gas flow characteristics. Thus the effect of gas rarefaction must be taken into consideration. The rarefied gas flow can be classified into three categories [10,11] on the basis of the Knudsen number K n , which is defined as the ratio of the gas molecular mean free path λ 0 to the characteristic length scale L, i.e., the slip region (0.001 < K n < 0.1), the transition region (0.1 < K n < 10), and the free molecular region (K n > 10). In the slip region, the velocity slip and temperature jump occur at the solid wall [12][13][14]. In the transition region, both the mutual collisions between the molecules and between gas molecule and boundary wall surface should be considered. The variation of microscopic particle velocity which results from collision can be ignored in the free molecular region and the Boltzmann equation without collision term and the Maxwellian distribution can be employed to study the gas flow in this region.
In most of the theoretical investigation of ultra-thin gas film lubrication, the Poiseuille flow rate coefficient, which is an important flow correction term, is incorporated into the compressible Reynolds equation to account for gaseous rarefaction effects. Burgdorfer [15] firstly introduced the first-order slip velocity boundary condition into the Reynolds equation for small Knudsen numbers in 1959. Other noted researchers, including Hsia [16], Mitsuya [17], Fukui, Kaneko [18], and Beskok [19], derived the modified Reynolds equations by incorporating their different Poiseuille flow factors for arbitrary Knudsen numbers. Subsequently, a number of approaches, based on the kinetic theory, to the practical analysis of rarefaction effects on micro fluidic systems have been reported that have created discussion. Hwang et al. [20] proposed a new modified Reynolds equation which includes three adjustable coefficients in the velocity boundary conditions, and the results showed that the model produces good approximation to the Boltzmann model for a wide range of Knudsen numbers. Kang et al. [21] extended the databases of Poiseuille flow rate and Couette flow rate by introducing new terms into system equations at track recording head/disk interface in hard disk drives by using the Bhatnagar-Gross-Krook kinetic equation. In 2007, Shen et al. [22] proposed a new first-order slip model which includes an additional term by utilizing the Chapman-Enskog solution of Boltzmann equation. In recent years, the problems connected with rarefaction effects for gas journal microbearings had also been investigated taking into account surface roughness, bearing pad deformation, turbulent flow, bearing thermal expansions and shaft misalignment effects. White [23] developed an analytical method that employs the perturbation method and multiple-scale analysis to discuss the influences of surface roughness and gas rarefaction on low clearance gas bearing with high wave number and high bearing number. Zhang et al. [24] adopted the finite volume method and atomic force microscope data to study the random surface roughness effect on the flow characteristics in gas microbearing. The calculations showed that the larger roughness exponent causes the bigger variation of the gas slip velocity at the bottom wall. Yang et al. [25] derived an extended slip velocity boundary condition, treating the slip coefficients as functions of wall accommodation coefficient and Knudsen number through the regularized 13 moment equations. The work showed that the solution of improved Reynolds equation has more accurate prediction of pressure distribution in the near transition flow regime. The effects of rarefaction and manufacturing defects on the linear and nonlinear bearing performances of micro flexure pivot tilting pad gas bearings were studied by Li et al. [26], they emphasized the importance of pad preload on the dynamic characteristics of power MEMS. Wang et al. [27] used the Greenwood-Williamson model and modified Reynolds equation to calculate the bearing force and asperity contact force, their studies concerned gas spiral groove thrust microbearing with different groove depth and standard deviation of asperity height. Research into the lubrication performance of textured journal bearings was stimulated by Ma et al. [28] who concluded that the concave spherical texture decreases the load capacity while the convex texture increases it. Wang and Ding [29] analyzed the stability characteristics of gas bearing-soft rotor system by two way coupled fluid structure interaction and derived the stiffness coefficients with the influence coefficients and harmonic excitation methods. They remarked that the stabilities of herringbone groove gas bearing and the plain gas bearing enhance as the bearing capacity increases. Zhang et al. [30] numerically and experimentally studied the characteristics of hybrid bearings including friction between silicon nitride and bearing steel by molecular dynamics methods. They showed that the effect of graphene on the lubrication of friction pair becomes greater when the pressure and shear velocity increases. However, probably due to the irregular oscillation and the drawback of poor convergence easily appeared in the finite difference methods for the dynamic analysis, a relatively small number of papers have focused on the dynamic behaviors of microbearings in the open literature. In addition, understanding the flow dynamics of rarefied gas flow is essential for the design of micro gas bearing and the improvement of its characteristics is important to guarantee the stability of micro-rotating machinery. Therefore, Appl. Sci. 2020, 10, 1634 3 of 18 a proper method with less complexity is necessary. Systematic investigations on the stiffness and damping coefficients of gas journal microbearings should be undertaken.
In this paper, the partial derivative method and iteration technique are applied to solve the modified Reynolds equations with different Poiseuille flow correctors, and the static and dynamic characteristics of micro gas journal bearings such as gas pressure distributions, dimensionless load capacity, friction coefficient and dynamic stiffness and damping coefficients are obtained for various Knudsen numbers. This work maybe provides new perspectives in understanding the dynamic properties of aerodynamic journal microbearings, which lays a foundation for the further stability of microbearing-rotor systems.

Gas rarefaction Models
In previous studies, several mathematical models describing rarefied gas flow in HDDs or MEMS applications were established. The first-order slip model of Burgdorfer, the second-order slip model of Hsia, the FK model of Fukui and Kaneko and the Boltzmann model are widely used in ultra-thin gas film lubrication. For small clearance spacing of the gas film, the Poiseuille flow rate correctors can be mainly expressed as follows: The first-order velocity slip model [15]: The second-order velocity slip model [16]: The modified Boltzmann model [31]: The Fukui-Kaneko (FK) model [18]: 2K n is the inverse Knudsen number, K n = λ 0 /L is the Knudsen number, and λ 0 is the mean free path of gas molecules, which takes the value of 65 nm. L is the characteristic length scale which is often replaced by gas-film thickness in the microbearing. The Poiseuille flow rate ratio Q, which is defined as the ratio of rarefied gas flow Poiseuille flow rate to continuum flow Poiseuille flow. For the above modified models that are used to predict the Poiseuille flow rate, the Boltzmann model and FK model, which are derived from the relationship between distribution function of gas molecules and the change rate of spatial location and time, are regarded as the accurate model for simulating ultra-thin film gas lubrication problems in engineering practice.

Governing Equations for Ultra-Thin Gas Film Lubrication
The schematic diagram of gas-lubricated journal microbearing is shown in Figure 1. For an isothermal, isoviscous, laminar and inertialess ideal gas, the compressible Reynolds equation which governs the gas pressure profiles considering the gaseous rarefaction effects is given below.
where p = p/p a is the dimensionless gas film pressure, H = h/C b is the dimensionless gas film thickness, p a and C b are the ambient pressure and the radius clearance respectively. ϕ and λ denotes the dimensionless circumferential angle coordinate and axial coordinate of micobearing normalized by the journal radius R. Λ = 6µωR 2 /p a C b 2 is the bearing number or compressibility parameter, µ is the dynamic viscosity of airflow, ω is the rotational speed of journal. Q represents the Poiseuille flow rate ratio and T is the dimensionless time.

Governing Equations for Ultra-Thin Gas Film Lubrication
The schematic diagram of gas-lubricated journal microbearing is shown in Figure 1. For an isothermal, isoviscous, laminar and inertialess ideal gas, the compressible Reynolds equation which governs the gas pressure profiles considering the gaseous rarefaction effects is given below. (5) where p = p/pa is the dimensionless gas film pressure, H = h/Cb is the dimensionless gas film thickness, pa and Cb are the ambient pressure and the radius clearance respectively. φ and λ denotes the dimensionless circumferential angle coordinate and axial coordinate of micobearing normalized by the journal radius R. Λ = 6μωR 2 /paCb 2 is the bearing number or compressibility parameter, μ is the dynamic viscosity of airflow, ω is the rotational speed of journal. Q represents the Poiseuille flow rate ratio and T is the dimensionless time. In the steady state equilibrium condition, the dimensionless film thickness and pressure distribution are time-independent, so the static modified Reynolds equation considering the gaseous rarefaction effect is simplified as follows. 3 3 () The dimensionless pressure boundary condition is depicted as: where B is the bearing width.
Letting PH = S, (PH) 2 = S 2 = Π [32], the governing Equation (6) can be transformed into the elliptic partial differential equation form -▽·(c▽u) + au = f as In the steady state equilibrium condition, the dimensionless film thickness and pressure distribution are time-independent, so the static modified Reynolds equation considering the gaseous rarefaction effect is simplified as follows.
The dimensionless pressure boundary condition is depicted as: where B is the bearing width.
Letting PH = S, (PH) 2 = S 2 = Π [32], the governing Equation (6) can be transformed into the elliptic partial differential equation form -·(c u) + au = f as The corresponding coefficients in Equation (8) are derived as The normalized load capacity C L of gas lubricated microbearing is given by integrating the gas-film pressure distribution.
where W is the dimensional load-carrying capacity in the vertical direction. The dimensionless circumferential friction on the journal surface is expressed as It is assumed that the journal whirl about its own static equilibrium position with the frequency ratio Ω, and positions of the eccentricity ε and attitude angle θ of journal consist of two parts: the static component and the dynamic component. The static equilibrium position is supposed to be (ε 0 , θ 0 ) and the dynamic small disturbance with regard to (ε 0 , θ 0 ) are indicated as E and Θ [33]. The eccentricity ratio ε and attitude angle θ of journal at a random position are given as: where E 0 is defined as the amplitude of perturbed eccentricity ratio in the complex field, Θ 0 is the small perturbation amplitude of the journal attitude angle with imaginary number domain. The dimensionless perturbation frequency Ω, which is defined as the ratio of journal disturbance frequency ν to journal rotational speed ω, i = √ −1. Therefore, the non-dimensional gas-film thickness and pressure between the bearing and journal surface can also be written as where H 0 = E 0 cos(ϕ − θ 0 ) + ε 0 Θ 0 sin(ϕ − θ 0 ). P 0 and H 0 are static gas-film thickness and pressure. Q gd and H gd are described as dynamic gas-film thickness and pressure, respectively. P 0 and H 0 are complex numbers which are the perturbation magnitudes of dynamic gas-film pressure and thickness. Substituting Equation (13) into Equation (1), and using Equation (6) for simplification with neglecting all higher order terms, the generalized dynamic Reynolds equation for ultra-thin gas-film lubrication appears as Dynamic stiffness and damping coefficients of ultra-thin gas film bearing under small journal disturbance can be calculated by solving Equation (14) using the partial derivative method. Setting Appl. Sci. 2020, 10, 1634 6 of 18 Differentiating Equation (14) with respect to E 0 and Θ 0 , the partial differential equations concerning the variables P E , P θ , H E and H θ are obtained through some mathematical transformation.
The boundary conditions for the above equations can be expressed as follows.
The above Equations (15)- (20) are solved simultaneously by using the finite element method and relaxed iterative algorithm, and P E and P θ are obtained. According to the coordinate system as shown in Figure 1, the stiffness coefficients K ij and damping coefficients D ij of gas microbearing are given by the following formula.
The dynamic coefficients in Cartesian coordinates are computed by the transformation matrix A.

Results and Discussion
The effects of gaseous rarefaction on the static and dynamic characteristics such as pressure distribution, load-carrying capacity as well as dynamic stiffness and damping coefficients with respect to different bearing parameters and perturbation frequencies are compared and analyzed by using the numerical scheme above.
To validate the method developed in this paper, the results are compared with Zhang et al. [34] based on the first order slip boundary condition. Figure 2 shows the dimensionless steady-state pressure profile along the bearing midplane with different Knudsen numbers, respectively. The geometry parameters of micro gas journal bearing are R = 2.0 mm, B = 0.4 mm, p a = 1.033 × 10 5 N/m, Λ = 2.4 and ε = 0.8. It can be seen that the present numerical results are in good agreement with those of Zhang et al. and errors are less than 3%. In the following discussions, the continuum flow model means that the effects of gas rarefaction effects is not considered.  Figure 3 demonstrates the effect of Knudsen number on the pressure profiles in the circumferential direction at the bearing mid-plane for different rarefaction models with ε = 0.6, Λ = 4 and B/D = 0.3. The normalized pressure distribution has similar changing tendency along the sliding direction. The gas film pressure decreases significantly as the Knudsen number increases, namely, the pronounced degree of the rarefied gas. It is found that the maximum gas film pressure of FK model is lowest, and the result of Boltzmann model is largest for the Kn = 0.0325 case. Nevertheless, when Kn = 0.325, by comparison with Boltzmann model, the first-order slip model overpredicts the dimensionless pressure profile while the second-order slip model underestimates it in the transition regime. The FK model is closer to the Boltzmann model than other rarefaction models. Figures 4 and 5 compare the dimensionless load capacities CL and attitude angles θ of Boltzmann model with eccentricity ratio ε for different values of Knudsen number. The CL is seen to increase with ε while the reverse trend is observed for the attitude angle. In comparison with the continuum flow case, the influence of gas rarefaction decrease the load carrying capacity. Therefore, the larger Knudsen number will result in the lower load capacity of microbearing.  Figure 3 demonstrates the effect of Knudsen number on the pressure profiles in the circumferential direction at the bearing mid-plane for different rarefaction models with ε = 0.6, Λ = 4 and B/D = 0.3. The normalized pressure distribution has similar changing tendency along the sliding direction. The gas film pressure decreases significantly as the Knudsen number increases, namely, the pronounced degree of the rarefied gas. It is found that the maximum gas film pressure of FK model is lowest, and the result of Boltzmann model is largest for the K n = 0.0325 case. Nevertheless, when K n = 0.325, by comparison with Boltzmann model, the first-order slip model overpredicts the dimensionless pressure profile while the second-order slip model underestimates it in the transition regime. The FK model is closer to the Boltzmann model than other rarefaction models. Figures 4 and 5 compare the dimensionless load capacities C L and attitude angles θ of Boltzmann model with eccentricity ratio ε for different values of Knudsen number. The C L is seen to increase with ε while the reverse trend is observed for the attitude angle. In comparison with the continuum flow case, the influence of gas rarefaction decrease the load carrying capacity. Therefore, the larger Knudsen number will result in the lower load capacity of microbearing. Appl. Sci. 2019, 9, Figure 6 indicates the relationship between the bearing number and non-dimensional load carrying capacity at different Knudsen numbers for ε = 0.6 and B/D = 1.4. It is noted that the load capacities obtained from all models increase when the bearing number Λ increases, whereas the magnitude of CL decreases dramatically with the increase of Knudsen number. The results with and without considering rarefied gas lubrication are approximately the same when Kn is equal to 0.0054. This is because the characteristic length scale of gas flow is less than 0.01 which the continuum hypothesis remains valid. As Kn increases to 0.0325 in the slip regime, the dimensionless load capacities of different modified models lie below that of the continuum model with nonslip boundary condition. The reason is that the gas flow is less restricted with the increasing slip velocity. Furthermore, the prediction of load carrying capacity for FK model is the smallest and for Boltzmann model is the largest. The curves of first-order and second-order slip model are situated between them. With the continuous increase of Kn, the dimensionless load capacity shows a linearly increasing trend when the bearing number becomes larger. In the meantime, the result of FK model is coincident with that of Boltzmann model, while the first-order slip model overestimates and the second-order slip model underpredicts in comparison with the Boltzmann model.    Figure 6 indicates the relationship between the bearing number and non-dimensional load carrying capacity at different Knudsen numbers for ε = 0.6 and B/D = 1.4. It is noted that the load capacities obtained from all models increase when the bearing number Λ increases, whereas the magnitude of C L decreases dramatically with the increase of Knudsen number. The results with and without considering rarefied gas lubrication are approximately the same when K n is equal to 0.0054. This is because the characteristic length scale of gas flow is less than 0.01 which the continuum hypothesis remains valid. As K n increases to 0.0325 in the slip regime, the dimensionless load capacities of different modified models lie below that of the continuum model with nonslip boundary condition. The reason is that the gas flow is less restricted with the increasing slip velocity. Furthermore, the prediction of load carrying capacity for FK model is the smallest and for Boltzmann model is the largest. The curves of first-order and second-order slip model are situated between them. With the continuous increase of K n , the dimensionless load capacity shows a linearly increasing trend when the bearing number becomes larger. In the meantime, the result of FK model is coincident with that of Boltzmann model, while the first-order slip model overestimates and the second-order slip model underpredicts in comparison with the Boltzmann model.   Figure 6 indicates the relationship between the bearing number and non-dimensional load carrying capacity at different Knudsen numbers for ε = 0.6 and B/D = 1.4. It is noted that the load capacities obtained from all models increase when the bearing number Λ increases, whereas the magnitude of CL decreases dramatically with the increase of Knudsen number. The results with and without considering rarefied gas lubrication are approximately the same when Kn is equal to 0.0054. This is because the characteristic length scale of gas flow is less than 0.01 which the continuum hypothesis remains valid. As Kn increases to 0.0325 in the slip regime, the dimensionless load capacities of different modified models lie below that of the continuum model with nonslip boundary condition. The reason is that the gas flow is less restricted with the increasing slip velocity. Furthermore, the prediction of load carrying capacity for FK model is the smallest and for Boltzmann model is the largest. The curves of first-order and second-order slip model are situated between them. With the continuous increase of Kn, the dimensionless load capacity shows a linearly increasing trend when the bearing number becomes larger. In the meantime, the result of FK model is coincident with that of Boltzmann model, while the first-order slip model overestimates and the second-order slip model underpredicts in comparison with the Boltzmann model.   The influence of aspect ratio and Knudsen number on load capacity and friction coefficient of Boltzmann model are shown in Figures 7 and 8 for ε = 0.6 and Λ = 5. As the value of aspect ratio increases, it is noted that both the load capacity and friction coefficient increase. The CL firstly increases and then subsequently asymptote to a constant value for a small Knudsen number regime. When the Knudsen number is further increased, the load capacity and friction coefficient are nearly proportional to the aspect ratio, and it causes a slight change in dimensionless load capacity at nanometer scale. This is because the lower clearance may weaken the shear stress of the rarefied gas flow and the extent of the flow in the sliding direction is quickly decreased.    The influence of aspect ratio and Knudsen number on load capacity and friction coefficient of Boltzmann model are shown in Figures 7 and 8 for ε = 0.6 and Λ = 5. As the value of aspect ratio increases, it is noted that both the load capacity and friction coefficient increase. The C L firstly increases and then subsequently asymptote to a constant value for a small Knudsen number regime. When the Knudsen number is further increased, the load capacity and friction coefficient are nearly proportional to the aspect ratio, and it causes a slight change in dimensionless load capacity at nanometer scale. This is because the lower clearance may weaken the shear stress of the rarefied gas flow and the extent of the flow in the sliding direction is quickly decreased. The influence of aspect ratio and Knudsen number on load capacity and friction coefficient of Boltzmann model are shown in Figures 7 and 8 for ε = 0.6 and Λ = 5. As the value of aspect ratio increases, it is noted that both the load capacity and friction coefficient increase. The CL firstly increases and then subsequently asymptote to a constant value for a small Knudsen number regime. When the Knudsen number is further increased, the load capacity and friction coefficient are nearly proportional to the aspect ratio, and it causes a slight change in dimensionless load capacity at nanometer scale. This is because the lower clearance may weaken the shear stress of the rarefied gas flow and the extent of the flow in the sliding direction is quickly decreased.      It is seen that the perturbation frequency has great influence on the dynamic coefficients. As the perturbation frequency increases, the dynamic stiffness and damping coefficients Kxx, Kxy, Kyy and Dyx increase. The increasing extent of Kyy is bigger than Kxx because of the gas film mainly supports the weight of high-speed rotor along the vertical direction, while Kyx, Dxx, Dxy and Dyy decrease quickly. With gas rarefaction effects taken into account, the direct stiffness coefficients Kxx and Kyy are linearly proportional to the perturbation frequency, but the direct terms of the damping coefficients Dxx and Dyy decrease slightly with the rise of frequency. Meanwhile, the cross coupled stiffness coefficients Kxy and Kyx gradually approach to the constant values with increasing Ω, yet the Dxy and Dyx are hardly affected by perturbation frequency. The dynamic stiffness coefficients for rarefaction models are less than those of continuum flow. When Ω > 2.0, the direct damping coefficients with the effects of gas rarefaction are higher than those of continuum flow. The reason might be that the energy transport by gas molecule collision is less, to some extent, which gives a larger value of fluid film stiffness and a larger gas film damping coefficient. The first-order slip model overestimates the stiffness coefficients of micro gas journal bearing compared with the simulation results of Boltzmann model, while the second-order velocity slip model underestimates it. Because the FK model is obtained based on the linearized Boltzmann equation, its results are in good agreement with those of the Boltzmann model. It is seen that the perturbation frequency has great influence on the dynamic coefficients. As the perturbation frequency increases, the dynamic stiffness and damping coefficients K xx , K xy , K yy and D yx increase. The increasing extent of K yy is bigger than K xx because of the gas film mainly supports the weight of high-speed rotor along the vertical direction, while K yx , D xx , D xy and D yy decrease quickly. With gas rarefaction effects taken into account, the direct stiffness coefficients K xx and K yy are linearly proportional to the perturbation frequency, but the direct terms of the damping coefficients D xx and D yy decrease slightly with the rise of frequency. Meanwhile, the cross coupled stiffness coefficients K xy and K yx gradually approach to the constant values with increasing Ω, yet the D xy and D yx are hardly affected by perturbation frequency. The dynamic stiffness coefficients for rarefaction models are less than those of continuum flow. When Ω > 2.0, the direct damping coefficients with the effects of gas rarefaction are higher than those of continuum flow. The reason might be that the energy transport by gas molecule collision is less, to some extent, which gives a larger value of fluid film stiffness and a larger gas film damping coefficient. The first-order slip model overestimates the stiffness coefficients of micro gas journal bearing compared with the simulation results of Boltzmann model, while the second-order velocity slip model underestimates it. Because the FK model is obtained based on the linearized Boltzmann equation, its results are in good agreement with those of the Boltzmann model. Appl. Sci. 2019, 9, Figures 11 and 12. It can be observed that the dynamic stiffness coefficients rise markedly with increasing eccentricity ratio. The reason is that the larger eccentricity ratio ε yields the smaller gas film thickness H of gas-lubricated microbearings, which leads to the increase in the Poiseuille flow rate. For small Knudsen number flows Kn < 0.325, the damping coefficients Dxy, Dyx and Dyy decrease while Dxx remains nearly constant as the eccentricity ratio increases. The increase in Knudsen number, which indicates that the rarefaction effect tends to be stronger. With the exception of Dyx, the dynamic coefficients gradually increase when Kn > 0.325. So the increase of direct damping coefficients result in a more stable state of the microbearing system.

K xx
Eccentricity ratio   To investigate the effect of eccentricity ratio on dynamical coefficients of Boltzmann model for Ω = 2, Λ = 3 and B/D = 0.5, comparisons of stiffness and damping coefficients at different values of Knudsen number from K n = 0.065 to K n = 2.1667 are shown in Figures 11 and 12. It can be observed that the dynamic stiffness coefficients rise markedly with increasing eccentricity ratio. The reason is that the larger eccentricity ratio ε yields the smaller gas film thickness H of gas-lubricated microbearings, which leads to the increase in the Poiseuille flow rate. For small Knudsen number flows K n < 0.325, the damping coefficients D xy , D yx and D yy decrease while D xx remains nearly constant as the eccentricity ratio increases. The increase in Knudsen number, which indicates that the rarefaction effect tends to be stronger. With the exception of D yx , the dynamic coefficients gradually increase when K n > 0.325. So the increase of direct damping coefficients result in a more stable state of the microbearing system.   Figures 11 and 12. It can be observed that the dynamic stiffness coefficients rise markedly with increasing eccentricity ratio. The reason is that the larger eccentricity ratio ε yields the smaller gas film thickness H of gas-lubricated microbearings, which leads to the increase in the Poiseuille flow rate. For small Knudsen number flows Kn < 0.325, the damping coefficients Dxy, Dyx and Dyy decrease while Dxx remains nearly constant as the eccentricity ratio increases. The increase in Knudsen number, which indicates that the rarefaction effect tends to be stronger. With the exception of Dyx, the dynamic coefficients gradually increase when Kn > 0.325. So the increase of direct damping coefficients result in a more stable state of the microbearing system.

K xx
Eccentricity ratio 

D xx
Eccentricity ratio    Figures 13 and 14 show the relation between the dynamic coefficients and the aspect ratio B/D for various Knudsen numbers with Ω = 3, ε = 0.6 and Λ = 5. It is found that all the stiffness coefficients increase significantly with the increase of B/D. The variations of damping coefficients is more complicated. The reason for this is the ultra-thin gas film is stiffened at higher aspect ratio and can reduce the dominant energy dissipation in such narrow spacing. For Kn < 0.325, the dimensionless damping coefficients Dxx, Dxy and Dyy first increase with increasing aspect ratio, reach a maximum and then begin to fall, while Dyx first decreases and then increases with increasing B/D. As the gaseous rarefaction effect is enhanced, the stiffness coefficients show a slightly more gradual changes at Kn =

D xx
Eccentricity ratio    Figures 13 and 14 show the relation between the dynamic coefficients and the aspect ratio B/D for various Knudsen numbers with Ω = 3, ε = 0.6 and Λ = 5. It is found that all the stiffness coefficients increase significantly with the increase of B/D. The variations of damping coefficients is more complicated. The reason for this is the ultra-thin gas film is stiffened at higher aspect ratio and can reduce the dominant energy dissipation in such narrow spacing. For Kn < 0.325, the dimensionless damping coefficients Dxx, Dxy and Dyy first increase with increasing aspect ratio, reach a maximum and then begin to fall, while Dyx first decreases and then increases with increasing B/D. As the gaseous rarefaction effect is enhanced, the stiffness coefficients show a slightly more gradual changes at Kn =  Figures 13 and 14 show the relation between the dynamic coefficients and the aspect ratio B/D for various Knudsen numbers with Ω = 3, ε = 0.6 and Λ = 5. It is found that all the stiffness coefficients increase significantly with the increase of B/D. The variations of damping coefficients is more complicated. The reason for this is the ultra-thin gas film is stiffened at higher aspect ratio and can reduce the dominant energy dissipation in such narrow spacing. For K n < 0.325, the dimensionless damping coefficients D xx , D xy and D yy first increase with increasing aspect ratio, reach a maximum and then begin to fall, while D yx first decreases and then increases with increasing B/D. As the gaseous rarefaction effect is enhanced, the stiffness coefficients show a slightly more gradual changes at K n = 2.1667 which means they are insensitive to the aspect ratio. The direct damping coefficients increase linearly with the aspect ratio up to B/D = 1 and D yx decreases marginally.

Conclusions
To study the lubrication mechanisms of rarefied gas flows in microbearings, the modified Reynolds equations incorporating different rarefaction models are derived and then solved by the partial derivative method and relaxed iterative algorithm. The effects of Knudsen number and bearing parameters on the static and dynamic characteristics of micro gas-lubricated bearings are investigated in detail.
The effect of gas rarefaction greatly decreases the load capacity and friction coefficient in comparison with the continuum flow case. In small Knudsen number regions, the maximum film pressure of FK model is lowest and the numerical result of Boltzmann model is largest. As the Knudsen number increases further, the first-order slip model overpredicts the microbearing performance while the second-order slip model underestimates it in the transition regime. The direct stiffness coefficients increase while the direct damping coefficients decrease with the increase of perturbation frequency for different rarefaction models. The stiffness coefficients of the Boltzmann model increase dramatically with increasing eccentricity ratio and aspect ratio. When Kn < 0.325, the damping coefficients Dxy, Dyx and Dyy decrease while Dxx has little change as the eccentricity ratio increases. Furthermore, Dxx, Dxy and Dyy first increase and then decrease with the growth of aspect ratio, while the reverse tendency occurs in Dyx. Thus, better understanding the gas rarefaction effect on ultra-thin film lubrication of microbearing is of practical importance for MEMS devices.

Conclusions
To study the lubrication mechanisms of rarefied gas flows in microbearings, the modified Reynolds equations incorporating different rarefaction models are derived and then solved by the partial derivative method and relaxed iterative algorithm. The effects of Knudsen number and bearing parameters on the static and dynamic characteristics of micro gas-lubricated bearings are investigated in detail.
The effect of gas rarefaction greatly decreases the load capacity and friction coefficient in comparison with the continuum flow case. In small Knudsen number regions, the maximum film pressure of FK model is lowest and the numerical result of Boltzmann model is largest. As the Knudsen number increases further, the first-order slip model overpredicts the microbearing performance while the second-order slip model underestimates it in the transition regime. The direct stiffness coefficients increase while the direct damping coefficients decrease with the increase of perturbation frequency for different rarefaction models. The stiffness coefficients of the Boltzmann model increase dramatically with increasing eccentricity ratio and aspect ratio. When K n < 0.325, the damping coefficients D xy , D yx and D yy decrease while D xx has little change as the eccentricity ratio increases. Furthermore, D xx , D xy and D yy first increase and then decrease with the growth of aspect ratio, while the reverse tendency occurs in D yx . Thus, better understanding the gas rarefaction effect on ultra-thin film lubrication of microbearing is of practical importance for MEMS devices.