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 Reference Bhatnagar, Gross and Krook1954), 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 Reference Sharipov and Graur2012):

(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 Reference Maxwell1879):

(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 Reference Chapman and Cowling1970):

(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 Reference Garcia-Colin, Velasco and Uribe2008).

Alternatively, by combining the moment method of Grad (Reference Grad1949) and the Chapman–Enskog expansion, the regularized 13-, 20- and 26-moment equations (Struchtrup & Torrihon Reference Struchtrup and Torrihon2003; Gu & Emerson Reference Gu and Emerson2009) can be derived from the Boltzmann equation to describe rarefied gas flows at different levels of rarefaction (Torrilhon Reference Torrilhon2016). 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 (Reference Gu and Emerson2009) 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. Reference Taheri, Rana, Torrilhon and Struchtrup2009; Gu, Emerson & Tang Reference Gu, Emerson and Tang2010) 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 Reference Gu and Emerson2009).

Macroscopic wall boundary conditions were obtained from the diffuse-specular boundary condition (Maxwell Reference Maxwell1879). 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 Reference Struchtrup2005), 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.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. Reference Chai, Lu, Shi and Guo2011):

(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 Reference Hadjiconstantinou2003), 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. (Reference Loyalka, Petrellis and Storvick1975) and Gibelli (Reference Gibelli2011). 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 Reference Chu1965; Aoki, Sone & Yamada Reference Aoki, Sone and Yamada1990) and the unified gas kinetic scheme (Huang, Xu & Yu Reference Huang, Xu and Yu2013). 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 (Reference Gu and Emerson2009). 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 Reference Wu, Reese and Zhang2014).

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. Reference Borner, Panerai and Mansour2017). 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. Reference Wu, White, Scanlon, Reese and Zhang2013); 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. (Reference Borner, Panerai and Mansour2017).

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. (Reference Lasseux, Valdes Parada and Porter2016), 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. Reference Lasseux, Valdes Parada and Porter2016):

(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 Reference Pan, Luo and Miller2006).

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 Reference Zhu, Zhong and Xu2016), 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. Reference Tang, Zhai, Tao, Gu and Emerson2013). 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. Reference Loyalka, Petrellis and Storvick1975) 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 Reference Takata and Funagane2011) 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 (Reference Klinkenberg1941) and Moghaddam & Jamiolahmady (Reference Moghaddam and Jamiolahmady2016).

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 (Reference Klinkenberg1941) 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.