On the apparent permeability of porous media in rarefied gas flows

TL;DR

This paper investigates the apparent permeability in porous media during rarefied gas flows, revealing key factors like tortuosity and gas-surface interactions that explain experimental observations and improve modeling accuracy.

Contribution

It identifies two main factors—tortuous flow paths and non-unitary tangential momentum accommodation—that explain Klinkenberg's experimental results, advancing understanding of gas flow in porous media.

Findings

The correction factor decreases with increasing Knudsen number due to tortuosity and surface interactions.

Klinkenberg's results are valid only when the permeability ratio is below approximately 30.

Numerical results serve as benchmarks for modeling rarefied gas flows in complex geometries.

Abstract

The apparent gas permeability of the 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 $\lesssim30$; 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 modeling/simulation of rarefied gas flows in complex geometries over a wide range of gas rarefaction.

Figures

Tables

Table 1. 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.

Geometry		circular cylinder in Fig. 3			complex geometry in Fig. 4
Spatial grid size		$801\times 401$			$3000\times 1500$
$Kn$	$K{n}^{\ast }$	iteration steps	time ($s$)	$K{n}^{\ast }$	iterations	time ($s$)
${10}^0$	$1.85\times {10}^0$	$2\times {10}^3$	$68$	$7.3\times {10}^1$	$4\times {10}^3$	$2076$
${10}^{-1}$	$1.85\times {10}^{-1}$	$2\times {10}^3$	$68$	$7.3\times {10}^0$	$3\times {10}^3$	$1558$
${10}^{-2}$	$1.85\times {10}^{-2}$	$7\times {10}^3$	$238$	$7.3\times {10}^{-1}$	$5\times {10}^3$	$2584$
${10}^{-3}$	$1.85\times {10}^{-3}$	$366\times {10}^3$	$12381$	$7.3\times {10}^{-2}$	$9\times {10}^3$	$4943$

Equations

$$\text{q}=-\frac{k_{\infty }}{\mu }\nabla p,$$

$$k_a=k_{\infty }\left(1+\frac{b}{\bar{p}}\right),$$

$$\begin{gathered} Kn=\frac{\lambda }{L}, \\ \text{and} \\ \lambda =\frac{\mu (T_0)}{\bar{p}}\sqrt{\frac{\pi R{T}_0}{2}}, \end{gathered}$$

$$\frac{k_a}{k_{\infty }}=\left[1+\frac{128}{15{\pi }^2}{\tan }^{-1}\left(4K{n}^{0.4}\right)Kn\right]\left(1+\frac{4Kn}{1+4Kn}\right),$$

$$\frac{\partial f}{\partial t}+v_1\frac{\partial f}{\partial x_1}+v_2\frac{\partial f}{\partial x_2}+v_3\frac{\partial f}{\partial x_3}=\mathcal{C}(f),$$

$$f_{eq}=\frac{\exp (-|\text{v}{|}^2)}{{\pi }^{3/2}},$$

$$\frac{\partial h}{\partial t}+v_1\frac{\partial h}{\partial x_1}+v_2\frac{\partial h}{\partial x_2}=\frac{\sqrt{\pi }}{2Kn}\left[\varrho +2u_1{v}_1+2u_2{v}_2+\tau \left(|\text{v}{|}^2-\frac{3}{2}\right)-h\right],$$

$$\begin{array}{r}\varrho =\int h{f}_{eq}d\text{v},(u_1,u_2)=\int (v_1,v_2)h{f}_{eq}d\text{v},\tau =\frac{2}{3}\int |\text{v}{|}^{2}h{f}_{eq}d\text{v}-\varrho .\end{array}$$

$$h\left(\mp 0.5,x_2,v_1,v_2,v_3\right)=\pm 1+h\left(\pm 0.5,x_2,v_1,v_2,v_3\right),\text{ when }{v}_1\gtrless 0,$$

$$h\left(x_1,0,v_1,v_2,v_3\right)=$$

$$h\left(x_1,0.5,v_1,v_2,v_3\right)=$$

$$h(\text{v})=\frac{2\alpha }{\pi }{\int }_{v_n^{\prime }<0}|v_n^{\prime }|h({\text{v}}^{\prime })\exp (-|{\text{v}}^{\prime }{|}^2)d{\text{v}}^{\prime }+(1-\alpha )h(\text{v}-2\text{n}{v}_n),$$

$$k_a=2\sqrt{\frac{1}{\pi }}Kn{G}_P,$$

$$\frac{\partial \rho }{\partial t}+\frac{\partial \rho u_i}{\partial x_i}=0,$$

$$\frac{\partial \rho u_i}{\partial t}+\frac{\partial \rho u_i{u}_j}{\partial x_j}+\frac{\partial {\sigma }_{ij}}{\partial x_j}=-\frac{\partial p}{\partial x_i},$$

$$\frac{\partial \rho T}{\partial t}+\frac{\partial \rho u_{i}T}{\partial x_i}+\frac{2}{3R}\frac{\partial q_i}{\partial x_i}=-\frac{2}{3R}\left(p\frac{\partial u_i}{\partial x_i}+{\sigma }_{ij}\frac{\partial u_j}{\partial x_i}\right).$$

$$f=f^{(0)}+Kn{f}^{(1)}+K{n}^2{f}^{(2)}+\cdots ,$$

$$f=f^{(0)}+Kn{f}^{(1)},$$

$$\begin{gathered} {\sigma }_{ij}=-2\mu \frac{\partial u_{<i}}{\partial x_{j>}} \\ \text{and} \\ {q}_i=-\frac{15}{4}R\mu \frac{\partial T}{\partial x_i}, \end{gathered}$$

$$\frac{\partial {\sigma }_{ij}}{\partial t}+\frac{\partial u_k{\sigma }_{ij}}{\partial x_k}+\frac{\partial m_{ijk}}{\partial x_k}=-\frac{p}{\mu }{\sigma }_{ij}-2p\frac{\partial u_{<i}}{\partial x_{j>}}-\frac{4}{5}\frac{\partial q_{<i}}{\partial x_{j>}}-2{\sigma }_{k<i}\frac{\partial u_{j>}}{\partial x_k},$$

$$\frac{\partial q_i}{\partial t}$$

$$\frac{\partial m_{ijk}}{\partial t}$$

$${\varphi }_{ijkl}=-\frac{4\mu }{C_1\rho }\frac{\partial m_{<ijk}}{\partial x_{l>}},R_{ij}=-\frac{24}{5}\frac{\mu }{p}RT\frac{\partial q_{<i}}{\partial x_{j>}},\Delta =-12\frac{\mu }{p}RT\frac{\partial q_k}{\partial x_k},$$

$$u_{\tau }=-\frac{2-\alpha }{\alpha }\sqrt{\frac{\pi RT}{2}}\frac{{\sigma }_{n\tau }}{p_{\alpha }}-\frac{5m_{nn\tau }+2q_{\tau }}{10p_{\alpha }},$$

$$RT-R{T}_w=-\frac{2-\alpha }{\alpha }\sqrt{\frac{\pi RT}{2}}\frac{q_n}{2p_{\alpha }}-\frac{RT{\sigma }_{nn}}{4p_{\alpha }}+\frac{u_{\tau }^2}{4}-\frac{75R_{nn}+28\Delta }{840p_{\alpha }}+\frac{{\varphi }_{nnnn}}{24p_{\alpha }},$$

$$k_{a}=\frac{1}{8\pi\left(1+2\xi{Kn}\sqrt{\frac{\pi}{\phi}}\right)}\bigg{[}-\ln\phi-\frac{3}{2}+2\phi-\frac{\phi^{2}}{2}+2\xi{Kn}\sqrt{\frac{\pi}{\phi}}\left(-\ln\phi-\frac{1}{2}+\frac{\phi^{2}}{2}\right)\bigg{]},$$

$$k_{a,\text{slip}}=\frac{1}{8\pi }\left[-\ln \varphi -\frac{3}{2}+2\varphi -\frac{{\varphi }^2}{2}+2\xi Kn\sqrt{\frac{\pi }{\varphi }}\left(1-2\varphi +{\varphi }^2\right)\right].$$

$$L^{\ast }=L\sqrt{\frac{12k_{\infty }}{ϵ}},$$

$$K{n}^{\ast }=\frac{\lambda }{L^{\ast }}=Kn\sqrt{\frac{ϵ}{12k_{\infty }}}.$$

$$\frac{k_a(i)-k_a(i-1000)}{k_a(i)}<10^{-10},$$

$$\frac{k_a}{k_{\infty }}=1+\frac{b}{\bar{p}}\equiv 1+b^{\prime }K{n}^{\ast },$$

$$\frac{k_a}{k_{\infty }}=1+6\xi K{n}^{\ast },$$

Full text

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 the 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 q (i.e. discharge per unit area, with units of length per time) is proportional to the pressure gradient $\nabla{p}$:

[TABLE]

where $\mu$ 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:

[TABLE]

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 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 $\lambda$ of gas molecules to the characteristic flow length $L$):

[TABLE]

where $\mu(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 regimes111Note 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 is 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 systematical experimental measurements, Klinkenberg found that the linear relation (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, 6, and 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 Fig. 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):

[TABLE]

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. Fig. 1 shows the numerical solutions of the Bhatnagar-Gross-Krook (BGK) kinetic equation (see § 2.1), as well as that of NSEs with velocity-slip boundary conditions. It is clear that both the BGK equation and NSEs with 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. NSEs with second-order velocity-slip boundary condition (Karniadakis et al., 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 about 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.

The accuracy of gas kinetic equations in modeling 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 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 et al. (2016) found that the AGP of the porous medium is not only a nonlinear function of $Kn$, but also has the 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.

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 analyze the limitation of the NSEs and FVBC. In § 3, numerical simulations are performed to analyze 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 analyzed, 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 Fig. 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,\textbf{x},\textbf{v})$ to describe the system state:

[TABLE]

where $\textbf{v}=(v_{1},v_{2},v_{3})$ is the three-dimensional molecular velocity normalized by the most probable speed $v_{m}=\sqrt{2RT_{0}}$, $\textbf{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 et al., 1954), resulting in the BGK equation.

When the porous medium is so long that the pressure gradient is small, namely, $|Ldp/pdx_{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

[TABLE]

and the perturbation $h$ is governed by:

[TABLE]

with macroscopic quantities such as the perturbed number density $\varrho$ of gas molecules, the velocity $u_{1}$ and $u_{2}$, and the perturbed temperature $\tau$ being calculated as

[TABLE]

The kinetic equation (7) 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):

[TABLE]

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

[TABLE]

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

[TABLE]

where n is the outer normal vector of the solid surface, and $\alpha$ 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-\alpha$, 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 $\alpha=1$ or $\alpha=0$, respectively. Note that most of the current studies, including the widely used equation (4), use the diffuse boundary condition.

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

[TABLE]

where $G_{p}=2\int_{0}^{1/2}u_{1}(x_{2})dx_{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 $\rho$, velocity $u_{i}$, and temperature $T$; and its dynamics is described by the Euler equations or 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 order of approximations.

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

[TABLE]

However, the above equations are not closed, since expressions for the shear stress $\sigma_{ij}$ and heat flux $q_{i}$ are not known. One way to close (13)-(15) is to use the Chapman-Enskog expansion, where the distribution function is expressed in the power series of $Kn$ (Chapman & Cowling, 1970):

[TABLE]

where $f^{(0)}$ is the equilibrium Maxwellian distribution function. When $f=f^{(0)}$, we have $\sigma_{ij}=0$ and $q_{i}=0$, and (13)-(15) reduce to the Euler equations. When the distribution function is truncated at the first-order of $Kn$, that is,

[TABLE]

we have

[TABLE]

where the angular brackets are used to denote the traceless part of a symmetric tensor, so that (13)-(15) reduce to NSEs. When $f=f^{(0)}+Knf^{(1)}+Kn^{2}f^{(2)}$, Burnett equations can be derived but are not used nowadays due to their intrinsic instabilities (Garcia-Colin et al., 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 R20 equations are used, which, in addition to (13)-(15), include governing equations for the high-order moments $\sigma_{ij}$, $q_{i}$, and $m_{ijk}$:

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

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 $\tau_{i}$ the tangential vector of the wall, the velocity-slip parallel to the wall $u_{\tau}$ and temperature-jump conditions are:

[TABLE]

where ${p_{\alpha}}=p+\sigma_{nn}/{2}-(30R_{nn}+7\Delta)/840RT-\phi_{nnnn}/24RT$. The rest of wall boundary conditions for higher-order moments are listed in Appendix A. Note that the velocity-slip boundary condition (23) is also of higher-order due to the appearance of the higher-order moment $m_{nn\tau}$.

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 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 some AGP of the porous medium as NSEs when $Kn\rightarrow 0$, and be more accurate than 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 et al., 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 Fig. 2. Using the NSEs and FVBC, when the porosity $\epsilon=1-\pi{D}^{2}/4L^{2}$ is large, the AGP can be obtained analytically (Chai et al., 2011):

[TABLE]

where $\phi=1-\epsilon$ is the solid fraction and $\xi$ is related to the viscous velocity slip coefficient. When the Maxwell’s diffuse boundary condition (11) is used (Hadjiconstantinou, 2003), we have $\xi=1.016\sqrt{{4}/{\pi}}=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 (7), 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 (Sharipov & Graur, 2012) and the unified gas kinetic scheme (Huang et al., 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 Fig. 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 “stair-case” 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 et al., 2014).

The results for a periodic porous medium with $\epsilon=0.8$ are shown in Fig. 3. It is seen that both the discrete velocity method on the “stair-case” geometry and the unified gas kinetic scheme on the curvilinear mesh have 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 re-normalized. Also, the mean free path defined in (3) is $15\pi/2(7-2\omega)(5-2\omega)$ times larger than that used in DSMC based on the variable hard sphere model, where $\omega=0.81$ for argon (Wu et al., 2013); so the Knudsen number should be re-scaled.

The accuracy of the permeability (25) 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 (25) 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 (17) adopted in the derivation of 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 NSEs with FVBC, higher-order moments in (19)-(21) and the higher-order velocity slip in (23), which are derived from the Boltzmann equation and the gas kinetic boundary condition (11) 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 (25), a nonlinear function of $Kn$, into the Taylor’s 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

[TABLE]

The slip-corrected AGP is plotted as the dotted line in Fig. 3(b). Surprisingly, this Taylor’s 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 about 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.

The conclusion that the NSEs with FVBC is accurate only to the first-order of $Kn$ not only holds for the simple porous medium in Fig. 2, but also applies to more complex porous media, for example, see the unit cell in Fig. 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-centered finite-volume discretization scheme, on unstructured grids. A body-fitted computational grid is generated using the OpenFOAM meshing tool, resulting in a mesh of about 600,000 cells of which the majority are hexahedra and the rest few close to the walls are prisms. Note that the results of 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):

[TABLE]

and the effective Knudsen number is

[TABLE]

From Fig. 4(b) we see that the AGP obtained from NSEs with FVBC increases linearly when $Kn\lesssim 0.002$ (or $Kn^{\ast}\lesssim 0.146$) and then quickly reaches to a constant. Again, the comparison with solution of the linearized BGK equation shows that the NSEs with FVBC are roughly accurate when the AGP is a linear function of $Kn$; in this region of $Kn$, the maximum AGP is only about 1.5 times larger than the intrinsic permeability $k_{\infty}$. Interestingly, like (26), 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 about 15% in a large range of the 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 Dell workstation (Precision Tower 7910) with Dual processors 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 roughly the same as the Gauss-Hermite quadrature, but the computational time is about 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

[TABLE]

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

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 about 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 et al., 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 350000 when $Kn=1$. The implicit unified gas kinetic scheme can speed up the convergence by hundreds of times (Zhu et al., 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 about an hour to obtain the converged solution for the complex porous medium case in Fig. 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 Fig. 3, single processor of Intel i7-4770 CPU at 3.4GHz 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, the 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 Fig. 1.

From the analytical solution (26), we find that, if we use a non-unitary TMAC (say, $\alpha=0.5$), then $\xi\approx 3.23$ (Loyalka et al., 1975), and the correction factor is about three times larger than that of $\alpha=1$ when $Kn$ is small. If at large Knudsen numbers the AGP of $\alpha=0.5$ is only slightly larger than that of $\alpha=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 Fig. 2 with squares. The only reason to do this is that the specular reflection in the boundary condition (11) can be accurately implemented.

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

[TABLE]

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.

Our numerical solutions for the linearized BGK equation, solved by the discrete velocity method, are shown in Fig. 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 Fig. 5(c,d) we see that, as the numerical results in § 3.1, when the TMAC is $\alpha=1$, the slip-corrected permeability agrees well with the numerical results of the linearized BGK equation, when the porosity is $\epsilon=0.4$ and $0.8$. Also, it is interesting to note that, when $Kn^{\ast}\rightarrow 0$, the correction factor $b^{\prime}$ is about 6 and 6.1 when the porosity is $\epsilon=0.4$ and 0.8, respectively. For Poiseuille flow between two parallel plates, the slip-corrected permeability is exactly

[TABLE]

where $\xi=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 $\alpha=0.5$ and 0.1, $\xi$ is respectively $3.23$ and 19.3, the correction factor $b^{\prime}$ from the NSEs with FVBC is about 17 and 100, respectively. From Fig. 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 $\alpha=0.5$ and 0.1, respectively; when the porosity is 0.4, the slip-corrected permeability overpredicts the AGP by 360% when $\alpha=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 Fig. 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 $\alpha=1$, it is seen that the correction factor $b^{\prime}$ increases with $Kn$, which clearly contradicts Klinkenberg’s experimental results. However, when $\alpha=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 (12) 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 behaviors for $G_{p}$ are observed for all the cases considered in this paper. According to (30), the equivalent correction factor $b^{\prime}$ at $Kn^{\ast}\rightarrow\infty$ can be calculated as follows:

[TABLE]

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 Fig. 5(e,f).

At small $Kn$, when $\alpha=0.5$ and 0.1, the correction factor decreases when $Kn$ increases. For instance, when the porosity is $\epsilon=0.8$, the correction factor $b^{\prime}$ with $\alpha=0.5$ drops from 17 when $Kn^{\ast}\rightarrow 0$ to the minimum value 9 at $Kn^{\ast}=0.4$, while $b^{\prime}$ with $\alpha=0.1$ drops from 100 at $Kn^{\ast}\rightarrow 0$ to the minimum value 13 at $Kn^{\ast}=0.6$. When the porosity is $\epsilon=0.4$, the correction factor $b^{\prime}$ with $\alpha=0.5$ drops from 17 at $Kn^{\ast}\rightarrow 0$ to the minimum value 12 at $Kn^{\ast}=0.5$, while $b^{\prime}$ with $\alpha=0.1$ drops from 100 when $Kn^{\ast}\rightarrow 0$ to the minimum value 27 at $Kn^{\ast}=1.1$. This behavior, is in qualitatively agreement222Note 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 Fig. 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 Fig. 6(c). Also, from the inset of Fig. 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 $\alpha=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 (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 Fig. 1) and the variation of porous media but with TMAC being one (see Figs. 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, 6, and 7 in Klinkenberg (1941) the maximum value of $k_{a}/\kappa_{\infty}$ is around 30. In this region of $k_{a}/\kappa_{\infty}$, from Fig. 5(c,e) and Fig. 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 the Klinkenberg’s famous experiment on the apparent gas permeability of the 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, Klinkenberg’s and others experimental results 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 modeling/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 $\alpha$ in the diffuse-specular boundary condition deviates from one. For example, when $\alpha=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 (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 systematical investigation.

Acknowledgements

LW acknowledges the support of an Early Career Researcher International Exchange Award from the Glasgow Research Partnership in Engineering (allowing him to visit the Hong Kong University of Science and Technology for one month) and the joint project from the Royal Society of Edinburgh and National Natural Science Foundation of China. LG thanks Lianhua Zhu for helpful discussions on the first-order velocity-slip boundary condition in OpenFoam. This work is also partly supported by the Engineering and Physical Sciences Research Council in the UK under grant EP/M021475/1.

Appendix A Wall boundary conditions for high-order moments

The wall boundary conditions for higher-order moments in the regularized 20-moment equations are given as follows:

[TABLE]

where ${\hat{u}_{\tau}}=u_{\tau}/\sqrt{RT}$ and ${\hat{T}_{w}}=T_{w}/T$.