On the apparent permeability of porous media in rarefied gas flows

Abstract

The apparent gas permeability of a porous medium is an important parameter in the prediction of unconventional gas production, which was first investigated systematically by Klinkenberg in 1941 and found to increase with the reciprocal mean gas pressure (or equivalently, the Knudsen number). Although the underlying rarefaction effects are well known, the reason that the correction factor in Klinkenberg’s famous equation decreases when the Knudsen number increases has not been fully understood. Most of the studies idealize the porous medium as a bundle of straight cylindrical tubes; however, according to the gas kinetic theory, this only results in an increase of the correction factor with the Knudsen number, which clearly contradicts Klinkenberg’s experimental observations. Here, by solving the Bhatnagar–Gross–Krook equation in simplified (but not simple) porous media, we identify, for the first time, two key factors that can explain Klinkenberg’s experimental results: the tortuous flow path and the non-unitary tangential momentum accommodation coefficient for the gas–surface interaction. Moreover, we find that Klinkenberg’s results can only be observed when the ratio between the apparent and intrinsic permeabilities is ${\lesssim}30$; at large ratios (or Knudsen numbers) the correction factor increases with the Knudsen number. Our numerical results could also serve as benchmarking cases to assess the accuracy of macroscopic models and/or numerical schemes for the modelling/simulation of rarefied gas flows in complex geometries over a wide range of gas rarefaction. Specifically, we point out that the Navier–Stokes equations with the first-order velocity-slip boundary condition are often misused to predict the apparent gas permeability of the porous medium; that is, any nonlinear dependence of the apparent gas permeability with the Knudsen number, predicted from the Navier–Stokes equations, is not reliable. Worse still, for some types of gas–surface interactions, even the ‘filtered’ linear dependence of the apparent gas permeability with the Knudsen number is of no practical use since, compared to the numerical solution of the Bhatnagar–Gross–Krook equation, it is only accurate when the ratio between the apparent and intrinsic permeabilities is ${\lesssim}1.5$.

1 Introduction

Unconventional gas reservoirs have recently received significant attention due to the shale gas revolution in North America (Wang et al. 2014), but the accurate prediction of unconventional gas production remains a grand research challenge. The permeability of a porous medium is an important parameter in the unconventional gas industry, as it describes how fast the gas can be extracted. The permeability at the representative elementary volume scale can be plugged into macroscopic upscaling equations to predict the gas production (Darabi et al. 2012; Lunati & Lee 2014). Therefore, it is necessary to investigate how the permeability changes with gas pressure, properties of the porous media and gas–surface interactions.

For laminar flows in highly permeable porous media, Darcy’s law states that the flux $\boldsymbol{q}$ (i.e. discharge per unit area, with units of length per time) is proportional to the pressure gradient $\unicode[STIX]{x1D735}p$ :

(1.1) $$\begin{eqnarray}\boldsymbol{q}=-\frac{k_{\infty }}{\unicode[STIX]{x1D707}}\unicode[STIX]{x1D735}p,\end{eqnarray}$$

where $\unicode[STIX]{x1D707}$ is the shear viscosity of the fluid and $k_{\infty }$ is the permeability of a porous medium that is independent of the fluid. For this reason, $k_{\infty }$ is known as the intrinsic permeability.

For gas flows in low permeable porous media, however, the permeability measured from experiments is larger than the intrinsic permeability and increases with the reciprocal mean gas pressure $1/\bar{p}$ (Klinkenberg 1941; Moghaddam & Jamiolahmady 2016). In order to distinguish it from the intrinsic permeability, the term of apparent gas permeability (AGP) $k_{a}$ has been introduced, which is expressed by Klinkenberg (1941) as:

(1.2) $$\begin{eqnarray}k_{a}=k_{\infty }\left(1+\frac{b}{\bar{p}}\right),\end{eqnarray}$$

where $b$ is the correction factor.

Darcy’s law can be derived from the Navier–Stokes equations (NSEs). The variation of the AGP with respect to the reciprocal mean gas pressure is due to the rarefaction effects, where infrequent collisions between gas molecules not only cause the gas slippage at the solid surface, but also modify the constitutive relation between the stress and strain rate (Struchtrup 2005), so that the NSEs have to be modified, or completely fail. The extent of rarefaction is characterized by the Knudsen number $Kn$ (i.e. the ratio of the mean free path $\unicode[STIX]{x1D706}$ of gas molecules to the characteristic flow length  $L$ ):

(1.3a,b ) $$\begin{eqnarray}Kn=\frac{\unicode[STIX]{x1D706}}{L},\quad \text{and}\quad \unicode[STIX]{x1D706}=\frac{\unicode[STIX]{x1D707}(T_{0})}{\bar{p}}\sqrt{\frac{\unicode[STIX]{x03C0}RT_{0}}{2}},\end{eqnarray}$$

where $\unicode[STIX]{x1D707}(T_{0})$ is the shear viscosity of the gas at a reference temperature $T_{0}$ and $R$ is the gas constant. Gas flows can be classified into four regimes (note that this partition of flow regime is roughly true for the gas flow between two parallel plates with a distance $L$ ; for gas flows in porous media, the region of $Kn$ for different flow regimes may change): continuum flow $(Kn\lesssim 0.001)$ in which NSEs are applicable; slip flow $(0.001<Kn\lesssim 0.1)$ where NSEs with appropriate velocity-slip/temperature-jump boundary conditions may be used; transition flow $(0.1<Kn\lesssim 10)$ and free-molecular flow $(Kn>10)$ , where NSEs break down and gas kinetic equations such as the fundamental Boltzmann equation are used to describe rarefied gas flows (Chapman & Cowling 1970).

Through systematic experimental measurements, Klinkenberg found that the linear relation (1.2) between the AGP and reciprocal mean gas pressure (or equivalently, the Knudsen number) is only an approximation, as it became ‘evident from tables 5 to 7, the value of constant $b$ increases with increasing pressure’, that is, $b$ decreases when $Kn$ increases. In a latest experiment, a stronger variation of $b$ with the pressure is observed, see figure 7 in Moghaddam & Jamiolahmady (2016).

Theoretically, by numerically solving the gas kinetic equation on the exact samples used by Klinkenberg, the nonlinear dependence between the AGP and Knudsen number can be recovered in the whole range of gas pressure. However, due to the complexity of numerical simulations, the porous medium is often simplified by a single straight cylindrical tube (or a bundle of straight tubes with the same radius), as initially suggested by Klinkenberg (1941). Then, using the Maxwell’s diffuse-specular boundary condition for the gas–surface interaction (Maxwell 1879), the AGP of a straight tube can be obtained numerically; the numerical solution is often fitted analytically, for example for the diffuse boundary condition, by the following expression (Beskok & Karniadakis 1999; Civan 2010):

(1.4) $$\begin{eqnarray}\frac{k_{a}}{k_{\infty }}=\left[1+\frac{128}{15\unicode[STIX]{x03C0}^{2}}\tan ^{-1}\left(4Kn^{0.4}\right)Kn\right]\left(1+\frac{4Kn}{1+4Kn}\right),\end{eqnarray}$$

to predict the unconventional gas production.

So far, although many attempts have been made to predict the AGP, none of them were able to explain the simple experiment by Klinkenberg 76 years ago, in a wide range of gas pressures. Figure 1 shows the numerical solutions of the Bhatnagar–Gross–Krook (BGK) kinetic equation (see § 2.1), as well as that of the NSEs with velocity-slip boundary conditions. It is clear that both the BGK equation and the NSEs with a first-order velocity-slip boundary condition (FVBC) cannot explain Klinkenberg’s experimental results that ‘the correction factor decreases when the Knudsen number increases’. In fact, NSEs give a constant correction factor, while the correction factor from the BGK equation increases with the Knudsen number. The NSEs with second-order velocity-slip boundary condition (Karniadakis, Beskok & Aluru 2005; Moghaddam & Jamiolahmady 2016) seem to produce the same trend as observed in Klinkenberg’s experiments at small $Kn$ , but they quickly predict a negative AGP as $Kn$ increases, which is definitely incorrect. In addition, in the experiment the AGP could increase approximately 30 times, but the use of second-order velocity-slip boundary condition predicts a maximum AGP that is only 1.5 times the intrinsic permeability.

Figure 1. The AGP versus the Knudsen number for a gas flow in a straight cylindrical tube, obtained from the numerical simulation of the Bhatnagar–Gross–Krook equation, where $\unicode[STIX]{x1D6FC}$ is the tangential momentum accommodation coefficient in the kinetic boundary condition for the gas–surface interaction (see § 2.1). Analytical solutions of the NSEs with first- and second-order velocity-slip boundary conditions are also shown, where the viscous velocity-slip coefficients are obtained from Loyalka, Petrellis & Storvick (1975) and Gibelli (2011).

The accuracy of gas kinetic equations in modelling rarefied gas flows has been tested extensively in various fields from aerodynamics to microfluidics. We believe that one of the reasons for the contradiction between the solution of the BGK equation and the experimental results is the use of oversimplified straight tubes, where the streamlines are straight instead of tortuous as in the experiment. The main objective of the present work is to conduct numerical simulations on simplified porous media (which are much simpler than the real samples to make the calculation tractable but more complicated than straight tubes) to pinpoint factors that lead to the same relation (trend) between the correction factor and Knudsen number, as initially observed by Klinkenberg (1941) and confirmed in many other experiments.

It should be noted that in a recent work, based on the NSEs with FVBC, Lasseux, Valdes Parada & Porter (2016) found that the AGP of the porous medium is not only a nonlinear function of $Kn$ , but also has a similar trend as observed in Klinkenberg (1941). This result, however, is highly questionable, because the NSEs were used beyond their validity. Through our theoretical analysis and numerical calculations, we show that the NSEs with FVBC can only predict the AGP of the porous medium to the first-order accuracy of $Kn$ . This forms the second objective of this paper, as the NSEs with FVBC are widely misused to predict the AGP of porous media.

Figure 2. A two-dimensional porous medium consisting of a periodic array of discs. A, B, C and D are the four corners of the unit rectangular cell (computational domain used below). The length of the side AB is $L$ , while that of AD is $L/2$ . Other complex porous media can be generated by adding random solids into the rectangle ABCD, and applying the periodic boundary condition at sides AD and BC but the symmetrical boundary condition along sides AB and CD, respectively.

The remainder of the paper is organized as follows: in § 2 we introduce the mesoscopic BGK kinetic model and the gas kinetic boundary condition, as well as the macroscopic regularized 20-moment (R20) equations to describe rarefied gas flows in porous media. We also analyse the limitation of the NSEs and FVBC. In § 3, numerical simulations are performed to analyse the variation of the AGP with the Knudsen number and to assess the validation of the NSEs with FVBC. The influence of the gas–surface interaction on the AGP is analysed, and possible factors that lead to the observation of Klinkenberg (1941) are identified in § 3.3. The paper closes with some final comments in § 4.

2 State of the problem

Consider a gas flowing through a periodic porous medium. Suppose the geometry along the $x_{3}$ -direction is uniform and infinite, the gas flow is effectively two-dimensional and can be studied in a unit rectangular cell ABCD, with appropriate governing equations and boundary conditions. One example of the porous medium consisting of a periodic array of discs is shown in figure 2; other more complex porous media can be generated by adding more solids randomly in the unit rectangular cell ABCD. We are interested in how the AGP varies with the Knudsen number, properties of porous media, and gas–surface interaction.

2.1 The mesoscopic description: the gas kinetic theory

The Boltzmann equation, which is fundamental in the study of rarefied gas dynamics from the continuum to the free-molecular flow regimes, uses the distribution function $f(t,\boldsymbol{x},\boldsymbol{v})$ to describe the system state:

(2.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}t}+v_{1}\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}x_{1}}+v_{2}\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}x_{2}}+v_{3}\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}x_{3}}={\mathcal{C}}(f),\end{eqnarray}$$

where $\boldsymbol{v}=(v_{1},v_{2},v_{3})$ is the three-dimensional molecular velocity normalized by the most probable speed $v_{m}=\sqrt{2RT_{0}}$ , $\boldsymbol{x}=(x_{1},x_{2},x_{3})$ is the spatial coordinate normalized by the length $L$ of the side AB, $t$ is the time normalized by $L/v_{m}$ , $f$ is normalized by $\bar{p}/v_{m}^{3}mRT_{0}$ with $m$ being the molecular mass, while ${\mathcal{C}}$ is the Boltzmann collision operator. In order to save the computational cost, ${\mathcal{C}}$ is usually replaced by the relaxation-time approximation (Bhatnagar, Gross & Krook 1954), resulting in the BGK equation.

When the porous medium is so long that the pressure gradient is small, namely, $|L\,\text{d}p/p\,\text{d}x_{1}|\ll 1$ with $p$ being the local gas pressure, the BGK equation can be linearized. The distribution function is expressed as $f=f_{eq}(1+h)$ , where the equilibrium distribution function is defined as

(2.2) $$\begin{eqnarray}f_{eq}=\frac{\exp (-|\boldsymbol{v}|^{2})}{\unicode[STIX]{x03C0}^{3/2}},\end{eqnarray}$$

and the perturbation $h$ is governed by:

(2.3) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}+v_{1}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x_{1}}+v_{2}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x_{2}}=\frac{\sqrt{\unicode[STIX]{x03C0}}}{2Kn}\left[\unicode[STIX]{x1D71A}+2u_{1}v_{1}+2u_{2}v_{2}+\unicode[STIX]{x1D70F}\left(|\boldsymbol{v}|^{2}-\frac{3}{2}\right)-h\right],\end{eqnarray}$$

with macroscopic quantities such as the perturbed number density $\unicode[STIX]{x1D71A}$ of gas molecules, the velocities $u_{1}$ and $u_{2}$ and the perturbed temperature $\unicode[STIX]{x1D70F}$ being calculated as

(2.4a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D71A}=\int hf_{eq}\,\text{d}\boldsymbol{v},\quad (u_{1},u_{2})=\int (v_{1},v_{2})hf_{eq}\,\text{d}\boldsymbol{v},\quad \unicode[STIX]{x1D70F}=\frac{2}{3}\int |\boldsymbol{v}|^{2}hf_{eq}\,\text{d}\boldsymbol{v}-\unicode[STIX]{x1D71A}.\end{eqnarray}$$

The kinetic equation (2.3) has to be supplied with the boundary condition. Suppose the pressure gradient is along the $x_{1}$ direction, on the inlet and outlet of the computational domain ABCD (the coordinates of the four corners A, B, C and D are $(-0.5,0),(0.5,0),(0.5,0.5)$ and $(-0.5,0.5)$ , respectively), the pressure gradient is applied and periodic condition for the flow velocity is used (Sharipov & Graur 2012):

(2.5) $$\begin{eqnarray}h(\mp 0.5,x_{2},v_{1},v_{2},v_{3})=\pm 1+h(\pm 0.5,x_{2},v_{1},v_{2},v_{3}),\quad \text{when }v_{1}\gtrless 0,\end{eqnarray}$$

at the lines AB and CD, the specular reflection boundary condition is used to account for the spatial symmetry:

(2.6) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}h\left(x_{1},0,v_{1},v_{2},v_{3}\right)=h\left(x_{1},0,v_{1},-v_{2},v_{3}\right),\quad \text{when }v_{2}>0,\\ h\left(x_{1},0.5,v_{1},v_{2},v_{3}\right)=h\left(x_{1},0.5,v_{1},-v_{2},v_{3}\right),\quad \text{when }v_{2}<0,\end{array}\right\} & & \displaystyle\end{eqnarray}$$

while at the solid surface, the diffuse-specular boundary condition is used (Maxwell 1879):

(2.7) $$\begin{eqnarray}h(\boldsymbol{v})=\frac{2\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x03C0}}\int _{v_{n}^{\prime }<0}|v_{n}^{\prime }|h(\boldsymbol{v}^{\prime })\exp (-|\boldsymbol{v}^{\prime }|^{2})\,\text{d}\boldsymbol{v}^{\prime }+(1-\unicode[STIX]{x1D6FC})h(\boldsymbol{v}-2\boldsymbol{n}v_{n}),\end{eqnarray}$$

where $\boldsymbol{n}$ is the outer normal vector of the solid surface and $\unicode[STIX]{x1D6FC}$ is the tangential momentum accommodation coefficient (TMAC). This boundary condition assumes that, after collision with the surface, a gas molecule is specularly reflected with probability $1-\unicode[STIX]{x1D6FC}$ , otherwise it is reflected diffusely (i.e. reflected towards every direction with equal probability, in a Maxwellian velocity distribution). Purely diffuse or specular reflections take place for $\unicode[STIX]{x1D6FC}=1$ or $\unicode[STIX]{x1D6FC}=0$ , respectively. Note that most of the current studies, including the widely used equation (1.4), use the diffuse boundary condition.

The AGP, which is normalized by $L^{2}$ , is calculated by

(2.8) $$\begin{eqnarray}k_{a}=2\sqrt{\frac{1}{\unicode[STIX]{x03C0}}}KnG_{P},\end{eqnarray}$$

where $G_{p}=2\int _{0}^{1/2}u_{1}(x_{2})\,\text{d}x_{2}$ is the dimensionless mass flow rate.

2.2 The macroscopic description: NSEs and moment equations

Historically, the state of a gas is first described by macroscopic quantities such as the density $\unicode[STIX]{x1D70C}$ , velocity $u_{i}$ and temperature $T$ ; and its dynamics is described by the Euler equations or the NSEs (based on the empirical Newton’s law for stress and Fourier’s law for heat flux). These equations, however, can be derived rigorously from the Boltzmann equation, at various orders of approximations.

By taking the velocity moments of the Boltzmann equation (2.1), the five macroscopic quantities are governed by the following equations:

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}u_{i}}{\unicode[STIX]{x2202}x_{i}}=0, & \displaystyle\end{eqnarray}$$

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}u_{i}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}u_{i}u_{j}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x2202}x_{j}}=-\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{i}}, & \displaystyle\end{eqnarray}$$

(2.11) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}T}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}u_{i}T}{\unicode[STIX]{x2202}x_{i}}+\frac{2}{3R}\frac{\unicode[STIX]{x2202}q_{i}}{\unicode[STIX]{x2202}x_{i}}=-\frac{2}{3R}\left(p\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{i}}+\unicode[STIX]{x1D70E}_{ij}\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{i}}\right). & \displaystyle\end{eqnarray}$$

However, the above equations are not closed, since expressions for the shear stress $\unicode[STIX]{x1D70E}_{ij}$ and heat flux $q_{i}$ are not known. One way to close (2.9)–(2.11) is to use the Chapman–Enskog expansion, where the distribution function is expressed in the power series of $Kn$ (Chapman & Cowling 1970):

(2.12) $$\begin{eqnarray}f=f^{(0)}+Knf^{(1)}+Kn^{2}f^{(2)}+\cdots \,,\end{eqnarray}$$

where $f^{(0)}$ is the equilibrium Maxwellian distribution function. When $f=f^{(0)}$ , we have $\unicode[STIX]{x1D70E}_{ij}=0$ and $q_{i}=0$ , and (2.9)–(2.11) reduce to the Euler equations. When the distribution function is truncated at the first order of $Kn$ , that is,

(2.13) $$\begin{eqnarray}f=f^{(0)}+Knf^{(1)},\end{eqnarray}$$

we have

(2.14a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D70E}_{ij}=-2\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}u_{{<}i}}{\unicode[STIX]{x2202}x_{j>}}\quad \text{and}\quad q_{i}=-\frac{15}{4}R\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}x_{i}},\end{eqnarray}$$

where the angular brackets are used to denote the traceless part of a symmetric tensor, so that (2.9)–(2.11) reduce to the NSEs. When $f=f^{(0)}+Knf^{(1)}+Kn^{2}f^{(2)}$ , the Burnett equations can be derived but are not used nowadays due to their intrinsic instabilities (Garcia-Colin, Velasco & Uribe 2008).

Alternatively, by combining the moment method of Grad (1949) and the Chapman–Enskog expansion, the regularized 13-, 20- and 26-moment equations (Struchtrup & Torrihon 2003; Gu & Emerson 2009) can be derived from the Boltzmann equation to describe rarefied gas flows at different levels of rarefaction (Torrilhon 2016). Here the R20 equations are used, which, in addition to (2.9)–(2.11), include governing equations for the high-order moments $\unicode[STIX]{x1D70E}_{ij}$ , $q_{i}$ , and $m_{ijk}$ :

(2.15) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}u_{k}\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x2202}x_{k}}+\frac{\unicode[STIX]{x2202}m_{ijk}}{\unicode[STIX]{x2202}x_{k}}=-\frac{p}{\unicode[STIX]{x1D707}}\unicode[STIX]{x1D70E}_{ij}-2p\frac{\unicode[STIX]{x2202}u_{{<}i}}{\unicode[STIX]{x2202}x_{j>}}-\frac{4}{5}\frac{\unicode[STIX]{x2202}q_{{<}i}}{\unicode[STIX]{x2202}x_{j>}}-2\unicode[STIX]{x1D70E}_{k<i}\frac{\unicode[STIX]{x2202}u_{j>}}{\unicode[STIX]{x2202}x_{k}},\end{eqnarray}$$

(2.16) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x2202}q_{i}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}u_{j}q_{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{1}{2}\frac{\unicode[STIX]{x2202}R_{ij}}{\unicode[STIX]{x2202}x_{j}}=-\frac{2}{3}\frac{p}{\unicode[STIX]{x1D707}}q_{i}-\frac{5}{2}p\frac{\unicode[STIX]{x2202}RT}{\unicode[STIX]{x2202}x_{i}}+\frac{\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x1D70C}}\left(\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{jk}}{\unicode[STIX]{x2202}x_{k}}\right)-RT\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x2202}x_{j}}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\frac{7}{2}\unicode[STIX]{x1D70E}_{ij}\frac{\unicode[STIX]{x2202}RT}{\unicode[STIX]{x2202}x_{j}}-\left(\frac{2}{5}q_{i}\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{j}}+\frac{7}{5}q_{j}\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{2}{5}q_{j}\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{i}}\right)-m_{ijk}\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{k}}-\frac{1}{6}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E5}}{\unicode[STIX]{x2202}x_{i}},\end{eqnarray}$$

(2.17) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}m_{ijk}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}u_{l}m_{ijk}}{\unicode[STIX]{x2202}x_{l}}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{ijkl}}{\unicode[STIX]{x2202}x_{l}} & = & \displaystyle -\frac{3}{2}\frac{p}{\unicode[STIX]{x1D707}}m_{ijk}-3\frac{\unicode[STIX]{x2202}RT\unicode[STIX]{x1D70E}_{{<}ij}}{\unicode[STIX]{x2202}x_{k>}}-\frac{12}{5}q_{{<}i}\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{k>}}\nonumber\\ \displaystyle & & \displaystyle +\,3\frac{\unicode[STIX]{x1D70E}_{{<}ij}}{\unicode[STIX]{x1D70C}}\left(\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{k>}}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{k>l}}{\unicode[STIX]{x2202}x_{l}}\right).\end{eqnarray}$$

The constitutive relationships between the unknown higher-order moments ( $R_{ij}$ , $\unicode[STIX]{x1D6E5}$ and $\unicode[STIX]{x1D719}_{ijkl}$ ) and lower-order moments were given by Structrup & Torrilhon (2003) and Gu & Emerson (2009) to close (2.9) to (2.17). For slow flows in porous media, it is adequate to use the gradient transport terms only (Taheri et al. 2009; Gu, Emerson & Tang 2010) and they are:

(2.18a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{ijkl}=-\frac{4\unicode[STIX]{x1D707}}{C_{1}\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}m_{{<}ijk}}{\unicode[STIX]{x2202}x_{l>}},\quad R_{ij}=-\frac{24}{5}\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}q_{{<}i}}{\unicode[STIX]{x2202}x_{j>}},\quad \unicode[STIX]{x1D6E5}=-12\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}q_{k}}{\unicode[STIX]{x2202}x_{k}},\end{eqnarray}$$

where the collision constant $C_{1}$ is 2.097 for Maxwell molecules (Gu & Emerson 2009).

Macroscopic wall boundary conditions were obtained from the diffuse-specular boundary condition (Maxwell 1879). In a frame where the coordinates are attached to the wall, with $n_{i}$ the normal vector of the wall pointing towards the gas and $\unicode[STIX]{x1D70F}_{i}$ the tangential vector of the wall, the velocity-slip ( $u_{\unicode[STIX]{x1D70F}}$ , parallel to the wall) and temperature-jump conditions are:

(2.19) $$\begin{eqnarray}\displaystyle & \displaystyle u_{\unicode[STIX]{x1D70F}}=-\frac{2-\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6FC}}\sqrt{\frac{\unicode[STIX]{x03C0}RT}{2}}\frac{\unicode[STIX]{x1D70E}_{n\unicode[STIX]{x1D70F}}}{p_{\unicode[STIX]{x1D6FC}}}-\frac{5m_{nn\unicode[STIX]{x1D70F}}+2q_{\unicode[STIX]{x1D70F}}}{10p_{\unicode[STIX]{x1D6FC}}},\quad & \displaystyle\end{eqnarray}$$

(2.20) $$\begin{eqnarray}\displaystyle & \displaystyle RT-RT_{w}=-\frac{2-\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6FC}}\sqrt{\frac{\unicode[STIX]{x03C0}RT}{2}}\frac{q_{n}}{2p_{\unicode[STIX]{x1D6FC}}}-\frac{RT\unicode[STIX]{x1D70E}_{nn}}{4p_{\unicode[STIX]{x1D6FC}}}+\frac{u_{\unicode[STIX]{x1D70F}}^{2}}{4}-\frac{75R_{nn}+28\unicode[STIX]{x1D6E5}}{840p_{\unicode[STIX]{x1D6FC}}}+\frac{\unicode[STIX]{x1D719}_{nnnn}}{24p_{\unicode[STIX]{x1D6FC}}},\quad & \displaystyle\end{eqnarray}$$

where $p_{\unicode[STIX]{x1D6FC}}=p+\unicode[STIX]{x1D70E}_{nn}/2-(30R_{nn}+7\unicode[STIX]{x1D6E5})/840RT-\unicode[STIX]{x1D719}_{nnnn}/24RT$ . The rest of wall boundary conditions for higher-order moments are listed in appendix A. Note that the velocity-slip boundary condition (2.19) is also of higher order due to the appearance of the higher-order moment $m_{nn\unicode[STIX]{x1D70F}}$ .

Following the above introduction, it is clear that the NSEs with FVBC are only accurate to the first order of $Kn$ ; therefore, any AGP predicted from the NSEs and showing the nonlinear dependence with $Kn$ is highly questionable. The R20 equations are accurate to the third order of $Kn$ (Struchtrup 2005), which should give the same AGP of the porous medium as the NSEs when $Kn\rightarrow 0$ , and be more accurate than the NSEs as $Kn$ increases. Numerical simulations are also performed in the following section to demonstrate this.

3 Numerical methods and results

In this section, we first present analytical/numerical solutions of the NSEs with FVBC, the linearized BGK equation, the R20 equation as well as the direct simulation Monte Carlo (DSMC) method (Borner, Panerai & Mansour 2017), for rarefied gas flows through porous media, to show the accuracy of the NSEs with FVBC. Then we give a possible explanation to understand Klinkenberg’s experimental results.

3.1 The accuracy of the NSEs and FVBC

We first investigate the rarefied gas flow through a periodic array of discs with diameter $D$ , as shown in figure 2. Using the NSEs and FVBC, when the porosity $\unicode[STIX]{x1D716}=1-\unicode[STIX]{x03C0}D^{2}/4L^{2}$ is large, the AGP can be obtained analytically (Chai et al. 2011):

(3.1) $$\begin{eqnarray}k_{a}=\frac{1}{8\unicode[STIX]{x03C0}\left(1+2\unicode[STIX]{x1D709}Kn\sqrt{{\displaystyle \frac{\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D719}}}}\right)}\bigg[-\ln \unicode[STIX]{x1D719}-\frac{3}{2}+2\unicode[STIX]{x1D719}-\frac{\unicode[STIX]{x1D719}^{2}}{2}+2\unicode[STIX]{x1D709}Kn\sqrt{\frac{\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D719}}}\left(-\ln \unicode[STIX]{x1D719}-\frac{1}{2}+\frac{\unicode[STIX]{x1D719}^{2}}{2}\right)\bigg],\end{eqnarray}$$

where $\unicode[STIX]{x1D719}=1-\unicode[STIX]{x1D716}$ is the solid fraction and $\unicode[STIX]{x1D709}$ is related to the viscous velocity-slip coefficient. When Maxwell’s diffuse boundary condition (2.7) is used (Hadjiconstantinou 2003), we have $\unicode[STIX]{x1D709}=1.016\sqrt{4/\unicode[STIX]{x03C0}}=1.15$ . The viscous velocity-slip coefficient for other values of TMAC can be found in Loyalka et al. (1975) and Gibelli (2011). The intrinsic permeability $k_{\infty }$ is obtained when $Kn=0$ .

For the linearized BGK equation (2.3), two reduced distribution functions were introduced to cast the three-dimensional molecular velocity space into a two-dimensional one, and the obtained two equations are solved numerically by the discrete velocity method (Chu 1965; Aoki, Sone & Yamada 1990) and the unified gas kinetic scheme (Huang, Xu & Yu 2013). In the unified gas kinetic scheme, a body-fitted structured curvilinear mesh is used, with 150 lines along the radial direction and 300 lines along the circumferential direction, see figure 3(a). In the discrete velocity method, a Cartesian grid with $801\times 401$ equally spaced points is used and the solid surface is approximated by the ‘staircase’ mesh. In solving the R20 equations, a body-fitted mesh with $201\times 101$ cells is used, and the detailed numerical method is given by Gu & Emerson (2009). The molecular velocity space in the BGK equation is also discretized: $v_{1}$ and $v_{2}$ are approximated by the $8\times 8$ Gauss–Hermite quadrature when $Kn$ is small ( $Kn<0.01$ in this case), and the Newton–Cotes quadrature with $22\times 22$ non-uniform discrete velocity points when $Kn$ is large (Wu, Reese & Zhang 2014).

The results for a periodic porous medium with $\unicode[STIX]{x1D716}=0.8$ are shown in figure 3. It is seen that both the discrete velocity method on the ‘staircase’ geometry and the unified gas kinetic scheme on the curvilinear mesh yield the same results. Also, although we use the linearized BGK model, our numerical results are quite close to those from the DSMC simulation (Borner et al. 2017). Note that in the DSMC simulation, the Knudsen number and AGP are normalized by the radius and radius square of the disc, respectively, so the data should be renormalized. Also, the mean free path defined in (1.3a,b ) is $15\unicode[STIX]{x03C0}/2(7-2\unicode[STIX]{x1D714})(5-2\unicode[STIX]{x1D714})$ times larger than that used in DSMC based on the variable hard sphere model, where $\unicode[STIX]{x1D714}=0.81$ for argon (Wu et al. 2013); so the Knudsen number should be rescaled.

Figure 3. (a) The body-fitted mesh used in the unified gas kinetic scheme, when the porosity of the porous medium in figure 2 is $\unicode[STIX]{x1D716}=0.8$ . For clarity, only $100\times 50$ cells are shown. (b) The AGP as a function of the Knudsen number when the diffuse boundary condition is used, i.e. $\unicode[STIX]{x1D6FC}=1$ in (2.7). The solid line and dots are numerical results of the linearized BGK solved by the discrete velocity method and the unified gas kinetic scheme, respectively. The dash-dotted and dashed lines are analytical solutions of the NSEs (3.1), with the no-slip and first-order velocity-slip boundary conditions, respectively, while the dotted line is the slip-corrected permeability obtained by expanding the analytical solution (3.1) to the first order of $Kn$ , see (3.2). The DSMC results are obtained from the recent simulation by Borner et al. (2017).

The accuracy of the permeability (3.1) is assessed by comparing to numerical solutions of the BGK and R20 equations. When $Kn\lesssim 0.15$ , our numerical simulations based on the linearized BGK equation and R20 equations agree with each other and the AGP is a linear function of $Kn$ . When $Kn\gtrsim 0.2$ , the R20 equations, although being accurate to the third order of $Kn$ , predict lower AGP than that of the BGK equation. For larger $Kn$ beyond the validity of the R20 equations and their boundary conditions, no converged numerical solutions can be obtained. The analytical permeability (3.1) increases linearly with $Kn$ only when $Kn\lesssim 0.02$ , and then quickly reaches to a maximum value when $Kn\gtrsim 0.2$ . This comparison clearly demonstrates that, the NSEs with FVBC are only accurate to the first order of $Kn$ . This result is in accordance with the approximation (2.13) adopted in the derivation of the NSEs from the Boltzmann equation. Although the ‘curvature of the solid–gas interface’, claimed by Lasseux et al. (2016), makes the AGP a concave function of $Kn$ in the framework of the NSEs with FVBC, higher-order moments in (2.15)–(2.17) and the higher-order velocity slip in (2.19), which are derived from the Boltzmann equation and the gas kinetic boundary condition (2.7) to the third-order accuracy of $Kn$ , restore linear dependence of the AGP on the Knudsen number when $Kn\lesssim 0.15$ .

Since theoretically the NSEs with FVBC are only accurate up to the first order of $Kn$ , we expand (3.1), a nonlinear function of $Kn$ , into a Taylor series of $Kn$ around $Kn=0$ , and retain only the zeroth- and first-order terms; the obtained permeability is called the slip-corrected AGP, and it reads

(3.2) $$\begin{eqnarray}k_{a,slip}=\frac{1}{8\unicode[STIX]{x03C0}}\left[-\ln \unicode[STIX]{x1D719}-\frac{3}{2}+2\unicode[STIX]{x1D719}-\frac{\unicode[STIX]{x1D719}^{2}}{2}+2\unicode[STIX]{x1D709}Kn\sqrt{\frac{\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D719}}}\left(1-2\unicode[STIX]{x1D719}+\unicode[STIX]{x1D719}^{2}\right)\right].\end{eqnarray}$$

The slip-corrected AGP is plotted as the dotted line in figure 3(b). Surprisingly, this Taylor expansion, which is valid at small $Kn$ , gives a good agreement with the numerical solution of the linearized BGK equation when $Kn\lesssim 0.35$ ; even in the transition flow regime where the Knudsen number is as high as one, this slip-corrected expression predicts an AGP only approximately 15 % smaller than the numerical result of the linearized BGK equation. We emphasis here that, however, this is only a coincidence; and we will explain this point in § 3.3 below.

Figure 4. (a) The geometry in a unit computational cell, when the porosity is $\unicode[STIX]{x1D716}=0.6$ . Solids of random size and position are shown in black. The periodic porous medium is generated by placing the whole computational domain inside the unit rectangular cell ABCD in figure 2. $3000\times 1500$ cells are used to discretize the spatial domain. (b) The ratio of the AGP to the intrinsic permeability $k_{\infty }=9.37\times 10^{-6}$ as a function of the Knudsen number. The solid and dashed lines are numerical results of the linearized BGK equation and the NSEs with FVBC, respectively. The dotted lines are the slip-corrected permeability, obtained by fitting the numerical solution of the NSEs with FVBC at small Knudsen numbers to the first order of $Kn$ . Note that here the effective Knudsen number $Kn^{\ast }$ defined in (3.4) is 73 times $Kn$ .

The conclusion that the NSEs with FVBC are accurate only to the first order of $Kn$ not only holds for the simple porous medium in figure 2, but also applies to more complex porous media, for example, see the unit cell in figure 4(a) where the porosity is 0.6. In this case, the linearized BGK equation is solved by the discrete velocity method, with a Cartesian mesh of $3000\times 1500$ cells; the grid convergence is verified, as using $6000\times 3000$ cells results in only 0.6 % increase of the AGP when $Kn=0.0001$ . The NSEs with FVBC are solved in OpenFOAM using the SIMPLE algorithm and a cell-centred finite-volume discretization scheme, on unstructured grids. A body-fitted computational grid is generated using the OpenFOAM meshing tool, resulting in a mesh of approximately 600 000 cells of which the majority are hexahedra and the rest few close to the walls are prisms. Note that the results of the R20 equation and the unified gas kinetic scheme are not available because currently the corresponding codes are not yet developed to deal with such a complex geometry.

Due to the small clearance between discs, the characteristic flow length is much smaller than the size of the computational domain $L$ . Thus, the effective characteristic flow length $L^{\ast }$ is often used instead, which is defined as (Lasseux et al. 2016):

(3.3) $$\begin{eqnarray}L^{\ast }=L\sqrt{\frac{12k_{\infty }}{\unicode[STIX]{x1D716}}},\end{eqnarray}$$

and the effective Knudsen number is

(3.4) $$\begin{eqnarray}Kn^{\ast }=\frac{\unicode[STIX]{x1D706}}{L^{\ast }}=Kn\sqrt{\frac{\unicode[STIX]{x1D716}}{12k_{\infty }}}.\end{eqnarray}$$

From figure 4(b) we see that the AGP obtained from the NSEs with FVBC increases linearly when $Kn\lesssim 0.002$ (or $Kn^{\ast }\lesssim 0.146$ ) and then quickly reaches a constant value. Again, the comparison with solution of the linearized BGK equation shows that the NSEs with FVBC are approximately accurate when the AGP is a linear function of $Kn$ ; in this region of $Kn$ , the maximum AGP is only approximately 1.5 times larger than the intrinsic permeability $k_{\infty }$ . Interestingly, like (3.2), when the numerical solution of the NSEs with FVBC is ‘filtered’ to keep only the zeroth- and first-order terms of $Kn$ , the resulting solution is only smaller than the numerical solution of the linearized BGK equation by approximately 15 % in a large range of gas rarefaction, i.e. $Kn^{\ast }\lesssim 73$ , where the AGP is larger than the intrinsic permeability by hundreds of times.

3.2 The computational time

It will be helpful to mention the computational time for the two examples in § 3.1. For the discrete velocity method, numerical simulations are performed on a Dell workstation (Precision Tower 7910) with Dual processor Intel Xeon CPU E5-2630 v3 2.40GHz. There are 32 cores in total but only 8 cores are used. For simplicity, the $8\times 8$ Gauss–Hermite quadrature is used here for all the Knudsen numbers shown in table 1. For the non-uniform velocity grids, the iteration numbers are approximately the same as the Gauss–Hermite quadrature, but the computational time is approximately 10 times higher due to the large number of discrete velocity points. The AGP is calculated every 1000 iteration steps, and our Fortran programme is terminated when the following convergence criterion

(3.5) $$\begin{eqnarray}\frac{k_{a}(i)-k_{a}(i-1000)}{k_{a}(i)}<10^{-10},\end{eqnarray}$$

is reached, where $i$ is the iteration step.

Table 1. Iteration steps and elapsed time when the linearized BGK equation is solved by the implicit discrete velocity method. The time is measured by wall clock time.

Generally speaking, when the implicit scheme is used to solve the gas kinetic equations, iteration steps increase when the Knudsen number decreases. From table 1 we find that this is true except for the complex geometry where the iteration steps at $Kn=0.1$ are less than that of $Kn=1$ . This is because we start the numerical simulation from $Kn=1$ ; and when the solution is converged, it is used as the initial guess in the computation of $Kn=0.1$ , and so on. Thus, it is possible that the solution for smaller $Kn$ requires fewer iteration steps, especially in complicated porous media. Note that for the circular cylinder case we need many iteration steps at $Kn=0.001$ ; this is because the flow is in the continuum regime where the collision dominates so that the exchange of information due to streaming is extremely slow and consequently the convergence is slow. In practical calculations, however, we do not have to calculate the AGP at such a small $Kn$ , because it is only approximately 2 % larger than the intrinsic permeability; we only have to calculate the intrinsic permeability by the Navier–Stokes solvers such as the OpenFOAM and the multiple relaxation-time lattice Boltzmann method for simulating flows in complex porous media (Pan, Luo & Miller 2006).

The unified gas kinetic scheme solves the linearized BGK equation explicitly, and hence the computational time is high. For example, the iteration steps are 350 000 when $Kn=1$ . The implicit unified gas kinetic scheme can speed up the convergence by hundreds of times (Zhu, Zhong & Xu 2016), but the code for porous media flow simulations is still under development. This scheme has the unique advantage that the spatial resolution can be coarser than that of the discrete velocity method to reach the same order of accuracy, so that the computational memory usage can be significantly reduced.

The OpenFOAM solver running on the same Dell workstation (using one core) takes approximately an hour to obtain the converged solution for the complex porous medium case in figure 4, at each Knudsen number.

The computational efforts to solve the moment equations are opposite to that for the BGK kinetic equation. When $Kn$ is small, a solution can be obtained quickly with large iteration relaxation factors (Tang et al. 2013). As the Knudsen number increases, the values of the relaxation factors have to be reduced and more iterations are required to reach to a converged solution. In the present study in figure 3, single processor of Intel Core i7-4770 CPU at 3.4 GHz is used and the computing times range from minutes to hours depending on the Knudsen number. Computer codes dealing with the complex geometry are still under development. We intend to couple deterministic solvers for moment equations and the linearized BGK equation, so that higher numerical accuracy and efficiency can be achieved simultaneously.

3.3 Possible explanation of the Klinkenberg’s experiment

From the numerical simulations in § 3.1, Klinkenberg’s finding that the correction factor decreases when $Kn$ increases is not observed, even when the porous media are more complicated than the straight cylindrical tube. Instead, we find that the correction factor increases with $Kn$ , which is exactly the same as the rarefied gas flowing through a straight cylindrical tube, see figure 1.

From the analytical solution (3.2), we find that, if we use a non-unitary TMAC (say, $\unicode[STIX]{x1D6FC}=0.5$ ), then $\unicode[STIX]{x1D709}\approx 3.23$ (Loyalka et al. 1975) and the correction factor is approximately three times larger than that of $\unicode[STIX]{x1D6FC}=1$ when $Kn$ is small. If at large Knudsen numbers the AGP of $\unicode[STIX]{x1D6FC}=0.5$ is only slightly larger than that of $\unicode[STIX]{x1D6FC}=1$ (unfortunately this is not the case for rarefied gas flows between two parallel plates or through a straight cylindrical tube), a decrease of the correction factor when $Kn$ increases may be observed in the numerical simulation of the linearized BGK equation in porous media.

To study the influence of the non-unitary TMAC on the AGP, we investigate the rarefied gas passing through an infinite square array of square cylinders, by replacing the discs in figure 2 with squares. The only reason to do this is that the specular reflection in the boundary condition (2.7) can be accurately implemented.

Since we use non-dimensional variables, we rewrite Klinkenberg’s famous equation (1.2) in the following form

(3.6) $$\begin{eqnarray}\frac{k_{a}}{k_{\infty }}=1+\frac{b}{\bar{p}}\equiv 1+b^{\prime }Kn^{\ast },\end{eqnarray}$$

and investigate the variation of the ‘equivalent correction factor’ $b^{\prime }$ with respect to the Knudsen number and TMAC. Note that $b^{\prime }$ is proportional to the correction factor $b$ and the proportionality constant is independent of the Knudsen number (reciprocal mean gas pressure) and TMAC.

Figure 5. (a,b) Streamlines in unit computational cells. (c,d) The ratio of the AGP to the intrinsic permeability and (e,f) the correction factor, as functions of $Kn^{\ast }$ . Note that the vertical axis in (e) and (f) is in the logarithmic scale. Triangles, squares and circles are numerical results of the linearized BGK equation (solved by the discrete velocity method), with the TMAC $\unicode[STIX]{x1D6FC}=1.0$ , 0.5 and 0.1, respectively. The black lines are the numerical solutions from the NSEs with FVBC, where only up to the first-order terms of $Kn$ are retained. In (a,c,e), the porosity is $\unicode[STIX]{x1D716}=0.4$ , $k_{\infty }=0.001$ and $Kn^{\ast }=5.69Kn$ , while in (b,d,f), $\unicode[STIX]{x1D716}=0.8$ , $k_{\infty }=0.018$ and $Kn^{\ast }=1.94Kn$ .

Our numerical solutions for the linearized BGK equation, solved by the discrete velocity method, are shown in figure 5. The slip-corrected permeability (obtained from the numerical solution of the NSEs with FVBC, but only keeping the linear dependence of the AGP with the Knudsen number at small $Kn$ ) is also shown for comparison. From figure 5(c,d) we see that, as the numerical results in § 3.1, when the TMAC is $\unicode[STIX]{x1D6FC}=1$ , the slip-corrected permeability agrees well with the numerical results of the linearized BGK equation, when the porosity is $\unicode[STIX]{x1D716}=0.4$ and $0.8$ . Also, it is interesting to note that, when $Kn^{\ast }\rightarrow 0$ , the correction factor $b^{\prime }$ is approximately 6 and 6.1 when the porosity is $\unicode[STIX]{x1D716}=0.4$ and 0.8, respectively. For Poiseuille flow between two parallel plates, the slip-corrected permeability is exactly

(3.7) $$\begin{eqnarray}\frac{k_{a}}{k_{\infty }}=1+6\unicode[STIX]{x1D709}Kn^{\ast },\end{eqnarray}$$

where $\unicode[STIX]{x1D709}=1.15$ for the diffuse boundary condition. This means that the correction factor is $b^{\prime }=6.9$ , which is close to the two values we obtained for complex porous media, when the gas–surface interaction is diffuse.

However, when the TMAC is small, the accuracy of the slip-corrected permeability is significantly reduced. Since at $\unicode[STIX]{x1D6FC}=0.5$ and 0.1, $\unicode[STIX]{x1D709}$ is respectively $3.23$ and 19.3, the correction factor $b^{\prime }$ from the NSEs with FVBC is approximately 17 and 100, respectively. From figure 5(c,d) we see that the slip-corrected permeability is only accurate when $Kn^{\ast }<0.1$ . At relative large Knudsen numbers (i.e. $Kn^{\ast }>0.5$ ), the accuracy of the slip-corrected permeability is greatly reduced. For instance, when the porosity is 0.8, the slip-corrected permeability overpredicts the AGP by 160 % and 670 %, when $\unicode[STIX]{x1D6FC}=0.5$ and 0.1, respectively; when the porosity is 0.4, the slip-corrected permeability overpredicts the AGP by 360 % when $\unicode[STIX]{x1D6FC}=0.1$ . These two examples clearly show that, when the TMAC deviates largely from one, the NSEs with FVBC should not be used to predict the AGP in porous media.

Now, we focus on the influence of the non-unitary TMAC on the correction factor. From figure 5(c,d) we see that, as expected, when $Kn$ is fixed, the AGP increases when TMAC decreases. However, at different $Kn$ the amount of increase is different, so the variation of $b^{\prime }$ with respect to $Kn$ is quite different among the three TMACs considered. When $\unicode[STIX]{x1D6FC}=1$ , it is seen from figure 5(e,f) that the correction factor $b^{\prime }$ increases with $Kn$ , which clearly contradicts Klinkenberg’s experimental results. However, when $\unicode[STIX]{x1D6FC}=0.5$ and 0.1, it is found that, when $Kn$ increases, $b^{\prime }$ first decreases and then increases. The reason that the correction factor increases with the Knudsen number when $Kn$ is large can be easily understood. It is well known that when $Kn\rightarrow \infty$ , the dimensionless mass flow rate $G_{p}$ in (2.8) increases with $Kn$ logarithmically for Poiseuille flow between two parallel plates (Takata & Funagane 2011) and approaches a constant for Poiseuille flow through the cylindrical tube. Similar behaviours for $G_{p}$ are observed for all the cases considered in this paper. According to (3.6), the equivalent correction factor $b^{\prime }$ at $Kn^{\ast }\rightarrow \infty$ can be calculated as follows:

(3.8) $$\begin{eqnarray}b^{\prime }=\frac{k_{a}/k_{\infty }-1}{Kn^{\ast }}=\frac{2KnG_{p}/\sqrt{\unicode[STIX]{x03C0}}k_{\infty }}{Kn^{\ast }}-\frac{1}{Kn^{\ast }}\equiv C-\frac{1}{Kn^{\ast }},\end{eqnarray}$$

where $C$ is proportional to $G_{p}$ . Depending on the structure of the porous medium, $C$ is either a constant or increases with $Kn$ when $Kn\rightarrow \infty$ . Therefore, the correction factor increases with $Kn$ . Also, since $G_{p}$ always increases when TMAC decreases, at large $Kn$ , the correction factor increases when TMAC decreases, see figure 5(e,f).

Figure 6. (a) The geometry in a unit computational cell, when the porosity is $\unicode[STIX]{x1D716}=0.84$ . Solids of random size and position are shown in black. The periodic porous medium is generated by placing the whole computational domain inside the unit rectangular cell ABCD in figure 2. $2000\times 1000$ cells are used to discretize the spatial domain. (b) The ratio of the AGP to the intrinsic permeability $k_{\infty }=4.77\times 10^{-5}$ and (c) the correction factor $b^{\prime }$ , as functions of $Kn^{\ast }$ . Triangles, squares and circles are numerical results of the linearized BGK equation, with the TMAC $\unicode[STIX]{x1D6FC}=1.0$ , 0.5 and 0.1, respectively. In the inset in (b), the solid, dashed, dotted lines are the slip-corrected permeability for $\unicode[STIX]{x1D6FC}=1.0$ , 0.5 and 0.1, respectively, where the slopes (correction factor $b^{\prime }$ ) are 6.45, 18.2 and 108.5, respectively. Note that $Kn^{\ast }=38.3Kn$ .

At small $Kn$ , when $\unicode[STIX]{x1D6FC}=0.5$ and 0.1, the correction factor decreases when $Kn$ increases. For instance, when the porosity is $\unicode[STIX]{x1D716}=0.8$ , the correction factor $b^{\prime }$ with $\unicode[STIX]{x1D6FC}=0.5$ drops from 17 when $Kn^{\ast }\rightarrow 0$ to the minimum value 9 at $Kn^{\ast }=0.4$ , while $b^{\prime }$ with $\unicode[STIX]{x1D6FC}=0.1$ drops from 100 at $Kn^{\ast }\rightarrow 0$ to the minimum value 13 at $Kn^{\ast }=0.6$ . When the porosity is $\unicode[STIX]{x1D716}=0.4$ , the correction factor $b^{\prime }$ with $\unicode[STIX]{x1D6FC}=0.5$ drops from 17 at $Kn^{\ast }\rightarrow 0$ to the minimum value 12 at $Kn^{\ast }=0.5$ , while $b^{\prime }$ with $\unicode[STIX]{x1D6FC}=0.1$ drops from 100 when $Kn^{\ast }\rightarrow 0$ to the minimum value 27 at $Kn^{\ast }=1.1$ . This behaviour is in qualitative agreement (note that a quantitative comparison between the experiment and our simulation is currently impossible because we do not have the detailed structure data of the porous media) with observations by Klinkenberg (1941) and Moghaddam & Jamiolahmady (2016).

Our findings also apply to complicated porous media. For example, we have investigated the rarefied gas flow through a porous medium consisting of squares of random size and position in figure 6(a). A similar trend in the variation of the correction factor with respect to the Knudsen number is observed for non-unitary TMACs, see figure 6(c). Also, from the inset of figure 6(b) we see that even the slip-corrected permeability, obtained from the NSEs with FVBC, is only accurate when $k_{a}/k_{\infty }\lesssim 1.5$ when $\unicode[STIX]{x1D6FC}=0.5$ and $0.1$ .

So far, we have identified two key factors that could lead to the observation that the correction factor $b$ in (1.2) decreases when $Kn$ (the reciprocal mean gas pressure) increases. The first factor is that the TMAC should be less than 1, and the second one is that flow streamlines should be tortuous. The two factors should be combined together to explain Klinkenberg’s observations, as we have shown that the variation of TMAC in the straight cylindrical tube (see figure 1) and the variation of porous media but with TMAC being one (see figures 3 and 4) cannot yield the result that the correction factor decreases when $Kn$ increases.

The reason that Klinkenberg did not observe the increase of the correction factor with the Knudsen number is probably because $Kn$ is not large enough in his experiments, as from tables 5 to 7 in Klinkenberg (1941) the maximum value of $k_{a}/\unicode[STIX]{x1D705}_{\infty }$ is approximately 30. In this region of $k_{a}/\unicode[STIX]{x1D705}_{\infty }$ , from figures 5(c,e) and 6(b,c) we see that, roughly, the correction factor $b^{\prime }$ decreases when $Kn^{\ast }$ increases.

4 Conclusions

The main purpose of this paper is devoted to understanding Klinkenberg’s famous experiment on the apparent gas permeability of porous media, i.e. to understand why in his experiments the correction factor decreases when the Knudsen number increases. By solving the linearized Bhatnagar–Gross–Krook equation in simplified but not simple porous media via the discrete velocity method, for the first time in the last 76 years we have pinpointed two key ingredients that lead to his discovery: the tortuous flow path and the coexistence of the diffuse and specular scatterings of the gas molecules when impinging on solid surfaces. Since in experiments the flow paths are always tortuous, the experimental results of Klinkenberg and others suggest that the tangential momentum accommodation coefficient for the gas–surface interaction is less than one.

We have also found that Klinkenberg’s results can only be observed when the ratio between the apparent and intrinsic permeabilities is ${\lesssim}30$ ; at large ratios the correction factor increases with the Knudsen number.

Our numerical solutions not only shed new light on the complex flow physics in porous media, but also could serve as benchmarking cases to check the accuracy of various macroscopic models and/or numerical methods for modelling/simulation of rarefied gas flows in complex geometries. Specifically, through theoretical and numerical analysis, we have shown that the Navier–Stokes equations with the first-order velocity-slip boundary condition can only predict the apparent gas permeability of porous media to the first-order accuracy of the Knudsen number. Although the slip-corrected expression, which only retains the linear dependence of the permeability and Knudsen number, works fine for the diffuse gas–surface boundary condition in a wide range of gas rarefaction, its accuracy is significantly reduced when the tangential momentum accommodation coefficient $\unicode[STIX]{x1D6FC}$ in the diffuse-specular boundary condition deviates from one. For example, when $\unicode[STIX]{x1D6FC}=0.5$ and 0.1, the slip-corrected expression is only accurate when the ratio between the apparent and intrinsic permeabilities is ${\lesssim}1.5$ . Thus, the slip-corrected expression is of no practical use. These issues must be properly taken into account since the Navier–Stokes equations are widely misused in rarefied gas flows to predict the apparent gas permeability.

Our work also implies that, all the currently widely used empirical solutions, such as (1.4), derived from straight cylindrical tube with the diffuse boundary condition should be reformulated to incorporate the tortuosity and diffuse-specular scattering, to predict the unconventional gas production by feeding the apparent gas permeability into upscaling models. This forms our future research directions. Also, it should be noted that all our conclusions are based on the simulation in two-dimensional geometry; whether the results hold or not for three-dimensional porous media requires future systematic investigation.