A mean free path approach to the micro/nanochannel gas flows

We investigate the gas flows near to solid surfaces in terms of the local spatial variation in the molecular mean free path (MFP). Molecular dynamics (MD) is the appropriate scientific tool for obtaining molecularly-accurate dynamic information in micro and nano-scale gas flows, and has been used to evaluate the molecular mean free path of gases. In the calibration procedure, the viscosity of a gas in the homogeneous case can be recovered in our MD simulations and reach good agreement with the theoretical prediction and data from NIST. In surface-bounded gas flows, if the collisions between gas molecules and walls are counted, a spatially-varying mean free path is presented, and for the first time we have observed that the distribution of the free paths deviates from the exponential one and spikes appear in their distributions at larger Kn, i.e. in the transition flow regime. Based on elementary kinetic theory, the effective viscosity of the gas derived from the mean free path has been incorporated into the framework of the continuum-fluid dynamics equations, and micro-Couette flows are performed to demonstrate this potential application.


Introduction
The Knudsen number (Kn) is usually adopted to quantify the deviation from the continuum behaviour in gas flows, i.e. the traditional index of the degree of rarefaction, and is defined as the ratio of the molecular mean free path (MFP) of gases to the characteristic dimension of the flow. The features of a gas flow are known to vary significantly depending on the value of Kn. For example, a gas flow at Kn = 1 may either refer to a low pressure condition and thus a large MFP in high altitude atmosphere (i.e. spaceship flying environments) or a higher pressure flow in very small channels (i.e. shale gas recovery [1,2]). Incidentally, the free molecular flow (i.e. Kn>10) in nano devices can be observed even at standard temperature and pressure (STP).
The physics of gas flows in rarefied regimes can significantly deviate from the continuum-based Navier-Stokes-Fourier (NSF) descriptions. Furthermore, fluid dynamics usually involve the motion of gases near the surface of a solid body, i.e. the Knudsen layer, in which the flow is in the local thermodynamic non-equilibrium state [3,4] and the standard NSF breaks down [5,6]. Thus Lockerby et al. [7] suggested to incorporate, rather than directly model, Knudsen layer effects within the conventional continuum-based fluid dynamics equations. This aim can be achieved by capturing the Knudsen layer structure within the NSF framework by scaling the stress/strain-rate relation in terms of correction functions. In the past decade, various researchers had proposed several different scaling functions for the rarefied gas flows: wall function [6], Kn-dependent functions [8], power-law scaling [9], and recently developed double power series [10,11]. These scaling functions come from the curve-fitting of kinetic theory or the direct simulation Monte Carlo (DSMC) calculations, and can provide the effective gas viscosity, which varies as a function of the normal wall distance. The primary advantage of the scaling constitutive relation is that it enables us efficiently to incorporate the rarefied effects within the NSF framework by paying less computational cost than the discrete velocity method (DVM) [12,13] and DSMC [14].
While the above mentioned scaling functions are derived from kinetic theory or DSMC data, it is still highly desirable to explore other case-specific scaling laws by counting gas-wall interactions (i.e. physical walls), allowing the scaling constitutive relation to be extended to complex flows such as the shale gas (i.e. polyatomic molecules) and the evaporation process of hydrocarbon fuel droplets (i.e. chain-like molecules) [15][16][17]. The novelty of this paper is to explore the potential of a molecular MFP-based scaling constitutive relation, which can be incorporated into the framework of continuum-based fluid dynamics equations to capture the effects that deviate from the equilibrium state such as the Knudsen layer in near-wall regions. In next section, we will introduce a procedure to evaluate the molecular MFP of gases bounded by solid surfaces in molecular dynamics simulations. Section 3 will present the spatial variation of MFP in micro/nano-channels, i.e. the effect of solid walls on the MFP, and illustrate the variation of the gas viscosity (i.e. effective viscosity) in near-wall regions. Importantly, a MFP-based scaling constitutive relation is applied in micro-Couette flows. Conclusions are drawn in Section 4.

Kinetic prediction
The molecular mean free path (MFP) is of great importance and fundamental in gas dynamics. Based on elementary kinetic theory arguments [18][19][20][21], the MFP for a gas composed of hard spheres is given by: or a general expression in terms of the gas viscosity μ given by [20]: where m is the molecular mass, d is the molecular diameter, ρ is the density of the gas, T is the gas temperature, and k B is the Boltzmann constant. Equation (1) assumes that the gas is dilute and spatially homogeneous. When the gas is bounded by solid walls, i.e. spatially-inhomogeneous cases, however, the expression in Eq. (1) needs to be revised to consider the collisions between gas molecules and solid walls. Stops [22] first derived a space-dependent expression for the MFP of gas molecules leaving a wall, distributed according to the diffusive Lambert's cosine law of reflection; three decades later, Guo et al. [23] extended the Stops' model for MFP of gases in microchannels; more recently, Abramov [24] has introduced a new expression for the MFP in confined geometries by replacing the cosine law with the equilibrium Maxwellian to give where y = z/λ 0 , z is the normal distance from the wall, λ 0 the nominal MFP, and E 1 (y) is the exponential integral function, i.e. E 1 (y) = ∞ y e −z /z dz. Beside its kinetictheoretical interest, a space-dependent expression for the MFP could be of practical relevance and has been used to propose scaling laws for the viscosity of a gas in near-wall regions [23][24][25][26].

Direct evaluation
Interestingly, the molecular MFP can also be calculated by averaging recorded individual free paths in molecular dynamics (MD) simulations [25][26][27][28][29]. MD simulations of the rarefied gas flows not only tackle the motion of gas molecules, but model the gas-surface interactions such as gas adsorption [30,31]. The distances travelled by gas molecules between every two successive collisions are calculated in MD simulations. In the present study, the monatomic gas argon (Ar) molecules are used. The molecules move according to Newton's second law, and the intermolecular interactions between molecules are modelled through the Lennard-Jones (LJ) potential [31,32]: where U is the intermolecular potential, r ij is the intermolecular distance, ε is the energy parameter, σ is the molecular diameter, and r c is the cutoff distance. For a micro-channel gas flow, we choose our wall atoms to be platinum (Pt), and the gas-wall (Ar-Pt) interaction is also modelled with the Lennard-Jones potential. The equations of motion are integrated using a velocity Verlet algorithm with a time step of t = 0.002τ , where the characteristic time unit is τ = σ Ar m Ar ε Ar = 2.15 × 10 −12 s. The parameters and their values that were used in our MD simulations are summarised in Table I of [33]. For the measurement of individual free paths in MD, as shown in Fig. 1, we set a condition to judge the occurrence of a collision event between gas molecules: if the distance between two gas molecules is equal to or less than the collision diameter, i.e. r ij ≤ σ col , where σ col is the collision diameter, the two molecules have collided and we stop recording the free paths of the involved molecules. Because the Lennard-Jones potential is used, the two molecules can still move closer to each other and the distance between them can be smaller than the collision diameter, i.e. r ij < σ col . So the finite time spent during collisions must be excluded when calculating the individual free paths. We therefore restart recording the free paths when the distance between the molecules is larger than the collision diameter, i.e. when r ij > σ col . In this procedure, the counter for recording the individual molecular free path is switched on after the last collision between two molecules and it is not switched off until the next collision. In addition, collisions between gas molecules and wall atoms are also counted: if the gas molecule approaches to the wall and has crossed a virtual plane that is placed σ col away from the walls, i.e. r z ≤ σ col or r z ≥ (H − σ col ), where r z and H are the coordinate of gas molecules in the z direction and distance between two Fig. 1 Schematic of the calculation of individual free paths in molecular dynamics (MD) simulations: the parameters l and v are the individual free path and velocity of molecules, respectively; r ij is the distance between two molecules, σ col is the collision diameter which is used to determine the collision event; t is the time step in MD, and n is the number of time steps between two successive collisions of the ith molecule. It should be noted that x and y are orthogonal coordinate directions (i.e. the walls are oriented along the xy-plane) and the distance between the walls is along z direction solid walls, respectively, we say it has collided with the wall (see Fig. 1) and stop recording its free path. The procedures of MD simulations and the treatment of bounding solid surfaces can be found in [33].

Recovering the viscosity of gases
In order to demonstrate the capacity of our MD technique, we first recover the viscosity of a gas in a homogeneous case using the equilibrium Green-Kubo approach [34,35] where V is the volume of the system, and the time integration is the ensemble average of the auto-correlation of the stress tensor P xz . Using this, the gas viscosity has been calculated for temperatures ranging from 273 K to 473 K at a pressure of 1 atm, and the results are shown in Table 1; the viscosities obtained in our MD simulations reach good agreement with the theoretical prediction of a Lennard-Jones (LJ) fluid [36,37] as well as with the data for argon from National Institute of Science and Technology (NIST, http://webbook.nist.gov/chemistry/fluid/). The data for the viscosity of argon for a Lennard-Jones (LJ) fluid [36,37] and from the National Institute of Science and Technology (NIST) are also shown for comparison It is worth noticing that, in principle, a temperature-dependent collision diameter could also be defined based on the Chapman-Enskog expression of the viscosity for a gas composed of hard spheres [36]: where d h is also called the molecular diameter of hard sphere in kinetic theory. The dependence of the collision diameter on the temperature predicted by Eq. (6) is shown in Fig. 2 by using the values of viscosity from NIST and the viscosity of Lennard-Jones (LJ) fluids [36,37]. In both cases, it can be seen that the relationship is roughly linear and the collision diameter also decreases as the gas temperature increases. However, the collision  [36,37] is also included for comparison. It can be seen that quasi-linear (i.e. MD simulations) or roughly linear (i.e. the hard-sphere gas) relationship exists in explored gas temperatures diameters predicted based on Eq. (6) significantly differ from the one obtained in MD: the change of collision diameter as the gas temperature is insignificant in MD simulations; but in both hard-sphere gas and LJ fluids the molecular diameter significantly depends on the gas temperature indicating a variable hard-sphere (VHS) model that is commonly used in kinetic theory.

Assessing the MFP in micro/nano-channels
We now consider a stationary gas confined in micro/nano-channels at STP, as shown in Fig. 1, and evaluate the MFPs using MD simulations. In the channel x and y are orthogonal coordinates, and the two solid walls are parallel and oriented along the xy-plane with periodic boundary conditions imposed along the x and y directions. We simulate cases with different distances between the two parallel walls in the z-direction (i.e. the gap across the channel), varying from 20λ 0 to λ 0 /2, where λ 0 is the nominal MFP given by Eq. (1). These correspond to Kn ranging from 0.05 to 2.0. For Kn = 0.05, there are 137,952 gas molecules in the micro-channel; as Kn increases (i.e. the channel height decreases), the number of gas molecules decreases accordingly. The spatial variation of the MFP is obtained by assigning the individual free paths to small constant-width bins that divide the space in the z-direction [26,27,29]. We have investigated the bin-independence of MFP results by varying the number of bins from 100 to 2000. It is found that the MFP profiles across the channel do not significantly change beyond 500 bins for Kn = 0.2, 1000 bins for Kn = 0.1, and 2000 bins for Kn = 0.05. These correspond to the same bin width in each case. The MD results reported below use these numbers of bins. We first investigate the influence of a solid wall on the MFP by analysing the probability distribution of the individual free paths across the channel. In elementary kinetic theory, the free path distribution of gas molecules is claimed to be exponential in a homogeneous system (i.e. the gas is non-bounded by the solid surfaces) and it had been validated in our previous study [33]. However, the free path of gas molecules can be terminated at the wall in confined geometries (i.e. an inhomogeneous case), and a power-law distribution of the free paths had been reported in rarefied gases [27]. Therefore, to consider the effect of solid walls, we measure the free path distribution of gas molecules in micro/nanochannels. Our previous study showed that the free path distributions are still exponential at small Kn (i.e. in the upper slip and early transition regimes) [33]. As Kn increases (∼0.5), as shown in Fig. 3, spikes are observed in distributions of the free paths, and move to the small free paths. The spikes may be caused by the gas molecules hitting the solid walls and get significant at higher Kn (i.e. in small channels narrower than λ 0 ). As reviewed in section 1, gas flows at large Kn originally refer to the low pressure conditions in high altitude atmosphere and large mean free paths are expected. On the other hand, the transition flow (i.e. 0.1 < Kn <10) in micro/nano confinements can also be observed at STP, and as shown in Fig. 3, the free path distribution of gas molecules is dominated by small free paths. Note that the distributions of free paths presented in Fig. 3 are different from the ones in [33], which measured the free path distributions of gas molecules in the central bins of channels.
It should be stressed that some studies argued that the MFP should be defined as the averaged distance travelled by gas molecules between two successive collisions with other gas molecules only, and that collisions between gas molecules and wall atoms or gas molecules adsorbed on the surface should not be taken into account when evaluating individual free paths [28]. With such a definition, a constant and isotropic MFP would be expected. Compared to the decrease of MFP at the walls when gas-wall collisions are counted, our previous MD predicted a constant MFP when gas-wall collisions were disregarded [33], whatever kind of walls (i.e. specular or diffusive) were utilised. Collisions between the gas and bounding surfaces can be encompassed in the concept of the MFP, and the resulting spatial variation of the MFP may help in understanding transport of gases in Knudsen layers near solid surfaces.
Apart from the spatial variation, the system MFP of a gas in the micro/nanochannel is also calculated in our MD simulations, and as shown in Fig. 4, it varies as a function of Kn. An exponential decay fit can be expressed as λ s (Kn)/λ 0 = 0.255 + 0.784 exp(-Kn/1.122), and we find that the system MFP deceases with increasing Kn (i.e. decreasing H) but approaches to λ 0 if Kn → 0. Perumanath Dharmapalan et al. [29,38] also observed the variation of the system MFP as a function of the channel height, but their MFPs in nanochannels (height varying from 6 nm to 16 nm) are only several nanometres, which are much smaller than ours. Barisik and Beskok [28] used a fixed Kn (∼1.14) by changing both the channel height and gas density, and their data reaches reasonable agreement with our calculated system MFPs.

Scaling constitutive relation in near-wall regions
It is well known that the standard fluid dynamics equations are no longer valid in the Knudsen layer. However, their range of validity can be extended by introducing the modified constitutive relation [23,25,39,40]. To summarise the scaling functions in the  Barisik and Beskok [28] (filled green diamond) and Perumanath Dharmapalan et al [29] (filled red squares) are also shown for comparison. Note that Kn was kept a constant by changing both the density of gas and channel height in [28], and the range of the channel height varied only from 6 nm to 16 nm in [29] literature, we propose a general form of the scaling law for the constitutive relation in onedimensional flow case (it can be easily extended to the three-dimensional case) as follows: where (z) is an expression of the higher order non-linear correction terms, which vary with the chosen curve-fitting function such as the wall function from kinetic theory [6] or power-law function from DSMC [9,41]; u is the velocity of gas flows, τ is the shear stress, and μ is the viscosity of gases. It should be noted that the higher order non-linear correction terms are only considered when the non-equilibrium effect becomes significant or plays an important role in the flow region such as the Knudsen layer. Otherwise, Eq. (7) reduces to the linear constitutive relation in the framework of continuum-fluid dynamics equations when (z) vanishes.
In the extended hydrodynamic models, the effective gas viscosity that varies with the normal distance to the nearest solid wall can be expressed by: where μ 0 is the nominal viscosity of the gas. We also find it convenient to recast the expression for the effective viscosity into an expression for the MFP of the gas. In elementary kinetic theory [18,36], the gas viscosity is proportional to the MFP. As a result, the spatial variation of MFP can be used to define a phenomenological scaling law for the viscosity of gases near any solid surfaces: For example, based on Eq. (3), Abramov [24,42] proposed a geometry-dependent viscosity; Dongari et al. [27] introduced a power-law based effective viscosity in rarefied gases; To et al. [26] also reported the local viscosity which varied with the MFP. According to the MFP-based scaling law, Guo et al. straightforwardly pointed out that the gas viscosity μ → 1 2 μ 0 at the wall (z = 0) if the distance between two plates becomes infinity (i.e. Kn → 0 in Kramers' problem); Fichman and Hetsroni [43] gave a fixed effective gas viscosity μ e = 1 2 μ 0 in the whole Knudsen layer by considering the interactions between gas molecules and the wall (i.e. momentum exchange events). Apart from the MFP-based scaling law, Lockerby et al. [6] predicted an effective gas viscosity of 0.59μ 0 at the solid surface when they used a wall-function expression in the constitutive relation; Reese et al. [44] generalised the scaling function to give the effective gas viscosity of 0.62μ 0 at the wall. However, Lilley and Sader [9] contradicted these claims, and demonstrated that the flow exhibits a striking power-law dependence on distance from the solid surface and the effective gas viscosity at the surface should be zero due to the singularity of velocity gradient. Later Lockerby and Reese [41] revisited their scaling law and proposed an exponential scaling function, which results in the zero-viscosity of the gas at the wall. All values of the gas viscosity at the wall predicted by different scaling laws are summarised in Table 2, and the expressions for different scaling functions can be found in Table I of [45].
Comparison of the gas viscosity in micro/nano-channels predicted by different scaling functions in Eq. (8) and by using Eq. (9) (as given by MD) illustrates that: all scaling functions have overpredicted the gas viscosity, particularly in the Knudsen layer; the generalised model proposed by Reese et al. [44] agrees well with the wall-function expression [6] except a larger gas viscosity at the wall; the power-law model introduced by Lilley and Sader [9] can only predict the gas viscosity within the Knudsen layer of one-MFP thickness, but it may fail when Kn varies. As can be seen in Fig. 5, the MFP-based gas viscosity also varies near the solid surface and decreases with increasing Kn. Our previous study has already applied the linear-response relationship (i.e. Newtonian's law) to calculate the gas viscosity and validate the MFP-based scaling function [33].
Although the Knudsen layer also affects the relationship between heat flux and temperature gradient, we restrict our attention to its effect on the stress/strain-rate relationship in the present study. Figure 6a shows the uniform distribution of the shear stress P xz across the channel but reduction at the wall due to the strong gas-wall interactions. It also illustrates that the shear stress P xz in the bulk increases with increasing Kn, which is Lockerby et al. [6] wall-function method 0.59 Reese et al. [44] constitutive scaling law 0.62 Lilley and Sader [9] power-law dependence 0 Lockerby and Reese [41] power-law scaling 0  [41,46]. We would like to point out that the greater value of the shear stress at the wall reported by Tu et al. [47] may be caused by counting the volume of gas within the bin close to the wall (i.e. part of wall atoms included) when sampling the properties of a gas across the channel in their MD simulations. Figure 6b indicates the distribution of the strain rate across the channel. It can be seen that the strain rate also increases as Kn increases and the larger gradient of the strain rate has been observed in near-wall regions. It indicates that the gas viscosity varies near a solid surface and decreases with increasing Kn, which was theoretically predicted by Fei et al. [48] in terms of the Green-Kubo relation. We also observe that the gas viscosity at the wall reduces almost to half value it has in the bulk, i.e. 49% for Kn = 0.05, 48% for Kn = 0.1, 44% for Kn = 0.2 and 48% for Kn = 0.5. It can be concluded that the MFP-based scaling law provides important insight into the use of the local viscosity of gas in near-wall regions, particularly in the Knudsen layer. Additionally, the thickness of the Knudsen layer is also measured and defined as the distance from the solid surface at which the stain rate reaches within 5% of the Navier-Stokes value [6]. Figure 7 presents the comparison of the Knudsen layer thickness between our MD simulations and other data: augmented Burnett approach gives the minimum thickness (∼ λ 0 ); regularised Burnett method estimates the maximum thickness (∼ 5λ 0 ); and our MD results reach good agreement with the data of BGK Burnett approach in Kramers' problem (∼ 2λ 0 ), i.e. Kn→0.

Application of MFP to non-equilibrium gas flows
In order to explore the potential application of the spatially-varying MFP, we discuss a MFP-based scaling solution for a non-equilibrium flow. According to the relation for the MFP given by Eq. (3), the viscosity of gases can be scaled in the same manner, i.e. Eq. (9). We assume that the temperature of the gas flow varies insignificantly, i.e. an isothermal and incompressible flow. Considering a variable viscosity, the Navier-Stokes velocity equation for the shear flow reduces to [24]: with where u 0 is the velocity of top and bottom walls which move in opposite directions along the x axis, and the local viscosity of gases across the channel μ(z) is given by Eq. (9). The  [50], R13 (filled downtriangles) [51], and Regularised Burnett (pink stars) [52]. The dashed line refers to an exponential fit for our MD measurement. Note that the thickness of the Knudsen layer is normalised by the nominal MFP λ 0 and is related to the Kramers' problem, i.e. Kn→0 presented in the literature solution for the shear flow (i.e. micro-Couette flow) derived from the MFP-based scaling law is illustrated in Fig. 8 with Kn ranging from 0.05 to 0.2. Without loss of generality, the wall-dependent MFP in Eq. (3) proposed by Abramov [24] has also been used to scale the gas viscosity. Apart from the scaling constitutive solution, our MD results for the profile of the crosschannel gas velocity have been presented in Fig. 8 for comparison. Note that the no-slip velocity boundary condition is used in the scaling constitutive solution and the possibly different velocity of the wall itself is disregarded for the sake of argument [24]. For the purpose of comparing the scaling constitutive solution and MD-computed velocity profiles, we remove the velocity slip at the solid wall in MD simulations by subtracting its zero-distance value from the velocity profiles. The Navier-Stokes solution with a constant viscosity, moreover, has been included for comparison. Figure 8 indicates that good agreement between two scaling constitutive solutions with variable viscosity and MD results has reached up to the early transition region (i.e. Kn∼0.2). However, the Navier-Stokes solution with a constant viscosity fails to capture the non-equilibrium gas flow in the Knudsen layer. The MFP-based scaling constitutive relation can accredit us to capture the non-equilibrium flow near any solid surfaces.
In the scaling or extended constitutive relation, the non-equilibrium gas flow in the Knudsen layer is much more sensitive to the slip coefficient at the solid surface and it  [24], respectively; open squares (green) refer to the full MD simulations results; the solid line refers to the Navier-Stokes solution with a constant viscosity. The data is plotted for only half the channel (i.e. above one of the surfaces) and the velocity is rescaled by half of the value for shearing walls has still been an open question. This tenet, together with the second-order constitutive and velocity slip and temperature jump models [40,45], may need further investigation in rarefied or micro/nano scale gas flows and it is beyond the scope of this paper.

Conclusion
The molecular mean free path (MFP) in a gas has been evaluated by using molecular dynamics (MD). In a spatially-homogeneous case, we obtained good agreement of the gas viscosities between our MD results and LJ fluids/NIST data. In micro-channels (spatiallyinhomogeneous cases) we found that the MFP at the walls reduces smoothly to half of its bulk value, as long as the collisions between the gas and the wall are accounted for. If it is used as a tool to scale transport properties near to bounding surfaces as a simple way to describe Knudsen layers, then a definition that counts gas-wall collisions should be preferred.
Evaluations of the molecular MFP of gases open up an exciting new direction: explicitly connecting the governing gas/gas and gas/surface interactions at a molecular scale to the continuum-fluid dynamics in the bulk. Based on elementary kinetic theory, the MFP-based scaling law has been proposed to obtain the local gas viscosity and the spatial variation of the gas viscosity as a function of the perpendicular distance from the solid surface has been implemented to the framework of continuum-fluid dynamics equations.
While we also presented the application of a spatially-varying MFP model to micro/nano-channel flows (i.e. micro-Couette flows), with the aim of developing greater understanding of transport behaviour in near-wall regions, our detailed investigations into the Knudsen layer near a solid surface have illustrated the capacity of this MFP-based scaling constitutive relation to capture the non-equilibrium gas flows. It is recommended that an extended MFP approach will be applicable to develop a higher order hydrodynamic model of the transport of gases in highly-confined geometries in the future.
Availability of supporting data Not applicable.