Interactive Effects of Rarefaction and Surface Roughness on Aerodynamic Lubrication of Microbearings

The aerodynamic lubrication performance of gas microbearing has a particularly critical impact on the stability of the bearing-rotor system in micromachines. Based on the Duwensee’s slip correction model and the fractal geometry theory, the interactive effects of gas rarefaction and surface roughness on the static and dynamic characteristics were investigated under various operation conditions and structure parameters. The modified Reynolds equation, which governs the gas film pressure distribution in rough bearing, is solved by employing the partial derivative method. The results show that high values of the eccentricity ratio and bearing number tend to increase the principal stiffness coefficients significantly, and the fractal roughness surface considerably affects the ultra-thin film damping characteristics compared to smooth surface bearing.


Introduction
The microfluidic devices are widely used in many applications such as ultra-precision machine tool spindles, inertial navigation system (INS), medical devices, and hard disk drives (HDDs) in micro-electro-mechanical systems (MEMS). Since the development of microsystems engineering technology, gas journal microbearing has been generally preferred over the electromagnetic bearing and rolling bearing owing to its advantages of simple structure, high rotary accuracy, high running speed, low friction power loss, and wide working temperature range [1][2][3]. The real surface of a mechanical part produced by various machining and finish operations is composed of a large number of distributed peaks and valleys. The lubricating film thickness between the surfaces of shaft and bearing has continually decreased, which results in the increase of the roughness heights that are of the same order of magnitude as the minimum clearance gap. Thus, the assumptions that the gas flow is typically treated as a continuum flow with no slip boundary condition and the bearing surface roughness is considered negligible in classical fluid mechanics have to be revised at microscales [4][5][6][7]. The interaction between rarefaction and surface roughness in microbearing can influence the reliability and operational efficiency of micro rotating machinery obviously, so the gas journal microbearing performance should be comprehensively analyzed.
In order to accurately predict the effects of surface roughness and rarefaction on the bearing characteristics, many researchers devoted numerous research efforts to study the complicated flow behaviors at very small clearances over the last few decades. In consideration of rarefaction effects in ultra-thin gas film lubrication, the Knudsen number K n is utilized to describe the rarefied gas flow, which is defined as the ratio of molecular mean free path λ 0 to the characteristic length of gas film thickness. Burgdorfer [8], Hsia et al. [9], and Mitsuya [10], based on the slip velocity boundary condition, derived the classical first order, second order, and 1.5 order slip models in the slider/disk interface for HDDs to take into account the effect of gas rarefaction. Fukui and Kaneko [11,12] developed a Poiseuille flow rate database for a wide Knudsen number range to modify the compressible Reynolds type equation including thermal creep flow and accommodation coefficient from linearized Boltzmann equation. In order to account for the effect of surface roughness, Christensen and Tønder [13,14] presented a stochastic model of hydrodynamic lubrication for finite width journal bearing in which they considered the lubricant film thickness as a stochastic process. The operating characteristics of bearing was theoretically analyzed with roughness pattern, nominal geometric features, and statistical properties by surface averaging techniques. Via linear transformation of random matrices, the Gaussian or non-Gaussian distribution of surface heights were generated by Patir [15] using the prescribed autocorrelation functions and frequency density functions. Patir and Cheng [16] further derived the average Reynolds equation suitable for various roughness structure and discussed the effect of roughness on mean hydrodynamic pressure, mean viscous friction, and mean bearing inflow in finite slider bearings. The average flow model of Patir and Cheng was extended by Tripp [17], in which the statistical expectation of flow factors were calculated with a perturbation expansion of the film pressure. The results showed that the flow factors are closely correlated with roughness parameters. White et al. [18] introduced the transverse sinusoidal roughness pattern to study the influence of surface roughness on steady-state pressure profiles of wedge bearing by variable grid implicit finite difference method and found that the load capacity could be decreased to a limiting value at higher bearing numbers. For the applications of perturbation technique and mapping function, Li et al. [19] studied the effects of roughness orientations and rarefaction on static performance of magnetic recording systems. The results demonstrated that the flow factors changed with the orientation angle and Peklenik number, and the effect of moving surface on surface characteristics is more significant than that of the stationary surface. Turaga et al. [20] proposed the stochastic finite method to solve Reynolds equation and obtained the static and dynamic performance of hydrodynamic journal bearings with the longitudinal, transverse and isotropic roughness pattern. Naduvinamani et al. [21] established the surface roughness by a stochastic random variable with nonzero mean, variance and skewness, and the average Reynolds equations were adopted to analyze the performance of porous step-slider bearings with Stokes couple stress fluid. Zhang et al. [22,23] presented the modified Reynolds equation by including fractal roughness effect and velocity slip boundary condition and concluded that the flow behaviors in gas-lubricated journal microbearings was appreciably affected by Knudsen number, bearing number and fractal dimension. The coupling effects of non-Newtonian micropolar fluids and roughness on the dynamic characteristics of plane slider bearings were investigated by Lin et al. [24] on the basis of the microcontinuum theory and Christensen stochastic roughness model. They indicated that the transverse roughness serves to somewhat increase bearing dynamic property, whereas the longitudinal roughness would tend to decrease the dynamic coefficients. Jao et al. [25] examined the influences of surface roughness and anisotropic slips on hydrodynamic lubrication of journal bearings. They described the lubricant flow in rough bearing surface by the product of flow factors and flow in nominal film thickness, and also identified that boundary slip reduced the effect of surface roughness. Kalavathi et al. [26] reported a generalized Reynolds equation for finite porous slider bearing with both longitudinal and transverse roughness. The authors showed the surface roughness enhanced the pressure distribution and load carrying capacity while the permeability parameter diluted the load. Quiñonez [27] utilized the linear superposition of perturbation method and Flourier transformation to provide a solution for the flow characteristics of wide exponential land slider bearings with rough surfaces. The results were in good agreement with the cases of sinusoidal and single Gaussian dent. The linear perturbation method was used by Wang et al. [28] to solve the unsteady Reynolds equation for rough aerostatic journal bearings during the iterative process, and the dynamic performance was obtained by taking into account the interactions of journal rotation and surface waviness. However, likely due to the nonlinear and complexity of dynamic flow behavior, previous papers were mainly focusing on steady-state characteristics in rough journal bearings, and the dynamic characteristics of hydrodynamic gas-lubricated microbearings were seldom reported in the research literature. Moreover, the statistical parameters (such as root mean square of asperity heights, surface slope, curvature, skewness, and kurtosis), which are conventionally applied to characterize surface roughness, vary with the sampling length and resolution of measuring equipment. A scale-invariant surface characterization should be considered. Hence, the analytical studies of surface roughness effect on dynamic characteristics of gas slider bearings with rarefaction coefficients in microfluidic engineering devices is motivated. In this paper, the Weierstrass-Mandelbrot (W-M) fractal function is used to characterize the homogeneous surface roughness, and the Boltzmann slip correction model is applied to represent the rheological behavior of compressible rarefied gas film. The generalized Reynolds-type equation considering gas rarefaction, as well as roughness effect, is mathematically derived and solved by the partial derivative method and relaxation iteration algorithm. Bearing performances (including the load-carrying capacity, friction coefficient and corresponding attitude angle, dynamic stiffness, and damping coefficients) are presented and discussed in comparison with smooth surface bearings. The work is expected to elucidate the performance characteristics of gas microbearings with Poiseuille flow and random asperities, which is conducive to understand the fluid mechanisms of very low clearance gas films for microfluidics devices.

Characterization of Fractal Rough Surface
The distribution of asperity heights on the bearing surface consisting of long narrow ridges and valleys in engineering practice. A rough surface is a random system and the fractal geometry is introduced to characterize the random and multiscale topographies. Mandelbrot [29] initially developed the fractal theory by researching the coastal geomorphology in 1967. He found that the most machined surfaces can be constructed using Weierstrass-Mandelbrot (W-M) function with randomness, multiscale nature, self-similarity, and self-affine property. Unlike the traditional characterization parameters of surface roughness, fractal roughness parameters are independent of scan lengths and provide all the surface topography information of rough profiles. The W-M function that explicitly expressed the homogeneous rough surface in self-acting gas-lubricated journal bearings is given by where h r (x,y) is the height of rough surface, x and y are the measure distances in the vertical and horizontal position, respectively. L is the sampling length of the profile of surface. D f is the fractal dimension, varying from 2 to 3 in three-dimensional surface topography. G is the scaling constant that relates to the roughness profile. γ is the scaling parameter (γ > 1), which determines the spectral density, γ is equal to 1.5 for a Gaussian and isotropic surface. M is the number of overlapped ridges on the surface, m and n are the frequency index, n max = int[log(L/L s )/logγ], ϕ m,n is the random phase, L s is the cut-off length that depends on cut-off wavelength of resolution in measuring machines. The values of asperity heights can be changed by the fractal dimension D f , comparisons of the distributions of asperity heights for different D f are illustrated in Figure 1. It is seen that the heights of rough surfaces increase as the fractal dimension decreases. In the three-dimensional isotropic roughness type, the schematic presentation of rough gas microbearing is plotted in Figure 2.

Numerical Model and Solution Method
In the current analysis, rarefied gas flow between the surfaces of bearing and journal is treated as isothermal, laminar flow with uniform viscosity. and elastic deformation of bearing surface is not considered. The Boltzmann slip correction is used to model the gas rarefaction effect at the gas-solid interface for arbitrary Knudsen numbers. The physical configuration of micro gas bearing with rough surface is shown in Figure 3. The modified Reynolds equation incorporating gas rarefaction and surface roughness effects in non-dimensional form appears as In the three-dimensional isotropic roughness type, the schematic presentation of rough gas microbearing is plotted in Figure 2. In the three-dimensional isotropic roughness type, the schematic presentation of rough gas microbearing is plotted in Figure 2.

Numerical Model and Solution Method
In the current analysis, rarefied gas flow between the surfaces of bearing and journal is treated as isothermal, laminar flow with uniform viscosity. and elastic deformation of bearing surface is not considered. The Boltzmann slip correction is used to model the gas rarefaction effect at the gas-solid interface for arbitrary Knudsen numbers. The physical configuration of micro gas bearing with rough surface is shown in Figure 3. The modified Reynolds equation incorporating gas rarefaction and surface roughness effects in non-dimensional form appears as

Numerical Model and Solution Method
In the current analysis, rarefied gas flow between the surfaces of bearing and journal is treated as isothermal, laminar flow with uniform viscosity. and elastic deformation of bearing surface is not considered. The Boltzmann slip correction is used to model the gas rarefaction effect at the gas-solid interface for arbitrary Knudsen numbers. The physical configuration of micro gas bearing with rough surface is shown in Figure 3. The modified Reynolds equation incorporating gas rarefaction and surface roughness effects in non-dimensional form appears as where P = p/p a , H = h/c b , ϕ = x/R, λ = z/R are the dimensionless gas film pressure, the dimensionless gas film thickness, and the coordinates in the circumferential and axial direction, p a is the ambient pressure, c b is the radius clearance, R is the radius of journal, p is the local gas pressure, h is the clearance spacing of ultra-thin gas film, ε is the eccentricity ratio and ε = e/c b , e is the eccentricity. Λ = 6µωR 2 /(p a c b 2 ) is the gas bearing number, µ is the viscosity coefficient, ω is the rotating angular velocity of journal, T is the dimensionless time.  In order to obtain the aerodynamic performance of gas microbearings, the modified Reynolds equation for homogeneous surface roughness and gas rarefaction should be solved numerically. Equation (2), governing the gas pressure distribution of the bearing, is a nonlinear two-dimensional partial differential equation (PDE). It is difficult to get its analytical solution. So, the partial derivative method [32,33], with a relaxation iteration algorithm, is employed to ensure a reasonable and efficient solution. By introducing the mathematical transformation PH = S, (PH) 2 = S 2 = Π, Equation (2) can be converted to the ellipse-type partial differential equation ( ) c u au f      in the following form: The pressure boundary conditions for the Reynolds equation are: where B is the bearing width.
The non-dimensional hydrodynamic gas film forces are obtained by integrating the film For the slip correction factors in rarefied gas flows under ultra-low spacing, the dimensionless Poiseuille flow rate between rotor and bearing surfaces with the gas thickness h is given as [30] The Poiseuille flow rate ratio Q appears as where the inverse Knudsen number D = √ π 2K n , the three adjustable coefficients a = 0.01807, b = 1.35355 and c = −1.17468. The Boltzmann Poiseuille flow rate ratio [31] Q is expressed as Here the lubricant film thickness h is the sum of nominal smooth film thickness h 0 and random roughness h r measured from the nominal smooth height.
As in the steady state, the transient term ∂(PH)/∂T of Equation (2) can be ignored, the static dimensionless modified Reynolds equation can be gained as In order to obtain the aerodynamic performance of gas microbearings, the modified Reynolds equation for homogeneous surface roughness and gas rarefaction should be solved numerically. Equation (2), governing the gas pressure distribution of the bearing, is a nonlinear two-dimensional partial differential equation (PDE). It is difficult to get its analytical solution. So, the partial derivative method [32,33], with a relaxation iteration algorithm, is employed to ensure a reasonable and efficient solution. By introducing the mathematical transformation PH = S, (PH) 2 = S 2 = Π, Equation (2) can be converted to the ellipse-type partial differential equation −∇ · (c∇u) + au = f in the following form: The pressure boundary conditions for the Reynolds equation are: where B is the bearing width. The non-dimensional hydrodynamic gas film forces are obtained by integrating the film pressure acting on the microbearing along the both horizontal and vertical directions.
The attitude angle θ of the gas-lubricated journal microbearing is calculated by The total non-dimensional load-carrying capacity is written as The non-dimensional skin friction coefficient on the journal surface can be computed by Suppose that the journal center whirls around its static equilibrium position with a small amplitude periodic motion under the perturbation frequency ratio Ω, the linear perturbation method is adopted for calculating the dynamic stiffness and damping coefficients. The steady-state position is indicated as (ε 0 ,θ 0 ), and its dynamic disturbance about (ε 0 ,θ 0 ) are denoted as E and Θ. The positions of eccentricity ratio ε and attitude angle θ of journal at a random position are represented by the static and dynamic components as follows: where E 0 and Θ 0 are the small perturbation amplitude of journal eccentricity ratio and attitude angle in the complex field. The dimensionless perturbation frequency Ω, which is defined as the ratio of journal disturbance frequency ν to rotating angular velocity ω of journal, i = √ −1. Thus, the non-dimensional gas-film pressure and gas-film thickness of the gas microbearing can be expressed as where , P 0 is the static gas-film pressure and H 0 is the static gas-film thickness. Q gd , H gd are the dynamic gas film pressure and gas film thickness, respectively. P 0 and H 0 are the perturbation magnitudes in terms of complex numbers for dynamic gas film pressure and gas film thickness. Substituting Equation (15) into Equation (2), the generalized dynamic lubrication equation for molecular model and surface roughness can be derived as: As mentioned in Reference [34], some variables are defined by After differentiating P 0 and H 0 in Equation (16) with respect to E 0 and Θ 0 , and combining with some mathematical transformation, the resulting dynamic PDE equations are obtained for rough surface microbearings concerning the variables P E , P θ , H E , and H θ .
According to the coordinate system at the bearing midplane as illustrated in Figure 3, by simultaneously solving the nonlinear Equations (18)-(21) using the partial derivative method with iteration procedure, the dynamic stiffness coefficients K ij and dynamic damping coefficients D ij of gas journal microbearing for fractal rough surface can be calculated by the following formula: The dynamic coefficients in the Cartesian coordinate system are given by the transformation matrix A.
It was confirmed that the dynamic characteristics of the gas journal bearing can greatly affect the critical speed, unbalance response, and instability threshold of a hydrodynamic gas-lubricated bearing-rotor system after long-term research and practice. The comprehensive analysis of dynamic coefficients is important in the design of gas microbearings in MEMS applications, the key findings presented in the next section will reveal some insights into the ultra-thin gas film lubrication problems.

Results and Discussion
On the basis of fractal geometry theory and the Boltzmann slip correction factor, the combined effects of gas rarefaction and surface roughness on the static and dynamic characteristics of ultra-thin film gas lubrication in journal microbearings are investigated in detail. The pressure distribution, load carrying capacity, friction coefficient and attitude angles of bearing, dynamic stiffness and damping properties are analyzed concerning the fractal dimension and bearing geometric parameters in this section, the primary design parameters shown in Figure 3 are R = 1mm, B = 200 µm, c = 1 µm, p a = 1.033 × 10 5 N/m 2 , and the aspect ratio is B/D = 0.1.

Steady-State Film Pressure
To verify the correctness of the developed theoretical model and program code employed in this paper, the non-dimensional pressure profile for the middle cross-section along the sliding direction obtained by the current solution is compared with Zhang et al. in Reference [22]. As indicated in Figure 4, the simulation results are in close agreement with the numerical predictions reported by Zhang et al. for ε = 0.7, B/D = 0.075, c = 12 µm, ω = 5 × 10 5 rpm. It can also be seen that the pressure randomly fluctuates when the effect of surface roughness is considered. The influence of bearing number Λ, eccentricity ratio ε, and surface roughness on pressure distributions P is shown in Figure 5. An important observation exhibited by Figure 4a,b is that the maximum gas film pressure becomes larger as Λ increases for fixed fractal parameters G and D. This phenomenon can be explained by the enhanced aerodynamic effect, which indicates that higher angular velocity of the journal dilutes the pressure diffusion. In Figure 4a,c, increasing the eccentricity ratio ε leads to the higher pressure distribution and magnitude of pressure fluctuations decrease. In comparison with the smooth surface case, the roughness effect increases the pressure profile for the rough surface at the same bearing number and eccentricity ratio, and the random surface roughness makes the pressure distribution across the entire lubricating film unpredictable. The contour plots show more clearly the detail and variation of pressure over the rough bearing surfaces. At greater bearing number Λ the influence of surface roughness decreases and the contour lines approach that of the smooth case, the fluctuations of the pressure contour plots are significant with the increase of fractal roughness in gas bearing. The influence of bearing number Λ, eccentricity ratio ε, and surface roughness on pressure distributions P is shown in Figure 5. An important observation exhibited by Figure 4a,b is that the maximum gas film pressure becomes larger as Λ increases for fixed fractal parameters G and D. This phenomenon can be explained by the enhanced aerodynamic effect, which indicates that higher angular velocity of the journal dilutes the pressure diffusion. In Figure 4a,c, increasing the eccentricity ratio ε leads to the higher pressure distribution and magnitude of pressure fluctuations decrease. In comparison with the smooth surface case, the roughness effect increases the pressure profile for the rough surface at the same bearing number and eccentricity ratio, and the random surface roughness makes the pressure distribution across the entire lubricating film unpredictable. The contour plots show more clearly the detail and variation of pressure over the rough bearing surfaces. At greater bearing number Λ the influence of surface roughness decreases and the contour lines approach that of the smooth case, the fluctuations of the pressure contour plots are significant with the increase of fractal roughness in gas bearing.  The influence of bearing number Λ, eccentricity ratio ε, and surface roughness on pressure distributions P is shown in Figure 5. An important observation exhibited by Figure 4a,b is that the maximum gas film pressure becomes larger as Λ increases for fixed fractal parameters G and D. This phenomenon can be explained by the enhanced aerodynamic effect, which indicates that higher angular velocity of the journal dilutes the pressure diffusion. In Figure 4a,c, increasing the eccentricity ratio ε leads to the higher pressure distribution and magnitude of pressure fluctuations decrease. In comparison with the smooth surface case, the roughness effect increases the pressure profile for the rough surface at the same bearing number and eccentricity ratio, and the random surface roughness makes the pressure distribution across the entire lubricating film unpredictable. The contour plots show more clearly the detail and variation of pressure over the rough bearing surfaces. At greater bearing number Λ the influence of surface roughness decreases and the contour lines approach that of the smooth case, the fluctuations of the pressure contour plots are significant with the increase of fractal roughness in gas bearing.   As the eccentricity ratio increases, which denotes the film thickness is thinner, the load capacity increases monotonically. The higher self-affine fractal dimensions yield the smaller roughness heights distribution on rough surfaces, and it is noted that the load carrying capacity is increased gradually when compared with the smooth bearing. The reason is that the increasing values of roughness heights reduces the sidewise leakage of airflow and the flow is restricted by the surface asperities. However, as the fractal dimension decreases further for Df = 2.2, the load carrying capacity tends to decline under this condition. This is because the minimum air film clearance between rotor and bearing may become too small so that the surface roughness effect which increases the load-carrying capacity is weaker than the gaseous rarefaction effects which reduces the dimensionless load capacity, thus causing a decrease in the bearing load capacity. The attitude angle θ is found to decrease with the growing eccentricity ratio. The decrease in attitude angle is more accentuated for a rougher surface as compared to a nominally smooth surface.  As the eccentricity ratio increases, which denotes the film thickness is thinner, the load capacity increases monotonically. The higher self-affine fractal dimensions yield the smaller roughness heights distribution on rough surfaces, and it is noted that the load carrying capacity is increased gradually when compared with the smooth bearing. The reason is that the increasing values of roughness heights reduces the sidewise leakage of airflow and the flow is restricted by the surface asperities. However, as the fractal dimension decreases further for D f = 2.2, the load carrying capacity tends to decline under this condition. This is because the minimum air film clearance between rotor and bearing may become too small so that the surface roughness effect which increases the load-carrying capacity is weaker than the gaseous rarefaction effects which reduces the dimensionless load capacity, thus causing a decrease in the bearing load capacity. The attitude angle θ is found to decrease with the growing eccentricity ratio. The decrease in attitude angle is more accentuated for a rougher surface as compared to a nominally smooth surface.

Load-Carry Capacity and Friction Coefficient
The variation of friction coefficient with the eccentricity ratio for different values of fractal dimension with fixed values of Λ = 20 and G = 1 × 10 −11 is depicted in Figure 8. It can be seen that friction coefficients monotonically increase as ε increases. Although the contact area between the rarefied gas flow and the surface asperities is larger when the homogeneous surface roughness becomes more and more obvious, the static friction coefficients show a slightly more gradual increase with fractal dimension. As the fractal dimension D f is equal to 2.2, the friction coefficient is lower than that of smooth surface for the same reason that the gaseous rarefaction effect is more pronounced and thus friction coefficient drops.  The variation of friction coefficient with the eccentricity ratio for different values of fractal dimension with fixed values of Λ = 20 and G = 1 × 10 −11 is depicted in Figure 8. It can be seen that friction coefficients monotonically increase as ε increases. Although the contact area between the rarefied gas flow and the surface asperities is larger when the homogeneous surface roughness becomes more and more obvious, the static friction coefficients show a slightly more gradual increase with fractal dimension. As the fractal dimension Df is equal to 2.2, the friction coefficient is lower than that of smooth surface for the same reason that the gaseous rarefaction effect is more pronounced and thus friction coefficient drops.  The variation of friction coefficient with the eccentricity ratio for different values of fractal dimension with fixed values of Λ = 20 and G = 1 × 10 −11 is depicted in Figure 8. It can be seen that friction coefficients monotonically increase as ε increases. Although the contact area between the rarefied gas flow and the surface asperities is larger when the homogeneous surface roughness becomes more and more obvious, the static friction coefficients show a slightly more gradual increase with fractal dimension. As the fractal dimension Df is equal to 2.2, the friction coefficient is lower than that of smooth surface for the same reason that the gaseous rarefaction effect is more pronounced and thus friction coefficient drops.    9 and 10 describe the effect of increasing bearing number on the load carrying capacity and attitude angle with various fractal dimensions for ε = 0.6 and G = 1 × 10 −10 . It is found that increasing the values of the bearing number from Λ = 3 up to Λ = 100 increases the carrying capacity and reduces the corresponding attitude angle of gas journal microbearing. The surface roughness effect in aerodynamic lubrication enhances the load capacity as compared with the smooth-bearing case, especially for the bearing operating at high bearing number. Figure 11 shows the comparison of the static friction coefficients with bearing number for different fractal dimensions. The friction coefficients exhibit a near-linear increasing trend with increasing Λ, while the increase extent in static friction coefficient is even higher at smaller D f values. Consequently, the strengthened gas-lubricated hydrodynamic effect and the bearing surface with roughness undulations have a significant influence on the skin friction at the interface.  Figure 9 and 10 describe the effect of increasing bearing number on the load carrying capacity and attitude angle with various fractal dimensions for ε = 0.6 and G = 1 × 10 −10 . It is found that increasing the values of the bearing number from Λ = 3 up to Λ = 100 increases the carrying capacity and reduces the corresponding attitude angle of gas journal microbearing. The surface roughness effect in aerodynamic lubrication enhances the load capacity as compared with the smooth-bearing case, especially for the bearing operating at high bearing number. Figure 11 shows the comparison of the static friction coefficients with bearing number for different fractal dimensions. The friction coefficients exhibit a near-linear increasing trend with increasing Λ, while the increase extent in static friction coefficient is even higher at smaller Df values. Consequently, the strengthened gas-lubricated hydrodynamic effect and the bearing surface with roughness undulations have a significant influence on the skin friction at the interface.

Dynamic Stiffness and Damping Coefficients
The variation of dynamic stiffness and damping coefficients with dimensionless perturbation frequency Ω for different values of fractal dimensions is plotted in Figure 12 and 13. The principal stiffness coefficients Kxx and Kyy increases as the perturbation frequency increases and the Kyy is much greater than Kxx because of the rarefied gas lubricating film supports the weight of journal in the vertical direction. The cross-couple stiffness Kxy increases at first, then decreases slightly with the growth of Ω, while Kyx decreases quickly at lower perturbation frequencies. Furthermore, the enhanced dynamic stiffness coefficients are seen for the rough bearing surface as compared to that of smooth bearing case. When Ω > 2, all the dynamic stiffness coefficients of micro gas-lubricated journal bearing at the fractal dimension Df = 2.3 are obviously greater than other rough surfaces. This is mainly due to the fact that the larger asperity heights lead to an increased Poiseuille flow component along the sliding direction and the side flow suffers the constriction resistance caused by homogeneous surface roughness. The influence of perturbation frequency on dynamic damping coefficients for different fractal roughness parameters can be observed from Figure 13. The principal

Dynamic Stiffness and Damping Coefficients
The variation of dynamic stiffness and damping coefficients with dimensionless perturbation frequency Ω for different values of fractal dimensions is plotted in Figures 12 and 13. The principal stiffness coefficients K xx and K yy increases as the perturbation frequency increases and the K yy is much greater than K xx because of the rarefied gas lubricating film supports the weight of journal in the vertical direction. The cross-couple stiffness K xy increases at first, then decreases slightly with the growth of Ω, while K yx decreases quickly at lower perturbation frequencies. Furthermore, the enhanced dynamic stiffness coefficients are seen for the rough bearing surface as compared to that of smooth bearing case. When Ω > 2, all the dynamic stiffness coefficients of micro gas-lubricated journal bearing at the fractal dimension D f = 2.3 are obviously greater than other rough surfaces. This is mainly due to the fact that the larger asperity heights lead to an increased Poiseuille flow component along the sliding direction and the side flow suffers the constriction resistance caused by homogeneous surface roughness. The influence of perturbation frequency on dynamic damping coefficients for different fractal roughness parameters can be observed from Figure 13. The principal damping coefficient D xx first increases with increasing dimensionless perturbation frequency, reaches a maximum, then starts to decline slowly. The absolute values of the cross-coupling terms of damping coefficients D xy and D yx decrease quickly at low Ω, then approaches to zero. The result show that the damping coefficients D xx , D yx increase with the increase in the isotropic and homogeneous roughness heights of gas slider bearing surface, and the D xy and D yy are first decreases and then dramatically increases as the fractal dimension D f decreases, whereas the difference in damping coefficients of rough and smooth cases appear to converge at higher values of Ω. Figures 14 and 15 display the relationship between the dynamic coefficients of gas microbearing and eccentricity ratio with various fractal dimensions. It is found that the dynamic stiffness coefficients increase as the eccentricity ratio increases for fixed values of Λ = 20 and G = 1 × 10 −11 . The higher ε corresponds to the thinner gas film thickness, which results in the increase of Knudsen number K n . With the growth of the fractal surface roughness, the stiffness coefficients first increase gradually, then decrease significantly at the same ε. As illustrated in Figure 15, the damping coefficients increase marginally as the eccentricity ratio increases and the effect of fractal dimension on damping coefficient D xy is negligible at lower eccentricity ratios. The damping coefficients become more sensitive to surface roughness at higher eccentricity ratios about ε > 0.7. It can be seen that increasing the random roughness heights increase the effect of gas rarefaction in small spacing. reaches a maximum, then starts to decline slowly. The absolute values of the cross-coupling terms of damping coefficients Dxy and Dyx decrease quickly at low Ω, then approaches to zero. The result show that the damping coefficients Dxx, Dyx increase with the increase in the isotropic and homogeneous roughness heights of gas slider bearing surface, and the Dxy and Dyy are first decreases and then dramatically increases as the fractal dimension Df decreases, whereas the difference in damping coefficients of rough and smooth cases appear to converge at higher values of Ω.  reaches a maximum, then starts to decline slowly. The absolute values of the cross-coupling terms of damping coefficients Dxy and Dyx decrease quickly at low Ω, then approaches to zero. The result show that the damping coefficients Dxx, Dyx increase with the increase in the isotropic and homogeneous roughness heights of gas slider bearing surface, and the Dxy and Dyy are first decreases and then dramatically increases as the fractal dimension Df decreases, whereas the difference in damping coefficients of rough and smooth cases appear to converge at higher values of Ω.   Figure 14 and 15 display the relationship between the dynamic coefficients of gas microbearing and eccentricity ratio with various fractal dimensions. It is found that the dynamic stiffness coefficients increase as the eccentricity ratio increases for fixed values of Λ = 20 and G = 1 × 10 −11 . The higher ε corresponds to the thinner gas film thickness, which results in the increase of Knudsen number Kn. With the growth of the fractal surface roughness, the stiffness coefficients first increase gradually, then decrease significantly at the same ε. As illustrated in Figure 15, the damping coefficients increase marginally as the eccentricity ratio increases and the effect of fractal dimension on damping coefficient Dxy is negligible at lower eccentricity ratios. The damping coefficients become more sensitive to surface roughness at higher eccentricity ratios about ε > 0.7. It can be seen that increasing the random roughness heights increase the effect of gas rarefaction in small spacing.  Figure 14 and 15 display the relationship between the dynamic coefficients of gas microbearing and eccentricity ratio with various fractal dimensions. It is found that the dynamic stiffness coefficients increase as the eccentricity ratio increases for fixed values of Λ = 20 and G = 1 × 10 −11 . The higher ε corresponds to the thinner gas film thickness, which results in the increase of Knudsen number Kn. With the growth of the fractal surface roughness, the stiffness coefficients first increase gradually, then decrease significantly at the same ε. As illustrated in Figure 15, the damping coefficients increase marginally as the eccentricity ratio increases and the effect of fractal dimension on damping coefficient Dxy is negligible at lower eccentricity ratios. The damping coefficients become more sensitive to surface roughness at higher eccentricity ratios about ε > 0.7. It can be seen that increasing the random roughness heights increase the effect of gas rarefaction in small spacing.       Kyy is near proportionally dependent on Λ and dynamic stiffness coefficients Kxx, Kxy and Kyx increase gradually with increasing bearing number. It is also observed that the increase in the dynamic stiffness coefficients is more accentuated for the fractal roughness surface as compared to the smooth surface bearing with the enhanced aerodynamic effect in gas journal bearings. The damping coefficients exhibit similar trends to the fractal dimension, namely the principal damping Dxx, Dyy and cross-couple damping Dxy for rough bearing surface become larger than the ones in the smooth bearing, whereas the roughness effect is rather marginal in the case of the cross-couple damping Dyx in Figure 17c. The damping coefficients Dxx, Dyy and Dxy increase quickly at first and then decreased, while the damping coefficient Dyx decreases with increasing bearing number. Therefore, the dynamic stiffness and damping characteristics of gas microbearings with isotropic and homogeneous roughness show relatively high values since the aerodynamic effect of the bearing is enhanced. Increment of the bearing number means the larger operating conditions. The principal stiffness K yy is near proportionally dependent on Λ and dynamic stiffness coefficients K xx , K xy and K yx increase gradually with increasing bearing number. It is also observed that the increase in the dynamic stiffness coefficients is more accentuated for the fractal roughness surface as compared to the smooth surface bearing with the enhanced aerodynamic effect in gas journal bearings. The damping coefficients exhibit similar trends to the fractal dimension, namely the principal damping D xx , D yy and cross-couple damping D xy for rough bearing surface become larger than the ones in the smooth bearing, whereas the roughness effect is rather marginal in the case of the cross-couple damping D yx in Figure 17c. The damping coefficients D xx , D yy and D xy increase quickly at first and then decreased, while the damping coefficient D yx decreases with increasing bearing number. Therefore, the dynamic stiffness and damping characteristics of gas microbearings with isotropic and homogeneous roughness show relatively high values since the aerodynamic effect of the bearing is enhanced.

Conclusions
The coupled effects of gas rarefaction and surface roughness on the static and dynamic characteristics of gas microbearings are studied using fractal geometry theory and Boltzmann model for Poiseuille flow. Based on the results mentioned in previous section, the following conclusions have been drawn: 1. Surface roughness and gaseous rarefaction effects are of great importance to the lubrication performance of gas journal microbearings. At small asperity height distributions, the gas film pressure distributions increase for higher bearing numbers and eccentricity ratios, consequently leads to an increment of load-carrying capacity.
2. The attitude angle increases with the larger fractal dimensions compared with smooth surface, whereas the skin friction coefficient yields a reversed trend because of the long narrow ridges and furrows impose a series of constrictions on gas lubricant in the sliding direction.
3. The principal stiffness coefficients increase with the increase in dimensionless perturbation frequency and all the stiffness coefficients become larger and larger at higher eccentricity ratios and bearing numbers. The damping coefficients increase as the eccentricity ratio increases and the principal damping can be first magnified and then diminished with increasing bearing number. 4. The stiffness coefficients increase as the fractal dimension decreases at the same perturbation frequency and bearing number, while the stiffness coefficients first increase and then decrease for increasing eccentricity ratio. The damping coefficients increase with decreasing fractal dimension except for the D yx at large bearing number values. 5. It seems that the fractal rough surface does not merely increases the steady-state and dynamic performance of gas microbearings. When the degree of rarefaction effect is more pronounced, the static load capacity and dynamic coefficients decrease quickly for larger values of distribution of asperity heights.