Comparative study of the discrete velocity and lattice Boltzmann methods for rarefied gas flows through irregular channels

TL;DR

This paper compares the discrete velocity and lattice Boltzmann methods for simulating rarefied gas flows in irregular channels, revealing accuracy limitations of the lattice Boltzmann approach in complex geometries.

Contribution

It provides a detailed comparison and benchmarking of the discrete velocity and lattice Boltzmann methods for complex, irregular geometries in rarefied gas flow simulations.

Findings

Lattice Boltzmann method's accuracy decreases significantly in rough and porous media.

Discrete velocity method offers more accurate results for complex geometries.

Benchmark cases established for future method development.

Abstract

Rooted from the gas kinetics, the lattice Boltzmann method is a powerful tool in modeling hydrodynamics. In the past decade, it has been extended to simulate the rarefied gas flow beyond the Navier-Stokes level, either by using the high-order Gauss-Hermite quadrature, or by introducing the relaxation time that is a function of the gas-wall distance. While the former method, with a limited number of discrete velocities (i.e. D2Q36), is accurate up to the early transition flow regime, the latter method, with the same discrete velocities as that used in simulating hydrodynamics (i.e. D2Q9), is accurate up to the free-molecular flow regime in the Poiseuille flow between two parallel plates. This is quite astonishing in the sense that more discrete velocities are less accurate. In this paper, by solving the Bhatnagar-Gross-Krook kinetic equation accurately via the discrete velocity method, we find that the accuracy of the lattice Boltzmann method is reduced significantly in the simulation of rarefied gas flows through the rough surface and porous media. Our simulation results could serve as benchmarking cases for future development of the lattice Boltzmann method for modeling and simulation of rarefied gas flows in complex geometries.

Figures

Equations

$$r(x_1/H)=A\sum _{n=n_1}^{\infty }\frac{\cos \left(2\pi {\gamma }^n{x}_1/H\right)}{{\gamma }^{(2-D)n}},$$

$$ϵ=\frac{\sigma }{H},$$

$$\frac{\partial f}{\partial t}+\text{v}\cdot \frac{\partial f}{\partial \text{x}}=\frac{\sqrt{\pi }}{2\text{Kn}}(F_{eq}-f),$$

$$F_{eq}=\frac{\rho }{\pi T}\exp (-\frac{|\text{v}-\text{u}{|}^2}{T}),$$

$$\lambda =\frac{\mu (T_0)}{\bar{p}}\sqrt{\frac{\pi R{T}_0}{2}}.$$

$$f=f_{eq}(1+h),f_{eq}=\frac{\exp (-|\text{v}{|}^2)}{\pi },$$

$$\frac{\partial h}{\partial t}+\text{v}\cdot \frac{\partial h}{\partial \text{x}}$$

$$\mathcal{L}(\varrho ,\text{u},\tau )$$

$$h(x_1,1,\text{v})=2\sqrt{\pi }{\int }_{v_2^{\prime }>0}{v}_2^{\prime }h(x_1,1,{\text{v}}^{\prime })f_{eq}({\text{v}}^{\prime })\mathrm{d}{\text{v}}^{\prime }+\frac{2U_0}{v_m}{v}_1,$$

$$h(x_1,x_2,\text{v})=2\sqrt{\pi }{\int }_{v_n^{\prime }<0}|v_n^{\prime }|h(x_1,x_2,{\text{v}}^{\prime })f_{eq}({\text{v}}^{\prime })\mathrm{d}{\text{v}}^{\prime },$$

$$u_1(x_1,x_2)=U_0\int v_{1}f(x_1,x_2,\text{v})\mathrm{d}\text{v},$$

$$Q=\frac{\iiint v_{1}f(x_1,x_2,\text{v})\mathrm{d}\text{v}\mathrm{d}{x}_1\mathrm{d}{x}_2}{5}.$$

$$v_{1,2}=\frac{4}{(N_v-1{)}^3}(-N_v+1,-N_v+3,\cdots ,N_v-1{)}^3,$$

$$\frac{\sqrt{\pi }}{2\text{Kn}}{h}^{(n+1)}+\text{v}\cdot \frac{\partial h^{(n+1)}}{\partial \text{x}}=\mathcal{L}({\varrho }^{(n)},{\text{u}}^{(n)},{\tau }^{(n)}),$$

$$Q^{\ast }=\frac{Q_r}{Q_s},$$

$$k_a=\sqrt{\frac{4}{\pi }}Q\text{Kn},$$

$$Q=\frac{\iiint v_{1}f(x_1,x_2,\text{v})\mathrm{d}\text{v}\mathrm{d}{x}_1\mathrm{d}{x}_2}{L},$$

$$H^{\ast }=H\sqrt{\frac{12k_{\infty }}{\varphi }},$$

$${\text{Kn}}^{\ast }=\frac{\lambda }{H^{\ast }}=\text{Kn}\frac{H}{H^{\ast }}=\text{Kn}\sqrt{\frac{\varphi }{12k_{\infty }}}.$$

Full text

text

I Introduction

When the mean free path $\lambda$ of gas molecules is comparable to or even larger than the characteristic flow length $H$, the Navier-Stokes (NS) equations derived from the continuum-fluid hypothesis fails, and the gas kinetic theory based on the fundamental Boltzmann equation for the velocity distribution function (VDF) of gaseous molecules is required to describe rarefied gas dynamics. According to the Chapman-Enskog expansion, NS equations are only the first-order approximation in the Knudsen number ($\text{Kn}=\lambda/H$) to the Boltzmann equation Chapman and Cowling (1970). Therefore, they are valid in the continuum flow regime where $\text{Kn}\lesssim 0.001$ Gad-el-Hak (1999). As the Knudsen number increases, higher-order (non-equilibrium) terms begin to dominate, and NS equations gradually loss validity; for instance, when $Kn$ is large, non-equilibrium terms embodied in the macroscopic mass, momentum and energy transport cannot be expressed in terms of the lower-order macroscopic quantities. The non-equilibrium effects not only cause noticeable velocity slip and temperature jump at solid surfaces in the slip flow regime ($0.001\lesssim\text{Kn}\lesssim 0.1$), but also modify the constitutive relations in the transition ($0.1\lesssim\text{Kn}\lesssim 10$) and free-molecular ($10\lesssim\text{Kn}$) flow regimes such that the Newton’s law for stress and strain and Fourier’s law for heat flux and temperature gradient do not hold anymore. This leads to a number of counterintuitive phenomena, including the thermal transpiration where gas molecules along solid surface move from the cold region to hot Reynolds (1879), the Knudsen paradox where the dimensionless mass flow rate in Poiseuille flow could increase when the gas pressure decreases Steckelmacher (1999), the temperature bimodality in the force-driven Poiseuille flow Mansour et al. (1997), the inverted velocity in cylindrical Couette flow Aoki et al. (2003), and the gas anti-resonance where the shear stress in an oscillating lid-driven cavity flow could be smaller than that of the one-dimensional (1D) Couette flow Wu et al. (2014a).

In the past decades, due to the rapid development of microelectromechanical systems Karniadakis et al. (2005) and the shale gas revolution in North America Wang et al. (2014), extensive works have been devoted to construct efficient numerical schemes to simulate gas flows at the micro scale, where the flow velocity is usually far smaller than the most probable speed of the gas molecules. In most of these applications, gas flows vary from the slip to the free-molecular flow regimes and the gas-surface interaction dominates flow behavior. High-fidelity numerical methods to solve the Boltzmann equation include the numerical kernel method Ohwada et al. (1989), the conservative projection-interpolation method Tcheremissine (2005), the low-variance direct simulation Monte Carlo method Homolle and Hadjiconstantinou (2007), and the fast spectral method Wu et al. (2014b), to name just a few. Due to the high computational cost, however, the Boltzmann equation is usually simplified by the Bhatnagar-Gross-Krook (BGK) equation under the single relaxation time (SRT) approximation Bhatnagar et al. (1954), which is often solved by the discrete velocity method (DVM) where the continuous molecular velocity space is represented by a few number of discrete velocities Chu (1965); Yang and Huang (1995); Sharipov (1996); Xu and Huang (2010); Guo et al. (2013). Generally speaking, rarefied gas flows with higher values of Kn need a larger number of discrete velocities to resolve the large variations/discontinuities in the VDF, for instance see the numerical examples in Refs. Wu et al. (2014b); S.Taguchi and Aoki (2012).

Rooted from the gas kinetics, the lattice Boltzmann method (LBM) is a popular and powerful tool in modeling the NS hydrodynamics and beyond. Historically, the SRT-LBM is firstly developed as an alternative solver for the NS equations. Since it uses a very limited but highly optimized number of discrete velocities, e.g. the D2Q9 scheme for 2D problems, the SRT-LBM can be viewed as a special type of DVM to solve the BGK equation He and Luo (1997). To capture the non-equilibrium effects beyond the NS hydrodynamics, a rigorous procedure for obtaining high-order approximations to the BGK equation is proposed Shan et al. (2006) and tested in various canonical problems Kim et al. (2008); Meng and Zhang (2011a, b); Shi et al. (2011); Ambruş and Sofonea (2012, 2014); Shi et al. (2015). By expanding the VDF into high-order Gauss-Hermite polynomials or Gauss quadrature in the spherical coordinate system, it is found that, with a not very large number of discrete lattice velocities, e.g. D2Q36 or D2Q64, high-order LBM schemes accurately describe rarefied gas flows up to the early transition flow regime, i.e. $Kn\lesssim 0.5$.

In addition to high-order quadrature, Zhang et al. Zhang et al. (2006) attempted to capture the Knudsen layer structure in rarefied gas flows within the framework of the SRT-LBM, but by introducing the gas kinetic boundary condition and the effective relaxation time that is a function of the gas-wall distance. In the pressure-driven flow, this “wall-scaling” approach provides a significant improvement for Knudsen numbers up to 0.5. Later, this scheme has also been applied to study the rarefied thermal Zhang et al. (2007) and oscillatory Couette flow Tang et al. (2008) in the early transition flow regime with great success.

Further improvement has been achieved by Guo, Zheng and Shi Guo et al. (2008), who developed a numerical scheme within the framework of multiple relaxation time (MRT) LBM. In addition to the “wall-scaling” of the relaxation time, the combined bounce-back and specular-reflection boundary condition is designed, such that the new scheme is equivalent to solve the NS equations with the second-order slip boundary condition in the slip flow regime. In the Poiseuille flow between two parallel plates, it is found that this MRT-LBM scheme is even able to predict the velocity profile and mass flow rate with good accuracy, in the whole transition flow regime. For this reason, the MRT-LBM has attracted significant attentions Guo et al. (2008); Li et al. (2011); Yuan and Rahman (2016); Deng et al. (2016); Zhao et al. (2016).

In the simulation of rarefied gas flows, it is quite astonishing that the MRT-LBM of Guo et al. with fewer number of discrete velocities Guo et al. (2008) are more accurate than the LBM of higher-order quadrature Shan et al. (2006); Kim et al. (2008); Meng and Zhang (2011a, b); Shi et al. (2011); Ambruş and Sofonea (2012, 2014); Shi et al. (2015). Although Guo et al. honestly and explicitly pointed out that their scheme is only designed and analyzed for plane walls Guo et al. (2008) and is hard to be extended to generalized gas-surface boundaries Guo et al. (2014), it has been applied to study the rarefied gas flows through the microchannel with rough surface and complex porous media Deng et al. (2016); Zhao et al. (2016). However, whether this scheme works on walls with curvatures and other complex geometries or not is not clear.

In this paper, aiming to address the accuracy of the MRT-LBM in simulating rarefied gas flows in complex geometries, the BGK equation will be solved by the DVM, and numerical results will be compared to those of MRT-LBM for the Couette flow through rough microchannels Deng et al. (2016) in Sec. II, and for the Poiseuille flow through porous media Zhao et al. (2016) in Sec. III.

II Couette flow in rough microchannels

As Ref. Deng et al. (2016), we first consider the Couette flow between two plates with a distance of $H$, see Fig. 1. The smooth top plate move in the $x_{1}$ direction with a speed $U_{0}$, while the stationary bottom plate has a roughness described by the following Weierstrass-Mandelbrot fractal function:

[TABLE]

where $D$ is the self-affine fractal dimension and $\gamma$ determines the frequency spectrum of the surface roughness. The aspect ratio of the microchannel is 5, and we choose $\gamma=1.5$ and $n_{1}=-4$ in Eq. (1). The scaling parameter $A$ is used to adjust the surface roughness $\epsilon$:

[TABLE]

where $\sigma$ is the root mean square of $r(x_{1})$.

We are interested in the horizontal velocity profile of the gas and the mass flow rate at different level of gas rarefaction, when the steady state is reached.

II.1 The gas kinetic theory and the DVM

The fundamental Boltzmann equation describes the evolution of the VDF in dependence of spatial position and time. Macroscopic quantities are obtained from the velocity moments of the VDF. In this paper, we adopt the BGK equation Bhatnagar et al. (1954) instead of the full Boltzmann collision operator. Numerical simulation for Poiseuille flow shows that the BGK equation can produce accurate mass flow rates Sharipov and Bertoldo (2009); Doi (2010).

The BGK equation for the 2D problem may be written in the following form

[TABLE]

where $f$ is the reduced two-velocity distribution function normalized by $\bar{p}/v_{m}^{2}mRT_{0}$ with $\bar{p}$ being the mean gas pressure at the reference temperature $T_{0}$ and $m$ the molecular mass, $\textbf{v}(v_{1},v_{2})$ is the molecular velocity normalized by the most probable speed $v_{m}=\sqrt{2RT_{0}}$ with $R$ being the gas constant, $\textbf{x}(x_{1},x_{2})$ is the spatial coordinate normalized by the height $H$, $t$ is the time normalized by $H/v_{m}$, and $F_{eq}$ is the equilibrium VDF defined as

[TABLE]

where $\rho$ is the number density of gas molecules normalized by $\bar{p}/mRT_{0}$, $\textbf{u}(u_{1},u_{2})$ is the flow velocity normalized by $v_{m}$, and $T$ is the temperature normalized by $T_{0}$. Note that the Knudsen number is defined as $\text{Kn}=\lambda/H$, where the mean free path of gas molecules is related to its shear viscosity $\mu$ as

[TABLE]

When the flow velocity is very small compared to $v_{m}$, we can linearize the VDF about the equilibrium state $f_{eq}$ as follows:

[TABLE]

and the perturbation VDF $h(t,\textbf{x},\textbf{v})$ is governed by the following linearized BGK equation Sharipov and Graur (2012):

[TABLE]

where $\varrho=\int{h}f_{eq}\mathrm{d}\textbf{v}$ is the perturbed number density and $\tau=\int|\textbf{v}|^{2}{h}f_{eq}\mathrm{d}\textbf{v}-\varrho$ is the perturbed temperature. The velocity u is also calculated from $h$, i.e. $u_{1}=\int v_{1}{h}f_{eq}\mathrm{d}\textbf{v}$ and $u_{2}=\int v_{2}{h}f_{eq}\mathrm{d}\textbf{v}$ respectively.

The kinetic equation (7) has to be solved together with the gas kinetic boundary condition which determines the VDF of the reflected gas molecules at the surface in terms of that of the incident molecules. At the inlet and outlet of the computational domain, the periodic boundary condition is used, while at the solid surface, the diffuse boundary condition is used Maxwell (1879). That is, at the moving smooth plate located at $x_{2}=1$, we have

[TABLE]

while at the bottom rough plate located,

[TABLE]

where $v_{n}$ is the normal velocity vector at the solid surface pointing into the fluid regime.

Note that the Weierstrass-Mandelbrot fractal surface is not differentiable, therefore, the rough surface is approximated by the “stair case” in the numerical simulation, as used in the MRT-LBM simulation Deng et al. (2016). Also, we take $U_{0}/v_{m}=1$, so the horizontal flow velocity is calculated as

[TABLE]

while the average mass flow rate, which is normalized by $\bar{p}HU_{0}/RT_{0}$, is given by

[TABLE]

In numerical simulations, the molecular velocity space v is represented by discrete velocities. To capture the discontinuities in the VDF at large Kn, $v_{1}$ and $v_{2}$ are represented by $N_{v}$ non-uniform points in each direction Wu et al. (2014b):

[TABLE]

where most of the discrete velocities are located near $v_{1,2}=0$. We choose $N_{v}=32$ in all our simulations.

The physical space is also discretized into Cartesian grids. In the $x_{1}$ direction, 1500 equidistant points are used. In the $x_{2}$ direction, the rough region $\min(r)\leq{}x_{2}<\max(r)$ is discretized by 60 equidistant points, while the rest is discretized by 100 equidistant points.

As we are interested in the steady-state solution, the derivative with respect to the time is omitted in Eq. (7), and the resulting equation is solved by the following iterative method:

[TABLE]

where the superscript $(n)$ and $(n+1)$ represent two consecutive iteration steps, and the spatial derivatives $\partial h/\partial x$ and $\partial h/\partial y$ are approximated by the second-order upwind finite difference. The iteration is terminated when the relative error in the horizontal flow velocity between two consecutive iteration steps is less than $10^{-5}$.

II.2 Numerical results

We first compare the velocity profiles between the numerical results obtained from the DVM and MRT-LBM, when $\text{Kn}=0.1$ and 1. When both plates are smooth, from Fig. 2 we find that MRT-LBM is accurate at $\text{Kn}=0.1$. However, at Kn=1, the slip velocity at the bottom surface from the MRT-LBM simulation is about 15 % smaller than that from the DVM. Nevertheless, the mass flow rates from both numerical schemes are quite close, as the MRT-LBM predicts a larger velocity close to the top plate. When the bottom surface is rough, the accuracy of the MTR-LBM is reduced significantly. For example, at $\text{Kn}=1$ and $x_{2}=0.1$ (i.e. the rough region), the average horizontal flow velocity from the DVM simulation is two and a half times larger than that from the MRT-LBM. Possible reasons are, (i) the combined bounce-back and specular-reflection boundary condition is only derived in a large flat plate, which may not work in a rough surface, and (ii) for rough surface, the “wall-scaling” of the relaxation time may not work properly. Actually, it has been observed that, the extend NS equations incorporating a wall-scaling model cannot give a correct prediction in velocity profiles even for a simple circular cylinder flow when Kn$>0.4$ Guo et al. (2014).

We then consider the influence of the surface roughness on the mass flow rate. As Ref. Deng et al. (2016), we study how the following ratio

[TABLE]

change with the Knudsen number. Here $Q_{r}$ and $Q_{s}$ are the average mass flow rates through the channels with rough and smooth surfaces, respectively. For the smooth channel with diffuse boundary condition, our numerical simulations show that $Q_{s}$ is always 0.5 in the whole flow regime. This can actually be proven in the continuum and slip flow regimes using the NS equations with slip boundary condition, and in the free-molecular flow regime using the gas kinetic theory. The numerical results for fractal surface with $D=1.5$ and $\epsilon=2\%$ are shown in Fig. 3. Clearly, the accuracy of MRT-LBM decreases as the Knudsen number increases. At $\text{Kn}=10$, $Q^{\ast}$ obtained from MRT-LBM is about 25% smaller than that from the DVM. This error, however, is in fact very large by taking into account that the value of $Q$ only decreases by 20% from the continuum flow regime to the free-molecular flow regime with $\text{Kn}=10$. In this sense, the MRT-LBM overestimates the decrease of $Q$ by twice.

We also checked the mass flow rate when the fractal dimension $D$ in Eq. (1) takes different values. Similar behaviors as that in Fig. 3 are observed (not shown). These examples show that the MRT-LBM is not reliable in simulating Couette flow in rough microchannels.

III Poiseuille flow in porous media

In this section, we further access the accuracy of the MRT-LBM in simulating rarefied Poiseuille flows through three porous medium (Medium 1 with irregular solids; Medium 2 with square solids; Medium 3 with circular solids), as shown in Fig. 4. The vertical size of each porous medium is assumed to be $H$. This kind of research becomes popular due to the shale gas revolution in North America, and LBM is widely applied to find the apparent permeability (a parameter describing how fast the gas can be extracted) of porous media. We will find the variation of the apparent permeability with respect to the Knudsen number using the DVM, and access the accuracy of MRT-LBM.

The DVM to solve the linearized BGK equation for this problem can be found in Ref. Wu et al. (2017) where the periodic boundary condition is used. We have also tested that the pressure boundary condition adopted in Ref. Zhao et al. (2016) leads to the same value of apparent permeability as long as the pressure difference is very small.

According to Ref. Wu et al. (2017), the apparent gas permeability $k_{a}$, which is normalized by $H^{2}$, is calculated by

[TABLE]

where $Q$ is the dimensionless mass flow rate

[TABLE]

with $L$ being length of the computational domain along the flow direction.

Since different porous medium has different distribution of pore radius, it will be useful to define the characteristic length as

[TABLE]

where $\phi$ is the porosity of the porous medium and $k_{\infty}$ is the intrinsic permeability, i.e. the apparent permeability when $\text{Kn}\rightarrow 0$. Hence the equivalent Knudsen number is defined as

[TABLE]

The comparison on the apparent gas permeability of the three porous media at different $Kn^{\ast}$ between the numerical results from the DVM and MRT-LBM is presented in Fig. 5. Our DVM simulations show that, for the porous media 1 and 3, the ratio of the apparent gas permeability to the intrinsic permeability, i.e. $k_{a}/k_{\infty}$, is nearly proportional to the equivalent Knudsen number $\text{Kn}^{\ast}$, while for the porous medium 2, $k_{a}/k_{\infty}$ is a convex function of $Kn^{\ast}$. For the three porous media, surprisingly, the results from the MRT-LBM simulations fall into the same line when $\text{Kn}^{\ast}>1$. It is interesting to note that, for the porous medium 1 both DVM and MRT-LBM results agree with each other very well, while for the porous medium 3, the apparent permeability from DVM is only larger than those from MRT-LBM by about 10 %. However, for the porous media 2, the apparent permeabilities predicted from DVM are nearly twice as large as those from MRT-LBM.

Therefore, in general, the MTR-LBM cannot describe the variation of the apparent permeability with respect to the shape and configuration of pores in porous media, and should not viewed as an accurate numerical method to simulate rarefied gas flows in porous media.

IV Conclusions

In summary, we provided a systematic assessment of the accuracy of the MRT-LBM method in simulating rarefied gas flows in rough microchannels and complex porous media, over a wide range of the Knudsen number. By solving the linearized Bhatnagar-Gross-Krook kinetic equation accurately via the discrete velocity method, we have shown that:

for Couette flow in microchannels, the MRT-LBM can predict the velocity profile with high accuracy in slip flow regime, when the two surfaces are smooth. However, as the Knudsen number increases or surface roughness emerges, the MRT-LBM losses accuracy in predicting the velocity slip and mass flow rate: it even becomes very inaccurate when the Knudsen number is one. 2. 2.

for Poiseuille flow through porous media, the MRT-LBM cannot describe the variation of the apparent permeability with respect to the shape and configuration of pores. In one case considered, it underestimates the apparent permeability by nearly twice.

Therefore, for rarefied flows through rough surface and complex porous media, the MRT-LBM should not be viewed as an accurate method. Possible reasons for the failure of the MRT-LBM include (i) the combined bounce-back and specular-reflection boundary condition is only derived in a large flat plate, which may not work in a rough surface, and (ii) for rough surface, the “wall-scaling” of the relaxation time may not work properly. Our simulation results from the discrete velocity method could serve as benchmarking cases for future development of the LBM for modeling and simulation of rarefied gas flows in complex geometry.

L.W. and H.L. acknowledge the financial support of the joint project from the Royal Society of Edinburgh and National Natural Science Foundation of China.