Challenges in Modeling Gas-Phase Flow in Microchannels: From Slip to Transition

It has long been recognized that the fluid mechanics of gas-phase microflows can differ significantly from the macroscopic world. Non-equilibrium effects such as rarefaction and gas-surface interactions need to be taken into account, and it is well known that the no-slip boundary condition of the Navier-Stokes equations is no longer valid. Following ideas proposed by Maxwell, it is generally accepted that the Navier-Stokes equations can be extended into the slip-flow regime, provided the Knudsen number is less than 10^{− 1}. Improvements in micro-fabrication techniques, however, are now enabling devices to be constructed with sub-micron feature sizes. At this scale, the flow will depart even further from equilibrium and will enter the transition regime. In recent years, there has been considerable success in the implementation of second-order slip-boundary conditions to extend the Navier-Stokes equations into the transition regime. Unfortunately, as yet, no consensus has been reached on the correct form of higher-order approach, with theoretical and experimental studies revealing large discrepancies in the magnitude of the second-order slip coefficient. It is believed that these discrepancies can be explained by the fact that continuum flow analyses neglect the Knudsen layer, which extends approximately one mean-free path from the channel wall. In addition, comparisons between kinetic and continuum slip-boundary formulations reveal another important source of error due to different definitions in the first-order slip coefficient. The paper explains how these discrepancies have arisen and describes future research directions that may help reconcile the different forms of higher-order approach.

INTRODUCTION

During the last decade, the development of precision batch-processing fabrication techniques for constructing microelectro-mechanical systems (MEMS) has led to the development of an increasing number of microfluidic technologies. The potential application areas for gas-phase microflows are numerous and include miniaturized heat exchangers for cooling integrated circuits, portable gas chromatography systems for the detection of air-borne pollutants, micro-reactors for generating small quantities of dangerous or expensive chemicals, and novel high-throughput gas flow cytometers. The emergence of MEMS, however, has resulted in the growing realization that the fluid mechanics in such small-scale devices is not necessarily the same as that experienced in the macroscopic world. This has led to numerous questions regarding the applicability of conventional analysis tools for gas-phase microflows [1, 2, 3].

The present paper highlights some of the challenges that need to be addressed when extending the Navier-Stokes equations into the transition regime. Initially, the paper focuses on the three fundamental assumptions that need to be satisfied for the Navier-Stokes equations to remain valid. We then describe the different length scales and dimensionless parameters that can be used to characterize the breakdown of the Navier-Stokes equations. Finally, the paper presents a summary of the main theoretical and experimental higher-order slip-velocity boundary conditions that can be used to extend the Navier-Stokes equations into the transition regime. Throughout the paper, we emphasize the discrepancies between the various higher-order slip formulations in an attempt to concentrate future research efforts into providing a clearer understanding of the boundary conditions within the transition regime.

NAVIER-STOKES ASSUMPTIONS

In conventional macroscale applications, the mass flow and heat transfer can be modeled using the principle of conservation of mass, momentum (Newton's second law), and energy (first law of thermodynamics). In addition, the heat transfer processes are constrained by the second law of thermodynamics. These principles are generally expressed using a continuum-based description of the fluid. The resulting mass flow and heat transfer equations, collectively known as the Navier-Stokes-Fourier equations, form a set of nonlinear coupled partial differential equations. However, as discussed by Gad-el-Hak [3], there are three fundamental assumptions that need to be satisfied for the equations to remain valid:

1.	Newtonian framework. The Navier-Stokes equations are based on a Newtonian framework that assumes the fluid motion is non-relativistic (i.e., the characteristic velocity has to be much smaller than the speed of light). As pointed out by Gad-el-Hak [3], the Newtonian framework is an excellent modeling tool for most flow problems (including MEMS), the only exception being the flow of very high-energy particles at the atomistic length scale or stars/galaxies at the other extreme of length scale.
2.	Continuum approximation. The Navier-Stokes equations are based on the assumption that the fluid is infinitely divisible. In other words, local flow properties such as density, pressure, velocity, and shear stress can be defined as averages computed over fluid elements that are sufficiently large compared to the microscopic structure of the fluid, but small enough in comparison to the macroscopic gradients to allow the use of differential calculus to describe the variation in the flow properties. Gad-el-Hak [3] has stated that the continuum approximation is almost always satisfied, the exception being flows where the spatial length scale is comparable to the mean distance between the molecules, such as in nanodevices and in shock waves that are extremely thin relative to the molecular separation distance. It should be noted that the continuum approximation leads to an indeterminate set of equations that can only be closed by expressing the functional relationship between the stresses and the rate of strain and between the heat flux and the temperature gradient. The continuum approximation by itself does not lead to the Navier-Stokes equations. For example, higher-order constitutive relations between the stresses and the rate of strain (as found in the Burnett equations) can be used to close the conservation equations (as seen in [4]).
3.	Thermodynamic equilibrium. The Navier-Stokes equations are also based on the assumption that all flow properties are in the local thermodynamic equilibrium, implying that the macroscopic quantities within the fluid have sufficient time to adjust to their surroundings. The rate at which a fluid can attain the equilibrium condition depends fundamentally on the time between molecular collisions. If the temporal and length scales of the molecular collisions are small in comparison to the macroscopic flow variations, then the fluid will be able to adjust quickly to the new conditions. The characteristic length scale describing the molecular collisions is known as the mean-free path, λ, and is defined as the average distance traveled by a gas molecule before colliding with another molecule. When the mean-free path is, for example, one order of magnitude smaller than the characteristic length scale of the device, the macroscopic flow quantities will vary almost linearly in space. The assumption of thermodynamic equilibrium, therefore, implies that the shear stress is linearly dependent on the rate of strain (Newton's law of viscosity), and the heat flux is linearly proportional to the temperature gradient (Fourier's law). In addition, thermodynamic equilibrium gives rise to the no-slip and no-temperature-jump boundary conditions found in all macroscale flow descriptions. The crucial issue in modeling gas-phase microflows is that the mean-free path of the gas molecules is often comparable to the length scale of the device, thereby invalidating the assumption of thermodynamic equilibrium. Under these circumstances, a discontinuity will occur in the velocity at the wall, and the gas will effectively slide over the solid surface.

GAS-PHASE MICROFLOWS

Length Scales and Non-Dimensionalized Parameters

The analysis of gas-phase microflows depends on a number of important characteristic length scales and non-dimensionalized parameters. At the molecular scale, an important parameter is the ratio of the mean molecular spacing, δ, and the mean molecular diameter, d. Gases which satisfy the following equation:

are often referred to as dilute gases [3, 5]. Conversely, if the condition in Eq. (1) is not satisfied, the gas can be described as being a dense gas.

In a dilute gas, the intermolecular forces can be neglected, and the molecules spend most of their time in free flight between successive collisions. The probability of more than two molecules colliding is extremely low, and thus most molecular interactions take the form of binary collisions.

The dilute gas approximation (binary collision assumption) leads to the kinetic theory of gases and the Boltzmann transport equation. For a simple hard-sphere gas in thermodynamic equilibrium, the mean-free path, λ, can be shown to be given by

where d is the mean molecular diameter and n is the number density (the number of molecules per unit volume), given by n = δ^{− 3}. In addition, it can be shown that the root mean square molecular speed c _rms is given by

where R is the specific gas constant and T is the temperature. The characteristic time between molecular collisions can thus be written as

The assumption of a thermodynamic equilibrium requires that the characteristic time between molecular collisions, t _c , is small in comparison to the characteristic timescale of the macroscopic variations. In practice, the validity of the thermodynamic equilibrium assumption in the Navier-Stokes equations can be related to the Knudsen number, Kn, which is the ratio of the mean-free path of the molecules to the characteristic length of the device, L:

The Navier-Stokes equations (in conjunction with a slip-velocity boundary condition) are applicable for Kn < 10^{− 1}, although the no-slip boundary condition requires a more stringent constraint of Kn < 10^{− 3}.

The continuum assumption in the Navier-Stokes equations is only valid when there are sufficient molecules within a given volume to achieve a stable estimate of the macroscopic flow properties. As detailed by Bird [5], the ratio of the characteristic length scale, L, to the mean molecular spacing, δ, should satisfy

to achieve a statistically stable estimate of the macroscopic properties. In other words, there should be at least 100 molecular spacings along each face of the sample volume, giving a total of 1 million molecules within the volume. If the region of interest is smaller than the constraint shown in Eq. (6), the macroscopic quantities obtained by averaging will suffer from excessive molecular fluctuations.

Classification of the Flow Regime

The earliest classification system to differentiate between the various flow regimes within a rarefied gas was proposed by Schaaf and Chambré [6]. Their classification system was based solely on the magnitude of the local Knudsen number:

•	For Kn < 10− 2, the continuum and thermodynamic equilibrium assumptions are appropriate, and the flow can be described by the Navier-Stokes equations using conventional no-slip boundary conditions, although Gad-el-Hak [2] has suggested that the breakdown in the thermodynamic equilibrium assumption is discernible at Knudsen numbers as low as Kn = 10− 3.
•	In the range, 10− 2 < Kn < 10− 1 (commonly referred to as the slip-flow regime), the Navier-Stokes equations remain valid, provided tangential slip-velocity and temperature-jump boundary conditions are implemented at the walls of the flow domain.
•	In the range, 10− 1 < Kn < 10 (the transition flow regime), the continuum and thermodynamic equilibrium assumptions of the Navier-Stokes equations begin to break down, and alternative methods of analysis using the Burnett equations [4] or particle-based direct simulation Monte Carlo (DSMC) approaches [5] must be employed. The difficulty in analyzing the transition regime from a continuum perspective arises from the fact that the stress-strain relationship for the fluid becomes non-linear within approximately one mean-free path of the wall (the so-called Knudsen layer).
•	Finally, for Kn > 10, the conditions can be described as being a free molecular flow. Under such conditions, the mean-free path of the molecules is far greater than the characteristic length scale, and consequently, molecules reflected from a solid surface travel on average many length scales before colliding with other molecules. The intermolecular collisions are thus negligible in comparison to the collisions between the gas molecules and the walls of the flow domain.

It should be emphasized that the limiting Knudsen numbers in the above classification system are somewhat empirical, and that the boundaries between the different flow regimes often depend upon the particular geometric details of the problem. This can be attributed to the fact that the choice of the characteristic length scale, L, is rarely unique [7]. As detailed elsewhere [2, 5, 7], it is preferable to define the characteristic length scale using the gradient of a macroscopic quantity, for example, the density, ρ:

Figure 1 presents a graphical representation of the flow regimes experienced by a range of gas-phase microfluidic components, as reported by Beskok [8]. It can be seen that most microsystems currently operate in the slip-flow or early transition regimes, with the only exception being the flow through hard disk drive reader heads, where the Knudsen number exceeds unity. However, rapid improvements in micro-fabrication techniques are enabling microsystems to be constructed with sub-micron feature dimensions. At this scale, the flow will depart even further from thermodynamic equilibrium and will enter the upper transition regime. It should be noted, of course, that microsystems used in low-pressure flows or vacuum systems will experience even higher Knudsen numbers.

Figure 1 Characteristic length scales of typical microfluidic components and the corresponding Knudsen number, Kn, at standard atmospheric conditions. Adapted from Beskok [8].

An alternative method of representing the continuum and thermodynamic equilibrium conditions is shown in Figure 2. The graph illustrates the conditions for a hard sphere gas with a molecular diameter, d = 4 × 10^{− 10} m, closely approximating the average molecular diameter of air. The left-hand ordinate represents the characteristic length scale, L, as defined in Eq. (7), while the right-hand ordinate shows the length scale normalized with the molecular diameter, L/d. The bottom abscissa represents the density normalized with a reference density, ρ/ρ₀, which is equivalent to the normalized number density, n/n ₀. As an aside, it should be noted that the density ratio is inversely proportional to (δ/d)³. Finally, the top abscissa represents the average distance between the molecules normalized with the molecular diameter, δ/d.

Figure 2 Limiting criteria for the application of the Navier-Stokes equations. Adapted from Bird [5].

Figure 2 shows the limits of applicability of the Navier-Stokes equations as defined by Eqs. (1, 5, and 6), using the values proposed by Bird [5] (δ/d = 7, L/δ = 100, and λ /L = 10^{− 1}). The vertical line in Figure 2 represents the boundary between a dilute and a dense gas. Air at standard ambient conditions has a pressure of 1.013 × 10⁵ N/m², a number density of 2.68666 × 10²⁵ m^{− 3}, a density ratio of unity (by definition) and δ/d = 8.5 (i.e., air at standard conditions can be assumed to be a dilute gas). However, it should be emphasized that the value of δ/d for standard air is close to the limit of the dilute gas assumption.

The L/δ = 100 line in Figure 2 represents the limit of molecular chaos. The continuum approximation is only valid if there are sufficient molecules within a given volume to achieve a statistically stable estimate of the macroscopic flow properties. Bird [5] estimated that there should be at least 100 molecular spacings along each face of the volume to achieve a reliable estimate of the macroscopic quantities. The continuum assumption is only valid above the L/δ = 100 line.

The remaining line (Kn = λ/L = 10^{− 1}) represents the limit of the validity of the thermodynamic equilibrium assumption. Because the mean-free path is inversely proportional to the number density, n, the gradient of the thermodynamic equilibrium line (on a logarithmic graph) is three times steeper than the molecular chaos line.

Inspection of Figure 2 reveals that for a dilute gas, as the flow dimensions are reduced in size, the thermodynamic equilibrium assumption of the Navier-Stokes equations fails first, followed by a failure in the continuum assumption. Conversely, for a dense gas, the continuum assumption breaks down first, followed by a failure in the thermodynamic equilibrium assumption. It is clear from Figure 2 that the continuum and thermodynamic equilibrium constraints are very different. Particular care should therefore be taken in the interpretation of the limiting conditions of the Navier-Stokes equations. It is also apparent that the usual flow classification system proposed by Schaaf and Chambré [6], based solely on the magnitude of the local Knudsen number, is unable to describe the system completely, as the Knudsen number is not the only parameter that needs to be taken into account.

PRESSURE-DRIVEN FLOWS IN MICROCHANNELS

Slip-Flow Regime

In traditional flow analyses, a no-slip velocity constraint is enforced along all solid walls. In practice, the no-slip condition is found to be appropriate up to a Knudsen number of approximately 10^{− 3}. If the Knudsen number is increased beyond this value, rarefaction effects start to influence the flow, and the molecular collision frequency per unit area becomes too small to ensure thermodynamic equilibrium. Under these conditions, a discontinuity will occur in the velocity at the wall, and the gas will effectively slide over the solid surface. The discontinuity in the velocity was first described by Maxwell [9], who showed that the tangential slip velocity at the wall could be written as

where u is the velocity, x and y are the streamwise and normal coordinates, λ is the mean-free path of the gas molecules, σ is the tangential momentum accommodation coefficient, μ is the viscosity, and ρ and T are the density and temperature of the gas at the wall. The second term in Eq. (8) is responsible for thermal creep (or thermal transpiration), which generates a slip velocity in the direction of increasing temperature. In the absence of any temperature variations along the wall (∂ T/∂ x = 0), Eq. (8) takes the form of a simple first-order slip-boundary condition in which the slip velocity is directly proportional to the velocity gradient normal to the wall.

The tangential momentum accommodation coefficient (TMAC) is introduced in Eq. (8) to account for the reduction in the momentum of gas molecules colliding with the wall. For an idealized surface (perfectly smooth at the molecular level), the angles of incidence and reflection are identical, and therefore the molecules conserve their tangential momentum. This is referred to as specular reflection and results in perfect slip at the boundary (σ → 0). Conversely, in the case of an extremely rough surface, the molecules are reflected at totally random angles and lose, on average, their entire tangential momentum, a situation referred to as diffusive reflection (σ = 1). For real walls, some molecules reflect specularly and some reflect diffusively, and therefore a proportion of the momentum of the incident molecules is lost to the wall. The tangential momentum accommodation coefficient is defined as the fraction of molecules reflected diffusively.

A number of experiments have attempted to measure the TMAC under various conditions [10, 11, 12]. The experiments indicate that the accommodation coefficient is a function of the molecular weight of the gas, the energy of the incoming molecules, the wall material, the temperature of the gas, and the condition of the surface. Experimental values for the TMAC are generally found to be between 0.2 and 1.0: the lower limit is associated with exceptionally smooth surfaces, while the upper limit is associated with very rough or highly oxidized surfaces. Recent experiments conducted by Arkilic et al. [13, 14, 15, 16], Maurer et al. [17], and Colin et al. [18] have measured the flow rates in silicon micro-machined channels and extracted tangential momentum accommodation coefficients ranging from 0.8 to 1.0, demonstrating that the highly polished and ordered crystalline surface of a silicon substrate exhibits sub-unity (or incomplete) momentum accommodation.

A similar expression to Eq. (8) was derived by Smoluchowski [19] to describe the temperature jump at the solid-gas interface:

where y is the distance normal to the wall, σ _T is the thermal accommodation coefficient, γ is the ratio of the specific heat capacities, and Pr is the Prandtl number.

Equation (8) is normally remembered as Maxwell's slip-velocity boundary condition, but closer scrutiny of the derivation indicates that the equation should only be applied to planar surfaces. As a consequence, Maxwell's boundary condition is commonly misapplied in practical situations involving surface curvature. Maxwell [9] initially derived the slip boundary condition in terms of the shear stress (Eq. 66, page 253) but then simplified the resulting expression, presumably believing that the variation in wall-normal velocity could be neglected. However, if the shear stress is retained, then the generalized version of Maxwell's slip-velocity boundary condition can be written as

where τ is the tangential shear stress at the wall, q is the tangential heat flux at the wall, and p is the pressure. As demonstrated by Barber et al. [20, 21] and Lockerby et al. [22], the application of Eq. (8) instead of Eq. (10) leads to the loss of crucial physics in the simulation of gas flows over curved or moving surfaces. A practical example where this might be important is the modeling of flow around serpentine bends. Another important application involves the modeling of gas flows around micro- and nano-particles.

Analytical Models: First-Order Formulations

Pressure-driven slip flow within ducts and channels has received considerable attention for many years, possibly as a result of the ability to formulate analytical or semi-analytical solutions. Early investigations of non-equilibrium gas flows in channels were conducted by researchers in the rarefied gas community who were primarily interested in low-density macroscale applications. Analytical solutions are available for a number of simple cross-sections: circular [23, 24, 25, 26], annular [27], and rectangular [27, 28]. In addition, high-aspect ratio rectangular channels can be approximated using the parallel-plate analogy, as demonstrated by Arkilic et al. [13, 14, 15, 16] and Maurer et al. [17]. Using the parallel plate assumption, Arkilic et al. [13] showed that the mass flow rate through a long microchannel under isothermal conditions is given by

where H, W, and L are the height, width, and length of the channel, respectively, R is the specific gas constant, T is the temperature, P ₁ and P ₀ are the pressures at the inlet and outlet of the channel, and Kn₀ is the Knudsen number at the outlet.

Analytical Models: Second-Order Formulations

Analytical models derived using the first-order slip boundary condition in Eq. (8) have been shown to be relatively accurate up to Knudsen numbers of approximately 10^{− 1} [16, 29]. For Kn > 10^{− 1}, however, experimental studies have shown that models based on the first-order boundary condition show considerable discrepancies against observed data [17, 18, 30]. This has prompted many authors to investigate the use of second-order slip-velocity boundary conditions to extend the validity of the slip-flow regime to higher Knudsen numbers. As yet, no consensus has been reached on the correct form of second-order formulation. For example, Deissler [31] proposed a second-order slip condition of the form

while Beskok et al. [32] defined the second-order slip boundary condition using a Taylor series expansion:

More recently, Beskok and Karniadakis [33] have proposed a “higher-order” boundary-condition of the form

where b is an empirical parameter dependent upon the geometric configuration and is the non-dimensionalized distance normal to the wall. Using a similar technique, Xue and Fan [34] have suggested an alternative higher-order boundary condition that again involves only the first derivative of the velocity gradient normal to the wall:

Other second-order formulations have been proposed to remove the difficulties associated with the evaluation of the second-order derivatives at the wall. For example, Jie et al. [35] utilized the following condition:

where and are the dimensionless streamwise and tangential coordinates normalized with respect to the channel height and Re is the local Reynolds number.

Recently, Lockerby et al. [22] have suggested that higher-order constitutive relations can be used to evaluate the shear stress and heat flux terms in Maxwell's generalized slip boundary condition, Eq. (10). For example, the Burnett equations, which are derived from terms up to the second-order in a series solution to the Boltzmann equation, can be used to provide a higher-order slip-velocity boundary condition. Substituting the Burnett stress tensor and heat flux vector into Eq. (10) and restricting the analysis to linear higher-order terms only, Lockerby et al. [22] proposed the following “linearized Maxwell-Burnett boundary condition”:

where u and v are the velocity components tangential and normal to the wall, respectively. In flows dominated by gas-surface interactions, there is some justification in adopting a more accurate model at the boundaries than in the interior of the flow.

It can be seen from Eqs. (12, 13, 14, 15, 16, 17) that there are numerous variations in the second-order slip boundary condition. Furthermore, there is no consensus as to whether the term involving the accommodation coefficient (i.e., (2−σ)/σ) should be associated with both the first- and second-order slip components or whether it should only be associated with the first-order term. Clearly, this is of little concern for fully-diffusive boundaries but may introduce considerable discrepancies in the case of incomplete (i.e., sub-unity) momentum accommodation.

To gain a better appreciation of the problems associated with the second-order slip boundary condition, it is useful to consider the simple case of Poiseuille flow between two parallel plates. The generalized second-order slip boundary condition then takes the form:

where A ₁ and A ₂ are the first- and second-order slip coefficients.

Table 1 presents a comparison of the theoretical values of A ₁ and A ₂ that have been proposed in the literature [9, 22, 31, 36, 37, 38, 39]. It can be seen that as yet, no agreement has been reached on the correct value of the second-order coefficient. As discussed by Lockerby et al. [22], there is also some evidence to suggest that the slip coefficients may be geometry-dependent.

Table 1 Theoretical first- and second-order slip coefficients proposed in the literature for a fully-diffusive boundary

Author	Date	A 1	A 2
Maxwell [9]	1879	1	0
Schamberg [36]	1947	1	5π /12
Cercignani and Daneri [37]	1963	1.1466	0.9756
Deissler [31]	1964	1	9/8
Hsia and Domoto [38]	1983	1	0.5
Karniadakis and Beskok [39]	2002	1	−0.5
Lockerby et al. [22]	2004	1	0.15–0.19*

The lack of a universally accepted second-order slip coefficient makes it difficult to extend slip-flow predictions into the transition regime. Moreover, it can be seen from Table 1 that second-order models based on a simple Taylor series expansion of Maxwell's boundary condition [39] predict a decrease in the slip velocity at the wall, whereas other second-order models predict an increase in the slip velocity. The disparity in the sign of the A ₂ coefficient has (until recently) led to considerable uncertainty and confusion when trying to extend slip-flow analyses to higher Knudsen numbers. However, experimental observations by Lalonde [40] and Aubert and Colin [41] have conclusively shown that first-order models underestimate the observed mass flow rate, demonstrating that the second-order contribution needs to increase the velocity at the wall.

Recent experimental studies [17, 18, 41] have shown that second-order slip models can be used with considerable success, provided the second-order coefficient is chosen with care. For example, Aubert and Colin [41] demonstrated that a second-order model based upon Deissler's boundary conditions [31] can give good agreement with experimental mass flow data up to an outlet Knudsen number of 0.25. Furthermore, experimental studies by Maurer et al. [17] have shown that second-order models can be extended to an outlet Knudsen number of 1.46, provided the second-order coefficient is suitably calibrated to fit the experimental data. However, the uncertainties over the correct value of the second-order coefficient will undoubtedly be a major problem in extending the Navier-Stokes equations into the transition flow regime.

ESTIMATION OF THE ACCOMMODATION COEFFICIENT

One of the most critical parameters affecting any slip-flow model is the tangential momentum accommodation coefficient, σ, accounting for the gas-surface interactions at the wall. The accommodation coefficient is known to be a function of the molecular weight of the gas, the wall material, the temperature and the condition (roughness) of the surface. For microchannel flows, the value of the accommodation coefficient is generally obtained from experimental mass flow rate versus applied pressure data. For example, recent experiments by Maurer et al. [17] and Colin et al. [18] have extracted tangential momentum accommodation coefficients for silicon micro-machined channels ranging from 0.87 to 0.93, as shown in Table 2. These values are similar to earlier experimental values determined by Arkilic et al. [15, 16]. Despite the apparent agreement between the different experimental and data-fitting techniques, it is informative to revisit the theoretical assumptions that are used to estimate the accommodation coefficient.

Table 2 Recent experimental estimates of the accommodation coefficient for silicon micro-machined channels

	Maurer et al. [17]		Colin et al. [18]
Author Gas	N2	He	N2	He
Min. Kn0	0.054	0.17	0.002	0.029
Max. Kn0	1.1	1.46	0.16	0.47
σ	0.87	0.91	0.93	0.93
Slip model	Planar (Poiseuille flow) second order		Rectangular second order (Deissler [31])

In the absence of thermal creep, Maxwell's first-order slip-velocity boundary condition can be written as

If the mean-free path of the gas molecules is defined as

as in Maxwell's original analysis, then the value of the coefficient α in Eq. (19) is unity. However, more rigorous kinetic analyses of the Boltzmann equation for planar flows [37, 42, 43] have shown that α needs to be modified to

The “exact” leading coefficient of 1.016191 in Eq. (21) was obtained numerically by Loyalka et al. [42] using a BGK model of the Boltzmann equation, whereas Wakabayashi et al. [43] solved the linearized Boltzmann equation and obtained a value of 0.98737. In practice, α has also been shown to depend upon the accommodation coefficient [42, 44]. Unfortunately, most of the recent experimental estimates of the accommodation coefficient for silicon microchannels have assumed α is unity in accordance with Maxwell's original analysis (e.g., [17, 18]). The disparity between the kinetic (Boltzmann equation) and hydrodynamic (Maxwell's) estimate of α has resulted in considerable confusion among the gas-phase microfluidic community.

Most experimental researchers appear to have adopted Maxwell's slip-velocity definition (α = 1) and have shown that the highly polished and ordered crystalline surface of a silicon substrate exhibits sub-unity (or incomplete) momentum accommodation, with σ typically varying between 0.87 and 0.93. However, implementing the kinetic theory estimate for α leads to a totally different interpretation of the gas-surface interactions at the wall. As shown in Table 3, applying the commonly accepted kinetic theory estimate of α = 1.1466 to the experimental data of Maurer et al. [17] and Colin et al. [18] leads to accommodation coefficients much closer to unity, implying that most of the gas molecules undergo diffusive reflections at the wall. It is clear that the disparity between the hydrodynamic and kinetic estimate of α has hindered the understanding of gas-phase microflows and, to a certain extent, has led to a polarization between the hydrodynamic- and kinetic-based research communities. The importance of standardizing the determination of the tangential momentum accommodation coefficient cannot be overemphasized. Without a universally accepted approach, gas-phase microfluidics will continue to use different definitions of the slip-velocity boundary condition, hampering the development of reliable modeling tools for the slip-flow regime.

Table 3 Comparison of recent experimental estimates of the accommodation coefficient for silicon micro-machined channels using hydrodynamic (α = 1) and kinetic (α = 1.1466) slip-formulations

		Maurer et al. [17]		Colin et al. [18]
Author Gas		N2	He	N2	He
α = 1.0	σ	0.87	0.91	0.93	0.93
α = 1.1466	σ	0.938	0.978	0.998	0.998

ESTIMATION OF THE SECOND-ORDER SLIP COEFFICIENT

Recent experimental work by Maurer et al. [17] on silicon microchannels has demonstrated that the second-order slip coefficient, A ₂, has a value of 0.26 ± 0.1 for nitrogen and 0.23 ± 0.1 for helium. On the other hand, earlier work by Sreekanth [30] on low-pressure pipe flow determined a somewhat lower value of 0.14 for nitrogen. It should be noted that there are considerable discrepancies between the experimentally-determined values of A ₂ and the theoretical values listed in Table 1. One possible explanation has been proposed by Hadjiconstantinou [45], who showed that considerable care is required in interpreting the experimental data, as the values of A ₂ are derived from mass flow rate measurements and not from direct observations of the slip-velocity. Hadjiconstantinou demonstrated that the effect of the Knudsen layers within approximately one mean-free path of the channel walls could result in a significant underestimation of the experimental value of A ₂ by as much as 0.3. The lack of consensus over the correct value of the second-order coefficient is a major issue in extending the Navier-Stokes equations into the transition regime.

KNUDSEN LAYER EFFECTS

It is interesting to note that both Arkilic et al. [16] and Maurer et al. [17] observed an apparent decrease in the value of the tangential momentum accommodation coefficient with increasing Knudsen number. In both these studies, however, the evidence to suggest that the TMAC varies with the Knudsen number was based, rather surprisingly, on an analytical model that assumes the accommodation coefficient is constant along the microchannel. Arkilic et al. [16] considered a first-order slip model and investigated flows up to an outlet Knudsen number of 0.4, while Maurer et al. [17] employed a second-order model and investigated flows up to an outlet Knudsen number of 1.46. We believe the apparent decrease in the accommodation coefficient is merely an indication of the breakdown in the thermodynamic equilibrium assumption in the Navier-Stokes equations and the growing importance of non-linear effects within the Knudsen layers. To support this view, we note that kinetic models are able to capture Knudsen layer phenomena without the need to modify the accommodation coefficient.

CONCLUSIONS

This paper has presented a review of currently available boundary conditions for modeling gas-phase microflows in the slip-flow regime and has highlighted some of the challenges that need to be addressed when extending the Navier-Stokes equations into the transition regime. It is apparent that there is no consensus on the theoretical value of the second-order slip coefficient, but it is also surprising to note that there are still significant discrepancies as to the correct value of the first-order slip coefficient. Determining the tangential momentum accommodation coefficient also remains problematic on account of the different first-order models. For example, experimental researchers appear to have adopted Maxwell's slip-velocity definition and have found that silicon microchannels exhibit incomplete momentum accommodation. On the other hand, boundary conditions derived from kinetic theory lead to accommodation coefficients much closer to unity, implying that most of the gas molecules undergo diffusive reflection at the wall. The lack of any real consensus on the nature of these gas-surface interactions emphasizes the challenges that need to be addressed in understanding gas-phase microflows. Furthermore, the large variations in the second-order slip coefficient and the exact form of the second-order component are likely to remain a major problem in developing reliable modeling tools based around the Navier-Stokes equations. It is apparent that there is an urgent need for precise experimental data that not only focus on the mass flow rate but also address the measurement of the velocity profile within the Knudsen layer. Experimental velocity data in the vicinity of the wall would be invaluable in helping to validate higher-order boundary conditions.

NOMENCLATURE

A 1	=	first-order slip coefficient
A 2	=	second-order slip coefficient
b	=	geometric parameter dependent on shape
c rms	=	root mean square molecular speed, m/s
d	=	mean molecular diameter, m
H	=	channel height, m
Kn	=	Knudsen number, λ/L
L	=	characteristic length scale, m
n	=	number density, 1/m3
p	=	pressure, N/m2
P 1	=	pressure at channel inlet, N/m2
P 0	=	pressure at channel outlet, N/m2
Pr	=	Prandtl number
q	=	heat flux, W/m2
R	=	specific gas constant, J/kg·K
Re	=	Reynolds number, ρ uL/μ
T	=	temperature, K
t c	=	characteristic collision time interval, s
u, v, w	=	Cartesian velocity components, m/s
W	=	channel width, m
x, y, z	=	Cartesian coordinates, m
Greek Symbols	=
α	=	coefficient in first-order slip equation
γ	=	ratio of specific heat capacities
δ	=	mean molecular spacing, m
λ	=	mean-free path, m
μ	=	viscosity, Ns/m2
ρ	=	density, kg/m3
σ	=	tangential momentum accommodation coefficient
σ T	=	thermal accommodation coefficient
τ	=	shear stress, N/m2

Acknowledgments

The authors would like to thank the UK Engineering and Physical Sciences Research Council (EPSRC) for its financial support of this research under grant GR/S77196/01. Additional support was provided by the EPSRC under the auspices of Collaborative Computational Project 12 (CCP12).

Notes

*A ₂ is dependent on the Prandtl number and the ratio of the specific heat capacities of the gas.