Analysis of micro-scale gas flows with pressure boundaries using direct simulation Monte Carlo method

Analysis of micro-scale gas ¯ows with pressure boundaries using direct simulation Monte Carlo method

J.-S. Wu

, K.-C. Tseng

Department of Mechanical Engineering, National Chiao-Tung University, 1001 Ta-Hsueh Road, Hsinchu 30050, Taiwan Received 3 January 2000; received in revised form 15 May 2000; accepted 14 August 2000

Abstract

The development of a two-dimensional direct simulation Monte Carlo program for pressure boundaries using unstructured cells and its applications to typical micro-scale gas ¯ows are described. For the mo- lecular collision kinetics, variable hard sphere molecular model and no time counter collision sampling scheme are used, while the cell-by-cell particle tracing technique is implemented for particle movement. The program has been veri®ed by comparison of simulated equilibrium collision frequency with theoretical value and by comparison of simulated non-equilibrium pro®les of one-dimensional normal shock with previous reported work. Applications to micro-scale gas ¯ows includes micro-manifold, micro-nozzle and slider air bearing. The aim is to further test the treatment of pressure boundaries, developed previously by the ®rst author, by particle ¯ux conservation for gas ¯ows involving many exits, complicated geometries and moving boundaries. For micro-manifold gas ¯ows, excellent mass ¯ow conservation between the inlet and two exits is obtained at low subsonic ¯ows. For micro-nozzle gas ¯ows, with ®xed inlet pressure, the mass ¯ow rate increases with decreasing pressure ratio (exit to inlet), but remains essentially the same at pressure ratios much lower than that obtained by continuum inviscid analysis. For higher speci®ed pressure ratios, the locations of maximum Mach number moves further downstream as the pressure ratio decreases;

while, for lower speci®ed pressure ratios, the Mach number increases all the way through the nozzle to the exit. Eventually, supersonic speed is observed at the exit for pressure ratios equal to or less than 0.143.

Finally, for slider air bearing gas ¯ows of the computer hard drive, the simulated gas pressures, at di€erent rotating speeds, agree very well with previous studies. However, there exists strong translational non- equilibrium in the gas ¯ows at the high rotating speeds. The applicability of the treatment of pressure boundaries using the equilibrium Maxwell±Boltzmann distribution function is discussed in terms of the magnitude of the local Knudsen number at the pressure boundary for micro-nozzles and slider air bearing applications. Ó 2001 Elsevier Science Ltd. All rights reserved.

*Corresponding author. Tel.: +886-3-5731693; fax: +886-3-5720634.

E-mail address: [email protected] (J.-S. Wu).

0045-7930/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved.

PII: S0045-7930(00)00029-3

Keywords: Direct simulation Monte Carlo; Knudsen number; Micro-nozzle; Micro-manifold; Slider air bearing;

Pressure boundary

1. Introduction

Micro-electro-mechanical systems (MEMS) are becoming prevalent both in commercial and industrial research due to the rapid growth and developments in the semiconductor industry.

MEMS are the ®eld with strong multi-disciplinary research needs. Applications for these devices include consumer products (e.g., airbag triggers, micro-mirror displays), industrial and medical tools (e.g., micro-valves, micro-motors), and instrumentation (e.g., micro-pressure sensors, micro- shear stress sensors) [1].

The shrinking size of MEMS (order of a micron) has several advantages over the corresponding macro-scale systems. First, they can be fabricated quite inexpensive and in large quantities, borrowing mature techniques developed in the semiconductor industry. Second, the miniature inertia (due to small size) allows them to respond very quickly to excitation, enabling the fabri- cation of actuators and sensors with frequency responses previously unimaginable for related mechanical systems [2].

In the past, most research was conducted concerning the electrical and mechanical properties of MEMS. Relatively few researches [15,20±22,28] focused on the thermal and ¯uid e€ects on MEMS performance, including rarefaction, thermal creep and compressibility. Hence, the un- derstanding of these physical e€ects is very important for the design optimization of these very small-scale devices.

Most research on the analysis of thermal and ¯uid problems in MEMS devices has relied on computational or theoretical approaches due to the experimental diculties presented in a micro- dimensional environment. The computational and analytical tools available are generally based on the Navier±Stokes equations. A crucial underlying assumption of these equations is that the

¯uid may be treated as a continuum rather than as a collection of discrete particles, as is the case with the more dicult to solve Boltzmann equation [3]. Based on the continuum assumptions,

¯ows can be described in terms of macroscopic variables such as temperature and pressure rather than microscopic variables such as molecular velocity distribution functions.

Unfortunately, continuum assumption begins to break down when the mean free path (k, average distance traveled by molecules before collision) becomes comparable to ¯ow character- istic length (L). The ratio of these two quantities is known as Knudsen number …Kn ˆ k=L† and is usually used to indicate the degree of rarefaction. Traditionally, ¯ows are divided into four re- gimes as follows [4]: Kn < 0:01 (continuum), 0:01 < Kn < 0:1 (slip ¯ow), 0:1 < Kn < 3 (transi- tional ¯ow) and Kn > 3 (free-molecular ¯ow). As Kn increases, the rarefaction becomes important and even dominates the ¯ow behavior. Hence, the N±S based computational ¯uid dynamics (CFD) techniques are often inappropriate for higher Kn ¯ows, such as slip ¯ow, transitional ¯ow and free-molecular ¯ow.

For most MEMS ¯ows, the Knudsen numbers are generally located in the slip-¯ow and transitional-¯ow regimes due to their micro-scales, even at atmospheric condition. An example of such a micro-scale device is a side-driven micro-motor [5,6]. The typical gap between the base and the rotor is 3 lm, the corresponding Kn equal to 0.05. Another example of large Kn is the modern Winchester-type hard disk drive mechanism, where the read/write head ¯oats about 60 nm or less above the surface of the spinning platter [7]. The corresponding Kn is equal to 1.1. Hence, the rarefaction e€ects, often neglected in traditional analysis, have to be carefully considered and modeled in micro-scale gas ¯ows.

The direct simulation Monte Carlo (DSMC) method, ®rst introduced by Bird [8,9], has been recognized as a defacto method for studying rare®ed gas dynamics. It has been applied very successfully to rare®ed hypersonic ¯ows [8,9] and other more fundamental scienti®c problems, such as ¯ow instabilities [10,11]. In addition to the space science applications, it has also been utilized in the analysis of ultra-high vacuum technology [12±14]. Very recently, it was applied to rare®ed internal gas ¯ow problems such as channel, pipe and duct ¯ows and the results agree very well with experiments [15]. Most importantly, DSMC is the only practical way to deal with ¯ows in the transitional regime, without resorting to the dicult Boltzmann equation, which requires modeling an integral±di€erential (collision) term. Hence, for the analysis of rare®ed gas dynamics, the DSMC method has proved to be a very powerful tool.

Most applications of the DSMC have used structuredgrids [8,9] to discretize the physical do- main. For problems with complicated geometry, multi-block meshing techniques were developed

®rst by Bird [9], which involved two steps: dividing the ¯ow ®eld into several blocks followed by discretizing each block into quadrilateral (two-dimensional) or cubic (three-dimensional) meshes.

Subsequent research has been directed to develop alternative meshing techniques such as the coordinate transformation method by Merkle [24], the body-®tted coordinate system by Shimada and Abe [25] and the trans®nite interpolation method by Olynick et al. [26]. However, all of these still used structured grids. It is much easier to program the code using structured grids, however, it requires tremendous problem speci®c modi®cation. To alleviate such restriction, an unstructured grid system is an alternative choice, although it might be computationally more expensive. Boyd's group [16±18] has applied such techniques to compute thruster plumes produced by spacecraft and found that the results are very satisfactory. Wilmoth et al. [19] have used two types of grid (unstructured tetrahedral and structured Cartesian grids) to compute the low-density, hypersonic

¯ows about reusable launch vehicle. Both methods were shown to give comparable results;

however, it was concluded that unstructured grid o€ers certain advantages for grid re®nement and structured grid appears to o€er greater overall advantages.

Wong and Harvey [27] have further employed remeshing adaptive grid techniques in combi- nation with unstructured meshes to study ¯ows with highly non-uniform density. However, the above-mentioned studies were mostly applied to high-speed gas ¯ows.

For low-speed gas ¯ow applications, Piekos and Breuer [20] used DSMC to analyze the gas dynamics of micro-mechanical devices using unstructured grids. In their research, they developed a special way to accommodate the in¯ow/out¯ow pressure boundary conditions. However, Ref.

[20]Õs boundary condition treatment is not as e€ective as they claimed when treating pressure boundaries; furthermore, it has not been tested in more complicated applications. Nance et al.

[22] adopted the theory of characteristics from continuum gas dynamics and incorporated the

procedure in the DSMC method to handle the pressure speci®ed micro-channel gas ¯ows at a subsonic exit.

Alexander et al. [28] used DSMC to simulate the nano-scale gas ¯ow of a computer hard drive slider air bearing. The results showed that the gas pressure inside the air bearing is not uniform.

The pressure generally increases from atmospheric level with distance in the disk rotating direc- tion, obtains its maximum value near the end of the write/read head and decreases rapidly to the atmospheric value. The location of the maximum pressure value moves further downstream with increasing disk speed. In treating these pressure boundary conditions, they have also designed a special way to update the inlet and exit speeds such that the inlet and exit pressures are the same as speci®ed. However, it is not clear how they precisely implemented the treatment.

Recently, Wu et al. [23] developed an e€ective way to process pressure boundaries with the DSMC method and successfully applied it to gas ¯ows with simple geometries such as micro- channel and backward-facing micro-step con®guration. The developed pressure boundary treat- ment is inherently mass conserved because particle ¯ux conservation is applied at each pressure boundary, assuming thermal equilibrium. Also it has been demonstrated [23] that the procedure requires fewer time steps to ``converge'' to the steady-state solution as compared with the ap- proach of Nance et al. [22]. However, their application to gas ¯ows with more complicated geo- metry or with many inlets and exits are yet to be veri®ed. In addition, it is also crucial to understand its limit of applicability in rare®ed or micro-scale gas ¯ows.

From the previous review, we see a need to develop an ecient and accurate simulation/

analysis tool for thermal/¯uid problems in micro-scale gas ¯ows with pressure boundaries. Thus, the objectives of the current research are to develop a DSMC code using an unstructured grid system and to apply the code to several micro-scale gas ¯ows, including micro-manifold, micro- nozzle and computer hard drive slider air bearing by using the pressure boundary treatment developed previously [23].

2. Numerical method

2.1. The direct simulation Monte Carlo method

Due to the expected rarefaction caused by the very small size of micro-scale or nano-scale devices, the current research is performed using the DSMC [8,9] method, which is a particle-based method. The basic idea of DSMC is to calculate practical gas ¯ows through the use of a method that has a physical rather than a mathematical foundation. The assumptions of molecular chaos and a dilute gas are required by both the Boltzmann formulation and the DSMC method [8]. The molecules move in the simulated physical domain so that the physical time is a parameter in the simulation and all ¯ows are computed as unsteady ¯ows. An important feature of DSMC is that the molecular motion and the intermolecular collisions are uncoupled over the time intervals that are much smaller than the mean collision time. Both the collision between molecules and the interaction between molecules and solid boundaries are computed on a probabilistic basis and, hence, this method makes extensive use of random numbers. In most practical applications, the number of simulated molecules is extremely small compared with the number of real molecules.

The details of the procedures, the consequences of the computational approximations can be found in Bird [8,9]. In the current study, the variable hard sphere (VHS) molecular model [8] and the no time counter (NTC) [8] collision sampling technique are used to simulate the molecular collision kinetics except for the computer hard disk slider air bearing problem. Note that the corresponding molecular data including reference diameter (dref), reference temperature (Tref), and the viscosity temperature exponent (x) for each species are taken from those listed in Ref. [8].

Solid walls for all cases considered in this study are assumed to be fully di€usive (100% thermal accommodation) and are equal to 300 Kunless otherwise speci®ed.

2.2. Pressure boundary condition treatment

In order to perform accurate simulation for in¯ow/out¯ow pressure boundaries, general pro- cedure for treating these conditions by using the concept of particle ¯ux conservation has been developed and incorporated into the basic DSMC algorithm [23]. The basic idea is to update the in¯ow and out¯ow velocities by applying particle ¯ux conservation, assuming thermal equilib- rium, at each pressure boundary, such that the mass ¯ow rate conserves automatically and the simulated boundary pressures coincide with the imposed values. Note that the assumption of thermal equilibrium at the pressure boundaries may not be justi®ed for very high local Knudsen number ¯ow, which can be seen clearly from the results presented later for speci®c test cases.

However, we will show that the thermal equilibrium is correct for most MEMS related subsonic

¯ows. For completeness, the implementation of updating the pressure boundary conditions is brie¯y described. For more details, see Ref. [23].

For a given mean speed (V) and temperature (T), the particle ¯ux across a boundary surface with area A in a particular direction h with the surface normal can be determined, assuming equilibrium Maxwell±Boltzmann distribution [8,9], as

_N

A ˆnVmpfexp… q²† ‡ 

pp

q‰1 ‡ erf…q†Šg 2 

pp ; …1a†

where q ˆ V

Vmpcosh; …1b†

Vmpˆ



2kT m r

; …1c†

m is the molecular mass, n is the number density, Vmp is the most probable speed and k is the Boltzmann constant.

At the in¯ow pressure boundary, in¯ow pressure pi and temperature Ti are both known and thus the in¯ow number density niis known. By applying the particle ¯ux conservation for each cell interface m (area Am) at the in¯ow pressure boundary, the updated in¯ow streamwise velocity at cell m (outside the ¯ow domain) is computed as

…ui†_mˆ _N‡ _N

niAm ; …2†

where _N‡ and _N are the computed particle ¯ux (Eqs. (1a)±(1c)) into and out of the ¯ow domain (cell m) at the in¯ow pressure boundary using the updated …ui†_mand sampled streamwise velocity at each in¯ow boundary cell m, respectively, with number density ni and temperature Ti. The value of …ui†_mand the sampled streamwise velocity for each in¯ow cell m will be varied during the simulation and eventually attain approximately the same value as steady state is reached.

At the out¯ow pressure boundary, all ¯uid properties except pressure are computed from the simulation. Again, we apply the concept of particle ¯ux conservation rather than the theory of characteristics [22] for each cell at the out¯ow pressure boundary. Then, a similar procedure treating in¯ow conditions is carried out for updating out¯ow streamwise velocity …ue†_m (outside the ¯ow domain) for each out¯ow boundary cell m. Note that the out¯ow temperature for each cell interface is not given in advance and is set to the temperature sampled inside for each out¯ow boundary cell m during simulation. Additionally, the out¯ow number density is computed using the equation of state.

By applying the above procedures at the in¯ow and out¯ow pressure boundaries, the simulated in¯ow and out¯ow pressures are found to be consistent with the speci®ed in¯ow and out¯ow pressures and mass conservation holds automatically. In Ref. [23], the results of a micro-channel

¯ow using the current approach agreed very well with those of Nance et al. [22] and it was shown that fewer samplings are required to ``converge'' to the steady-state solution. Typical pressure pro®le, in a micro-channel, as a function of streamwise location from the current pressure boundary treatment along with that from Ref. [22] is shown in Fig. 1 with the following ¯ow conditions: argon gas (VHS gas), inlet temperature Ti ˆ 300 K, pressure ratio Pi=Peˆ 3 (Piˆ 160:839 kPa) and aspect ratio L=H ˆ 5 (L ˆ 5 lm). From the excellent agreement between these results, we are con®dent that the results of the current approach applying particle ¯ux conservation can at least be consistent with that of Nance et al. [22] using the method of char- acteristics concept.

Fig. 1. Comparison of the simulated pressure pro®les for the pressure boundary treatment by Nance et al. [22] and the current approach.

2.3. Other concerns

The macroscopic quantities such as mean velocities and temperature are sampled and averaged from cells, which may be triangular, quadrilateral or hybrid with both. Because of the unstruc- tured meshes used, cell-by-cell particle tracing technique similar to Piekos and Breuer [20,21] is adopted to locate the ®nal particle position at the end of each time step during computation. This might increase the computational load in locating the cell number of the ®nal particle position during each time step, however, it o€ers the greatest ¯exibility of handling di€erent types of boundary conditions and required much less problem speci®c change in programming. In prac- tice, the computational time required using a particle tracing technique is comparable to those of structured mesh with careful selection of time step and cell sizes using the same number of sim- ulated particles.

3. Results and discussion 3.1. Code veri®cation

For code veri®cation, we have computed the equilibrium collision frequency within an adia- batic enclosure, non-equilibrium pro®les within a one-dimensional (1-D) normal shock, and sensitivity test by using di€erent mesh types. Results from these investigations are brie¯y de- scribed.

First, we computed the equilibrium collision frequency of argon and helium, respectively, in a square enclosure with adiabatic wall conditions. The time step used is about 1/6 of the mean collision time and the cell size is about 1/2 of the mean free path. This rule of thumb of selecting simulation time step applies to all cases studied in this research. Further re®nement of cell size does not improve the results much. The ¯ow was run for 1000 time steps to allow any initial condition e€ects to die out, then sampled every two time steps until 2000 samples were obtained.

The results, insensitive for argon and helium, show that the computed values are within 1%, 1.9%, 3% and 5% of theoretical values, when the average number of particle per cell equals 100, 50, 20 and 10, respectively. These ®ndings are comparable to those of Piekos [21].

Second, a 1-D Mach 8 argon normal shock is simulated using the stabilization technique proposed by Bird [8]. The upstream velocity, temperature and number density are equal to 2548.6 m/s, 293 Kand 1E20 particles/m³, respectively. Typical results are illustrated in Fig. 2 along with those obtained by Bird [8]. Computed temperature pro®les illustrate that there exists strongly non- equilibrium between the normal and parallel translational degrees of freedom at such high Mach number. The present results are in good agreement with Bird [8].

Finally, the code is used to predict the argon gas ¯ows in micro-channels using triangular, quadrilateral and hybrid (triangular and quadrilateral) meshes. The purpose is to test the sensi- tivity of the results to the type of meshes used and the particle tracing technique for handling the di€erent types of meshes. Results of centerline streamwise velocities are presented in Fig. 3, where non-dimensional location is de®ned as x=L with the length of the channel L ˆ 22:5 lm and the channel height is 4.5 lm. The inlet pressure and temperature are set as 7.146 kPa and 300 K, respectively, with inlet-to-exit pressure ratio equal to 3. For the three meshes tested, the simulated

centerline streamwise velocities coincide very well with each other. This agreement applies to other

¯uid properties as well. This demonstrates that the current code is mesh-type insensitive. Thus, the current DSMC code is veri®ed to be accurate in tracing particles in unstructured meshes and in predicting both the equilibrium and non-equilibrium phenomena in rare®ed gas dynamics. For clarity of discussion in the following, all subscripts with ``i'' and ``e'' stand for inlet and exit of the

¯ows, respectively.

Fig. 3. Micro-channel centerline u-velocity pro®les calculated with di€erent meshes (L ˆ 22:5 lm).

Fig. 2. Computed temperatures in a Ma ˆ 8 argon shock, including parallel and normal components.

3.2. T-shape micro-manifold

A micro-manifold is chosen as the ®rst application case for the purpose of testing the general treatment of pressure boundaries with multiple exits, a device with potential applications in medical technology. The sketch of a T-shaped micro-manifold is shown in Fig. 4 with one inlet and two exits, where exit 1 is at the upper right hand side (horizontal channel) and exit 2 at the lower bottom side (vertical channel). For the current simulation, L=hiˆ 7, H=hiˆ 3, hi ˆ he1ˆ he2, and hiˆ 1 lm. The inlet pressure is ®xed and the pressures at exits 1 and 2 are varied such that the pressure ratios between the inlet and exits are 3:1:1 (case I), 3:1:2 (case II) and 3:2:1 (case III), respectively. Argon is the working gas with the inlet value (uniform) of temperature and pressure equal to 300 Kand 26.79 kPa, respectively. The resulting Knudsen number based on the inlet conditions is 0.2. The initial conditions inside the micro-manifold are set as follows: zero mean velocities, uniform temperatures as Ti and uniform number densities as ni. The simulation time step used is 8:33  10 ¹¹ s.

Table 1 shows the speci®ed, simulated pressures and their di€erence at the inlet and exits for all the three cases considered. Note that the simulated pressures are computed from the average of inlet or exit boundary cells with almost negligible gradients across the channel. From the results, the simulated pressures are generally within 1.2% of those speci®ed at the inlet and exits. This implies that the treatment of pressure boundaries by applying ¯ux conservation is successful in handling the micro-manifold with single inlet and multiple exits at low subsonic ¯ows.

Typical normalized mass ¯ow rates as a function of time are presented in Fig. 5 for the case of speci®ed pressures, Pi:Pe1:Pe2ˆ 3:1:1 (case I). This illustration monitors how the mass ¯ow rate at each pressure boundary converges to its steady value and how the total mass ¯ow conserves as well. Before the sampling starts, the mass ¯ow rates at the inlet and exits ¯uctuates violently (sometimes negative, not shown in the ®gure); however, after the sampling starts, the corre- sponding mass low rates begins to converge to the steady values. The sum of the exit ¯ow rates is greater than the inlet ¯ow rate before reaching steady state (for time less than approximately 0.01 s in Fig. 5) because of the initial conditions imposed. The initial sum of the exit mass ¯ow rates is larger than the inlet mass ¯ow rate because of the larger di€erence between the initial ¯ow ®eld density (the same as ni) and the initial boundary density at the exits. We have tried di€erent initial

Fig. 4. Sketch of T-shape micro-manifold (hiˆ he1ˆ he2ˆ 1 lm, L=H ˆ 7=3, L=hiˆ 7, Kn_iˆ 0:2).

conditions by varying the initial uniform density in the ¯ow ®eld from vacuum (lowest) to inlet value (highest). As expected, the ®nal steady-state mass ¯ow rates are the same for all cases. By applying the particle ¯ux conservation at each pressure boundary, the mass ¯ow rate summation of exits 1 and 2 quickly approaches the mass ¯ow rate of the inlet, within 0.6% in this case. Also about 54% of the mass ¯ows out of exit 1 (horizontal channel). This behavior is mainly due to the larger pressure drop required for the gas particles to turn into the vertical exit 2. Consequently, the mass ¯ow rate is less at exit 2, as illustrated in Fig. 5. Table 2 records the normalized steady mass ¯ow rates into or out of each pressure boundary and the total mass ¯ow rate for the three cases considered in the present study. As tabulated, the summation of simulated mass ¯ow rates at exit 1 and exit 2 are within 0.9% of the inlet mass ¯ow rates for the three cases tested. Mass ¯ow rate decreases with increasing exit pressures for either exit 1 or 2 as expected. However, as speci®ed pressure at exit 2 is higher than that at exit 1 (case II), most gas particles (more than 97%)

¯ow directly into exit 1 (horizontal channel), which e€ectively blocks exit 2 (vertical channel).

General ¯ow properties for the T-shaped micro-manifold, including velocity vector and u- velocity, v-velocity, Mach number, temperature and pressure contours, are illustrated in Fig. 6a±f

Fig. 5. The normalized mass ¯ow rate for a T-shape micro-manifold with Pi:Pe1:Pe2ˆ 3:1:1 and Piˆ 26:79 kPa.

Table 1

Comparison of simulated and speci®ed I/O pressures for T-shape micro-manifolds Case Pi:Pe1:-

Pe2

Inlet^a Exit 1^a Exit 2^a

Pspc:b Psim:c …Psim: Pspc:†=

Pspc: (%) Pspc: Psim: …Psim: Pspc:†=

Pspc: (%) Pspc: Psim: …Psim: Pspc:†=

Pspc: (%)

I 3:1:1 1.0 1.001 0.10 0.333 0.335 0.60 0.333 0.336 0.90

II 3:1:2 1.0 1.001 0.10 0.333 0.337 1.20 0.667 0.668 0.15

III 3:2:1 1.0 1.004 0.40 0.667 0.668 0.15 0.333 0.336 0.90

aAll pressure values are normalized to inlet speci®ed pressure value, Piˆ 26:79 kPa.

bSubscript ``spc.'' represents the speci®ed value.

cSubscript ``sim.'' represents the simulated value.

with Pi:Pe1:Pe2ˆ 3:1:1. From Fig. 6a±c, we can see that gas ¯ow accelerates gradually in the horizontal channel and achieves maximum speed at about 0.17 times the most probable speed based on inlet conditions, cmpi, at the turning corner and then decelerates rapidly across the corner region and ®nally accelerates gradually again further downstream in the horizontal channel with a centerline exit speed at about 0.17cmpi. Similarly, the gas ¯ow from the inlet of the horizontal channel accelerates across the corner region down to the exit 2 in the vertical channel with a smaller centerline exit speed (about 0.15cmpi) than that of exit 1. However, at this ¯ow condition with the same pressures at both the exits, more gas particles ¯ow into exit 1 than exit 2 (about 8%

more, as discussed earlier) since a higher pressure drop for the gas ¯owing into exit 2 is required at the turning corner (larger pressure gradient), as can be seen from Fig. 6f. Therefore, the Mach number at the center part of exit 1 (0.20) is higher than that (0.17) at the center part of exit 2 (Fig. 6d). Also appreciable slip velocities at the solid wall, in the range of 0.06±0.12cmpi, are observed at the inlet Knudsen number of 0.2. The simulated temperature ®eld (Fig. 6e) is noisier than the velocity ®eld, which is expected due to the low-speed ¯ows involved and the large sta- tistical variance of the instantaneous particle velocities used in obtaining the temperature data.

The resulting temperature is almost uniform at the inlet temperature, except for the 3±5% drop near the exit centers caused by gas expansion. Finally, the pressure contours, as shown in Fig. 6f, illustrate the larger pressure gradient in the ®rst half of the horizontal channel, while a much smaller value is observed in the second half of the horizontal channel and vertical channel. That is, most of the pressure drop occurs before the corner region in the horizontal channel and this also explains that the small subsonic speeds obtained at both exits 1 and 2.

Although there is no direct experimental data or theoretical prediction available to support the simulation conducted in the T-shaped micro-manifold, the physical trends qualitatively coincide with physical expectation. Thus, the generality of the pressure boundary treatment, using particle

¯ux conservation, for multiple inlet and exits is at least qualitatively established.

3.3. Micro-nozzle

A micro-nozzle is selected as the second tested case because of its complex geometry and its potential for application as a MEMS device ± micro-thrusters on spacecraft, for example. Ad- ditionally, this might serve as a test case to identify the limit of applicability of pressure boundary treatment [23], assuming thermal equilibrium, and results are presented later that clearly identify the problems at large Knudsen number ¯ows. A micro-nozzle with area ratio (ratio of inlet to throat area; AR, hereafter) of 2 is considered. The dimensions are shown in Fig. 7. L=h ˆ 6:0, where L and h are the nozzle length and throat height, respectively, and h ˆ 1 lm. The inlet pressure, Pi, remains ®xed at 80.4 kPa, while the speci®ed exit pressure, Pe, is reduced such that

Table 2

Comparison of simulated mass ¯ux of inlet and exits for T-shape micro-manifolds

Case Pi:Pe1:Pe2 Inlet Exit 1 Exit 2 Summation Error

… _mi= _mi;I† … _mi= _mi† … _me1= _mi† … _me2= _mi† … _me1‡ _me2†= _mi … _me1‡ _me2 _mi†= _mi (%)

I 3:1:1 1.00 1.00 0.541 0.465 1.006 0.60

II 3:1:2 0.755 1.00 0.981 0.028 1.009 0.90

II 3:2:1 0.751 1.00 0.098 0.910 1.008 0.80

the resulting pressure ratio of exit to inlet ranges from 0.667 to 0. The ``0'' case refers to vacuum speci®ed at the nozzle exit. Again, argon is used as the working gas. The resulting Knudsen number based on the inlet conditions is 0.067, while that based on exit conditions (simulated) is in the range of 0.098±0.537 (Table 3). For the vacuum case, any particles crossing the exit boundary during simulation are removed and no particles are introduced from outside this boundary, while at inlet pressure boundary particle ¯ux conservation is imposed.

Fig. 6. The vector and contour plots of micro-manifold with Pi:Pe1:Pe2ˆ 3:1:1 (a) vector plot (b) u-velocity (c) v-velocity (d) Mach number (e) temperature (f) pressure.

As shown in Table 3, the simulated centerline pressures at the cells adjacent to the nozzle inlet and exit are within 1% of those speci®ed at the inlet and exit for Pe=PiP 0:143. However, the deviation tends to deteriorate as the pressure ratio decreases. For example, for Pe=Pi ˆ 0:067, the simulated gas pressure is 37.3% higher than that speci®ed at the nozzle exit. The reasons for this over-prediction are explained as follows. The simulated pressure ratio for the vacuum case is 0.086, which is much lower (about 50%) than the value of 0.177 as predicted from isentropic supersonic expansion with an inviscid continuum analysis. For the vacuum case, supersonic speed at the exit center is also observed; however, the value is smaller than that from the continuum analysis. This result should be attributed to the viscous dissipation and two-dimensional phe- nomena caused by the cross-sectional area variation along the nozzle, which is neglected in the quasi-1-D inviscid continuum analysis. In addition, we can ®nd that, for Pe=Piˆ 0:067, the mass

¯ow rate is nearly the same as that for vacuum case. This suggests that the predicted exit pressure should be the same as that for the vacuum case, 0.086, if the speci®ed pressure is lower than this value. However, the predicted pressure is 0.092, which is 7% higher than 0.086. Hence, the de- viation between predicted and speci®ed pressures at the exit should be caused by the high local Knudsen number at the exit, de®ned as KnUx ˆ k=jU=…oU=ox†j, where U is either temperature or u-velocity. This can be seen clearly from Table 3, which lists the local Knudsen numbers, Knuxand KnTx at the center of the inlet and exit boundaries, respectively, at di€erent speci®ed pressure ratios. For example, the Knux at the inlet (Table 3) are generally on the order of 0.01 for all pressure ratios, while those at the exit are larger due to the non-negligible u-velocity gradient and larger mean free path. The value of Knux at the exit increases with decreasing speci®ed pressure ratio. Similar trends are also found for KnTxand Knqxwhere both numbers are included in Table 3.

In fact, it can be shown that the deviation between the particle number ¯uxes, based on the Chapmann±Enskog distribution function and the Maxwell±Boltzmann distribution function, is a function of both KnTx and Knux [8]. If we compare these results with pressure data in Table 3, we can ®nd that the predicted exit pressure begins to deviate more from the speci®ed value as the local Knudsen number, both KnTx and Knux, exceed approximately 0.05 (Pe=Pi6 0:067). This discrepancy increases with local exit Knudsen number (hence, decreasing speci®ed pressure ratio) as explained previously. This high local Knudsen number makes the equilibrium Maxwell±

Boltzmann distribution function problematic for introducing particles into the ¯ow ®eld. This de®nitely requires future study by using the Chapmann±Enskog distribution function, instead of the Maxwell±Boltzmann distribution function, for higher local Knudsen number pressure boundaries. Anyway, the current results con®rm that the pressure boundary treatment used herein is successful for gas ¯ows with local Knudsen numbers (both KnTx and Knux) less than 0.05 at the pressure boundary. Again the mass conservation at steady state between inlet and exit is excellent,

Fig. 7. Sketch of micro-nozzle (AR ˆ 2, L=h ˆ 6, h ˆ 1 lm).

Table3 Comparisonofsimulatedandspeci®edI/Opressuresformicro-nozzleswithARˆ2 Pe=PiInletaExitaError Pspc:bPsim:c…Psim:Pspc:†= Pspc:(%)KnTxKnuxPspc:Psim:…Psim:Pspc:†= Pspc:(%)KnTxKnuxKnd…_mi=_mi;vac†…_me=_mi;vac†…_me_mi†=_mi (%) 0.6671.01.001‡0.100.002450.022320.6670.6660.150.003550.036370.0980.5460.5400.67 0.5001.00.9990.100.001980.012530.5000.4980.400.007350.025330.1300.7440.7410.41 0.3331.00.9990.100.001990.014040.3330.3320.300.006250.036140.1910.8930.895‡0.23 0.2001.01.000‡0.000.000980.010000.2000.1981.010.021650.024540.2960.9670.968‡0.20 0.1431.01.002‡0.200.000360.014790.1430.144‡0.700.039540.042320.3780.9840.988‡0.37 0.0671.01.000‡0.000.003520.013060.0670.092‡37.310.145900.160740.5121.0101.014‡0.40 0.0001.00.9990.100.001460.01504±0.086±0.188170.181600.5371.0001.003‡0.26 a Allpressurevaluesarenormalizedtoinletspeci®edpressurevalue,Piˆ80:4kPa;Kniˆ0:067. bSubscript``spc.''representsthespeci®edvalue. cSubscript``sim.''representsthesimulatedvalue. dKnudsennumberbasedonlocalaveragetemperatureandpressure.

which shows (Table 3) that the maximum di€erence is less than 0.7%, although the data for Pe=Piˆ 0:067 might be problematic. In addition, the mass ¯ow rate increases with reduced pressure ratio and approaches a maximum value for a pressure ratio less than 0.143, as shown in both Fig. 8 and Table 3. However, this choked pressure ratio is much lower than that, 0.81, as predicted by 1-D, inviscid continuum analysis.

Figs. 9±11 illustrate the pressure ratio e€ects on the Mach number, v-velocity and temperature with AR ˆ 2 at Kniˆ 0:067 and piˆ 80:4 kPa. They will be explored in detail in the following.

Fig. 9a±f illustrate the Mach number contours in micro-nozzle in the order of decreasing exit pressure with the same inlet pressure. It is clearly that the exit Mach number increases with de- creasing pressure ratio as expected. From Fig. 9a±c, the centerline Mach number increases along the nozzle, a maximum value is obtained at the downstream side of the nozzle throat and then decreases further downstream in the diverging section of the nozzle. Note that the position of maximum Mach number moves further downstream as the exit pressure decreases. The exit Mach number is still subsonic for Pe=PiP 0:2. This is obviously di€erent from the continuum analysis, assuming a quasi-1-D inviscid ¯ow, due to the viscous e€ects and strong rarefaction [18]. The most signi®cant di€erence in results for di€erent pressure ratios is that the viscous layer grows to the point where it completely ®lls the nozzle exit plane for Pe=PiP 0:33. In Fig. 9d, with Pe=Piˆ 0:20, the centerline Mach number increases all the way from the converging section to the diverging section, and a substantial Mach number, Meˆ 0:74, is observed at the exit center. In the case of a vacuum exit (Fig. 9f), supersonic speed with Mach number of approximately 1.25 is obtained; however, as expected, it is much lower than the Mach number of 1.73 predicted from continuum analysis. Note that the data at Pe=Pi ˆ 0:20 is essentially the same as that for the vacuum case, as can be seen clearly from Fig. 9e and f. The data presented in Fig. 9a ¯uctuate more than the other due to the lower speed involved and the inherent statistics used in the DSMC method, as explained earlier. The situation improves as pressure ratio decreases (hence, the speed

Fig. 8. Variations of normalized mass ¯ow ratio with pressure ratio …Pe=Pi† of micro-nozzle (AR ˆ 2).

increases). Appreciable slip velocities are observed along the nozzle walls. For the vacuum exit condition, Mach numbers adjacent to the wall on the order of 0.70, close to the exit, are observed (Fig. 9f). The most striking features with micro-nozzle ¯ows is that there is no shock predicted inside the diverging part of the nozzle, as would be expected from 1-D inviscid continuum

Fig. 9. Mach number contours of micro-nozzle at di€erent pressure ratios with Pe=Piequal to (a) 0.667 (b) 0.50 (c) 0.33 (d) 0.20 (e) 0.067 (f) 0 (Kniˆ 0:067, AR ˆ 2).

analysis, for all the exit pressures considered. The absence of shocks might be due to the strong viscous and rarefaction e€ects, which deserves further study.

Fig. 10a±f present the v-velocity contour in the order of decreasing exit pressure. Obviously, the v-velocities are mainly caused by the geometry and ¯ow acceleration e€ects. At higher pressure ratios such as those in Fig. 10a and b, the v-velocities are about one to two orders of magnitude less than the u-velocities; hence, they are essentially zero (less than 2% of inlet most probable

Fig. 10. v-velocity contours of micro-nozzle at di€erent pressure ratios with Pe=Piequal to (a) 0.667 (b) 0.50 (c) 0.33 (d) 0.20 (e) 0.067 (f) 0 (Kniˆ 0:067; AR ˆ 2).

speed) considering the statistical uncertainties in the samplings. As the pressure ratio decreases, the v-velocities remain small in the converging section and increases to some appreciable amount

Fig. 11. Temperature contours of micro-nozzles at di€erent pressure ratios with Pe=Piequal to (a) 0.667 (b) 0.50 (c) 033 (d) 0.20 (e) 0.067 (f) 0 (Kniˆ 0:067, AR ˆ 2).

in the diverging section. For example, the v-velocity at the exit can be as high as 0.18cmpias shown in Fig. 10f with a vacuum exit.

Temperature contours in the order of decreasing pressure ratio are shown in Fig. 11a±f. The data (Fig.11a and b) exhibit a high degree of scatter for the lower pressure ratio case, the reasons are the same as explained earlier in T-shaped micro-manifold. Temperature distribution obviously is not uniform in both spatial directions, especially in the diverging part of the nozzle. Temper- ature tends to decrease along the nozzle due to ¯ow acceleration (gas expansion) as shown clearly in Fig. 11c±f. In addition, the exit temperature decreases with decreasing pressure ratio and can be as low as 0.66 times the inlet temperature, as is the case for the vacuum exit condition. Tem- perature jump in the converging part of the nozzle is negligible; while it becomes substantial in the diverging part of the nozzle due to the high local Knudsen number close to the walls. For ex- ample, temperature jump can grow to 0.11 times that of the inlet temperature with a vacuum exit, as shown in Fig. 11f.

3.4. Slider air bearing

The main purpose of choosing a slider air bearing as the ®nal tested case is because of the availability of similar simulation data available in the literature and the opportunity to test pressure boundary treatment for a moving wall. In addition, this inherently high Kn number test case, also serves to identify the applicability limits of the pressure boundary treatment. The test con®guration, as shown in Fig. 12, is the same as that in Alexander et al. [28] and is brie¯y described as follows. A hard sphere argon gas (dref ˆ 0:366 nm; x ˆ 0:5) with Tiˆ 273 K; ®xed height ratio, hi=heˆ 2…L=heˆ 100†, and Ma ˆ 0:08 (Kni ˆ 0:625; heˆ 50 nm), 0.50 (Kniˆ 2:084;

heˆ 15 nm) and 1.0 (Kni ˆ 0:625; heˆ 50 nm). Note that Tiis the inlet gas temperature, hiis the inlet height of air bearing, heis the exit height of air bearing, L is the length of the platter, Ma is the Mach number de®ned by the platter speed and inlet acoustic speed, and Kniˆ k=hi. Although the inlet and exit pressures are both at atmospheric level (P0), it is highly rare®ed due to the nano- scale gap between the write/read head and the rotating disk.

Present simulated centerline gas pressure pro®les are presented in Fig. 13 along with those obtained by Ref. [28] for comparison, where x=L is the non-dimensional location in the platter moving direction. Note that the data illustrated as lines are present simulated centerline gas pressures, while the solid symbols are present simulated wall pressures, obtained by the sum- mation of incident and re¯ective momentum of particles on the write/read head surface, and the

Fig. 12. Sketch of slider air bearing (L=heˆ 100, hi=heˆ 2, Piˆ Peˆ 1 atm, where heˆ 0:05 lm for case I and case III, and heˆ 0:015 lm for case II).

open symbols are gas pressures obtained in Ref. [28]. Note that gas pressure is obtained using ideal gas law, where total average temperature is used. As can be seen, the present simulated gas pressure pro®les are in good agreement with those in Ref. [28] for all three cases considered. Note that only gas pressures are reported in Ref. [28]. All the predicted pressures increase at ®rst along the air bearing, achieving the maximum values, and then decreasing rapidly at the end of the air bearing. The initial rate of increase of the gas pressures is essentially the same for Ma P 0:5, while the rate of decrease after the maximum pressure and the maximum value both increase with Mach number. Table 4 lists the speci®ed, present and previous predicted gas pressures [28] and nor- malized mass ¯ow rates at the inlet and exit for all cases tested. At low disk rotating speed, Ma ˆ 0:08, both the predicted gas pressures at the inlet and exit coincide with speci®ed atmo- spheric pressure very well. However, both the predicted gas pressures at the exit exceed those speci®ed for high disk rotating speed cases, as Ma P 0:5. It seems the current prediction performs better than Ref. [28] considering the pressure di€erence between predicted and speci®ed at the exit.

However, it is hard to compare both, since we do not really know how they process the pressure boundaries. Similar to micro-nozzle ¯ows, the inlet and exit local Knudsen numbers based on u-velocity and temperature gradients are also listed in Table 4. We can see that as the Knuxat the exit increases beyond 0.05 for the Ma ˆ 0:5 and 1.0 cases, the gas pressure is overpredicted by about 2.3% and 9.3% at the exit, respectively. The boundary pressure deviation trend is similar to that observed for the micro-nozzle ¯ow, although the critical value of Knux and KnTx is di€erent.

Comparing case II …Kni ˆ 2:08, Ma ˆ 0:5) with case III (Kniˆ 0:63, Ma ˆ 1:0), we ®nd that, even the Knux of case II is higher or about the same as those of case III, the deviation of boundary pressure prediction is lower. In addition, in case I, local Knudsen numbers are on the order of 0.05, however, the deviation is negligible. This seems to suggest that this non-equilibrium caused by high local Knudsen numbers is in¯uenced by local velocity at the pressure boundary. As

Fig. 13. Gas pressure pro®les for slider air bearing at di€erent Knudsen numbers and Mach numbers (L ˆ 5 lm for case I and case III, and L ˆ 1:5 lm for case II).

Table4 Comparisonofsimulatedandspeci®edI/Opressureforthin-®lmairbearings CaseMaInletaExitaError Pspc:bPsim:c…Psim:Pspc†= Pspc:(%)KnKnTxKnuxPspc:Psim:…Psim:Pspc:†= Pspc:(%)KndKnTxKnux…_mi=_mi;I†…_mi=_mi†…_me=_mi†…_me_mi†=_mi (%) I0.081.00.9990.100.62500.00610.04101.00.9990.100.960.03200.04331.001.000.9950.50 II0.501.00.9980.202.08350.03760.00611.01.022‡2.203.010.27800.85312.101.000.9920.80 III1.001.01.003‡0.300.62500.01370.01771.01.093‡9.300.870.29590.656213.231.000.9950.50 a Allpressurevaluesarenormalizedtoinletspeci®edpressurevalue,Piˆ101.3kPa. bSubscript``spc.''representsthespeci®edvalue. cSubscript``sim.''representsthesimulatedvalue. dSubscriptrepresentstheKnudsennumberbasedonlocalaveragetemperatureattheexitplane.

pointed out earlier, more research is needed to clearly identify the criterion to resolve this problem.

For a low rotating disk speed, Ma ˆ 0:08, the present predicted wall pressures agree well with predicted gas pressures. This means that at low rotating disk speed, the gas particles are in thermal equilibrium, and hence the gas pressures are isotropic, resulting in the same pressure distribution both at wall and in the gas. However, at Ma ˆ 0:5 and Kniˆ 2:08, obvious discrepancies are observed between the present predicted gas pressures and the wall pressures. The present pre- dicted wall pressures are generally lower than the present predicted gas pressures, especially those at the inlet region. This discrepancy increases as Ma increases. The deviations between the gas pressures and the wall pressures are due to the fast moving platter speed and the highly rarefaction of the ¯ow. Particles entering from the inlet are entrained to collide with the moving plate at the bottom such that substantial downward mean velocities (v-velocities) are formed at the inlet re- gion. Therefore, there are comparatively fewer particles from outside the inlet colliding with the upper wall, especially near the inlet. This trend of deviation is more clear for the higher Mach number case, Ma ˆ 1:0. This di€erence can be explained partly next by the non-equilibrium be- tween di€erent translational temperatures in the gas.

Fig. 14 presents the predicted ratios of translational temperature, Tx=Ty, as a function of non- dimensional location, x=L, at Ma ˆ 0:08, 0.5 and 1.0. For low disk rotating speed, Tx=Ty is nearly unity, which means the ¯ow is in thermal equilibrium. As the disk rotating speed increases up to Ma ˆ 0:5, Tx=Ty attains a value of approximately 1.06, which represents a slightly non-equilibrium gas ¯ow. A stronger non-equilibrium gas ¯ow occurs as Ma ˆ 1:0 with Tx=Ty equal to 1.15. These non-equilibrium gas ¯ows at high rotating disk speeds explains why the wall pressures are less

Fig. 14. Ratio of temperature components, Tx=Ty, as a function of x-position in slider air bearing for di€erent Mach numbers (L ˆ 5 lm for case I, and case III, and L ˆ 1:5 lm for case II).

than the gas pressures inside the air bearing since the gas pressures are obtained by averaging the translational temperatures (hence pressures) in three directions.

4. Conclusions

This investigation has successfully developed a DSMC code using unstructured meshes. The code is veri®ed by comparison with theoretical equilibrium collision frequency and with previous results for the temperature pro®les of a 1-D normal shock. The code is then used to test the general pressure boundary treatment developed previously [23] by applying it to three di€erent typical micro-scale (or nano-scale) gas ¯ows, including a T-shaped micro-manifold, a micro- nozzle and a slider air bearing of computer hard drive. The major conclusions of the study are as follows.

1. The developed code is highly ¯exible in handling ¯ow ®elds with complicated geometry without much program modi®cation.

2. The pressure boundary treatment of Ref. [23], assuming thermal equilibrium and applying particle ¯ux conservation, is proved successful in computing gas ¯ows with complicated geometry, multi-exits and moving boundary at low local Knudsen numbers. However, a clear criterion for the applicability of the equilibrium Maxwell±Boltzmann distribution function requires further study.

3. Simulated results of a T-shaped micro-manifold ¯ow show that the ¯ow accelerates in the horizontal channel, decelerates rapidly across the corner region, and ®nally accelerates gradually to the horizontal exit. It is found that as the exit pressure of the vertical channel is twice that of the horizontal channel, the ¯ow into the vertical channel is almost blocked with less than 2% of total mass ¯ow rate.

4. Simulated results of a micro-nozzle in the transitional regime show that the ¯ow accelerates at ®rst and decelerates to the exit for larger pressure ratio (inlet to exit). However, the ¯ow always accelerates along the nozzle for smaller pressure ratio. For the vacuum exit case, supersonic ¯ow is observed; however, with smaller value compared with that of the continuum case.

5. Simulated results of a nano-scale slider air bearing show that predicted gas pressures are in good agreement with those of Ref. [28]. In general, the gas pressure increases gradually, attains a maximum value in the later portion of the air bearing and then decreases rapidly to the atmo- spheric level at the exit. The non-equilibrium becomes stronger as the speed (Mach number) of the moving disk increases, which makes the predicted gas pressures to deviate more from the pre- dicted wall pressures.

6. For a micro-nozzle ¯ow, the over-prediction of exit pressure becomes appreciable as the Knux

and KnTxexceed approximately 0.05; however, the deviation is negligible at similar values of Knux

and KnTx for a slider air bearing. This suggests that the applicability of the current pressure boundary treatment depends not only on the Knuxand KnTxbut also on the magnitude of interface speed at the pressure boundary.

This study raises the question of using the equilibrium Maxwell±Boltzmann distribution function in highly non-equilibrium gas ¯ows. Also, in the current research all the particle-wall interactions are assumed to be fully di€usive, which may be suspicious for the silicone-made

micro-scale devices as shown by the experimental work of Arkilic [29]. Thus, more studies towards resolving the above two concerns are worthwhile and are currently in progress.

Acknowledgements

This investigation was supported by the National Science Council of Taiwan, Grant Nos. NSC- 88-2212-E-009-019 and NSC-89-2212-E-009-034.

References

[1] Scott WB. Micromachines hold promise for aerospace. Aviation Week and Space Technology 1989;138:1±5.

[2] Trimmer WSN. Microrobots and micromechanical system. Sensors and Actuators 1989;19:267±87.

[3] Gombosi TI. Gaskinetic theory. UK: Cambridge University Press; 1994 (Chapter 5).

[4] Scha€ S, Chambre P. Fundamentals of gas dynamics Princeton, NJ: Princeton University Press; 1958 (Chapter H).

[5] Mehregany M, Nagarkar P, Senturia SD, Lang JH. Operation of microfabricated harmonic and ordinary side- driven motor. IEEE Micro Electro Mechanical System Workshop, Napa Valley, CA, February 1990.

[6] Tai YC, Fan LS, Muller RS. IC-processed micro-motors: design, technology, and testing. IEEE Micro Electro Mechanical System Workshop, Salt Lake, February 1989.

[7] Tagawa N. State of the art for ¯ying head slider mechanisms in magnetic recording disk storage. Wear 1993;168:43±7.

[8] Bird GA. Molecular gas dynamics and the direct simulation of gas ¯ows, New York: Oxford University Press;

1994.

[9] Bird GA. Molecular gas dynamics, Oxford, England: Clarendon Press; 1976.

[10] Stefanov S, Cercignani C. Monte Carlo simulation of the Taylor-Couette ¯ow of a rare®ed gas. J Fluid Mech 1993;256:199±213.

[11] Golshtein E, Elperin T. Convective instabilities in rare®ed by direct simulation Monte Carlo method. J Thermophys Heat Transf 1996;10:250±6.

[12] Nanbu K, Kubota H, Igarashi S, Urano C, Enosawa H. Application of DSMC to the design of spiral grooves on a rotor of turbomolecular pumps. Trans Jpn Soc Mech Engrs 1991;B57:172±8.

[13] Lee YK, Lee JW. Direct simulation of compression characteristics for a simple drag pump model. Vacuum 1996;47:807±9.

[14] Lee YK, Lee JW. Direct simulation of pumping characteristics for a model di€usion pump. Vacuum 1996;47:297±

[15] Beskok A. A model for rare®ed internal gas ¯ows. Microscale Thermophys Engng 2000, in press.306.

[16] Boyd ID, Jafry Y, Beukel JW. Particle simulation of helium microthruster ¯ows. J Spacecraft Rockets 1994;31:271±81.

[17] Boyd ID, Penko PF, Meissner DL, DeWitt KJ. Experimental and numerical investigations of low-density nozzle plume ¯ows of nitrogen. AIAA J 1992;30:2453±61.

[18] Kannenberg KC. Computational method for the direct simulation Monte Carlo technique with application to plume impingement. PhD Thesis, Cornell University, Ithaca, NY, 1998.

[19] Wilmoth RG, Lebeau GJ, Carlson AB. DSMC grid methodologies for computing low-density. Hypersonic Flow About Reusable Launch Vehicles, AIAA paper no. 96-1812, 1996.

[20] Piekos ES. Breuer KS. Numerical modeling of micromechanical devices using the direct simulation Monte Carlo method. Trans ASME J Fluids Engng 1996;118:464±9.

[21] Piekos ES. Numerical modeling of micromechanical devices using direct simulation Monte Carlo method. MS Thesis. MIT, Cambridge, USA, 1996.

[22] Nance RP, Hash DB, Hassan HA. Role of boundary conditions in Monte Carlo simulation of microelectrome- chanical systems. J Thermophys Heat Transf 1998;12:447±9 (Technical notes).

[23] Wu JS, Lee WS, Lee F, Wong SC. Pressure boundary treatment in internal gas ¯ows at subsonic speed using the DSMC method. 22nd International Symposium on Rare®ed Gas Dynamics, Sydney, Australia, July 9±14, 2000.

[24] Merkle CL. New possibilities and applications of Monte Carlo methods. In: Belotserkovsk II, editor. Rare®ed gas dynamics, 13th International Symposium, 1985. p. 333±48.

[25] Shimada T, Abe T. Applicability of the direct simulation Monte Carlo method in a body-®tted coordinate system.

In: Muntz et al., editors. Rare®ed gas dynamics, Progress in Astronautics and Aeronautics, AIAA 1989. pp. 258±

70.

[26] Olynick DP, Moss JN, Hassan HA. Grid generation and application for the direct simulation Monte Carlo method to the full shuttle geometry, AIAA, paper 90-1692.

[27] Wang L, Harvey JK. The application of adaptive unstructured grid technique to the computation of rare®ed hypersonic ¯ows using the DSMC method. In: Harvey J, Lord G, editors. Rare®ed Gas Dynamics. 19th International Symposium, 1994. p. 843±9.

[28] Alexander FJ, Garcia AL, Alder BJ. Direct simulation Monte Carlo for thin-®lm bearings. Phys Fluids A 1994;6:3854±60.

[29] Arkilic EB. Measurement of the mass ¯ow and tangential momentum accommodation coecient in silicon micromachined channels. PhD Thesis. MIT, Cambridge, USA, 1997.