Discrete effect on the anti-bounce-back boundary condition of lattice Bhatnagar-Gross-Krook model for convection-diffusion equations

TL;DR

This paper investigates the discrete effects of boundary conditions in the lattice Boltzmann method for convection-diffusion equations, demonstrating how to eliminate numerical slip by adjusting relaxation time and distance ratio.

Contribution

It introduces a non-halfway anti-bounce-back boundary condition analysis within the BGK model, showing how to eliminate numerical slip by tuning parameters, extending applicability to other LB models.

Findings

Numerical slip depends on relaxation time and distance ratio.

Relaxation time can be freely adjusted to eliminate slip.

Results validated through simulations with straight and curved boundaries.

Abstract

The discrete effect on the boundary condition has been a fundamental topic for the lattice Boltzmann method in simulating heat and mass transfer problems. In previous works based on the halfway anti-bounce-back (ABB) boundary condition for convection-diffusion equations (CDEs), it is reported that the discrete effect cannot be commonly removed in the Bhatnagar-Gross-Krook (BGK) model except for a special value of relaxation time. Targeting this point in the present paper, we still proceed within the framework of BGK model for two-dimensional CDEs, and analyze the discrete effect on a non-halfway ABB boundary condition which incorporates the effect of the distance ratio. By analyzing an unidirectional diffusion problem with a parabolic distribution, the theoretical derivations with three different discrete velocity models show that the numerical slip is a combined function of the relaxation time and the distance ratio. Different from previous works, we definitely find that the relaxation time can be freely adjusted by the distance ratio in a proper range to eliminate the numerical slip. Some numerical simulations are carried out to validate the theoretical derivations, and the numerical results for the cases of straight and curved boundaries confirm our theoretical analysis. Finally, it should be noted that the present analysis can be extended from the BGK model to other lattice Boltzmann (LB) collision models for CDEs, which can broaden the parameter range of the relaxation time to approach 0.5.

Tables

Table 1. Table 1: Numerical slip ϕ s subscript italic-ϕ 𝑠 \phi_{s} and relaxation parameter corresponding to ϕ s = 0 subscript italic-ϕ 𝑠 0 \phi_{s}=0 in the present work within the BGK model and those ShuC16 within the MRT model. The listed results are based on the D2Q4, D2Q5 and D2Q9 lattice models. The relaxation parameter s 1 = 1 / τ ϕ subscript 𝑠 1 1 subscript 𝜏 italic-ϕ s_{1}=1/\tau_{\phi} is related to the diffusion coefficient in ShuC16 , while s 2 subscript 𝑠 2 s_{2} is served as the free relaxation parameter.

Discrete lattice model	${\varphi }_s$
	BGK(Present)	MRT(Ref. ShuC16 )
D2Q4	$\Delta \varphi \frac{{\delta }_x^2}{H^2}\frac{4{\tau }_{\varphi }^2+4\gamma \left(2{\tau }_{\varphi }-1\right)-8{\tau }_{\varphi }-4{\gamma }^2+3}{4}$	$\Delta \varphi \frac{{\delta }_x^2}{H^2}\frac{2-s_1-s_2}{2s_1{s}_2}$
D2Q5	$\Delta \varphi \frac{{\delta }_x^2}{H^2}\frac{6{\tau }_{\varphi }^2+5\gamma \left(2{\tau }_{\varphi }-1\right)-11{\tau }_{\varphi }-5{\gamma }^2+4}{5}$	$\Delta \varphi \frac{{\delta }_x^2}{H^2}\frac{s_1{s}_2-12\left(s_1+s_2\right)+24}{20s_1{s}_2}$
D2Q9	$\Delta \varphi \frac{{\delta }_x^2}{H^2}\frac{8{\tau }_{\varphi }^2+6\gamma \left(2{\tau }_{\varphi }-1\right)-14{\tau }_{\varphi }-6{\gamma }^2+5}{6}$	$\Delta \varphi \frac{{\delta }_x^2}{H^2}\frac{s_1{s}_2-8\left(s_1+s_2\right)+16}{12s_1{s}_2}$
Discrete lattice model	Relaxation parameter(${\varphi }_s=0$)
	BGK(Present)	MRT(Ref. ShuC16 )
D2Q4	${\tau }_{\varphi }=\frac{2\left(1-\gamma \right)+\sqrt{8{\gamma }^2-4\gamma +1}}{2}$	$s_2=2-s_1$
D2Q5	${\tau }_{\varphi }=\frac{11-10\gamma +\sqrt{5\left(44{\gamma }^2-20\gamma +5\right)}}{12}$	$s_2=\frac{12\left(s_1-2\right)}{s_1-12}$
D2Q9	${\tau }_{\varphi }=\frac{7-6\gamma +\sqrt{3\left(28{\gamma }^2-12\gamma +3\right)}}{8}$	$s_2=\frac{8\left(s_1-2\right)}{s_1-8}$

Table 2. Table 2: Relaxation time τ ϕ subscript 𝜏 italic-ϕ \tau_{\phi} versus distance ratio γ 𝛾 \gamma for ϕ s = 0 subscript italic-ϕ 𝑠 0 \phi_{s}=0 with the D2Q4, D2Q5 and D2Q9 lattice models.

$\gamma$	${\tau }_{\varphi }$
	D2Q4	D2Q5	D2Q9
0.2	$\frac{8+\sqrt{13}}{10}\approx 1.1606$	$\frac{45+\sqrt{345}}{60}\approx 1.0596$	$\frac{29+\sqrt{129}}{40}\approx 1.0090$
0.5	$1.0000$	$\frac{6+\sqrt{30}}{12}\approx 0.9564$	$\frac{2+\sqrt{3}}{4}\approx 0.9330$
0.8	$\frac{2+\sqrt{73}}{10}\approx 1.0544$	$\frac{15+\sqrt{2145}}{60}\approx 1.0219$	$\frac{11+\sqrt{849}}{40}\approx 1.0034$
1.2	$\frac{-2+\sqrt{193}}{10}\approx 1.1892$	$\frac{-5+\sqrt{5545}}{60}\approx 1.1578$	$\frac{-1+3\sqrt{241}}{40}\approx 1.1393$
1.5	$\frac{\sqrt{13}-1}{2}\approx 1.3028$	$\frac{\sqrt{370}-4}{12}\approx 1.2696$	$1.2500$

Equations

$${\partial }_t\varphi +\nabla \cdot (\varphi \mathbf{u})=\nabla \cdot (D\nabla \varphi )+R,$$

$$f_i(\mathbf{x}+{\mathbf{c}}_i{\delta }_t,t+{\delta }_t)-f_i(\mathbf{x},t)=-\frac{1}{{\tau }_{\varphi }}\left[f_i(\mathbf{x},t)-f_i^{(eq)}(\mathbf{x},t)\right]+{\delta }_t(1-\frac{1}{2{\tau }_{\varphi }})R_i(\mathbf{x},t),$$

$$f_i^{(eq)}(\mathbf{x},t)$$

$$R_i(\mathbf{x},t)$$

$$\varphi (\mathbf{x},t)=\sum _i{f}_i(\mathbf{x},t)+\frac{{\delta }_t}{2}R(\mathbf{x},t).$$

$$\text{Collision}:f_i^{\ast }(\mathbf{x},t)=f_i(\mathbf{x},t)-\frac{1}{{\tau }_{\varphi }}\left[f_i(\mathbf{x},t)-f_i^{(eq)}(\mathbf{x},t)\right]+{\delta }_t(1-\frac{1}{2{\tau }_{\varphi }})R_i(\mathbf{x},t),$$

$$\text{Streaming}:f_i(\mathbf{x}+{\mathbf{c}}_i{\delta }_t,t+{\delta }_t)=f_i^{\ast }(\mathbf{x},t),$$

$$D\frac{{\partial }^2\varphi }{\partial y^2}+R=0,$$

$$\varphi (x,y=0)={\varphi }_B,\varphi (x,y=H)={\varphi }_T,$$

$$R=2D\Delta \varphi /H^2,\Delta \varphi ={\varphi }_T-{\varphi }_B,$$

$$\varphi (y)={\varphi }_B+\Delta \varphi \frac{y}{H}(2-\frac{y}{H}).$$

$$f_2=\left(1-\frac{1}{2\gamma }\right)f_2^{\ast }-\frac{1}{2\gamma }{f}_4^{\ast }+\frac{{\omega }_2}{\gamma }{\varphi }_B,$$

$$f_5=\left(1-\frac{1}{2\gamma }\right)f_5^{\ast }-\frac{1}{2\gamma }{f}_7^{\ast }+\frac{{\omega }_5}{\gamma }{\varphi }_B,$$

$$f_6=\left(1-\frac{1}{2\gamma }\right)f_6^{\ast }-\frac{1}{2\gamma }{f}_8^{\ast }+\frac{{\omega }_6}{\gamma }{\varphi }_B.$$

$${\varphi }_1=\frac{\gamma +1}{\gamma }{\varphi }_0-\frac{{\varphi }_B}{\gamma }-\frac{4{\tau }_{\varphi }^2+8(\gamma -1){\tau }_{\varphi }+3}{2\gamma (2{\tau }_{\varphi }-1)}{\delta }_{t}R$$

$${\varphi }_1=\frac{\gamma +1}{\gamma }{\varphi }_0-\frac{{\varphi }_B}{\gamma }-\frac{6{\tau }_{\varphi }^2+(10\gamma -11){\tau }_{\varphi }+4}{2\gamma (2{\tau }_{\varphi }-1)}{\delta }_{t}R$$

$${\varphi }_1=\frac{\gamma +1}{\gamma }{\varphi }_0-\frac{{\varphi }_B}{\gamma }-\frac{8{\tau }_{\varphi }^2+2(6\gamma -7){\tau }_{\varphi }+5}{2\gamma (2{\tau }_{\varphi }-1)}{\delta }_{t}R$$

$${\varphi }_j={\varphi }_B+\Delta \varphi \frac{y_j}{H}\left(2-\frac{y_j}{H}\right)+{\varphi }_s,$$

$${\varphi }_s=\frac{\Delta \varphi }{4}\frac{{\delta }_x^2}{H^2}\left[4{\tau }_{\varphi }^2+4\gamma (2{\tau }_{\varphi }-1)-8{\tau }_{\varphi }-4{\gamma }^2+3\right],$$

$$\gamma =\frac{1}{2}\Rightarrow {\varphi }_s=\frac{\Delta \varphi }{4}\frac{{\delta }_x^2}{H^2}\left(4{\tau }_{\varphi }^2-4{\tau }_{\varphi }\right).$$

$${\varphi }_s=\frac{\Delta \varphi }{5}\frac{{\delta }_x^2}{H^2}\left[6{\tau }_{\varphi }^2+5\gamma (2{\tau }_{\varphi }-1)-11{\tau }_{\varphi }-5{\gamma }^2+4\right],$$

$$\gamma =\frac{1}{2}\Rightarrow {\varphi }_s=\frac{\Delta \varphi }{5}\frac{{\delta }_x^2}{H^2}\left(6{\tau }_{\varphi }^2-6{\tau }_{\varphi }+\frac{1}{4}\right).$$

$${\varphi }_s=\frac{\Delta \varphi }{6}\frac{{\delta }_x^2}{H^2}\left[8{\tau }_{\varphi }^2+6\gamma (2{\tau }_{\varphi }-1)-14{\tau }_{\varphi }-6{\gamma }^2+5\right],$$

$$\gamma =\frac{1}{2}\Rightarrow {\varphi }_s=\frac{\Delta \varphi }{6}\frac{{\delta }_x^2}{H^2}\left(8{\tau }_{\varphi }^2-8{\tau }_{\varphi }+\frac{1}{2}\right).$$

$${\varphi }_s=0\Rightarrow {\tau }_{\varphi }=\frac{2(1-\gamma )+\sqrt{8{\gamma }^2-4\gamma +1}}{2},$$

$$\gamma =\frac{1}{2},{\varphi }_s=0\Rightarrow {\tau }_{\varphi }=1.$$

$${\varphi }_s=0\Rightarrow {\tau }_{\varphi }=\frac{11-10\gamma +\sqrt{5(44{\gamma }^2-20\gamma +5)}}{12},$$

$$\gamma =\frac{1}{2},{\varphi }_s=0\Rightarrow {\tau }_{\varphi }=\frac{6+\sqrt{30}}{12}.$$

$${\varphi }_s=0\Rightarrow {\tau }_{\varphi }=\frac{7-6\gamma +\sqrt{3(28{\gamma }^2-12\gamma +3)}}{8},$$

$$\gamma =\frac{1}{2},{\varphi }_s=0\Rightarrow {\tau }_{\varphi }=\frac{2+\sqrt{3}}{4}.$$

$${\delta }_t=\eta {\delta }_x^2,\eta =\frac{\chi ({\tau }_{\varphi }-\frac{1}{2})}{D},$$

$$\varphi (r)=\frac{{\varphi }_a\text{ln}(R_b/r)+{\varphi }_b\text{ln}(r/R_a)}{\text{ln}(R_b/R_a)},R_a\le r\le R_b.$$

Full text

††preprint: APS/123-QED

I Introduction

In the past couples of decades, the lattice Boltzmann method (LBM) has been gradually developed as an effective and powerful technique for a wide range of application areas Succi01 ; Guob13 , such as single-phase flows, multiphase flows, microgaseous flows and porous flows ChenSY98 ; Aidun10 ; HuangHB15 ; ZhangJ11 ; Liang19 . Unlike the conventional computational methods, the LBM solves the discrete Boltzmann equation instead of the macroscopic continuum equations. The kinetic nature of the LBM possesses several attractive features in flow simulations, such as simple program, intrinsically parallel computation, and easy boundary treatment. Among the other successful extensions, the LBM has also been adapted to solve convection-diffusion equations (CDEs), which are commonly encountered in studying heat and mass transfer associated with fluid flows. So far, there have been many LB models proposed for CDEs Sman00 ; Huang14 ; WangL18 . More detailed reviews about these works can be found in Refs. ShiB09 ; Yoshida10 ; Chai13 ; ZhangSub .

To completely solve CDEs by the LBM, apart from the numerical algorithm for the LB equation (LBE), the boundary condition should also be specified for the unknown distribution functions at boundary nodes (i.e., lattice nodes nearest to the physical boundary). It is a critical issue and has attracted increasing researchers’ efforts towards accurate boundary treatments. In several recent publications ZhangT12 ; ChenQ13 ; HuangJ15 ; Kruger17 , the reader can trace some existing LBM boundary conditions such as the ABB scheme Ginzburg05 ; Li17 and the non-equilibrium extrapolation scheme Chai16 ; GuoZ02 . The terminology of ABB is in contrast to the bounce-back (BB) scheme for fluid flows. That is, the outgoing population reflects back in the opposite direction with the BB scheme, while it changes its sign with the ABB scheme Ginzburg17 . As have recognized in the boundary conditions of LBM for flow simulations, it is known that the discrete effect on the boundary condition also must be minimized to derive correct results for CDEs ZhangT12 ; Dubois10 ; ShuC16 . However, there has not been extensive investigations on this topic as those for the fluid flow simulations. Based on the developed Taylor expansion method Dubois07 ; Dubois08 , Dubois et al. Dubois10 analyzed the ABB boundary condition within the framework of multiple-relaxation-time (MRT) model for one-dimensional diffusion equation with the Dirichlet boundary condition. They demonstrated that the halfway ABB (HABB) boundary condition can be accurate up to order two in space under a specific combination of the relaxation rates. Within the BGK model framework, Zhang et al. ZhangT12 proposed a HABB boundary condition for CDEs, and also analyzed the discrete effect of their boundary condition. For the diffusion in Couette flow with wall injection, they derived mathematically that the concentration jump or the numerical slip is related with the relaxation time and has a second-order dependence with the lattice spacing. It is also shown that the numerical slip cannot commonly be removed in the BGK model. As for the discrete effect of the HABB boundary condition, Cui et al. ShuC16 revisit this topic based on the MRT model with three discrete lattice models, and derive the numerical slip relating with two relaxation rates and the square of lattice spacing. Their theoretical analysis and numerical results show that the discrete effect on the HABB boundary condition can be removed owing to the free relaxation parameter $s_{2}$ in the MRT model, while it cannot be eliminated except for a special value of the relaxation time in the BGK model. However, we note that the boundary condition in the above works is concentrated to the halfway boundary scheme, which intrinsically disregards the possible degree of freedom from the wall arrangements between lattice nodes.

Actually, in the boundary conditions for CDEs, the wall can be located between two lattice nodes with an arbitrary but not only halfway intersection distance Li17 ; HuangJ15 ; ChenQ13 ; Dubois19 . This means that if the wall location is embodied in the boundary condition, it may appear as a free parameter besides the relaxation time in the derived numerical slip. Therefore, it naturally brings out a fundamental question about the discrete effect of the boundary condition: whether the numerical slip can be eliminated while not limited at a special relaxation time in the BGK model. To our knowledge, no publications have been reported on this topic. In this work, we will analyze the discrete effect on the non-halfway ABB (NHABB) boundary condition for CDEs within the framework of BGK model. The boundary condition proposed in Ref. HuangJ15 is adopted here for its locality and ability to adjust the wall location arbitrarily between lattice nodes. And importantly, we will show how to choose the relaxation time to eliminate the discrete effect freely by tuning the parameter of wall location. From this point, in addition to resorting to other LBE models (e.g., MRT model) for more degrees of freedom, the present work reveals another way to eliminate the discrete effect on boundary condition of the BGK model for CDEs.

The paper is organized as follows. In Sec. II, the BGK-LBE for the CDE with a source term is presented. Sec. III is devoted to analyzing the discrete effect of the halfway and non-halfway ABB boundary conditions. In Sec. IV, some numerical experiments and discussions are given, and followed by some conclusions finally presented in Sec. V.

II Lattice Bhatnagar-Gross-Krook model for convection-diffusion equations

In this work, our analyses are specially focused on the BGK model for the convection-diffusion equation. For the two-dimensional case, the CDE with a source term reads

[TABLE]

where $\phi$ is the scalar variable as a function of time and space, $D$ is the diffusion coefficient, $\bm{u}=(u_{x},u_{y})^{T}$ is the convection velocity with $T$ denoting the transposition operator, and $R(\bm{x},t)$ is the source term. The BGK-LBE to sovle the CDE (1) is written as follows

[TABLE]

where $\{f_{i}(\bm{x},t):i=0,1,\cdots,b-1\}$ are the distribution functions associated with the discrete velocities $\{\bm{c}_{i}:i=0,1,\cdots,b-1\}$ at position $\bm{x}$ and time $t$, $\tau_{\phi}$ is the relaxation time, $\delta_{t}$ is the evolution time increment; $f_{i}^{(eq)}(\bm{x},t)$ is the equilibrium distribution function, and $R_{i}(\bm{x},t)$ is the discrete source term, which are respectively defined as

[TABLE]

where $\omega_{i}$ is the weight coefficient, and $c_{s}$ is the sound speed.

The discrete velocity set $\bm{c}_{i}$ is subjected to the $\text{D}n\text{Q}b$ ($\text{D}n\text{Q}b$ denotes $b$ velocity directions in $n$D space) lattice models reported in the literature Qian92 . In this work, the discrete effect of the ABB boundary condition is inspected with three discrete lattice models. As adopted in Ref. ShuC16 for the MRT collision model, the D2Q4, D2Q5 and D2Q9 models are also considered for subsequent analysis connected with the BGK model. The corresponding parameters for these three models are given as follows: for the D2Q4 model, $\{\bm{c}_{i}:i=1,2,3,4\}=\{(\pm 1,0)c,(0,\pm 1)c\}$, $\omega_{1-4}=\frac{1}{4}$, and $c_{s}^{2}=\frac{1}{2}c^{2}$; for the D2Q5 model, $\{\bm{c}_{i}:i=0,1,2,3,4\}=\{(0,0)c,(\pm 1,0)c,(0,\pm 1)c\}$, $\omega_{0-4}=\frac{1}{5}$, and $c_{s}^{2}=\frac{2}{5}c^{2}$; for the D2Q9 model, $\{\bm{c}_{i}:i=0,1,\cdots,8\}=\{(0,0)c,(\pm 1,0)c,(0,\pm 1)c,(\pm 1,\pm 1)c\}$, $\omega_{0}=\frac{4}{9}$, $\omega_{1-4}=\frac{1}{9}$, $\omega_{5-8}=\frac{1}{36}$, and $c_{s}^{2}=\frac{1}{3}c^{2}$; where $c=\delta_{x}/\delta_{t}$ is the lattice speed with $\delta_{x}$ the lattice spacing.

The macroscopic variable $\phi$ is determined by the distribution functions as

[TABLE]

With this definition, the CDE with a source term, Eq. (1), can be recovered from the BGK model through the Chapman-Enskog analysis WangL18 ; Chai13 . Also, the diffusion coefficient can be derived and determined by the relaxation time $\tau_{\phi}$ as $D=c_{s}^{2}(\tau_{\phi}-\frac{1}{2})\delta_{t}$.

Numerically, the evolution of the BGK-LBE (2) is implemented via two steps, i.e., the collision and streaming step:

[TABLE]

where $f_{i}^{*}(\bm{x},t)$ is the postcollision distribution function.

III Discrete effect on the anti-bounce-back boundary condition of BGK model for CDE

As mentioned previously, the discrete effect on the ABB boundary condition for CDEs is concentrated on the halfway scheme in the existing analysis ZhangT12 ; ShuC16 . However, the discrete effect is unclear when it is affected by the wall location between lattice nodes. To resolve this gap, we will restrict within the framework of BGK model to analyze the discrete effect of the NHABB boundary condition. For clarity of illustration, our analysis is based on the problem used in Ref. ShuC16 within the MRT framework. The considered problem is an unidirectional and time-independent diffusion in a straight channel (see Fig. 1) in which $u_{x}$ is constant, $u_{y}=0$, and $\partial_{x}\phi=0$ for any scalar variable $\phi$.

For the constant $\phi_{B}$ and $\phi_{T}$ corresponding to the bottom and top walls (i.e., the Dirichlet boundary condition), the problem can be described by the following equations

[TABLE]

where $H$ is the height of the channel. As the source term $R$ is further defined by

[TABLE]

we can obtain the analytical solution to this simple problem

[TABLE]

When the BGK model (2) is implemented to solve the above diffusion problem, after a time step $\delta_{t}$, the unknown distribution functions (see Fig. 2) should be specified by proper boundary conditions of the LBM.

To this end, we adopt the NHABB boundary condition proposed by Huang and Yong HuangJ15 using the asymptotic analysis technique. As displayed in Fig. 2, the wall boundary, say the bottom wall, is located away from its nearest inner lattice nodes with the distance of $\gamma\delta_{x}$. When the location of bottom wall is adjustable, the distance ratio $\gamma$ can be considered as a free parameter, which is conventionally used in the range of $0\leq\gamma\leq 1$ ZhangSub ; ChenQ13 ; HuangJ15 . However, as will be shown later, we would like to note that the distance ratio $\gamma$ would not be intuitively limited in $0\leq\gamma\leq 1$, but can be larger than unity to derive accurate results. With this point, the unknown distribution functions at the layer $j=0$ are then determined by the following equations.

D2Q4 or D2Q5 lattice model:

[TABLE]

D2Q9 lattice model:

[TABLE]

One can see that if the parameter $\gamma=\frac{1}{2}$, this boundary condition will reduce to the HABB scheme ZhangT12 ; ShuC16 . It should be noted that the above boundary condition is a local scheme and has second-order accuracy for the case of straight walls HuangJ15 . Additionally, we would like to point out that following the procedures presented in Ref. Zhao17b for Dirichlet boundary condition of the Navier-Stokes equations, the above boundary condition can also be obtained by the Maxwell iteration method Yong16 ; Zhao17a with the diffusive scaling $\delta_{t}=\eta\delta_{x}^{2}$ and an adjustable parameter $\eta$.

Based on the adopted boundary schemes and the assumptions for the diffusion problem, one can follow the derivations in Refs. ZhangT12 ; ShuC16 to derive that

[TABLE]

for the D2Q4 lattice model, and

[TABLE]

for the D2Q5 lattice model, and

[TABLE]

for the D2Q9 lattice model. Here, $\phi_{0}$ and $\phi_{1}$ are the scalar variables at the layer of $j=0$ and $j=1$.

During the above derivations, we can also deduce that the numerical scalar variable $\phi_{j}$ satisfies $D(\phi_{j+1}-2\phi_{j}+\phi_{j-1})=-R\delta_{x}^{2}$, where the diffusion coefficient $D$ is also given by $D=c_{s}^{2}(\tau_{\phi}-\frac{1}{2})\delta_{t}$. Clearly, this is the central finite-difference discretization of Eq. (7a), meaning that the BGK-LBE is an equivalent solver for the CDEs. However, due to the discrete effect from the boundary condition, the LB results will deviate from the analytical solution to the problem [Eq. (9)]. As a result, the solution of the BGK model with the NHABB boundary condition can be expressed as

[TABLE]

where $y_{j}=(j+\gamma)\delta_{x}$, and $\phi_{s}$ is the numerical slip originated from the discrete effect of the boundary condition. By substituting $\phi_{0}$ and $\phi_{1}$ from Eq. (15) respectively into Eqs. (12), (13) and (14), we can obtain the numerical slips $\phi_{s}$ from the D2Q4, D2Q5 and D2Q9 lattice models.

D2Q4 lattice model:

[TABLE]

D2Q5 lattice model:

[TABLE]

D2Q9 lattice model:

[TABLE]

From each of the above equations, one can find that the HABB and NHABB boundary conditions generate a nonzero numerical slip $\phi_{s}$, which has second-order accuracy in space owing to the term of $\delta_{x}^{2}/H^{2}$. It is noted that the results of $\phi_{s}$ for the HABB boundary condition ($\gamma=1/2$) [Eqs. (16b), (17b) and (18b)] here are identical to those of the BGK model given in Ref. ShuC16 . However, the numerical slip $\phi_{s}$ of the halfway boundary condition is not available for the non-halfway boundary condition. Due to the fixed location of wall with $\gamma=1/2$, $\phi_{s}$ of the HABB boundary condition is only related with $\tau_{\phi}$. This indicates that the discrete effect of the HABB boundary condition always exists in the BGK model unless a special relaxation time is used ZhangT12 ; ShuC16 . In contrast, owing to the adjustable distance ratio $\gamma$ as revealed above, $\phi_{s}$ of the NHABB boundary condition [Eqs. (10) and (11)] is dependent with the relaxation time $\tau_{\phi}$ and the distance ratio $\gamma$. Thus, the relaxation time $\tau_{\phi}$ has more degree of freedom to minimize the discrete effect on the boundary condition. The above results inspire us that the numerical slip $\phi_{s}$ of the BGK model could be eliminated freely by the relaxation time $\tau_{\phi}$ with the help of the free parameter $\gamma$.

Now let us focus on how to choose the relaxation time $\tau_{\phi}$ tuned by the distance ratio $\gamma$ to guarantee $\phi_{s}=0$. Mathematically, this can be done by solving the quadratic equation $\phi_{s}=0$ from Eqs. (16a), (17a) and (18a) respectively for the D2Q4, D2Q5 and D2Q9 discrete lattice model. Because of the stability condition as well as the positivity of diffusivity, there is only one root of $\tau_{\phi}$ that is determined by $\gamma$ from Eqs. (19a), (20a) and (21a). For the HABB boundary scheme ($\gamma=1/2$), the corresponding relaxation time $\tau_{\phi}$ is obtained by Eqs. (19b), (20b) and (21b).

D2Q4 lattice model:

[TABLE]

D2Q5 lattice model:

[TABLE]

D2Q9 lattice model:

[TABLE]

For the case of $\gamma=1/2$, the corresponding values of $\tau_{\phi}$ are fixed and identical to those reported in Ref. ShuC16 . However, with the adopted NHABB boundary condition, the relaxation time $\tau_{\phi}$ is related with the distance ratio $\gamma$, and hence can be freely tuned by $\gamma$ to ensure $\phi_{s}=0$. To see this more clearly, the dependence of $\tau_{\phi}$ on $\gamma$ as given above for $\phi_{s}=0$ is shown in Fig 3.

Take the D2Q4 model as an example. It is seen that as the distance ratio $\gamma$ increases, $\tau_{\phi}$ can change continuously to fulfil $\phi_{s}=0$, while as $\gamma=1/2$ for the HABB boundary condition, $\phi_{s}=0$ determines the merely fixed $\tau_{\phi}=1$. When the distance ratio $\gamma$ varies in the region of $0\leq\gamma\leq 1$, the relaxation time $\tau_{\phi}$ will take values limitedly between 1 and 1.5. To achieve a wider parameter range for $\tau_{\phi}$, the distance ratio $\gamma$ should not be confined to the region of $0\leq\gamma\leq 1$. This point will be examined in the subsequent numerical examples, and the computations therein reveal that reasonable results can be also obtained as $\gamma$ is beyond 1. From the figure, it is also find that there are two values of $\gamma$ with $\gamma\leq 2$ corresponding to the same $\tau_{\phi}$ in the range of $1\leq\tau_{\phi}\leq 1.5$, while there is only one $\gamma$ with $\gamma>2$ corresponding to a certain $\tau_{\phi}>1.5$. Similar results stored in Eqs. (20a) and (21a) are also observed for the D2Q5 and D2Q9 lattice models. It should be noted that as $\gamma\rightarrow 0$, the boundary conditions (10) and (11) may lose the numerical stability since the included term $\frac{1}{\gamma}$ will become very large. This will be also affirmed in the subsequent numerical examples. Therefore, the relaxation time $\tau_{\phi}$ should be chosen carefully to avoid very small values of $\gamma$ in the computations. However, we would note that such limitation of parameter range in $\gamma$ may be remedied through recomposing the distribution functions and their coefficients in the adopted boundary condition ZhangSub .

From the above derivations, it is clear that due to the distance ratio $\gamma$, the numerical slip $\phi_{s}$ can be theoretically eliminated within the framework of BGK model. As for the MRT model, it has been commonly recognized that the numerical slip $\phi_{s}$ can be overcome owing to its multiple relaxation parameters ShuC16 . For an explicit comparison of the two model frameworks, Table 1 presents the numerical slip $\phi_{s}$ and the relaxation parameter corresponding to $\phi_{s}=0$ derived in this work together with those deduced with the MRT model in Ref. ShuC16 .

One can find that the relaxation parameter under $\phi_{s}=0$ is related with another relaxation rate in the MRT model, while it is related with the distance ratio $\gamma$ here in the BGK model. Based on this, we note that the elimination of numerical slip in the MRT model is ascribed to the degree of freedom from the relaxation parameter of the evolution equation, while in the BGK model here, the degree of freedom is from the wall location of the NHABB boundary condition.

IV Numerical results and discussions

To examine the above theoretical analysis, the BGK-LBE with the halfway and non-halfway ABB boundary conditions HuangJ15 are executed in the numerical simulations. The diffusion problems considered here are the same as those adopted in Ref. ShuC16 . In the following simulations, the lattice spacing $\delta_{x}$ is determined by $\delta_{x}=H/(M+2\gamma)$ with $\gamma$ the distance ratio, where $M$ is the grid number between the walls in the vertical direction. The distance ratio $\gamma$ is set as an input variable, and other related parameters are given by

[TABLE]

where $\chi$ is a model-dependent constant defined by $c_{s}^{2}=\chi c^{2}$, and equals to $1/2,2/5,1/3$ respectively for the D2Q4, D2Q5 and D2Q9 lattice model.

IV.1 Unidirectional diffusion in a straight channel

The first problem is shown in Fig. 1, where $H=1$, $u_{x}=0.1$, and the diffusion coefficient $D=0.1$. The periodic boundary condition is applied to the inlet and outlet of the channel, and the ABB boundary condition with the distance ratio $\gamma$ is applied to the top and bottom walls. In Fig. 4, the simulated results of numerical slip, normalized by the results at the case of $\gamma=1/2$, are presented as a function of $\gamma$ at $\tau_{\phi}=1.2$ and $M=15$. The normalized theoretical results [Eqs. (16a), (17a) and (18a)] are also included for comparison. Clearly, the numerical predictions are well consistent with the theoretical derivations for the three discrete lattice models. In particular, the unambiguous agreement between such two results is observed when $\gamma$ is greater than $1$ up to $2$, and even at $\gamma=3$ (the results are not shown here).

Moreover, we find that the computations will break down as $\gamma$ decreases to 0.1. These twofold results verify the aforementioned statements about the choice of distance ratio $\gamma$ in the boundary conditions. Additionally, as the distance ratio $\gamma$ increases, it is observed that the numerical slip varies increasingly first and then decreasingly after one certain $\gamma$ due to its quadratic function as derived above. It is noted that similar results as shown in Fig. 4 can also be obtained at other relaxation times.

The relations between $\gamma$ and $\tau_{\phi}$ are next examined especially for the numerical slip $\phi_{s}=0$. To this end, simulations with different grid sizes are carried out for two different values of $\tau_{\phi}$ at each of two distance ratios $\gamma=0.6$ and $\gamma=1.5$. One relaxation time is given by $\gamma$ to satisfy $\phi_{s}=0$ as derived above, while the other relaxation time (e.g., $\tau_{\phi}=3.0$) is not the case. For $\gamma=0.6$ and $\gamma=1.5$ considered here, the corresponding relaxation times $\tau_{\phi}$ to ensure $\phi_{s}=0$ can be obtained from Eq. (19a) as $\tau_{\phi}=\frac{4+\sqrt{37}}{10}$ and $\tau_{\phi}=\frac{\sqrt{13}-1}{2}$ for the D2Q4 model, Eq. (20a) as $\tau_{\phi}=\frac{25+\sqrt{1105}}{60}$ and $\tau_{\phi}=\frac{\sqrt{370}-4}{12}$ for the D2Q5 model, and Eq. (21a) as $\tau_{\phi}=\frac{19}{20}$ and $\tau_{\phi}=\frac{5}{4}$ for the D2Q9 model. Figs. 5, 6 and 7 respectively present the simulated results of the D2Q4, D2Q5 and D2Q9 lattice models.

As clearly shown in the figures, only when the relaxation time $\tau_{\phi}$ is determined by $\gamma$ while guaranteeing $\phi_{s}=0$, the results of the BGK model agree well with the analytical solution even with four grid points. However, if this requirement is not satisfied, clear discrepancies between the LBE results and the analytical solution can be observed even at $M=15$. Furthermore, to quantify the differences between such two results, the relative errors in $L^{1}$ norm are evaluated under different grid sizes. In Fig. 8, the relative errors $E(\phi)$ of the scalar variable $\phi$ are plotted against the grid size $M$.

It is clearly shown that as compared with the case of $\tau_{\phi}=3.0$, the obviously large errors are significantly reduced near zero when $\tau_{\phi}$ is given by $\gamma$ to ensure $\phi_{s}=0$. This further strengthens and supports our theoretical derivations.

The results exhibited in Fig. 4 has shown that the numerical slip of the NHABB boundary condition is different from that of the HABB boundary condition ($\gamma=1/2$). This indicates that the relaxation time $\tau_{\phi}$ derived from $\phi_{s}=0$ for $\gamma=1/2$ [Eqs. (19b), (20b) and (21b)] must be amended for the NHABB boundary condition to derive accurate results. In what follows, the values of relaxation time $\tau_{\phi}$ are inspected versus different values of $\gamma$ under the numerical slip $\phi_{s}=0$. In Tab. 2, the approximations of the calculated values of $\tau_{\phi}$ from Eqs. (19), (20) and (21) are listed against $\gamma$.

As seen from the table, the relaxation times $\tau_{\phi}$ are approximate to 1.0 at $\gamma=0.5$ in the D2Q5 and D2Q9 lattice model. Thus, as have pointed out in Ref. ShuC16 , satisfactory results can be usually obtained as $\tau_{\phi}=1.0$ even if the numerical slip is not strictly removed. However, we should note that this result is derived and valid for the halfway boundary scheme (i.e., $\gamma=0.5$). In fact, when the distance ratio $\gamma$ deviates away from 0.5, e.g., $\gamma=1.5$, the relaxation time $\tau_{\phi}$ is definitely larger than 1.0, which is also reflected in Fig. 3. This clearly indicates that when $\gamma$ varies away from 0.5, accurate results cannot be achieved any longer if $\tau_{\phi}$ still remains at 1.0. In other words, to derive accurate results ($\phi_{s}=0$), the relaxation time $\tau_{\phi}$ must be adjusted with $\gamma$.

IV.2 Diffusion between two concentric cylinders

In this section, we investigate a more complex problem, i.e., the steady diffusion between two concentric circular cylinders, as shown in Fig. 9.

The inner cylinder has radius of $R_{a}$ and boundary value $\phi_{a}$, while the outer cylinder has radius of $R_{b}$ and boundary value $\phi_{b}$. The outer cylinder boundary is separated from the square region with a distance of $\gamma_{0}\delta_{x}$. There is no source term for diffusions between the two concentric cylinders. From Eq. (1) in polar coordinates, the analytical solution to this problem can be solved and read as Carslaw13

[TABLE]

In the simulations, the two cylinders are positioned at the center of a square region with length $L=1.0$. The radius ratio of the two cylinders is $R_{a}:R_{b}=1:2$, the diffusion coefficient is set to $D=0.001$ and the boundary values are $\phi_{a}=0.0,\phi_{b}=1.0$. Unlike the previous problem with straight walls, the curved boundary geometries herein may bring different distance ratios, denoted by $\gamma_{in}$ and $\gamma_{out}$ for the boundary nodes respectively of the inner and outer cylinders. To have an unique relaxation time $\tau_{\phi}$ in ensuring $\phi_{s}=0$ (see Eqs. (19a), (20a) and (21a)), we approximate the distance ratio $\gamma$ by the average values of all $\gamma_{in}$ and $\gamma_{out}$ at a given $\gamma_{0}$ and grid number $M$.

Two cases of $\gamma_{0}$, i.e., $\gamma_{0}=0.2,0.8$, are simulated with two relaxation times $\tau_{\phi}$. As done in the above problem, one $\tau_{\phi}$ is given by the average $\gamma$ from each $\gamma_{0}$ to meet $\phi_{s}=0$. The distributions of $\phi$ along the centerline are predicted by the D2Q4, D2Q5 and D2Q9 lattice models. Figs. 10, 11 and 12 delineate the profiles of $\phi$ between the two cylinders under different grid sizes $M$.

Take the D2Q4 lattice model as an example. When the relaxation time $\tau_{\phi}$ obeys Eq. (19a) even through the average $\gamma$, the predictions are all much more accurate than the cases where $\tau_{\phi}$ dissatisfies $\phi_{s}=0$ (e.g., $\tau_{\phi}=8.0$). The similar results in Figs. 11 and 12 from the D2Q5 and D2Q9 lattice models again demonstrate our analysis in this work. In addition, with a careful look at Fig. 5 and Fig. 10 for $\phi_{s}=0$ under a small grid number, the agreement of the present results with the analytical solution is found not as closely as those achieved in the previous problem. This result can be expected since the approximated average distance ratio is used for the curved boundary of the cylinders.

V Conclusions

In this work, the discrete effect on the ABB boundary condition has been analyzed in the framework of BGK model for the CDE. Different from previous works on the HABB boundary condition, the boundary scheme adopted in this paper incorporates the distance ratio $\gamma$ of boundary nodes as a free parameter HuangJ15 . The theoretical derivations clearly shows that unlike the HABB boundary scheme ($\gamma=1/2$), the numerical slip $\phi_{s}$ of the NHABB boundary condition can be relieved from only relating with the relaxation time $\tau_{\phi}$ but together with $\gamma$. Therefore, as the numerical slip $\phi_{s}=0$ is guaranteed, the relaxation time $\tau_{\phi}$ can be freely adjusted as a function of the distance ratio $\gamma$, which cannot be realized for the HABB boundary condition. Concretely, for the distance ratio $\gamma$ varying in a proper range, if the relaxation time $\tau_{\phi}$ changing with $\gamma$ conforms to Eq. (19a) in the D2Q4 lattice model, Eq. (20a) in the D2Q5 lattice model, or Eq. (21a) in the D2Q9 lattice model, the discrete effect of the NHABB boundary condition can be eliminated within the framework of BGK model, while in the HABB boundary condition, the discrete effect always exists except for a special value of the relaxation time $\tau_{\phi}$. On the basis of the BGK model, the non-halfway and halfway ABB boundary conditions are both implemented to validate the theoretical analysis. For the unidirectional diffusion with a parabolic distribution in a straight channel, the numerical results show that owing to the free parameter of $\gamma$, a much wider range of the relaxation time $\tau_{\phi}$ can be achieved to produce accurate results. For the diffusion between two concentric circular cylinders, satisfactory agreements between the numerical results and the analytical solution can be obtained even with the average distance ratio.

We would like to point out that due to the quadratic dependence on $\gamma$ in $\phi_{s}=0$, the minimum relaxation time $\tau_{\phi}$ can reach only around 1 while not near 0.5, as shown in Fig. 3. However, we also note that this limitation can be improved by adding more degree of freedom in determining $\tau_{\phi}$ from the numerical slip $\phi_{s}$. One straightforward strategy for this is to extend the present analysis from the framework of BGK model to the two-relaxation-time (TRT) or the multiple-relaxation-time (MRT) model. This topic will be investigated in our forthcoming work.

Acknowledgements.

This work is financially supported by the National Natural Science Foundation of China (No. 51776068, No. 51606064 and No. 11602075) and the Fundamental Research Funds for the Central Universities (No. 2018MS060). L. Wang would like to thank Profs. Wen-An Yong, Zhaoli Guo and Dr. Weifeng Zhao for their fruitful discussions and advices.