A Multicoefficient Slip-Corrected Reynolds Equation for Micro-Thin Film Gas Lubrication

This work investigates and analyzes the performance of conventional slip models among various regimes of Knudsen number and developes a new multicoefficient slip-velocity model, by using Taguchi quality control techniques and numerical analysis. A modified Reynolds equation is also derived based on the new slip-flow model. The multicoefficient slip model and its slip-corrected Reynolds equation are suitable to a wide Knudsen range from slip to transition regime. In comparison with other conventional slip models, it is found that the current results have a better agreement with the solution obtained from the linearized Boltzmann equation and direct simulation of Monte Carlo method (DSMC).


INTRODUCTION
The Reynolds equation for gas film lubrication has been widely accepted for the description and analytical investigations of the microflow in magnetic storage system.However, in the modern computer magnetic storage device, the increasing demand for smaller, more compact disk storage requires minimum film thickness.When the spacing between flying head and the rotating disk approximates 10 3 of molecular mean-free path or less, the gas dynamics cannot be described from the continuum transport theory directly [1,2], and the gaseous rarefaction effects must be taken into account.Recently, the minimum spacing under the slider is within the order of 50 nm or less.This rarefaction effect can be incorporated by introducing phenomenological slip flow at the solid surfaces to the continuum momentum transport equation.
The rarefied flow is typically classified into four flow regimes according to their Knudsen number, defined as the ratio of molecular mean-free path λ to characteristic length, or the minimum spacing h, in a gas film lubrication prob-This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.lem [2,3].They are continuum flow (Kn ≤ 10 −2 ), slip flow (10 −2 < Kn < 0.1), transition flow (0.1 ≤ Kn < 10), and molecular flow (Kn ≥ 10).Commonly, the continuum flow can be described by the Navier-Stokes equation and the free molecular flow is described by the collisionless Boltzmann equation.
Traditionally, in slip-flow or transition-flow regimes, the flow is predicted by the Reynolds equation with slip corrections (slip-corrected Reynolds equation), or by a molecular generalized lubrication equation [4].
Slip flow may give a limited correction in slip-flow regime by using conventional slip-flow models.Such as the wellknown first-order slip model of Maxwell [5], and the firstorder-corrected Reynolds equation of Burgdorfer [6]; the second-order model of Hsia and Domoto [7], and the 1.5order slip of Mitsuya [8]; the Beskok and Karniadakis' model [9].On the contrary, in transition-flow regime, especially when the Knudsen number approaches or is larger than one, the gas rarefaction effect becomes very sensitive to the slipflow model.The conventional models are limited to a narrow Kn range, and would result in some deviations over the transition-flow regime [8,9].Therefore, to develop a universal slip model suitable to both slip and transition regime is meaningful for the study and development of modern computer storage system [10].
The geometry of gas film lubrication.
In the current study, the authors devote themselves to develop an improved slip model suitable to both slip and transition regimes by means of the variable analysis using Taguchi quality control method [11].Based on this novel model, a slip-corrected Reynolds equation is finally derived and solved numerically.The results confirm that the multicoefficient slip model and the slip-corrected Reynolds equation perform better than the previous models for ultrathin film gas lubrication problem.

THEORETICAL ANALYSIS OF SLIP MODEL
From the momentum transfer between the gas and the plates as shown in Figure 1, the well-known slip-flow models are given by (i) first order: (ii) second order: (iii) 1.5th order: where u s is slip velocity; σ is an accommodation coefficient and y is in the bottom wall-normal direction.
It can be noticed that only the first order is derived using the physical momentum transfer model.Therefore, we use only one parameter A for the part of first order, but using two variables B and C in the other part to make it more flexible.The 1.5th-order slip model presented in (3) lacks physical sense.It is included here solely for comparison purpose.Strictly speaking, the order of the slip boundary conditions should be integers, and they should correspond to the order of the Chapman-Enskog expansion of the Boltzmann equations.For example, the first-order expansion in Kn results in the Navier-Stokes equations, while the second-order expansion gives the Burnett equations, and so forth.The slip boundary condition (slip velocity) can be written accordingly as a universal form: For ease in analysis, the accommodation coefficient is taken as 1.According to (1), (2), and (3), the parameters A, B, and C can be obtained as shown in Table 1.
Using this slip-flow model, velocity boundary conditions for the lubricating film in Figure 1 have the form of Here, h is the height of channel in stream-wise location, and h 0 is the height at the channel trailing edge (slider exit); L is the length of channel; p 0 is ambient pressure; U is bottom wall velocity.
The velocity distribution is obtained by solving the momentum equation without the inertia term subject to the above boundary conditions (5), and this avoid difficulty in computing the second-derivative term Here, µ is absolute viscosity of gas.The Poiseuille flow rate with the universal slip model can also be derived as Here, D is inverse Knudsen number.Compared to the conventional slip models, there are three parameters in the universal model, which can be adjusted to give suitable values after investigating their effects.

PARAMETRIC STUDY USING TAGUCHI METHOD
Taguchi method [11] is a very efficient tool for analyzing and developing high quality products at a low cost.Instead of the real experiment, which could be simulated by using various analytical tools or numerical techniques, the Taguchi method  is implemented using a series of combinations from the different factors and levels involved.The grouping of these iterations is called an orthogonal array.Each parameter is referred to as a factor and each different value for each factor is termed as a level.
In the present parametric study, we solve (7) by using three-factor and three-level orthogonal arrays L 27 (3 3 ) to identify the factor/parameter that has the most influence on the Poiseuille flow rate.This procedure will make the following numerical analysis easier and more accurate.We can then determine the value of parameters A, B, and C better and faster.The three factors are the parameters A, B, and C. The three levels are shown in Table 2, which are chosen to be around their values from the conventional models (Table 1).
Before applying the Taguchi method, we need to set several monitor points so that to produce the calculated results.
Here the three points are selected for the wider possible range of flow rate and defined as shown in Figure 2. In the first case (case 1), three monitoring points are set at Kn = 0.8, Kn = 8.0, and Kn = 80.The relative parameters of the three points (monitor locations of ∆ on mass flow rate) P ∆1 , P ∆2 , and P ∆3 are included in Table 3.The associated three values ∆1, ∆2, and ∆3 are the differences between the current result and the Boltzmann result of mass flow rate, respectively.The averages of the values ∆1, ∆2, and ∆3 are expressed by ∆1, ∆2, and ∆3 for the three-level parameters A, B, and C, respectively, as shown in Table 4.The difference between the biggest value and the smallest value for each parameter is expressed as R(∆), which provides the tabulated results of the effects of parameters upon the mass flow rate.The higher the R(∆) value, the greater the effect of the parameter has on the difference of mass flow rate.For the first location, parameter A has the highest effect on the difference of mass flow rate ∆1(Kn = 0.8) because there is a largest value of R(∆1) = 0.35449076 for A, whereas B has a little smaller (R(∆1) = 0.30342839) effect on ∆1, and C with R(∆1) = 0.02707428 has an even smaller effect on ∆1.It is also seen that B has the largest effect on ∆2(Kn = 8.0), and C has the biggest effect on ∆3(Kn = 80.0), as tabulated in Table 4.
In the second test example (case 2), the parametric study is investigated around the critical Kn = 1.0.The three monitor points are included in Table 5.The results in Table 6 indicate that the parameter A has the highest effect on the difference of mass flow rate ∆1(Kn = 0.1), and both parameters A and B have the same effect on ∆2(Kn = 1.0), whilst parameter C has negligible effect as compared to A and B. Finally, parameterB has the highest effect on ∆3(Kn = 10.0).
In test case 3, the parametric study is investigated around the critical Kn = 4.5.The three monitor points are included in Table 7.The results in Table 8 at monitor location P ∆1 , where Kn = 0.45, shows that factor A had the largest influence on the mass flow rate as R(∆1) was 0.354490757.Factor B would then affect the mass flow rate more than factor C, as their R(∆1) values were 0.204426944 and 0.065018550, respectively.At monitor location P ∆2 , (Kn = 4.5), factor B would affect the mass flow rate most followed by factors C and A. At monitor location P ∆3 , (Kn = 45.0),factor C had the largest influence on the mass flow rate followed by factors B and C.

DETERMINATION OF MODEL COEFFICIENTS
Using the Taguchi results in case 1, the parameter C is adjusted first because it has a significant effect on the mass flow rate.From the ∆3 values in Table 4, the smallest C(= −0.5) results in the smallest average value of the difference of mass flow rate between the current and Boltzmann results.Therefore, a series of results are produced using smaller C(≤ −0.5) as shown in Figure 3.The results with C = −0.6 and C = −0.7 have the similar trends with Boltzmann result.Therefore, in the analysis of other parameters, C is taken as −0.65.Also in case 1, for the result of ∆2, parameter B has a bigger effect (Table 4).The calculated results for different B are shown in Figure 4. Apparently, the result with B = 0.25 is preferred.
According to the investigated results of Taguchi method, parameter A has a dominating effect on ∆1 in case 1 (Kn = 0.8).After the determination of parameters B and C, the parameter A can be determined and then is adjusted in detail according to the results of case 2. Table 4 indicates that the larger A is preferable; therefore the calculation is processed with A ≥ 1.0.
From the results of case 2 in Table 6, when A is taken around 1.2, a smaller value of ∆1 (the average value of difference of mass flow rate between the current and Boltzmann results at Kn = 0.8) is obtained.Thus an investigation with value of A around 1.2 is performed in detail in the regime of 0.1 ≤ Kn ≤ 10.The results in Figure 5 suggest that A should be taken between 1.1 and 1.2.For the slip model used in the regime around Kn = 1, a value of A = 1.15 is found reasonable.Thereafter, a new slip model is constructed.It can be expressed as This new slip model is a multicoefficient slip model with three adjustable coefficients and the corresponding velocity distribution is (1.15)λh + 2(0.25)Kn (−0.65) λ 2 + hy − y 2 ∂p ∂x

THE MODIFIED REYNOLDS EQUATION
By using the newly derived slip model ( 5) and the velocity distribution ( 6), as well as the nondimensional parameters and for an isothermal process, one can obtain twodimensional Reynolds equation for an ultrathin gas bearing lubrication using the newly proposed multicoefficient slip model as shown in the following.(The details of its derivation can refer to [16]).This equation is termed as multicoefficient slip-corrected Reynolds equation with here Λ is bearing number, Kn = Kn0p/P 0 , where subscript 0 shows the reference conditions at the slider exit, and This is the two-dimensional Reynolds equation.The Poiseuille flow rate would be   [8] is included here for better comparison.In a wide regime, where the inverse Knudsen number varies from D = 100 to D = 0.01, the first-order slip-flow approximation underestimates the flow rate from the Boltzmann equation, while second-order slip-flow and 1.5th-order slip-flow approximations overestimate it.Using the multicoefficient slip model, the flow rate produces a good approximation that is closer to the flow rate from the linear Boltzmann solution.Note that the current effort gets the coefficients of the arbitrary slip model by matching the flow rate.This of course would ensure the correct flow rate, but not necessarily the correct velocity profile.The boundary conditions should be used to solve for the velocity profile, which will in return integrate to the correct volumetric flow rate.Flow rate is an integral value and infinitely many parabolic velocity profiles with slip will give the desired flow rate.This is the basic inconsistency in most high-order slip models, which model the entire Knudsen regime with slip corrections.It would thus be interesting to compare the velocity profiles from the present model with the Boltzmann solutions of pressure-driven channel flow from Sone et al. [12].Also, in the 2D Boltzmann equation solutions, the volumetric flow rate increases logarithmically with Kn in the free molecular flow regime, which is physically not possible.It is because there are intermolecular collisions existing in the spanwise direction for degenerating of 2D geometry.This model, originally by Fukui and Kaneko [13], is usually evaluated by a lookup table.Other curve fit in common use is also known [14].Thus the current method reduces to a sophisticated interpolation of experimental data, where the fitting parameters are introduced in the boundary conditions.
The integral Reynolds equation ( 11) is solved using a fourth-order Runge-Kutta method.The pressure profiles obtained from the multicoefficient slip-corrected Reynolds equation and the previous Reynolds equation are shown in Figures 7 and 8.In Figure 7, the Kn number is 1.24 with the bearing number of 61.6, and the current result DSMC result of Alexander et al. [14] and Reynolds equation results with the first-order slip and the second-order slip models are included.In Figure 8, the pressure is calculated with the Kn number of 1.24 but the bearing number of 132.2, and the DSMC result is obtained from Ng et al. [15,16].It is obvious that the Reynolds equation result using multicoefficient slip model is in good agreement with the DSMC computation.The differences between the DSMC prediction and the conventional approximations using previous models are significant.The current results indicated that the slip-corrected Reynolds equation with the newly proposed multicoefficient slip model performs better both in slip and transition regimes than the previous models.It is thus warranting to further compare the current model with Beskok and Karniadakis' [9] new slip model for various applications and test cases of microflows.

CONCLUSION
The authors have developed an improved slip model, termed multicoefficient slip model, particularly suitable to the transition regime where the conventional slip models fail to describe the slippage phenomena accurately.With the current new slip model, a slip-corrected Reynolds equation is further derived.The predicted results of flow rate for plane Poiseuille flow, and the results of pressure distributions for slope flow are obtained using the new slip model and the modified Reynolds equation, and are benchmarked with those available in the literatures.The present effort yields better results than conventional models in a wider range of Knudsen numbers that may encounter in the modern computer magnetic storage application.

Figure 2 :
Figure 2: The definition of the flow rate differences between slip model and Boltzmann result.

Figure 6
Figure6compares the flow rates of Poiseuille case for the current and the previous slip-flow models.The flow rate from the Boltzmann equation[8] is included here for better comparison.In a wide regime, where the inverse Knudsen number varies from D = 100 to D = 0.01, the first-order slip-flow approximation underestimates the flow rate from the Boltzmann equation, while second-order slip-flow and 1.5th-order slip-flow approximations overestimate it.Using the multicoefficient slip model, the flow rate produces a good approximation that is closer to the flow rate from the linear Boltzmann solution.Note that the current effort gets the coefficients of the arbitrary slip model by matching the flow rate.This of course would ensure the correct flow rate, but not necessarily the correct velocity profile.The boundary conditions should be used to solve for the velocity profile, which will in return integrate to the correct volumetric flow rate.Flow rate is an integral value and infinitely many parabolic velocity profiles with slip will give the desired flow rate.This is the basic inconsistency in most high-order slip models, which model the entire Knudsen regime with slip corrections.It would thus be interesting to compare the velocity profiles from the present model with the Boltzmann solutions of pressure-driven channel flow from Sone et al.[12].Also, in the 2D Boltzmann equation solutions, the volumetric flow rate increases logarithmically with Kn in the free molecular flow regime, which is physically not possible.It is because there are intermolecular collisions existing in the spanwise direction for degenerating of 2D geometry.This model, originally by Fukui and Kaneko[13], is usually evaluated by a lookup table.Other curve fit in common use is also known[14].Thus the current method reduces

Table 1 :
The factors/parameters in three slip-flow models.

Table 2 :
The factors and levels used in Taguchi method.

Table 3 :
The locations of three monitor points for case 1.

Table 4 :
Results for effects of parameters on case 1.

Table 5 :
The locations of three monitor points for case 2.

Table 6 :
Results for effects of parameters on case 2.

Table 7 :
The locations of three monitor points for case 3.

Table 8 :
Results for effects of parameters on case 3.