A new extended Reynolds equation for gas bearing lubrication based on the method of moments

An extended Reynolds equation based on the regularised 13 moment equations and lubrication theory is derived for gas slider bearings operating in the early to upper transition regime. The new formulation performs well beyond the capability of the conventional Reynolds equation modified with simple velocity-slip models. Both load capacity and pressure distribution can be reliably predicted by the extended Reynolds equation and are in good agreement with available direct simulation Monte Carlo data. In addition, the equations are able to provide both velocity and stress information, which is not conveniently recovered from available kinetic models. Tests indicate that the equations can provide accurate data for Knudsen numbers up to unity but, as expected, begin to deteriorate afterwards.


Introduction
In modern computer hard disk drives, the gas between the read-write head and the surface of a spinning disk forms an air slider bearing that supports the head floating above the disk. Traditionally, the classical Reynolds lubrication 1 3 23 Page 2 of 12 2010). Modification of the slip boundary condition alone is not sufficient to compensate for the substantial deviation from the continuum assumption. The Reynolds lubrication equation with a first-order velocity-slip boundary condition will over-predict the pressure rise (Burgdorfer 1959), while a second-order boundary condition consistently under-predicts the pressure rise (Hsia and Domoto 1983). To account for the non-equilibrium effects from the wall and the Knudsen layer, Fukui and Kaneko (1988) derived a generalised lubrication equation (FK model) for arbitrary Knudsen number based on the linearised Boltzmann equation with the Bhatnagar-Gross-Krook (BGK) collision model (Bhatnagar et al. 1954). Their derivation focused on the Poiseuille flow rate and no analytical expressions for the velocity profile and shear stress were given, although they can be recovered from the solution of the pressure equation. Due to the complexity of the original expression of the FK model, a piecewise curve fit formula was provided for its practical use (Fukui and Kaneko 1990). However, Cercignani et al. (2004) pointed out that the database given by Fukui and Kaneko (1990) is inaccurate. Furthermore, as noted by Chen and Bogy (2010), there is a contact pressure singularity in the FK model. Cercignani et al. (2007) derived a Reynolds equation on the basis of the BGK and ellipsoidal statistical BGK Boltzmann equation. Its solution of pressure requires an accurate evaluation of the Abramowitz functions (Abramowitz and Stegun 1970), and their results show that the ellipsoidal statistical BGK model slightly under-predicts the pressure rise compared to the BGK model. In both kinetic approaches, an accurate solution can be obtained from the linearised Boltzmann equation, but each requires significant effort and diligent programming.
To study the physics of gas slider bearings in the transition regime, the direct simulation Monte Carlo (DSMC) approach has generally been employed to simulate the flow details (Alexander et al. 1994;Liu and Ng 2001;Jiang et al. 2005;John and Damodaran 2009). Gas velocity profiles between the disk and slider, the shear stress on the wall, as well as the pressure distribution, can be obtained by averaging particle movement statistically. Data obtained using DSMC provide an accurate and reliable way to validate our lubrication theory in the transition regime. However, the method is computationally expensive and is generally inappropriate for engineering design and analysis. A lubrication equation that can capture Knudsen layer effects offers great potential for design engineers. In the present study, we have used the method of moments to derive an extended Reynolds equation that can be used in the transition regime and can provide design information beyond the FK pressure distribution and without substantial effort.
The method of moments (Grad 1949) provides an approximation to the Boltzmann equation. In this approach, the Boltzmann equation is satisfied in a certain average sense rather at the molecular level. Macroscopic quantities, such as pressure, p, gas velocity, u i , stress, σ ij , and heat flux, q i , are moments of the molecular distribution function. As the constitutive relationships for σ ij and q i in the Navier-Stokes-Fourier (NSF) equations are no longer valid in the transition regime, governing equations for σ ij and q i are used, which can be derived from the Boltzmann equation. However, the fluxes of σ ij and q i appear in a set of governing equations, which are necessary to close the system of equations. In the original approach of the moment method, the fluxes of σ ij and q i are calculated from a truncated distribution function of Hermite polynomials. The resultant set of equations is referred to as Grad's 13 moment equations (G13) (Grad 1949). Applying a Chapman-Enskoglike expansion, a set of regularised 13 moment equations (R13) was developed (Struchtrup and Torrilhon 2003), which significantly improves the capability of the G13 equations, particularly in confined geometries. The regularisation procedure can be readily extended to incorporate higher moments such as the R26 equations . Recently, a velocity-slip boundary condition constructed from the R13 solution of a half-space problem was adopted to modify the traditional NS-based Reynolds equation (Yang, Zhang and Liu 2015). In the present work, the extended Reynolds equation is derived from kinetic theory and the R13 equations, rather than the NS equations.
In small-scale devices, such as micro-electro-mechanical systems (MEMS), Reynolds numbers are usually very small. Consequently, a linearised set of moment equations is generally adequate to capture the necessary flow features for many practical purposes. Macroscopic quantities in the Knudsen layer can be expressed as superpositions of exponentials of different widths obtained from the linearised moment equations (Struchtrup 2008;Gu et al. 2010;Gu and Emerson 2014). The linearised R26 equations have been shown to provide a more accurate description of the Knudsen layer than the R13 equations. However, the linearised R13 equations are adopted in the present study to derive an extended Reynolds equation because of their relative simplicity. As the pressure variation in the bearing system is significant, it is inappropriate to linearise the pressure-related terms for thin-film gas flow between two plates. To maintain the accuracy of the R13 equations, the pressure-related nonlinear terms need to be retained. In this way, the R13 equations are semi-linearised. Indeed, we will show that this new system of equations achieves a good level of accuracy and is suitable for engineering design and analysis.
In the next section, the air bearing problem is specified and the semi-linearised R13 equations are presented. We will then apply lubrication theory to the moment equations to reduce the complexity of the equation set. Using this approximation, the solutions derived from the extended Reynolds equations are presented and are compared with available DSMC data (Alexander et al. 1994;Liu and Ng 2001;Ng et al. 2002).

Problem specification and mathematical formulation
A schematic of a three-dimensional slider bearing geometry sitting in a Cartesian coordinate system, x i = x, y, z , is shown in Fig. 1. The bottom plate, located at z = 0, moves with a constant velocity, u w , in the x-direction and has a length, L, and a width, W. The upper plate is fixed and has a constant pitch angle, θ . The gap between the plates in the z-direction, H, varies and has a minimum value of H 0 at the end of the bearing. The gap is much smaller than the length and width of the plates so that ε = H 0 L ≪ 1, H 0 W ≪ 1 , with sin θ ≃ 0 and cos θ ≃ 1. The gas flow between the plates is considered to be isothermal and at a constant temperature, T. Outside the domain of the slider bearing, the pressure is at atmospheric conditions,p 0 , and the mean free path of gas molecules, , can be estimated by where µ is the gas viscosity and R the gas constant, respectively. Based on the minimum gap height, H 0 , the Knudsen number at the exit, Kn, is defined by It is further assumed that the gas obeys the ideal gas law, p = ρRT . The gas flow inside the slider bearing obeys the mass conservation law: where ρ is the density of the gas and u i = [u, v, w] is the gas velocity. The momentum equation is required to describe the gas movement. In terms of stress, σ ij , and pressure, p, the momentum equation reads Instead of using a constitutive relationship for velocity and stress to close Eqs. (3) and (4), the governing equations for the stress, σ ij , and the heat flux, q i , are derived from the moments of the Boltzmann equation to capture any nonequilibrium effects. The semi-linearised equations of σ ij and q i for low-speed isothermal flows, in terms of primitive variables, are: and in which the linear term of the higher moments m ijk , R ij , and for Maxwell molecules (MM) is approximated by (Struchtrup 2005(Struchtrup , 2008: and Equations (3)-(9) constitute a semi-linearised set of regularised 13 moment equations where the underlined terms in Eqs. (5) and (6) indicate the pressure-related nonlinear components.
Recently, Gupta and Torrilhon (2012) calculated the hard sphere (HS) collision coefficients for the R13 moment equation. The only difference between HS and MM collision models for the R13 moment equation is the approximation for the higher moment, R ij . Instead of Eq. (8), the expression for R ij for the HS model is: Comparison of Eqs. (8) and (10) shows that the value of the gradient coefficient is reduced from 4.8, for MM, to 3.8 for the HS model. The remaining coefficients for the other moments of MM and HS models are identical. As a consequence, both MM and HS models produce results that are

Fig. 1
Schematic of a three-dimensional slider geometry similar, which will be shown in Sect. 4. For this paper, the coefficients of the MM model are employed to derive the Reynolds equation for the slider bearing.

Lubrication approximation and extended Reynolds equation
Lubrication theory is an approximate description of flows at low Reynolds numbers through narrow geometries with slow changes in curvature. Under the condition that ε = H 0 L ≪ 1, lubrication theory scales the variables in the following way: Inserting Eq. (11) into Eqs. (4)-(6) and neglecting inertial terms in the momentum equations under lubrication theory, the resulting momentum equations for a three-dimensional thin gas flow are expressed in non-dimensional form as: and The relationship, σ xx +σ yy +σ zz = 0, has been used in Eq. (14). The equations for σ ij and q i in scaled non-dimensional form in Cartesian coordinates are written as: By using lubrication theory and neglecting all terms of O ε 2 and higher, the momentum equations can be simplified significantly. In particular, Eq. (14) indicates that the pressure gradient in the z-direction can be approximated as zero so that the pressure can be regarded as a two-dimensional variable. For convenience and ready comparison with available kinetic data, it is written as: Without the terms of O ε 2 and higher, Eqs. (12)-(18) can be grouped into two independent sets: one for flow in the x-direction, the other in the y-direction. For flow in the x-direction, they are: Similarly, the equations for flow in the y-direction can be written as, and in which Equation (19) indicates that P is independent of z, and the underlined nonlinear terms in Eqs. (21), (22), (26) and (27) reduce to linear terms. There were no explicit wall boundary conditions available for the moment equations, which hamper the use of the moment method in a wide range of practical applications, although the principle for constructing wall boundary conditions was suggested by Grad (1949). The first set of explicit wall boundary conditions (Gu and Emerson 2007) for the regularised moment equations is constructed from Maxwell's kinetic boundary condition (Maxwell 1879). The boundary conditions are further improved for the R13 and R26 moment equations . For Eqs. (20)-(29), the linearised boundary conditions for a solid wall with normal vector, n i , and tangential vector, τ i , are expressed by ): (33) u n = 0, with where α is the accommodation coefficient, which indicates that a fraction (1 − α) of gas molecules will undergo "specular" reflection, while the remaining fraction α will be "diffusely" reflected with a Maxwellian distribution function at the wall temperature (Maxwell 1879). Here, ū τ =ū i τ i and ū n =ū i n i , are the normalised tangential velocity slip at the wall and the velocity normal towards the wall, respectively, and σ nτ =σ ij n i τ j , q τ =q i τ i , m nnτ =m ijk n i n j τ k and R nτ =R ij n i τ j are the tangential components of σ ij , q i , m ijk , R ij relative to the wall.
Using the ideal gas law, p = ρRT , the mass continuity Eq. (3) is rewritten as: Integrating Eq. (52) using the no-mass penetration condition given by Eq. (33) leads to a depth-averaged continuity equation: (43) (45) q y = C 1y e −az + C 2y e az + 3 2 If we substitute Eqs. (37) and (47) into Eq. (53), we obtain the governing equation for pressure within the slider bearing as: in which and The underlined terms in Eq. (54) correspond to the Navier-Stokes solution with a first-order velocity-slip boundary, and the non-equilibrium behaviour introduced through the Knudsen layer is represented by F P and F C in Eqs. (55) and (56). It should be pointed that Eqs. (55) and (56) are merely algebraic expressions. In the Reynolds equation based on the Boltzmann equation (Fukui and Kaneko 1988;Cercignani et al. 2007), a significant number of terms are expressed in Abramowitz functions, which need to be evaluated accurately. Equation (54) is a two-dimensional, time-dependent second-order partial differential equation, and it can be used to study the dynamic response of the gas between two walls. The numerical and computational effort in solving the extended Reynolds equations is negligible in comparison with the DSMC method. The DSMC approach needs to simulate the whole geometry of the slider bearing, while the extended Reynolds equation is one dimension less than the geometry. Furthermore, the DSMC method requires a sufficient number of particles to travel and collide for a certain time to ensure statistical accuracy. For example, it takes hours to simulate the DSMC validation cases referenced (Liu and Ng 2001) in the next section, but it takes less than one second for the extended Reynolds equation to obtain the solution. In practical time-dependent engineering applications, the computational advantage of the present approach over the DSMC is very significant.
In the next section, this equation will be validated against available DSMC data for slider bearing configurations. However, its application is beyond that. For example, it can readily be applied to study squeeze film air damping in MEMS (Bao and Yang 2007).

Results and discussions
Knowledge of the pressure distribution within the slider bearing is essential for bearing design. The DSMC data of Alexander, Garcia and Alder (1994) and Liu and Ng (2001) are used to evaluate the accuracy of the extended Reynolds Eq. (54) obtained from the R13 moment equations. In their DSMC studies, a steady two-dimensional slider bearing configuration was considered. For the simulations, the length of the channel was L = 5μm and the minimum channel height was fixed at, H 0 = 50 nm, which results in an aspect ratio ε = H 0 L = 0.01. The angle of inclination, θ , varied from 0.002 to 0.016 rad, so that the entrance channel height, H 1 , ranges from 60 to 130 nm. The working fluid in the DSMC study was argon gas with viscosity, µ = 2.08 × 10 −5 Pa s and ambient conditions set to be T = 273 K and p 0 = 1 atm (101,325 Pa). This gives the Knudsen number, Kn, a value of 1.24. The walls are treated as fully diffusive, i.e. α = 1. In the simulation of Alexander et al. (1994), hard sphere particles with a diameter of 3.66 × 10 −10 m were used, while in the study of Liu and Ng (2001), a variable hard sphere collision model with a diameter of 4.17 × 10 −10 m was employed.
The extended Reynolds Eq. (54) for two-dimensional steady-state flows can be reduced to: Equation (57) can be solved numerically by the finite difference method, and we used 100 grid points uniformly across the slider bearing length. It takes less than a second using a personal computer. The computed pressure distribution obtained from Eq. (57) with two different molecular collision models is shown in Fig. 2 and compared with the DSMC data (Alexander et al. 1994;Liu and Ng 2001) for the case of u w = 25 m/ s, a bearing number, Λ = 61.6, and pitch angle, θ = 0.01 rad. As the Maxwell molecules and hard sphere model have similar collision constants, their results are close to each other. In addition, the pressure distributions for the MM and HS models are in good agreement with the DSMC data, although they are closer to the data of Alexander et al. (1994) than the results of Liu and Ng (2001). The discrepancy between the two sets of DSMC (57) data is mainly due to different particle size and collision models employed. Based on Fig. 2, only the results using the extended Reynolds equation with the MM model will be presented in the following sections. The improvement of the extended Reynolds equation over the conventional Reynolds equation derived with a velocity-slip boundary condition is presented in Fig. 3. It is clear that the first-order slip boundary condition overpredicts the pressure rise, while the second-order boundary condition under-predicts the pressure change. The discrepancies between the conventional Reynolds equation with slip models and the DSMC data are so significant that they are unacceptable for engineering analysis. The extended Reynolds equation based on the R13 equations is also in good agreement with the FK model derived from the linearised Boltzmann equation (Fukui and Kaneko 1988;1990), and both results are in close agreement with the data of Alexander et al. (1994).

Effect of pitch angle θ
During its operation, the slider floats above the disk with an averaged pitching angle. The dynamic response of the pressure distribution to the variation of θ can be readily computed from Eq. (54). However, it is not easy to perform DSMC simulations with dynamic variation of θ . Instead, a series of steady-state flow cases with different pitching angles were simulated with DSMC (Liu and Ng 2001). The pressure distributions predicted by Eq. (57) in the slider bearing with Kn = 1.24, Λ = 61.6 and θ = 0.002, 0.006, 0.01 and 0.016 rad, respectively, are presented in Fig. 4 and compared with the DSMC data (Liu and Ng 2001). At a low pitching angle, θ = 0.002 and 0.006 rad, the extended Reynolds equation reproduces the DSMC data fairly accurately. When the pitching angle is increased to θ = 0.016 rad, Eq. (57) underestimates the peak pressure by about 1.5 %. As we continue to increase the value of θ , the assumptions that underpin lubrication theory start to break down and the accuracy of the predictions using the extended Reynolds equation begin to deteriorate. The pitching angle, θ , influences the load capacity of the slider bearing, which is defined by the bearing force, W , as The predicted bearing force, W , from the solution of the extended Reynolds equation, plotted against pitch angle, θ , is shown in Fig. 5 and compared with available DSMC data (Liu and Ng 2001). Agreement is good, and both results indicate that as the value of θ increases, the value of W will increase, which is consistent with the trends in Fig. 4. The load centre, x c , is the focal point of the resultant pressure on the slider surface, and it can be calculated from the pressure distribution by The position of the load centre moves towards the narrow end of the slider bearing as the pitch angle increases, as indicated in Fig. 5. However, the change of the centre position is only slight.

Effect of bearing number Λ and Knudsen number Kn
It can be seen from Eqs. (54) or (57) that the bearing number, Λ, is one of the important parameters and can influence the pressure distribution significantly. For a fixed value of Kn = 1.24 and θ = 0.01 rad, the DSMC data of three different bearing numbers are available to validate the extended Reynolds equation. Liu and Ng (2001) simulated the operating conditions of low and medium bearing numbers, Λ = 61.6, 123.2, and 221.8, respectively. The predicted pressure distributions from the extended Reynolds equation agree with the DSMC data quite well. For a high bearing number, Λ = 758, the predicted pressure distribution from the macroscopic equation is in excellent agreement with the DSMC data of Alexander, Garcia and Alder  1994). In this case, the corresponding Mach number is 1.0, yet the extended Reynolds equation still performs well. The pressure distributions at different values of Kn for a slider bearing with Λ = 61.6 and θ = 0.01 rad are shown in Fig. 7. At a transitional Knudsen number of 0.5, the bearing generates more than a 55 % pressure rise. The pressure profile predicted by the extended Reynolds equation is close to the FK model (1988,1990) with a slight over-prediction of the peak pressure. As the Knudsen number increases, with the bearing number fixed, the disk velocity is reduced so that the magnitude of the pressure rise reduces. At Kn = 1.24, the pressure rise is reduced by about 30 %. The FK model (1988, 1990, DSMC data (Alexander et al. 1994), and Eq. (57) are all in good agreement. As the Knudsen number increases further, the extended Reynolds equation gradually starts to under-predict the pressure rise predicted by the FK model. This is due to the fact that a set of 13 moments is not sufficient to predict correctly the flow rate for Poiseuille flow at such large Knudsen numbers . It should be noted that the pressure in the extended Reynolds equation is an averaged quantity across the gap. The good agreement of pressure between the DSMC data and the extended Reynolds equation at a Knudsen number around unity does not imply the R13 equations can accurately predict the flow field but they will provide qualitatively reliable information.
Comparing Figs. 6 and 7 indicates that Λ and Kn have opposite effects on the pressure rise, which will affect the load capacity of the slider bearing. For a slider bearing with θ = 0.01 rad and ε = 0.01, the bearing force diagram constructed from the solution of Eq. (57) is shown in Fig. 8 for different values of Kn and Λ. For a fixed value of Λ, the bearing force, W, decreases as the Knudsen number increases. The lower the bearing number, the more rapidly the bearing force drops as Kn increases. For a fixed Knudsen number, the bearing force increases as the bearing number increases.

Predicted velocity profiles and velocity at the wall
Using the extended Reynolds equation, the velocity profile between the moving disk and the fixed slider can be calculated from Eq. (37). It is worth noting that the tangential heat flux, q x , is involved in Eq. (37), even if the conditions are assumed to be isothermal. This non-equilibrium effect is clearly captured by the exponentials in Eq. (35) and contributes to the characteristic Knudsen layer description. The predicted velocity profiles from Eq. (37) at x = 0, 0.5, and 1 are shown in Fig. 9 for a slider bearing with ε = 0.01, Kn = 1.24 and Λ = 61.6. In general, the macroscopic predictions qualitatively follow the DSMC data (Liu and Ng 2001) at the three locations, although they cannot capture the Knudsen layer accurately. As indicated in a previous study of Kramers problem (Gu et al. , 2010, higherorder moment equations, such the R26 equations , are necessary to describe the Knudsen layer effect accurately. Nevertheless, the present extended Reynolds equation and its associated velocity expression are useful in bearing analysis and design. The velocity predicted for the walls indicate a higher level of slip on both the static upper wall and the moving lower wall, as shown in Fig. 9. Figure 10  a slider bearing with ε = 0.01, θ = 0.01 rad and Kn = 1.24, in comparison with DSMC data (Liu and Ng 2001). It can be seen that the velocity slip is over-predicted for all three cases, particularly for the case with a low bearing number, Λ = 61.6. The over-prediction of the slip velocity is related to the poor description of the Knudsen layer close to the wall. As the bearing number increases, the prediction of the slip velocity gradually begins to improve.

Shear stress on the moving lower wall
The shear stress on the wall can be as important as the pressure distribution in the design of a slider bearing. In the present approach, the shear stress can be calculated from Eq. (36) once the solution of the extended Reynolds equation has been obtained. An extensive DSMC study of the shear stress on the lower wall was carried out by Ng et al. (2002). For a slider bearing with Kn = 1.24 and θ = 0.01 rad, three aspect ratios ε = 0.005, 0.01 and 0.02 were set up with three wall velocities, u w = 25, 50 and 100 m/s, respectively, to achieve a bearing number, Λ = 123.2. The shear stress for the three cases is shown in Fig. 11, with the lines representing the solution from Eq. (36) and symbols the DSMC data. The agreement between the present study, using Eq. (36), and the DSMC data is good with the extended Reynolds equation predicting slightly higher values of stress. From the values of shear stress, the drag force on the walls can be calculated. Using lubrication theory, we have derived an extended Reynolds equation based on the regularised 13 moment equations. This new set of equations is able to capture all of the important non-equilibrium effects associated with gaseous transport between narrow gaps and can readily be used in the design of new slider bearings. Results obtained with this new set of equations have been compared with available direct simulation Monte Carlo data and are shown to be significantly better than current Navier-Stokes models that incorporate velocity slip. Moreover, the equations are also able to predict shear stress distributions and velocity fields, although the latter indicates that resolution of the Knudsen layer is restricted with the 13 moment equations and a higher-order set of equations is necessary if a more detailed analysis of the Knudsen layer is required. In general, however, this new set of equations performs very well, is straightforward to implement for robust engineering design, and will provide a useful advance in modelling gaseous transport in small gaps.

Appendix
The terms in the integration constants of Eqs. (40)- (43) and (48)    Lines: stress Eq. (36). Symbols: DSMC data (Ng et al. 2002) the financial support from the National Natural Science Foundation of China (Grant No. 11102071) and China Scholarship Council (Grant No.201208330295) and would like to thank STFC Daresbury Laboratory and Professor Y.H. Zhang in University of Strathclyde to host his academic visit in UK.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.