## Abstrac

The force-driven Poiseuille flow of dense gases between two parallel plates is investigated through the numerical solution of the generalized Enskog equation for two-dimensional hard discs. We focus on the competing effects of the mean free path λ, the channel width L and the disc diameter σ. For elastic collisions between hard discs, the normalized mass flow rate in the hydrodynamic limit increases with L/σ for a fixed Knudsen number (defined as Kn = λ/L), but is always smaller than that predicted by the Boltzmann equation.

#### Authors

##### Lei Wu^{1}, Haihu Liu^{2}, Jason M. Reese^{3} and Yonghao Zhang^{1}

^{1}James Weir Fluids Laboratory, Department of Mechanical and Aerospace Engineering, University of Strathclyde, Glasgow, UK. ^{2}School of Energy and Power Engineering, Xi’an Jiaotong University, 28 West Xianning Road, China. ^{3}School of Engineering, University of Edinburgh, UK.

Received 6 December 2015; revised 25 January 2016; accepted 27 February 2016; first published online 30 March 2016.

© Cambridge University Press 2016

Under ultra-tight confinement, the famous Knudsen minimum disappears, and the mass flow rate increases with Kn, and is larger than that predicted by the Boltzmann equation in the free-molecular flow regime; for a fixed Kn, the smaller L/σ is, the larger the mass flow rate. In the transitional flow regime, however, the variation of the mass flow rate with L/σ is not monotonic for a fixed Kn: the minimum mass flow rate occurs at L/σ ≈ 2–3. For inelastic collisions, the energy dissipation between the hard discs always enhances the mass flow rate. Anomalous slip velocity is also found, which decreases with increasing Knudsen number. The mechanism for these exotic behaviours is analysed.

## 1. Introduction

Gas transportation along micro/nanoscale channels exhibits unusual behaviour that cannot be captured by conventional computational fluid dynamics (Karniadakis, Beskok & Aluru 2005; Holt et al. 2006). One well-known phenomenon is the Knudsen paradox (Knudsen 1909): when a constant pressure difference drives flow along a narrow channel, the mass flow rate (MFR) displays a characteristic minimum as the inlet pressure is reduced, although the conventional Navier–Stokes equations simply predict a monotonic decreasing MFR. As the ratio of the gas molecular mean free path to the channel width (the Knudsen number, Kn) becomes appreciable in narrow channels, the Boltzmann equation is commonly used to describe the collective motion of a dilute gas from the hydrodynamic to the free-molecular flow regimes (Chapman & Cowling 1970). From numerical analysis of the Boltzmann equation (Ohwada, Sone & Aoki 1989a), the Knudsen minimum may be understood as a consequence of the competition between the gas velocity slip at the channel wall and the effective shear viscosity, which both increase with the Knudsen number: while larger velocity slip increases the MFR, higher effective viscosity decreases the MFR by making the velocity profile flatter. The minimum MFR occurs at Kn ∼ 1. Poiseuille flow with non-negligible Knudsen numbers has received much attention in recent years, due to the development of micro/nano-electromechanical systems (Karniadakis et al. 2005) and the shale gas revolution (Wang et al. 2014). The MFR obtained from the Boltzmann equation is a key parameter in the generalized Reynolds equation to describe gas-film lubrication problems such as the gas slide bearing (Fukui & Kaneko 1987, 1990; Cercignani, Lampis & Lorenzani 2007), as well as in the upscaling methods used to predict the gas permeability of ultra-tight shale strata (Darabi et al. 2012; Mehmani, Prodanovic & Javadpour 2013; Lunati & Lee 2014; Ma et al. 2014).

While the Boltzmann equation is a fundamental model for gas dynamics at low pressures, its applicability to a typical shale gas production scenario, with gas pressures of the order of 10^{6} –10^{7} Pa, is problematic (Ma et al. 2014). This is because the Boltzmann equation is derived under the assumption that the mean free path is much larger than the molecular dimension, but this is not the case for large gas pressures. Instead, the Enskog equation, which takes into account the finite size of gas molecules and non-local collisions, has been developed for dense gases (Chapman & Cowling 1970).

Although the original equation is for hard-sphere molecules, long-range intermolecular attractive forces can be included through a mean-field theory (Grmela 1971; Karkheck & Stell 1981), so that the van der Waals’ equation of state can be recovered. The extended Enskog equation has successfully described gas–liquid phase transitions, net condensation/evaporation and liquid menisci between two solid surfaces (Frezzotti, Gibelli & Lorenzani 2005; Kon, Kobayashi & Watanabe 2014; Barbante, Frezzotti & Gibelli 2015), with numerical results in good agreement with molecular dynamics simulations (Frezzotti 1998; Frezzotti et al. 2005).

For the Poiseuille flow of dense gases, there are three characteristic length scales: the mean free path, the channel width and the molecular diameter. The purpose of this paper is to investigate the influence of these competing length scales on the MFR in force-driven flow, through numerical solution of the Enskog equation. Our investigation is not, however, limited to elastic collisions between gas molecules, but is extended to dissipative collisions in dense granular gases, i.e. agitated ensembles of macroscopic grains modelled as inelastic hard discs or spheres (Brilliantov & Pöschel 2004; Aranson & Tsimring 2006). For simplicity, we consider the dense flow of hard discs in this paper, which is modelled by the generalized Enskog equation in § 2. In § 3, numerical results of the Enskog equation are presented, and analytical models are proposed to explain the unexpected behaviour of the MFR and the slip velocity. Finally, we draw conclusions in § 4.

## 2. The generalized Enskog equation for dense granular gases

Consider a granular gas composed of smooth hard discs of diameter σ and mass m, subject to an external acceleration a in the x direction, and confined between two infinite parallel plates of temperature T_{w} at y=±(L + σ )/2. Since the hard discs cannot fully occupy the region less than half a disc diameter away from a plate, L is regarded as the effective channel width.

The impact of two hard discs conserves mass and momentum but not kinetic energy. Suppose **v** and **v**^{∗} are, respectively, the precollision velocities of a first and a second hard disc. The corresponding post-collision velocities **v** and **v**^{∗} are given by

where **u** = **v** − **v**_{∗} is the relative collision velocity,** k** is the unit vector from the centre of the first disc to the centre of the second at the time of their impact, and α is the restitution coefficient characterizing the loss of the normal relative velocity, i.e. **k** **·** (**v’** − **v’**** _{∗}** ) = −α

**k · u**. The collision is elastic when α = 1 and inelastic when 0 ≤α < 1. When the disc number density n is small, so that the mean free path (1/2√(2)nσ χ (n), where χ is the pair correlation function defined below) is much larger than the disc diameter, the dynamics of the dilute granular gas is modelled by the generalized Boltzmann equation (Goldstein & Shapiro 1995). However, as the number density increases, the mean free path could become much smaller than the disc diameter, in which case the localized collision assumption in the Boltzmann equation is invalid and the finite area of the hard discs becomes important. The generalized Enskog equation (Esteban & Perthame 1991; Brey, Dufty & Santos 1997; Bobylev, Carrillo & Gamba 2000) describes the flow of dense granular gases, where the velocity distribution function f (y,

**v**), which is a function of the hard-disc velocity

**v**= (v

_{x}, v

_{y}) and spatial location y, is governed by

in which time is omitted because we are only interested in the steady-state profiles. Here F = m_{a}L/2k_{B}T_{w} is the dimensionless acceleration, k_{B} is the Boltzmann constant, n_{0} is the average number density, and **J** is the Jacobian of the transformation from (**v ^{˜}** ,

**v**

^{˜}_{∗}) into (

**v**,

**v**

_{∗}), where

**v˜**and

**v˜**

**are the pre-collision velocities producing post-collision velocities**

_{∗}**v**and

**v**

_{∗}. For hard-disc gases, collisions are instantaneous, and multiple encounters can be neglected (Chapman & Cowling 1970). The binary collision is non-local, as the distance between the mass centres of the two colliding hard discs is the hard-disc diameter. Such non-local collisions are characterized by the product of the distribution function and a pair correlation function at different locations, i.e.

Here k_{y} is the projection of **k** in the y direction, and the pair correlation function χ is given by (Sanchez 1994)

which, compared to other expressions (Henderson 1975; Baus & Colot 1987), is accurate for the solid fraction η= πσ^{2}nn_{0}/4 up to the square-packing density η = π/4. Note that in writing (2.2) the location y is normalized by the effective channel width L, the number density n is normalized by n_{0}, velocity is normalized by the most probable speed v_{m}=√2k_{B}T_{w}/m, temperature T is normalized by T_{w}, and the velocity distribution function f is normalized by n_{0}/v^{2}_{m}.

When the velocity distribution function is known, the normalized number density, flow velocity and temperature are given as follows:

The MFR, which is normalized by n_{0} mL^{2} a/v_{m} , is defined as

where the integration range of the normalized y is from −1/2 to 1/2 because the centres of the hard discs are confined to this region such that macroscopic quantities outside of this region of y vanish.

## 3. Numerical and analytical results

We solve the generalized Enskog equation (2.2) by using the fast spectral method (Wu, Zhang & Reese 2015) under the constraint

, for various values of L/σ and the global Knudsen number

together with the diffuse boundary condition (Galvin, Hrenya & Wildman 2007)

where

Note that the collisions in the Boltzmann equation are local, and the Boltzmann equation is recovered if we choose g(y, **k**,** v**) in (2.3) as g(y, **k**, **v**) = f (y, **v**) and Kn=1/2√(2)n_{0}σL.

**3.1. Elastic collisions: ****α**** = 1**

We consider the case where the external acceleration is very small (F = 0.0001). The Enskog equation (2.2) is solved from η ≈ 0, until the maximum solid fraction reaches the square-packing limit. The MFR is shown in figure 1 as a function of the Knudsen number at several values of L/σ . For the Boltzmann equation, the Knudsen minimum is clearly seen at Kn ∼ 1, and when Kn < 1 (> 1), the MFR decreases (increases) monotonically with Kn. For the Enskog equation, the MFR is significantly influenced by the tightness of the wall confinement. When the channel width is large (L/σ≥50), the MFR at Kn ≥ 1 is the same as that from the Boltzmann equation, while the MFR for Kn < 1 is smaller than that from the Boltzmann equation, and the difference increases as Kn decreases.

As L/σ decreases, the MFR becomes smaller, and the difference from the Boltzmann result gradually extends to large Knudsen numbers. For L/σ > 5, the MFR is not a monotonically decreasing function of Kn in the slip flow regime (Kn < 0.1). Instead, there exists a specific value of Kn at which the MFR is locally maximum, i.e. Kn = 0.02 when L/σ = 30, Kn = 0.05 when L/σ = 15, and Kn = 0.08 when L/σ = 10. For L/σ≤5, the Knudsen minimum actually disappears, and the MFR only increases with the Knudsen number. Also, in the free-molecular regime (Kn > 10), the MFR is larger than the Boltzmann prediction.

*FI G U R E 1**. (Colour online) The mass flow rate as a function of the Knudsen number, for various values of L/**σ** , when the collision between the hard discs is elastic. Dashed lines: numerical results of the Navier–Stokes equation (3.8) with the first-order velocity slip boundary condition (3.10). Dash-dotted lines: the asymptotic analytical solutions (3.5) as Kn → ∞; from top to bottom, L/**σ** = 1, 2 and 10.*

We now examine the influence of the channel width on the MFR in the free-molecular regime. This regime occurs at small solid fractions. As with the Boltzmann equation, the normalized number density is unitary across the channel; however, the non-local collision in the Enskog equation makes the collision frequency spatially varying, which consequently affects the MFR. For small external accelerations, the distribution function can be linearized around the global equilibrium distribution function f_{eq}(v) = exp(−**v**^{2} )/π as f (y,**v**) = f_{eq}(**v**)[1+h(y,**v**)]. The equilibrium collision frequency reads (as the pair correlation function is approximately 1)

where H(y) is unitary for |y|≤1/2 and zero otherwise. The mean collision frequency ν^{˜}(y)= ∫ƒeq(**v**)ν(y,**v**)d**v** is plotted in figure 2 for various values of L/σ . For the Boltzmann equation with localized collisions, the mean collision frequency is a straight line with the value 2.507. For the Enskog equation, however, the non-local collision reduces the collision frequency in the region within one hard-disc diameter from the wall. For instance, the mean collision frequency at the wall for L/σ≥1 is half that of the Boltzmann equation because half of the collision is shielded by the wall.