An Alternative SimulationMethod for Calculation of Microgas Flows under Flying Head Sliders

The precise knowledge of the force and moment generated by the air squeezed under the read-write slider by the rotating disc is an engineering necessity in designing the air bearing surface slider. This paper reviews methods addressing the thin gas film bearings problem. It firstly reviews briefly the relatively well-known two methods of calculations of the microgas flows under flying head sliders, the generalized Reynolds equation, having given a number of useful results of slider design, and the DSMC method, which is precise and appropriate for the flow of complex configurations but is restricted to miniature (∼micrometer) size sliders. The main purpose of the paper is to introduce to the reader an alternative method, the information preservation (IP) method, for use in simulation of the flows under air bearing surfaces. Some recent results of IP simulation of slider flows published on conference proceedings are introduced here.


Introduction
The problem of thin gas film bearings in the gap between the flying head slider and the magnetic disc now has an increasing interest among scientists concerning computation of actual problems of micro gas flows.To calculate precisely the force and moment generated by the air squeezed by the rotating disc is essential in the design of the head slider.The typical slider length in disk drives is about 1 mm and the width is the same order of the length.The size of the clearances between the slider and the disc is much smaller and is constantly decreasing to increase the recording density.The flying height in the early stage of the disc recording head was of the order of 8∼10 micrometer [1] and 1 micrometer [2] while the drives areal data densities of the order 100 Gbit/in 2 and 1 Tbit/in 2 in consideration require sliders to have a flying height of 5∼10 nm and less (see e.g., [3]).The flow regime in the modern flying sliders is definitely beyond the slip flow regime, and application of the tools of computation in transitional flow regime is inevitable.
The thin film air flow between the slider air bearing surface and the disc is most appropriately described by the Reynolds equation which is a differential equation relating the pressure p, density ρ, platter velocity U, and the height h of the gap, firstly developed by Reynolds for continuum fluid [4].The derivation of it can be found in [5].Burgdorfer [6] introduced the velocity slip correction to the Reynolds equation.Fukui and Kaneko [7,8] developed a generalization of the equation suitable for the transitional regime.The work in [9] provides a simple and enlightening derivation showing that the equation is the expression of balance of the flow rate of the Poiseuille flow and the flow rate of the Couette flow, and it also shows how the equation can be extended to slip regime and transitional regime.The equation modified to include slip and transitional effects is still called Reynolds equation.It is essentially a mass conservation relation applied to the cross sections of the squeezed air flow and is obtained from the continuity equation by integrating it over the vertical direction with the employment of the momentum equation, and now the generalized Reynolds equation is in routine use to calculate the air bearing parameters in sliders with complicated air bearing surfaces (see e.g., [3,[10][11][12][13]).The direct simulation Monte-Carlo (DSMC) method [14] is appropriate for simulation of thin film flows in sliders.But the low information-to-noise ratio for low Mach number flows makes the computational process very time-consuming and only simulation results of miniature (micrometer, not authentic millimeter) sizes are available [15][16][17].The information preservation (IP) method [18,19] was proposed to overcome the difficulty of low information-to-noise ratio of DSMC method.It was successfully applied to many low-speed microflows [20][21][22][23][24][25][26] and the IP results for channel flows compared well with DSMC results, experimental data, and exact kinetic theory [27].The present paper firstly reviews briefly the successful results obtained by the Reynolds equation and the DSMC method, relatively well-known in the literature.But its main purpose is to introduce the IP method and some recent development of the method in calculating the thin film flow under the air bearing surface not familiarized to the readers and show its feasibility as an alternative tool in calculating the micro gas flows under flying head sliders.

The Reynolds Equation
The Reynolds equation was first used in the continuum regime, the derivation is enlightening and with its essential assumption clearly revealed, and its framework can be easily used in the slip flow and transitional flow cases; so we give a simplified version of derivation of it.For simplicity the twodimensional case is considered (see Figure 1), assuming that the head width D is larger than the height h o ; so the spanwise (z direction) motion can be neglected.Writing the continuity equation in the form ∂ρv ∂y and integrating it over y across the whole flow region yields The left-hand side of (3) vanishes, as there is no fluid flown into or out of the walls.Interchanging the integration and differentiation gives For thin film flow with the inertial terms neglected the steady momentum Navier-Stokes equation has the following form: Integration across the gap with the nonslip boundary conditions yields the solution of the streamwise velocity component u: Substituting ( 7) into (4) and accomplishing integration over y, the following equation is attained (we emphasis here that the assumption of ∂p/∂y = 0 is essentially necessary in the derivation): This is the general form of the Reynolds equation for the twodimensional case.By introducing X = x/L, H = h/h o , P = p/ p o , and the bearing number: Equation ( 8) for steady and two-dimensional case can be written in the normalized form [15]: The first term of ( 7) is the slip-less solution of the velocity in the Couette flow when the upper plate is stationary and the lower plate moves towards the right with velocity U (see [9] compare with its (5.63) for ζ = 0), and the second term is the slip-less solution of the velocity in the Poiseuille flow when the x axis is aligned along the lower plate (see [9], the second term of ( 7) can be obtained from (5.69) in [9] by a simple translation of the ordinate y).Equation (10) shows that the flow rate across any cross section is the sum of the flow rate of the Couette flow and the Poiseuille flow and this rate does not change from one cross section to another in steady flow.
In [9] (see (5.73) there) one sees that the flow rate of the Poiseuille flow with slip boundary condition surpasses that of the slip-less case by a factor: As for the Couette flow the flow rates have a specific feature, and in the slip and the transitional flow cases and the slipless case they are identical and have the following value independent of the Knudsen number owing to the symmetry of the flow (see Figure 2): From the flow rate expressions (11) and (12) for Poiseuille and Couette flows in the slip flow case one can conclude that in the slip flow regime the following Reynolds equation is obtained in place of (10): where Kn = λ/h is local Knudsen number.When the slip boundary conditions instead of the nonslip boundary condition (6) is employed in solving the momentum equation ( 5), and the resulted velocity profile is substituted into the mass conservation relation (4), one would arrive at the same slip corrected Reynolds equation ( 13) [6,15].Fukui and Kaneko [7] showed that the solution of the linearized Boltzmann equation for the thin film bearing problem can be decomposed into the solutions of the plane Couette flow and the plane Poiseuille flow.On this basis they derived the generalized Reynolds equation for the thin film air bearing problem by employing the flow rates of the fundamental Poiseuille and Couette flows solved by the linearized Boltzmann equation.This generalized Reynolds equation in the isothermal case can be written as [7], where Q P,TR (Kn) is the flow rate in transitional regime (normalized by the slip-less value Q P,C ) calculated from the linearized Boltzmann equation for Poiseuille flow.A tabled database of the calculated values of Q P,TR (Kn) for σ = 1, σ = 0.9, σ = 0.8, and σ = 0.7 is provided in [8], and a fitted formula approximation for diffuse reflection (σ = 1) by Robert is recorded in [15] (there the second term on the right-hand side is misprinted as 6A √ πKn): where A = 1.318889 and B = 0.387361.Alexander et al. [15] used the DSMC method to simulate the short head length air bearing problem (L = 5 μm, h o = 50 nm = 0.05 μm, U = 25 m/s, σ = 1) and found excellent agreement of the DSMC simulation with the generalized Reynolds equations ( 15) and ( 16).But at that time they described the latter as continuum hydrodynamic Reynolds equation corrected for slip (in fact only ( 13) is such an equation); now we have cognized that the generalized Reynolds equation is a global mass conservation equation, but at each cross section its flow is governed by exact kinetic theoretical equation which is appropriate for transitional regime.The comparison made in [15] for the cases (the ratio of the inlet to outlet heights is kept as 2 : 1) showed good agreement of the generalized Reynolds equation with the results of DSMC simulation; this just confirms that the generalized Reynolds equation is a mass conservation equation in form (although some framework of the N-S equation has been used) but in fact it balances the flow rates of Poiseuille flow and Couette flow calculated from the kinetic theory; so it has the kinetic theoretical merit and can be used to solve the air bearing problem in the entire transitional flow regime.By analogy with the above derivation of Reynolds equation for the two-dimensional case, it is a simple matter to derive the unsteady and three-dimensional Reynolds equation in the following form (cf., e.g., the 3D Reynolds equation shown in [3]): where T is time normalized by 1/ω, ω being angular velocity of the rotating disc, Z, and W is the dimensionless coordinate and velocity in the head width direction z.Fukui and Kaneko [7] used the generalized Reynolds equation and obtained lubrication characteristic results valid for large Knudsen numbers.By using Fukui-Kaneko's generalized Reynolds equations some useful results of the slider design were obtained: Hu et al. [10] investigated partial contact air bearing characteristics of tripad sliders.Hu et al. [11] investigated the air bearing dynamics of two configurations of authentic sized sub-ambient pressure sliders and found the way to ensure the reliability of the unloading performance of the types of sliders considered.Tagawa et al. [12] computed the slider's responses to nanotextured disc surfaces.The computation techniques are developed for solving the slider air bearing problem [3] and new first-and second-order slip models are introduced to the generalized Reynolds equations [13] by Wu and Bogy.
In deriving the Reynolds equation in the form of ( 17) we have seen that it is essential to use the assumption ∂p/∂y = 0 (see the note in the brackets before ( 8)).So the Reynolds equation is applicable undoubtedly for cases, where the assumption ∂p/∂y = 0 is valid, for example, for sliders with smooth thin gas layers between the slider and the disc.
In applying the generalized Reynolds equation to the modern authentic sliders having rails on the peripheries (sometimes with thin terraces on them) and on the center part of the slider and fully recessed regions (see, e.g., the NSIC (National Storage Industry Consortium) slider cited in [3]) one should have caution.The depth of the air bearing surface varies drastically at the vertical wall profile regions.This makes the pressure change drastically along the vertical direction and causes perceivable gas flow along the vertical wall direction.The condition ∂p/∂y = 0 exists no more.The extension of the Reynolds equation in the form of (17) to problems of modern configuration air bearing surfaces is questionable and needs further validation.Some comparisons with DSMC or other touch stone results or experimental data are desirable.For comparison with DSMC results even computation of small sized slider with complex configurations is of value.

DSMC Method
The direct simulation Monte-Carlo (DSMC) method [14] has been used for simulation thin film flows in sliders.Alexander et al. [15] studied the gas flow under a flat head surface above the moving drive platter and found excellent agreement with Reynolds equation, but their statement that the Fukui Kaneko's generalized Reynolds equation was a continuum hydrodynamic equation corrected for slip was wrong, and their simulations were restricted to small length sliders (L = 5 μm).Huang et al. [16] extended the DSMC simulation from two-dimensional to three-dimensional cases and again compared with Fukui-Kaneko's generalized Reynolds equation (they called it MGL model) and found that, overall, the two solutions agree well with each other and the agreement is better when the spacing decreases to about 5 nm, that is, when the Knudsen number is large.We stress here that the Fukui-Kaneko's generalized Reynolds equation is a global mass conservation equation, but the detailed flow field is calculated using the exact solution of the linearized Boltzmann equation.
So the agreement between the DSMC results and the results of generalized Reynolds equation is not surprising.The size simulated is again very small (L = 4 μm) in comparison with the authentic size of the modern slider.Wang and Liu [17] by using the DSMC method studied the pressure distribution in head/disc interface at the same position of different disc speeds and at different positions of the same disc speed and the lift force change percentage when air flow velocity changes or slider flying height changes.The size studied was again of the order of 4 μm, far smaller than the authentic size of 1 mm.It is the low information-to-noise ratio for low Mach number flows that makes the computational process of DSMC very time-consuming and only simulation results of miniature (not authentic) sizes are available.But the DSMC method is valuable in the microgas flows under flying sliders, even not being able to simulate authentic size, as it has the testing merit, it can be used for small sized but complicated configurations to testify other methods.

Information Preservation (IP) Method and Test of It by Channel Flows
The information preservation (IP) method is proposed in [18,19] to treat the problems encountered by the DSMC method of the huge ratio of the noise to the useful information and the demand of extremely large sample size.This is a method imbedded in the DSMC method in which each simulated molecule is assigned two velocities: thermal velocity c and information velocity u i .The former is just the molecular velocity c in the DSMC method and is used to calculate the motion, collision, and the reflection of molecules at the surfaces following the same algorithms and models as the DSMC method.Besides c we suppose that each molecule carries the so-called information velocity (IP velocity) u i to record the collective velocity of the enormous number of real molecules represented by each simulated molecule.The IP velocities do not produce any influence on the motion of molecules and are used only for summation to obtain the macroscopic velocities; the primitive information is taken from the oncoming flow and the body surface.When the molecules reflect from the surface, collide with each other, experience force action, and enter from boundary, the IP velocities attain new values.
The IP method also assigns each cell the IP velocity U i and IP density ρ.The comprehensive elucidation of the rules governing the renewing of the information velocities u i and U i and ρ can be found in [9] and [28].In general 2D and 3D cases the mass and momentum conservation equations for U i and ρ are The integrals are taken over the whole volume and surface of the cell concerned, respectively.After a time step Δt the cell IP density attains increment according to (21): from where the density and pressure are also renewed: p = nkT.The increment of the IP velocity u i is, according to (22), The IP velocity U i of the cell is obtained as the average of u i in the cell.
Omitting description of IP simulation of many different flow fields we give here only some detailed simulation of micro channels, for the geometrical forms and flow patterns of two cases of channel and disc flow are similar.Besides there have been two means to test the calculation of the channel flow to show the validity of the IP method.Firstly abundant experimental results of pressure distribution and flow rated through the channel for Knudsen numbers in the transitional regime.And the degenerated Reynolds equation suited for channel flows has been suggested to be derived easily for the generalized Reynolds equation originally derived for the flows under flying sliders and has the merit of kinetic theoretical test stone.So we use both the experimental data and strict theoretical computations to show the validity of the IP method for calculation of the microflows.
The channel flow of fluids has a long history beginning with the famous Poiseuiulle flow, that is, the fluid flow driven by constant pressure gradient.But for gas because of the global mass conservation the pressure distribution along the channel could not be linear and the gas Poiseuille flow with constant pressure gradient is only a hypothetic flow.With the emergence of the MEMS (micro-electromechanical systems) and the technique of manufacture of them, LIGA (abbreviation of German words Lithographie Galvanoformung Abformung) and EDM (electrodischarge machining), fine experiments were carried in new micromachined channels (see, e. g., [29][30][31][32][33][34]). Microchannel/pressure sensor systems with microsized pressure sensors integrated with the flow system were fabricated and exact pressure distributions were obtained.The flow rates were also measured exactly.Nitrogen, argon, and helium have been used as the fluid medium.The Knudsen numbers of ∼0.5 to 2.5 and larger are obtained and provide database to test the transitional flow.Distinguished nonlinear pressure distributions were also obtained.In [30] the ratio of the pressure difference of the inlet to the outlet is approximately 0.34∼1.7.In [31] for helium pressure difference ratio is approximately 0.59∼ 1.29.In [34] the pressure difference ratio is even higher; for argon it is 0.94∼2.96.The experimental work provides a solid foundation for testing various methods intended to calculate the microflows.
The generalized Reynolds equation originally is derived for application in the thin film air bearing problem with the lower plate moving with a velocity U and the upper plate tilted.Shen suggests to use this Reynolds equation to solve the microchannel flow problem [27] in which the lower plate is stationary and the upper plate is parallel to the lower one.
Owing to the steadiness of the lower plate the right-hand side term vanishes; as U = 0 and Λ = 0, there is no any contribution of the Couette flow.Owing to the parallelity of the two plates the value H is a constant and also can be dropped from the equation.So the generalized Reynolds equation for application to the microchannel problems is degenerated to the form: The values of P on the inlet and outlet of the channel are to be specified to make the microchannel problem solvable.For the case of diffuse reflection, the fitted formula approximation of Q P,TR (Kn), ( 16), can be used, and the degenerated Reynolds equation attains the form: For the ease of integration the local Knudsen number Kn is most conveniently expressed through P; for example, for hard sphere model it can be written as where for we have for hard sphere (see [9, (2.222)]).p 0 is the pressure at the outlet, T 0 is the temperature of the gas, and μ is the viscosity of the gas at T 0 .The constant C has the physical meaning of the Knudsen number Kn out at the outlet of the channel (see (25), at outlet P = 1).Substituting ( 25) into (24), one arrives at where C m,TR,CFR is an unspecified constant to be determined from the integration and has the physical meaning of the flow rate across the channel normalized by the slipless flow rate value.Equations ( 23) and ( 24) are in fact the generalized Reynolds equation, degenerated for the microchannel flows from where pressure distributions can be calculated.The flow rate across the channel C m,TR,CFR has also been integrated in closed [35] C m,TR,CFR = 1 2   [26] with experimental data of [31].Kn o = 0.16.The pressure values indicated are the inlet pressures.

And the dimensional corresponding mass flow rate is
We emphasize here, that as the assumption ∂p/∂y = 0 is strictly satisfied here for flat channel, the degenerated Reynolds equation has a strict kinetic theoretical foundation.
The microchannel flow problem was solved by IP method in [26].The approach, where the pressure p is fixed at the ends as the same of the prescribed (experimental) condition and U i is allowed to change continuously in the process of computation and finally to reach the steady solution, is used here.In the practice of general IP method the IP process has no reverse influence on the DSMC process.In the solution of channel flow, where the velocities on the boundaries are to be regulated during the simulation, the varying IP velocities on the boundaries are used to continuously adjust the boundary conditions of the DSMC-IP procedure.This enables the DSMC much quickly to have the correct value on the boundaries.To make the calculation process convergent it is essential to use the conservative form of the mass conservation equation.And a superrelaxation technique is employed to speed up the convergence process.For more detailed elucidation of the solution of the channel flow problem by the IP method see [9,26].Here only the comparison of some results of pressure distribution given by the IP method and the experimental data and the degenerated Reynolds equation are given (see Figures 3, 4,  and 5).
As shown in the above comparisons the IP method gives pressure and mass flow results of microchannels in excellent agreement with the degenerated Reynolds equation and also experimental data; so it can be considered as a reliable tool in dealing with the microscale slow gas flows including the micro flows under flying sliders.The detailed IP method      simulation of internal rarefied gas flows can be found in [28] and [9].

Microflows under Flying Head Slider Solved by Information Preservation Method
The IP method was applied to simulate the 2D flat slider problem [36].The conservative scheme of continuum equation and the superrelaxation scheme suggested in [26] are also used here to accelerate the converging process.As the flat air bearing surface is inclined, except rectangular cells, some incised cells appear, and they are combined with the lower adjacent regular cells to avoid extreme small number of simulated molecules in one cell.At the same time a simple weighted averaging method is introduced to smoothen the no-physical oscillation of the cell density: where c 1 is a weighting factor.The pressure distributions for different slider sizes have obtained by the IP method.The L = 5 μm result is in good agreement with DSMC method and the Reynolds equation.For L = 25 μm good agreement is obtained between the IP and Reynolds equation results.
It is remarkable that while DSMC methods only provided results for L = 5 μm, the IP simulation result for authentic size head slider L = 1 mm is in excellent agreement with Reynolds equation solution (see Figure 6).This is the only comparison of the Reynolds equation with other methods for authentic-sized slider.
In [37] the IP method was used to solve the flows under 3D flat slider and a miniature head slider of authentic    complicated configuration.As flows under flat surfaces have been compared between the DSMC method and the Reynolds equation in [16], here only the comparisons of the IP method result with that of the DSMC method are cited from [37]. Figure 7 shows the pressure under the flat slider obtained by the two methods for flying height of 10 nm and the speed of U = 25 m/s.The length of the slider is 4 μm and its width is 3.3 μm.
Note that the DSMC sample process starts and lasts 0.2 million time steps but for IP sample only 200 time steps are used.The total computation time (the time of convergence process plus the time of sampling process) of the DSMC method is about 100 times longer than the IP method.When the slider surface speed is lower, for example, when U = 2 m/s, the DSMC results with 0.6 million sample time steps still have remarkable fluctuations (see Figure 8).

Nomenclature
A: Numerical constant B: Numerical constant c: Molecular thermal velocity D: The width of the slider along z direction h: Height of the gap between the head slider and the magnetic disc h 0 : Minimum height of the gap between the head slider and the magnetic disc H: Dimensionless height h/h 0 L: The length of the slider along x direction p: Pressure p 0 : Pressure at the front edge of the slider Q: Dimensionless time tω (sometimes for temperature: p = nkT) u: Gas velocity along x direction in Reynolds equation u i : Molecular information velocity in IP method U: Boundary velocity of the magnetic plate along x direction in Reynolds equation U i : Information velocity of cells in IP method v: Gas velocity along y direction in Reynolds equation w: Gas velocity along z direction (span wise direction) in Reynolds equation W: Boundary velocity of the magnetic plate along z direction (span wise direction) in Reynolds equation x: Length direction along the moving disc surface X: Dimensionless coordinate x/L y: Direction vertical to the moving disc surface z: Width direction along the moving disc surface Z: Dimensionless coordinate z/D ρ: Density of the gas μ: Viscosity of the gas σ: Reflectioncoefficient λ: Mean free path of the gas Λ: Bearing number Λ 1 : Bearing number calculated by U B Λ 3 : Bearing number calculated by W B

UFigure 1 :
Figure 1: A schematic model of the thin film air bearing flow.

Figure 2 :
Figure 2: Velocity profiles and the flow rates of the slip-less and slip Couette flows (the transitional case also is symmetric relative to the central line but is not shown here).

Figure 3 :
Figure3: Comparison of streamwise pressure distributions of helium flow given by IP[26] with experimental data of[31].Kn o = 0.16.The pressure values indicated are the inlet pressures.

Figure 4 :
Figure4: Comparison of the pressure distribution in a 1.2 × 40 × 4000 μm 3 microchannel for helium from the solution of the degenerated (24) (solid line) with experimental data[31] (symbols) and the IP method[26] (dashed lines almost obstructed by the solid lines).Kn o=0.15579.
Flow rate Q P : Flow rate of the Poiseuiulle flow Q P,C : Flow rate of the Poiseuiulle flow for continuum flow case Q P,SL : Flow rate of the Poiseuiulle flow for slip flow case Q P,TR : Flow rate of the Poiseuiulle flow for transitional flow case Q P,TR : Dimensionless flow rate of the Poiseuiulle flow for transitional flow case Q P,TR /Q P,C t: Time T: