On the numerical simulation of rarefied gas flows in micro-channels

Classical continuum methods break down for rarefied gas flows in micro-channels. Two non-continuum hydrodynamic models are investigated here: compressible Korteweg fluid-like and Bi-Velocity (Volume Diffusion) hydrodynamic models. The pressure driven rarefied gas flow in a rectangular micro-channel in the whole range of Knudsen number is numerically considered. The Korteweg model shows improved results in comparison with standard Navier–Stokes up to Knudsen number of unity. The Bi-Velocity method is found to allow a match between experiments and numerical results over the full range of Knudsen number.


Introduction
Rarefied gas flows are modelled by solving the Boltzmann equation which is valid for the full range of Knudsen number [1]. Application of the Boltzmann equation and associated kinetic methods such as the Direct Simulation Monte Carlo (DSMC) remains computationally expensive and their deployment to solving complex engineering cases (such as those involving multiphase flows) are still to be explored [2]. The inappropriateness of the traditional continuum flow model, namely the Navier-Stokes equations in rarefied gases and the need for additional contributions to its constitutive equations was discussed early by Maxwell [3]. Van der Waals theory of capillarity led to the observation that a liquid-vapour interface actually represents a rapid but smooth transition of physical quantities between two fluids [4]. Based on this, Korteweg obtained constitutive equations for the fluid flow equations within which phase transition phenomena occur by incorporating in the stress tensor contributions from the gradient of density [5].
The development of extended continuum models in gases has shifted since. Higher order hydrodynamics equations are derived using Chapman-Enskog expansion [6]. Grad's moments method [7] and associated regularizations (e.g., R13 equations) are the second route usually explored [8]. Several issues related to these higher order (Burnett type) equations are now well-known and detailed in the literature [9,10]. García-Colín et al claim that problems with Burnett equations are in reality limitations of those methods and highlighted the need of alternative approaches [11]. Other routes to ameliorate predictions by the Navier-Stokes at large Knudsen numbers by the use of slip boundary condition corrections have also been extensively explored [12,13].
The inappropriateness of Navier-Stokes system to describe temperature field in a stationary gas in the hydrodynamic limit has been shown by Sone et al [14]. These authors proved that corrections to the standard constitutive equations to involve Korteweg diffuse interface type stress tensor are required to acquire appropriate description of heat transfer processes even in the vanishing Knudsen number limit [14]. An alternative approach to extend the continuum fluid model appeared more recently in the literature based on the concept of two different conceptual velocities [15,16]. According to this theory, several fluid velocities can be distinguished: (a) the mass velocity U m which is linked to the classical notion of mass motion and (b) the volume velocity U v associated with changes in the local volume of the fluid material. This model is termed Bi-Velocity (or Volume Diffusion) hydrodynamics. In classical Navier-Stokes and other continuum flow descriptions one deals with one velocity (U m ). Bi-Velocity model as an extension to the classical Navier-Stokes system deviates from the earlier mentioned extended gas flow equations that are for example based on solving the Boltzmann equation using series expansion in Knudsen number [11]. Similarities with Korteweg equations are however found [17]. In a recent review by Gorban and Karlin with an extensive analysis of the Korteweg model it was noted that it is the first post Navier-Stokes equation which remains within continuum mechanics while capturing non-equilibrium kinetic level phenomena. They concluded that Korteweg equations 'give us a hint of how post Navier-Stokes equations may look like' [18].
Bi-Velocity model was tested in a mixed convection problem and was found to provide reasonable results compared to Navier-Stokes for high Knudsen number flows [19].
Ewart et al provided experimental data for the mass flow rate of rarefied gas flows in a micro-channel for the whole range of Knudsen number [20]. Dadzie and Brenner derived an associated analytical solution for the rectangular micro-channel helium gas flow using a model of first-order slip boundary conditions [21]. Their solution (which was based on Bi-Velocity) showed agreement with the experimental data up to Knudsen number 5. Following this method, Lv et al incorporated an effective volume diffusivity coefficient to show an agreement of the Bi-Velocity equation with the experimental data to a Knudsen number of 50 [22].
The present paper expands on these investigations. The Korteweg model of compressible flows and the Bi-Velocity model are considered separately for the numerical simulations of the gas flow in a rectangular microchannel.
The paper is organized as follows: In section 2 we present the two new non-continuum hydrodynamic models as alternative to Navier-Stokes. In section 3 we describe their numerical implementations and associated boundary conditions. Section 4 presents the numerical results, comparison with experimental data and discussion. The paper conclusion is drawn in section 5.

The two non-continuum hydrodynamics models
It is well accepted that Navier-Stokes equations fail to describe rarefied gas flows. Among other continuum fluid models that are being investigated are Korteweg fluid-like and Bi-Velocity models [5,16]. These models have been reported to account for non-equilibrium effects in other configurations [23].

Compressible Korteweg fluid-like model
From Van der Waals theory of capillarity authors noted that near a liquid-vapour phase transition a thermodynamic formulation of a fluid flow equation includes contributions from a density gradient-energy term. This formulation introduces a capillarity stress tensor in the momentum equation [24]. This stress tensor, also referred to as Korteweg stress tensor, may be obtained as [25]: where ρ is the fluid density, I is the identity tensor and the term α(∇ρ ⊗ ∇ρ) is also sometimes referred to as 'the Korteweg tensor'. Phenomenological coefficients α, β, γ are material functions dependent upon ρ and called capillary coefficients. Their exact value and expression are still elusive [26]. From a dimensional analysis, they may be represented as: where α * is a scalar and μ is the dynamic viscosity. The standard Navier-Stokes set may be modified by adding the Korteweg stress (equation (2.1)) to its viscous stress tensor. A compressible Korteweg fluid flow system of equations is then written as:

Conservation of momentum
where the shear stress term Π is given by, with U denoting the transpose tensor of ∇U and I the identity matrix. The energetic heat flux is then expressed as [26]: where k is the thermal conductivity and K E is the gradient (internal) energy coefficient which is assumed to be constant and T the fluid temperature [26]. In the above equations, U is the fluid mass velocity, p is the scalar pressure and e in is the fluid specific internal energy.

Bi-Velocity model
The kinetic level derivation of the Volume Diffusion continuum equations is based on the inclusion of an additional transport term originating from particle molecular level spatial diffusion [16]. Two different macroscopic fluid velocities can then be derived based on different averaging methods. Bi-Velocity method recognizes existence of two (conceptually) different velocities. In one hand, mass velocity U m is proportional to the mass flux and is found in the continuity equation. On the one hand, volume velocity, U v , accounts for variation in volume occupied by the fluid mass. A volume flux density relates the two: with k m an additional transport coefficient: the molecular (or volume) diffusivity coefficient. Figure 1 depicts the two different concepts of macroscopic fluid velocity. Flux J c is viewed as a molecular level diffusive flux associated with gas molecule concentration and −J c is the gas volume diffusive flux. U m and U v are both macroscopic concepts by definition. In the description of Volume Diffusion, a fluid particle is no longer perceived as a point mass like in classical mechanics but with a varying volume element associated with it. The diffusive flux −J c is oriented in the direction from high density to low density in the case of compression and the reverse in the case of expansion. With no density variation J c vanishes.
Bi-Velocity hydrodynamics set of equations can be written as [19,27]:

Conservation of momentum
where the shear stress term Π is now given by, and the energetic heat flux, J u , is given by, The specific internal energy of the fluid is e in = (3/2)RT, with R being the specific gas constant. Note that in the set (2.9) the mass velocity U m is now denoted U as it is equivalent to the traditional mass velocity like in a classical system (2.3). The Fourier's law now applies to the entropic heat flux as [16]: Here, the volume diffusivity coefficient is related to dynamic viscosity, μ, as: where our investigation led to α * as: This value of the coefficient α * was previously observed in other classical rarefied gas flow configurations [19,27]. The mean free path is calculated using G. Bird formula, λ = (k b T)/((2πd 2 p) 1/2 ),where k b is the Boltzmann constant and d is the molecular diameter [28].
By substituting equations (2.14) and (2.15) into (2.8) we have the final expression of the molecular level volume (actually, concentration) flux as:

Numerical implementation
Based on the two hydrodynamic models described above we propose two numerical solvers developed on OpenFOAM platform.

Modification to rhoCentralFOAM in OpenFOAM
Numerical implementations of the two sets of hydrodynamics equations are done on the OpenFOAM platform. These are adaptations of the solver rhoCentralFOAM which comprises of a finite volume (FV) discretization using semi-discrete, non-staggered central schemes for collocated variables prescribed on a mesh of polyhedral cells that have an arbitrary number of faces and solved on a three-dimensional unstructured mesh of polygonal cells. rhoCentralFOAM was initially developed to simulate compressible flows with better shock capturing [29].
The solver is density-based solver and uses a physical form of the equations described in section 2. It initially starts with solving the inviscid version of the set of conservative equations in an iterative order. Density-weighted fields are calculated first. The momentum density is calculated from U  = ρU and the total energy density where the total energy is E e .
in U 2 2 = + | | The set of equations are consequently solved for ρ, U  and E.  Based on these the temperature, T, is evaluated by the subtraction of kinetic energy from the total energy as: where c v is the specific heat capacity at constant volume. The momentum and heat diffusion are introduced to the solvers by the inclusion of the appropriate diffusive terms in the governing constitutive equations. U and T are evaluated explicitly since the momentum and energy equations are solved explicitly as noted in the inviscid version above.
The solver starts with the solution of the momentum equation by solving for U  Mass velocity is then updated by solving U U . r =  This step is done before solving a diffusion correction for U For Bi-Velocity method the explicit component of the stress tensor is expressed: In equation (3.3) we note the appearance of the term ∇ · (μ∇J c ) which is neglected for the Korteweg solver. The Laplacian terms from the deviatoric tensor ∇ · (μ∇U) and ∇ · (μ∇J c ) are implemented implicitly and they form coefficients within the solution matrices, rather than values in the source vectors. Similar procedure is followed to implement the Korteweg compressible fluid model where the additional explicit component of the stress tensor is now: To obtain the solution for the energy equation, similar procedure is followed by solving first for E  as: The temperature, T is updated from ρ, U and E  from equation (3.1) before solving a diffusion correction equation for T via: We note, however, that the experimental rectangular micro-channel case considered here is isothermal and therefore solutions to the energy equations are not as important as the momentum equations.

Boundary conditions
Maxwell type first-order slip boundary conditions are adopted to accompany the different hydrodynamics equations. In more detail, we used the compressible slip boundary conditions as presented below [30]: where U w is the wall velocity, σ u is the tangential momentum accommodation coefficient and τ v is the new tangential shear stress, τ v = S · (n · Π), extracted from the new modified shear stress in (2.4) and (2.10). The tensor S = I−nn with j = J · S based on the new energetic heat flux in equations (2.6) and (2.12). The unit normal vector n is defined as positive in the direction of the flow domain. The specific heat ratio is named γ.

The micro-channel configuration
Dimensions of the channel were accurately set to the original from Ewart et al [20]. The purpose of their experiments was to complete the database of mass flow rate measurements, obtained for a gas flow in a single micro-channel (see figure 2) ranging from the hydrodynamics to the near free molecular regime. The average pressure in the experiments was between 67 000 Pa and 30 Pa. The experimental method that was used to measure the mass flow rate through the micro-channel involves the use of two large constant volume tanks that are much larger than the volume of the micro-channel. More details of the experimental procedure can be found in Ewart et al [20]. Their data (mass flow rate results) are used as a benchmark for our simulations. Fluid properties and physical coefficients that are used in the simulations are listed in table 1. The mean Knudsen number is defined: . It can be noted that the length (L) and the width (W) of the micro-channel are extensively greater than its height (h).
Tangential momentum accommodation coefficient was set to 0.9 for all simulations along the whole range of Knudsen number in this paper.
To perform cell independency tests on the three numerical methods developed during the present investigations, three grids composed of 100 × 50 × 50, 200 × 100 × 100 and 300 × 150 × 250 cells are considered. figure 3 shows the corresponding normalized velocity profiles along the vertical center line of the channel for all methods. It is observed that numerical results are equivalent for 200 × 100 × 100 and 300 × 150 × 250 grids. As a result, the grid containing 200 × 100 × 100 cells were selected for the reported results.

Results and discussion
We compare numerical predictions by Korteweg, Bi-Velocity, Navier-Stokes and a previously derived analytical model in [21] with the experimental data of Ewart et al [20]. The following non-dimensional mass flow rate equation is used for the analysis: The mass flow rate, M  through the micro-channel is given by: Figure 4 presents the comparison between the models. Compressible Korteweg fluid-like model with its additional capillarity contribution to the shear stress gives agreement with experiments up to Knudsen number of about unity with minor deviation observed around Kn mean between 0.2-0.6. Bi-Velocity model based on volume diffusivity coefficient as expressed in equations (2.14) and (2.15) appears to fit the experimental data over the full flow regime, including excellent agreement in the free molecular regime (10 < Kn mean < 100). Dadzie-Brenner solution which was based also on the Bi-Velocity hydrodynamics model, but with a form of a slip condition and another expression of the volume diffusivity coefficient agrees with the experimental data for Kn mean < 5. Knudsen studied gas flows through tubes in the transition and free molecular regimes. From his results, the normalized volumetric flow rate showed a minimum at a Knudsen number near unity [31]: The Knudsen paradox. The present numerical solution by the   Knudsen paradox is usually difficult to capture using a Navier-Stokes with slip conditions [13].
In general, figure 4 reveals that the two continuum hydrodynamics models which involved density gradient expressions in the shear stress constitutive equations are better in predicting this micro-channel gas flow. This corroborates recent observations by Gorban and Karlin [18].  Figure 5 presents the normalized pressure distribution along the streamwise direction for the Bi-Velocity, Korteweg and Navier-Stokes. The convex curvature is captured by the three methods identically for the low Knudsen number. Difference between Bi-Velocity and Korteweg is insignificant across the Knudsen number range. This may be explained by the fact that both method constitutive equations for the shear stresses contain the density gradient capillarity effects. For a higher Knudsen number that corresponds to the free molecular regime (Kn = 80), Navier-Stokes profile may be described as unphysical compared to the other two models that appear to predict a rather monotonic variation in the profiles from the low to the higher Knudsen number. Figure 6 shows the normalized velocity profiles in a cross-section of the streamwise direction in a position where the flow is fully developed for all methods. U 0 is the average velocity taken at the channel entrance. For Kn = 0.043 the mass velocity as derived from the Navier-Stokes has a small difference from Bi-Velocity and Korteweg mass velocity. The velocity profile differs from the parabolic shape expected from a pure hydrodynamic regime pressure-driven flow. The mass velocity decreases as the Knudsen number increases and vanishes for the Navier-Stokes at higher Knudsen number (figure 6(c)) which is consistent with the no mass flow predicted by this model at higher Knudsen number. Bi-Velocity predicts the highest mass velocity at the highest Knudsen number which is in agreement with the mass flow rate plotted in figure 4.

Conclusion
Two extended hydrodynamic models based on Korteweg and Bi-Velocity theories have been investigated for the simulation of micro-channel gas flows. The simulation results are compared with existing experimental data. The two methods agree well for pressure and velocity profiles. Regarding the Korteweg model, it has been shown that it can capture non-equilibrium effects and match the experimental data up to Kn = 1. Korteweg model diverges from solutions obtained with the Bi-Velocity model at higher Knudsen numbers. The Bi-Velocity model shows the excellent agreement for the mass flow rate over the whole range of Knudsen number. The investigation reveals new insights into the expression of the volume diffusivity coefficient. The capability of the Bi-Velocity model to capture non-equilibrium effects in rarefied gas flows was demonstrated previously. The present study reinforces this. The results provide further evidences of the possibility of simulating rarefied gas flow phenomena with relatively high Knudsen number by a pseudo-continuum fluid flow model. Volume Diffusion model can agree well with the experimental data over the whole range of Knudsen number up into the free molecular regime. These results trigger further future investigations into the mathematical structure of both the Korteweg and the Volume Diffusion type of continuum flow equations in describing non-local-equilibrium flows beyond the classical model. The study not only provides a simple, accurate and reliable numerical method for simulating rarefied gas flows in micro-channels but also shows how post Navier-Stokes equations may look like.