Discrete effects on some boundary schemes of multiple-relaxation-time lattice Boltzmann model for convection-diffusion equations

TL;DR

This paper analyzes the discrete effects of boundary schemes in the multiple-relaxation-time lattice Boltzmann model for convection-diffusion equations, providing relations to eliminate numerical slip and extending results to complex-valued equations.

Contribution

It offers a general analysis of boundary scheme discrete effects, introduces a flexible parameter adjustment method, and extends the analysis to complex-valued convection-diffusion equations.

Findings

Discrete effects are independent of the transform matrix choice.

Relations are derived to eliminate numerical slip in ABB and BB schemes.

Numerical examples confirm the effectiveness of the proposed relations.

Abstract

In this paper, we perform a more general analysis on the discrete effects of some boundary schemes of the popular one- to three-dimensional DnQq multiple-relaxation-time lattice Boltzmann model for convection-diffusion equation (CDE). Investigated boundary schemes include anti-bounce-back(ABB) boundary scheme, bounce-back(BB) boundary scheme and non-equilibrium extrapolation(NEE) boundary scheme. In the analysis, we adopt a transform matrix $\textbf{M}$ constructed by natural moments in the evolution equation, and the result of ABB boundary scheme is consistent with the existing work of orthogonal matrix $\textbf{M}$. We also find that the discrete effect does not rely on the choice of transform matrix, and obtain a relation to determine some of the relaxation-time parameters which can be used to eliminate the numerical slip completely under some assumptions. In this relation, the weight coefficient is considered as an adjustable parameter which makes the parameter adjustment more flexible. The relaxation factors associated with second moments can be used to eliminate the numerical slip of ABB boundary scheme and BB boundary scheme while the numerical slip can not be eliminated of NEE boundary scheme. Furthermore, we extend the relations to complex-valued CDE, several numerical examples are used to test the relations.

Tables

Table 1. Table 1: The a 0 subscript 𝑎 0 a_{0} and a 1 subscript 𝑎 1 a_{1} in different lattice models.

	Different models		$a_0$		$a_1$
	$D1Q2$		$0$		${\omega }_1$
	$D1Q3$		${\omega }_0$		${\omega }_1$
	$D2Q4$		$2{\omega }_1$		${\omega }_1$
	$D2Q5$		${\omega }_0+2{\omega }_1$		${\omega }_1$
	$D2Q9$		${\omega }_0+2{\omega }_1$		${\omega }_1+2{\omega }_5$
	$D3Q7$		${\omega }_0+4{\omega }_1$		${\omega }_1$
	$D3Q13$		${\omega }_0+4{\omega }_1$		$4{\omega }_1$
	$D3Q15$		${\omega }_0+4{\omega }_1$		${\omega }_1+4{\omega }_7$
	$D3Q19$		${\omega }_0+4{\omega }_1+4{\omega }_7$		${\omega }_1+4{\omega }_7$
	$D3Q27$		${\omega }_0+4{\omega }_1+4{\omega }_7$		${\omega }_1+4{\omega }_7+4{\omega }_{19}$

Table 2. Table 2: The GREs of D2Q9 MRT model with ABB boundary scheme and different relaxation parameters ( w 0 = 4 / 9 subscript 𝑤 0 4 9 w_{0}={4}/{9} , w 1 = 1 / 9 subscript 𝑤 1 1 9 w_{1}={1}/{9} , w 5 = 1 / 36 subscript 𝑤 5 1 36 w_{5}={1}/{36} ).

Different values				$N=5$		$N=9$		$N=17$
$s_1=0.1$		$s_3=0.0$		$1.6143\times {10}^{-14}$		$1.1575\times {10}^{-14}$		$4.6266\times {10}^{-15}$
		$s_3=1.0$		$9.1778\times {10}^{-16}$		$4.5187\times {10}^{-16}$		$3.2051\times {10}^{-16}$
		$s_3=s_1$		$6.4495\times {10}^{-16}$		$6.5046\times {10}^{-16}$		$2.8975\times {10}^{-16}$
		$s_3=s_2$		$7.0977\times {10}^{-16}$		$5.5918\times {10}^{-16}$		$6.5757\times {10}^{-16}$
$s_1=0.6$		$s_3=0.0$		$1.4288\times {10}^{-14}$		$9.2039\times {10}^{-15}$		$2.1330\times {10}^{-8}$
		$s_3=1.0$		$4.8550\times {10}^{-16}$		$2.4793\times {10}^{-15}$		$2.1372\times {10}^{-8}$
		$s_3=s_1$		$2.6330\times {10}^{-16}$		$1.5328\times {10}^{-15}$		$2.1355\times {10}^{-8}$
		$s_3=s_2$		$4.5732\times {10}^{-16}$		$4.3549\times {10}^{-15}$		$2.1393\times {10}^{-8}$
$s_1=1.071797$		$s_3=0.0$		$2.1428\times {10}^{-14}$		$1.8222\times {10}^{-8}$		$1.1939\times {10}^{-7}$
		$s_3=1.0$		$2.5713\times {10}^{-15}$		$1.8272\times {10}^{-8}$		$1.1947\times {10}^{-7}$
		$s_3=s_1$		$2.3383\times {10}^{-15}$		$1.8275\times {10}^{-8}$		$1.1948\times {10}^{-7}$
		$s_3=s_2$		$2.2578\times {10}^{-15}$		$1.8275\times {10}^{-8}$		$1.1948\times {10}^{-7}$
$s_1=1.9$		$s_3=0.0$		$2.2912\times {10}^{-7}$		$8.1846\times {10}^{-7}$		$3.0786\times {10}^{-6}$
		$s_3=1.0$		$2.2926\times {10}^{-7}$		$8.1861\times {10}^{-7}$		$3.0787\times {10}^{-6}$
		$s_3=s_1$		$2.2938\times {10}^{-7}$		$8.1873\times {10}^{-7}$		$3.0789\times {10}^{-6}$
		$s_3=s_2$		$2.2914\times {10}^{-7}$		$8.1849\times {10}^{-7}$		$3.0786\times {10}^{-6}$

Table 3. Table 3: The E m ​ a ​ x subscript 𝐸 𝑚 𝑎 𝑥 E_{max} of D2Q9 MRT model with ABB boundary scheme and different relaxation parameters ( w 0 = 4 / 9 subscript 𝑤 0 4 9 w_{0}={4}/{9} , w 1 = 1 / 9 subscript 𝑤 1 1 9 w_{1}={1}/{9} , w 5 = 1 / 36 subscript 𝑤 5 1 36 w_{5}={1}/{36} ).

Different values				$N=5$		$N=9$		$N=17$
$s_1=0.1$		$s_3=0.0$		$2.1427\times {10}^{-14}$		$1.5488\times {10}^{-14}$		$4.6266\times {10}^{-15}$
		$s_3=1.0$		$1.3322\times {10}^{-16}$		$5.6899\times {10}^{-16}$		$4.4409\times {10}^{-16}$
		$s_3=s_1$		$7.7716\times {10}^{-16}$		$7.7716\times {10}^{-16}$		$4.4409\times {10}^{-16}$
		$s_3=s_2$		$7.7716\times {10}^{-16}$		$7.7716\times {10}^{-16}$		$4.4409\times {10}^{-16}$
$s_1=0.6$		$s_3=0.0$		$1.7097\times {10}^{-14}$		$1.0214\times {10}^{-15}$		$2.2029\times {10}^{-8}$
		$s_3=1.0$		$4.4409\times {10}^{-16}$		$2.5535\times {10}^{-15}$		$2.2073\times {10}^{-8}$
		$s_3=s_1$		$3.3307\times {10}^{-16}$		$1.6653\times {10}^{-15}$		$2.2055\times {10}^{-8}$
		$s_3=s_2$		$4.4409\times {10}^{-16}$		$4.6629\times {10}^{-15}$		$2.2095\times {10}^{-8}$
$s_1=1.071797$		$s_3=0.0$		$2.2759\times {10}^{-14}$		$1.8820\times {10}^{-8}$		$1.2330\times {10}^{-7}$
		$s_3=1.0$		$2.4425\times {10}^{-15}$		$1.8871\times {10}^{-8}$		$1.2339\times {10}^{-7}$
		$s_3=s_1$		$2.2204\times {10}^{-15}$		$1.8875\times {10}^{-8}$		$1.2339\times {10}^{-7}$
		$s_3=s_2$		$2.2204\times {10}^{-15}$		$1.8875\times {10}^{-8}$		$1.2339\times {10}^{-7}$
$s_1=1.9$		$s_3=0.0$		$2.3665\times {10}^{-7}$		$8.4530\times {10}^{-7}$		$3.1796\times {10}^{-6}$
		$s_3=1.0$		$2.3679\times {10}^{-7}$		$8.4546\times {10}^{-7}$		$3.1797\times {10}^{-6}$
		$s_3=s_1$		$2.3691\times {10}^{-7}$		$8.4558\times {10}^{-7}$		$3.1798\times {10}^{-6}$
		$s_3=s_2$		$2.3667\times {10}^{-7}$		$8.4534\times {10}^{-7}$		$3.1796\times {10}^{-6}$

Table 4. Table 4: The GREs of D2Q5 and D2Q9 MRT models with ABB boundary scheme and different parameters.

Different models				$N=5$		$N=9$		$N=17$
$s_1=0.1$		$D2Q9,{\omega }_0=\frac{4}{9},{\omega }_1=\frac{1}{9},{\omega }_5=\frac{1}{36}$		$9.1778\times {10}^{-16}$		$4.5187\times {10}^{-16}$		$3.2051\times {10}^{-16}$
		$D2Q5,{\omega }_0=\frac{1}{5},{\omega }_1=\frac{1}{5}$		$5.7786\times {10}^{-16}$		$5.2053\times {10}^{-16}$		$3.3281\times {10}^{-16}$
$s_1=0.6$		$D2Q9,{\omega }_0=\frac{4}{9},{\omega }_1=\frac{1}{9},{\omega }_5=\frac{1}{36}$		$4.8550\times {10}^{-16}$		$2.4793\times {10}^{-15}$		$2.1372\times {10}^{-8}$
		$D2Q5,{\omega }_0=\frac{1}{5},{\omega }_1=\frac{1}{5}$		$2.8491\times {10}^{-16}$		$1.5632\times {10}^{-16}$		$1.5599\times {10}^{-8}$
$s_1=1.9$		$D2Q9,{\omega }_0=\frac{4}{9},{\omega }_1=\frac{1}{9},{\omega }_5=\frac{1}{36}$		$2.2926\times {10}^{-7}$		$8.1861\times {10}^{-7}$		$3.0787\times {10}^{-6}$
		$D2Q5,{\omega }_0=\frac{1}{5},{\omega }_1=\frac{1}{5}$		$1.4640\times {10}^{-7}$		$6.7983\times {10}^{-7}$		$2.6060\times {10}^{-6}$

Table 5. Table 5: The E m ​ a ​ x subscript 𝐸 𝑚 𝑎 𝑥 E_{max} of D2Q5 and D2Q9 MRT models with ABB boundary scheme and different parameters.

Different models				$N=5$		$N=9$		$N=17$
$s_1=0.1$		$D2Q9,{\omega }_0=\frac{4}{9},{\omega }_1=\frac{1}{9},{\omega }_5=\frac{1}{36}$		$1.3322\times {10}^{-15}$		$5.6899\times {10}^{-16}$		$4.4409\times {10}^{-16}$
		$D2Q5,{\omega }_0=\frac{1}{5},{\omega }_1=\frac{1}{5}$		$6.1062\times {10}^{-16}$		$5.5511\times {10}^{-16}$		$4.4409\times {10}^{-16}$
$s_1=0.6$		$D2Q9,{\omega }_0=\frac{4}{9},{\omega }_1=\frac{1}{9},{\omega }_5=\frac{1}{36}$		$4.4409\times {10}^{-16}$		$2.5535\times {10}^{-15}$		$2.2073\times {10}^{-8}$
		$D2Q5,{\omega }_0=\frac{1}{5},{\omega }_1=\frac{1}{5}$		$3.3307\times {10}^{-16}$		$2.2204\times {10}^{-16}$		$1.6110\times {10}^{-8}$
$s_1=1.9$		$D2Q9,{\omega }_0=\frac{4}{9},{\omega }_1=\frac{1}{9},{\omega }_5=\frac{1}{36}$		$2.3679\times {10}^{-7}$		$8.4546\times {10}^{-7}$		$3.1797\times {10}^{-6}$
		$D2Q5,{\omega }_0=\frac{1}{5},{\omega }_1=\frac{1}{5}$		$1.5121\times {10}^{-7}$		$7.0213\times {10}^{-7}$		$2.5966\times {10}^{-6}$

Table 6. Table 6: The GREs of D3Q7 and D3Q19 MRT models with ABB boundary scheme and different parameters.

Different models				$N=5$		$N=9$		$N=17$
$s_1=0.1$		$D3Q19,{\omega }_0=\frac{1}{3},{\omega }_1=\frac{1}{18},{\omega }_7=\frac{1}{36}$		$3.0474\times {10}^{-10}$		$1.7407\times {10}^{-10}$		$1.9787\times {10}^{-10}$
		$D3Q7,{\omega }_0=\frac{1}{3},{\omega }_1=\frac{1}{9}$		$1.5758\times {10}^{-10}$		$7.4854\times {10}^{-11}$		$5.6273\times {10}^{-11}$
$s_1=0.6$		$D3Q19,{\omega }_0=\frac{1}{3},{\omega }_1=\frac{1}{18},{\omega }_7=\frac{1}{36}$		$3.2280\times {10}^{-9}$		$1.7372\times {10}^{-9}$		$2.1372\times {10}^{-8}$
		$D3Q7,{\omega }_0=\frac{1}{3},{\omega }_1=\frac{1}{9}$		$4.4101\times {10}^{-11}$		$2.4137\times {10}^{-9}$		$8.7858\times {10}^{-8}$
$s_1=1.9$		$D3Q19,{\omega }_0=\frac{1}{3},{\omega }_1=\frac{1}{18},{\omega }_7=\frac{1}{36}$		$2.2926\times {10}^{-7}$		$8.1861\times {10}^{-7}$		$3.0787\times {10}^{-6}$
		$D3Q7,{\omega }_0=\frac{1}{3},{\omega }_1=\frac{1}{9}$		$3.1045\times {10}^{-7}$		$1.2000\times {10}^{-6}$		$4.6762\times {10}^{-6}$

Table 7. Table 7: The E m ​ a ​ x subscript 𝐸 𝑚 𝑎 𝑥 E_{max} of D3Q7 and D3Q19 MRT models with ABB boundary scheme and different parameters.

Different models				$N=5$		$N=9$		$N=17$
$s_1=0.1$		$D3Q19,{\omega }_0=\frac{1}{3},{\omega }_1=\frac{1}{18},{\omega }_7=\frac{1}{36}$		$3.7204\times {10}^{-10}$		$1.7549\times {10}^{-10}$		$3.1573\times {10}^{-10}$
		$D3Q7,{\omega }_0=\frac{1}{3},{\omega }_1=\frac{1}{9}$		$1.2115\times {10}^{-10}$		$1.1272\times {10}^{-11}$		$6.4008\times {10}^{-11}$
$s_1=0.6$		$D3Q19,{\omega }_0=\frac{1}{3},{\omega }_1=\frac{1}{18},{\omega }_7=\frac{1}{36}$		$3.3340\times {10}^{-9}$		$1.7942\times {10}^{-9}$		$2.2073\times {10}^{-8}$
		$D3Q7,{\omega }_0=\frac{1}{3},{\omega }_1=\frac{1}{9}$		$4.5550\times {10}^{-11}$		$2.4929\times {10}^{-9}$		$9.0739\times {10}^{-8}$
$s_1=1.9$		$D3Q19,{\omega }_0=\frac{1}{3},{\omega }_1=\frac{1}{18},{\omega }_7=\frac{1}{36}$		$2.3679\times {10}^{-7}$		$8.4546\times {10}^{-7}$		$3.1797\times {10}^{-6}$
		$D3Q7,{\omega }_0=\frac{1}{3},{\omega }_1=\frac{1}{9}$		$3.2065\times {10}^{-7}$		$1.2394\times {10}^{-6}$		$4.8296\times {10}^{-6}$

Table 8. Table 8: The GREs of D2Q5 MRT model for the complex cases with ABB boundary scheme ( ω 0 = 1 / 3 , ω 1 = 1 / 6 formulae-sequence subscript 𝜔 0 1 3 subscript 𝜔 1 1 6 \omega_{0}={1}/{3},\omega_{1}={1}/{6} ).

Different models				$N=5$		$N=9$		$N=17$
$MRT$		${\tau }_r=1.0,{\tau }_i=0.5$		$1.2775\times {10}^{-16}$		$2.0708\times {10}^{-9}$		$1.1467\times {10}^{-7}$
		${\tau }_r=10.0,{\tau }_i=9.5$		$4.1977\times {10}^{-16}$		$1.2100\times {10}^{-14}$		$1.9386\times {10}^{-10}$
		${\tau }_r=0.501,{\tau }_i=0.001$		$4.4416\times {10}^{-6}$		$1.5648\times {10}^{-5}$		$5.8706\times {10}^{-5}$

Table 9. Table 9: The E m ​ a ​ x subscript 𝐸 𝑚 𝑎 𝑥 E_{max} of D2Q5 MRT model for the complex cases with ABB boundary scheme ( ω 0 = 1 / 3 , ω 1 = 1 / 6 formulae-sequence subscript 𝜔 0 1 3 subscript 𝜔 1 1 6 \omega_{0}={1}/{3},\omega_{1}={1}/{6} ).

Different models				$N=5$		$N=9$		$N=17$
$MRT$		${\tau }_r=1.0,{\tau }_i=0.5$		$1.1102\times {10}^{-16}$		$2.5946\times {10}^{-9}$		$1.3973\times {10}^{-7}$
		${\tau }_r=10.0,{\tau }_i=9.5$		$7.2165\times {10}^{-16}$		$1.7431\times {10}^{-14}$		$2.3966\times {10}^{-10}$
		${\tau }_r=0.501,{\tau }_i=0.001$		$6.4009\times {10}^{-6}$		$2.2385\times {10}^{-5}$		$8.4811\times {10}^{-5}$

Table 10. Table 10: The GREs of D2Q5 MRT model for the complex cases with BB boundary scheme ( ω 0 = 1 / 3 , ω 1 = 1 / 6 formulae-sequence subscript 𝜔 0 1 3 subscript 𝜔 1 1 6 \omega_{0}={1}/{3},\omega_{1}={1}/{6} ).

Different models				$N=5$		$N=9$		$N=17$
$MRT$		${\tau }_r=1.0,{\tau }_i=0.5$		$3.2814\times {10}^{-16}$		$1.4988\times {10}^{-11}$		$1.0332\times {10}^{-7}$
		${\tau }_r=10.0,{\tau }_i=9.5$		$1.6129\times {10}^{-11}$		$5.6665\times {10}^{-11}$		$1.0665\times {10}^{-9}$
		${\tau }_r=0.501,{\tau }_i=0.001$		$2.8179\times {10}^{-6}$		$1.2367\times {10}^{-5}$		$5.2043\times {10}^{-5}$

Table 11. Table 11: The E m ​ a ​ x subscript 𝐸 𝑚 𝑎 𝑥 E_{max} of D2Q5 MRT model for the complex cases with BB boundary scheme ( ω 0 = 1 / 3 , ω 1 = 1 / 6 formulae-sequence subscript 𝜔 0 1 3 subscript 𝜔 1 1 6 \omega_{0}={1}/{3},\omega_{1}={1}/{6} ).

Different models				$N=5$		$N=9$		$N=17$
$MRT$		${\tau }_r=1.0,{\tau }_i=0.5$		$4.4409\times {10}^{-16}$		$2.3118\times {10}^{-11}$		$1.4393\times {10}^{-7}$
		${\tau }_r=10.0,{\tau }_i=9.5$		$2.2474\times {10}^{-11}$		$7.3081\times {10}^{-11}$		$1.6065\times {10}^{-9}$
		${\tau }_r=0.501,{\tau }_i=0.001$		$4.4960\times {10}^{-6}$		$1.8838\times {10}^{-5}$		$7.7969\times {10}^{-5}$

Equations

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

$$Collision:$$

$$Streaming:$$

$$f_i^{eq}(\text{x},t)=w_i\varphi (1+\frac{{\text{c}}_i\cdot \text{u}}{c_s^2}),$$

$$R_i={\omega }_{i}R.$$

$$\textbf{M}_{0}=\left(\begin{array}[]{ccc}1&1&1\\ -1&0&1\\ 1&0&1\\ \end{array}\right).$$

$${\text{C}}_d=diag(1,c,c^2),$$

$$\text{S}=diag(s_0,s_1,s_2).$$

$$\textbf{c}=\left(\begin{array}[]{ccccccccc}0&1&0&-1&0&1&-1&-1&1\\ 0&0&1&0&-1&1&1&-1&-1\\ \end{array}\right)c,$$

$$\textbf{M}_{0}=\left(\begin{array}[]{ccccccccc}1&1&1&1&1&1&1&1&1\\ 0&1&0&-1&0&1&-1&-1&1\\ 0&0&1&0&-1&1&1&-1&-1\\ 0&1&0&1&0&1&1&1&1\\ 0&0&1&0&1&1&1&1&1\\ 0&0&0&0&0&1&-1&1&-1\\ 0&0&0&0&0&1&1&-1&-1\\ 0&0&0&0&0&1&-1&-1&1\\ 0&0&0&0&0&1&1&1&1\\ \end{array}\right).$$

$${\text{C}}_d=diag(1,c,c,c^2,c^2,c^2,c^3,c^3,c^4),$$

$$\text{S}=diag(s_0,s_1,s_1,s_2,s_2,s_2,s_3,s_3,s_4).$$

$$\varphi =\sum _i{f}_i+\frac{\theta R}{2}\delta t.$$

$$f_i^j=\{\begin{array}{rlr}{f}_i^{j,+}& ,& i=0\\ f_i^{j-1,+}& ,& i=1\\ f_i^{j+1,+}& ,& i=-1\end{array}$$

$$f_i^{j,+}=f_i(x_j,t)-({\text{M}}^{-1}\text{S}\text{M}{)}_{ik}(f_k(x_j,t)-f_k^{eq}(x_j,t))+{\delta }_t[{\text{M}}^{-1}(\text{I}-\frac{\theta \text{S}}{2})\text{M}{]}_{ik}{R}_k,i=0,1,-1.$$

$$\frac{{\varphi }_{k+1}{u}_{k+1}-{\varphi }_{k-1}{u}_{k-1}}{2\delta x}=D\frac{{\varphi }_{k+1}-2{\varphi }_k+{\varphi }_{k-1}}{\delta x^2}+R,$$

$$f_i^j=\{\begin{array}{rlr}{f}_i^{j,+}& ,& i=0,1,3\\ f_i^{j-1,+}& ,& i=2,5,6\\ f_i^{j+1,+}& ,& i=4,7,8\end{array}$$

$$f_{013}^k=$$

$$+[({\omega }_1+2{\omega }_5)\theta (s_2-s_0)+({\omega }_0+2{\omega }_1)(1-\frac{\theta s_0}{2})]{\delta }_{t}R,$$

$$f_{256}^{k+1}=f_{256}^k-(\frac{s_2}{2}-\frac{s_1}{2})(f_{478}^k-f_{478}^{k,eq})-(\frac{s_2}{2}+\frac{s_1}{2})(f_{256}^k-f_{256}^{k,eq})+({\omega }_1+2{\omega }_5)(1-\frac{\theta s_2}{2}){\delta }_{t}R,$$

$$f_{478}^{k-1}=f_{478}^k-(\frac{s_1}{2}+\frac{s_2}{2})(f_{478}^k-f_{478}^{k,eq})-(\frac{s_2}{2}-\frac{s_1}{2})(f_{256}^k-f_{256}^{k,eq})+({\omega }_1+2{\omega }_5)(1-\frac{\theta s_2}{2}){\delta }_{t}R,$$

$$f_{013}^k={\varphi }_k-f_{256}^k-f_{478}^k-\frac{\theta R}{2}\delta t.$$

$$f_{256}^k+f_{478}^k={\varphi }_k-f_{013}^{k,eq}+A{\delta }_{t}R,$$

$$f_{256}^{k+1}=(1-s_1)f_{256}^k+s_1{f}_{256}^{k,eq}+B{\delta }_{t}R,$$

$$f_{478}^{k-1}=(1-s_1)f_{478}^k+s_1{f}_{478}^{k,eq}+B{\delta }_{t}R,$$

$$f_{256}^k=(1-s_1)f_{256}^{k-1}+s_1{f}_{256}^{k-1,eq}+B{\delta }_{t}R,$$

$$f_{478}^k=(1-s_1)f_{478}^{k+1}+s_1{f}_{478}^{k+1,eq}+B{\delta }_{t}R.$$

$$f_{256}^k=(1-s_1)({\varphi }_{k-1}-f_{478}^{k-1}-f_{013}^{k-1,eq}+A{\delta }_{t}R)+s_1{f}_{256}^{k-1,eq}+B{\delta }_{t}R,$$

$$f_{478}^k=(1-s_1)({\varphi }_{k+1}-f_{256}^{k+1}-f_{013}^{k+1,eq}+A{\delta }_{t}R)+s_1{f}_{478}^{k+1,eq}+B{\delta }_{t}R.$$

$$a_1\frac{s_1-2}{s_1}({\varphi }_{k+1}+{\varphi }_{k-1}-2{\varphi }_k)=\frac{{\varphi }_{k+1}{u}_{y,k+1}-{\varphi }_{k-1}{u}_{y,k-1}}{2c}+{\delta }_{t}R,$$

$$\frac{{\varphi }_{k+1}{u}_{y,k+1}-{\varphi }_{k-1}{u}_{y,k-1}}{2\delta x}=D\frac{{\varphi }_{k+1}-2{\varphi }_k+{\varphi }_{k-1}}{\delta x^2}+R,$$

$$f_{478}^1=(1-s_1)f_{478}^2+s_1{f}_{478}^{2,eq}+B{\delta }_{t}R.$$

$$f_{478}^1=(1-s_1)({\varphi }_2-f_{013}^{2,eq}+A{\delta }_{t}R-f_{256}^2)+s_1{f}_{478}^{2,eq}+B{\delta }_{t}R.$$

$$f_{478}^1=(1-s_1)(2a_1{\varphi }_2+A{\delta }_{t}R-(1-s_1)(2a_1{\varphi }_1-f_{478}^1+A{\delta }_{t}R)+2a_1{s}_1{\varphi }_1+B{\delta }_{t}R)+2a_1{s}_1{\varphi }_2+B{\delta }_{t}R.$$

Full text

††journal: Computers & Mathematics with Applications

1 introduction

In recent years, the lattice Boltzmann method (LBM) has gained much attention, and has been wildly used in many fields [1, 2, 3, 4]. The LBM has some distinct advantages over traditional methods in dealing with Navier-Stokes equations [5, 6, 7, 8, 9] and convection-diffusion equations (CDEs) [10, 11, 12, 13, 14, 15]. One of the advantages of LBM is dealing with the complex boundary conditions in porous media [16, 17, 18, 19, 20, 21]. When we solve the macroscopic partial differential equation, there are always discrete errors in the numerical scheme, the boundary discrete effect exists between the real boundary condition and the numerical solution in the boundary point.

To our knowledge, the discrete effect of the bounce-back(BB) scheme was first discussed for the Poiseuille flow. Ginzburg and Adler [22] first performed a boundary condition analysis for the face-centered-hypercubic lattice Boltzmann (LB) model applied to the Poiseuille flow and a plane stagnation flow. After that, He $et$ $al.$ [23] analyzed the discrete effect of BB boundary scheme in the Bhatnagar-Gross-Krook (BGK) model, and found that the relaxation time $\tau$ has a significant influence on the BB scheme for the no-slip boundary condition. In a similar way, Guo $et$ $al.$ [24] studied the existing discrete effect of the discrete Maxwell’s diffuse-reflection (DMDR) scheme and the combined bounce-back/specular-reflection (CBBSR) scheme. Then, they simulated the Poiseuille flow in the slip flow regime with the multiple-relaxation-time (MRT) LB model, and found that the BGK model cannot yield correct results in this regime owing to the discrete effect [25]. Due to find that the boundary schemes considered in Refs. [24, 25] are nonlocal, they are not suitable for fluid flows in complex geometries, Chai $et$ $al.$ [26] developed a local scheme combined halfway bounce-back boundary condition and full diffusive boundary condition for microscale gas flows in complex geometries, and illustrated that to realize the exact slip boundary condition, the discrete effect must be included and corrected. Lu $et$ $al.$ [27] proposed an immerse boundary MRT LB model, and presented a special relaxation between two relaxation time parameters in which can reduce the numerical boundary slip effectively. Recently, Ren $et$ $al.$ [28] analyzed the discrete effects in the DMDR and CBBSR schemes for the rectangular LBE, and presented a reasonable approach to overcome these discrete effects in these two schemes.

We noted that all of above works focus on the discrete effect of BB boundary scheme for fluid flows. Subsequently, there are also some works on the discrete effect of anti-bounce-back (ABB) boundary schemes for CDEs. Zhang $et$ $al.$ [29] presented a general ABB boundary scheme of the BGK model for CDEs. They performed an analysis on the discrete effect of the ABB boundary condition, and suggested that there is a numerical slip related to the lattice size in the diffusion of Couette flow between solid walls, which cannot be eliminated in the BGK model. Then, Cui $et$ $al.$ [30] analyzed the ABB boundary condition of the MRT model for CDEs. They presented a theoretical analysis on the discrete effect of the ABB boundary scheme for the simple problems with a parabolic distribution in one direction, and observed that the numerical slip can be eliminated in the MRT LB model by choosing the free relaxation parameters properly. However, the analysis is limited to some special MRT LB models, e.g., D2Q4, D2Q5, and D2Q9 model. Recently, based on the two-relaxation-times(TRT) model, Ginzburg $et$ $al.$ [31] presented a more general relation between the two relaxation factors through equating the set of closure relations of the given boundary scheme to the Taylor expansion. In this work, based on the existing works [30], we firstly conduct the discrete effect on the ABB boundary scheme of the more general MRT model composed of the natural moments for CDEs, and then derived a relation with four parameters : the weight coefficient, the relaxation factors $s_{1}$ and $s_{2}$ associated with first and second moments and a model parameter $\theta$ for adjustment to elimate the numerical slip. After that, we conduct the discrete effect on BB boundary condition and non-equilibrium extrapolation(NEE) boundary condition, observed that the discrete effect can be elimated when $s_{1}+s_{2}=2$ on BB boundary condition and can not be elimated on NEE boundary condition. Furthermore, we observed that the relations is applicable to both real- and complex-valued problems, and has a general expression from one to three dimensions.

The paper is organized as follows. In Sec. II, we introduce the MRT model composed of natural moments. Then we derive the equivalent finite-difference scheme of the MRT model for CDEs, and discuss the discrete effects on the ABB, BB, NEE boundary conditions in Sec. III. Numerical tests are performed in Sec. IV. Finally, we give a brief summary in Sec. V.

2 MRT LB model for convection-diffusion equation

Firstly, we introduce the MRT model composed of the natural moments for CDEs. The n-dimensional (nD) CDEs can be written as

[TABLE]

where $\phi$ is a scalar function of position x and time t, $\mathbf{\nabla}$ is the gradient operator with respect to the position x in n dimensions. $D$ is the diffusion coefficient, u is the convection velocity and $R$ is the source term.

The evolution equation of the MRT model with D$n$Q$q$ lattice for the CDE can be written as

[TABLE]

where $\delta_{t}$ is time step, I is the identity matrix, and S is a diagonal relaxation matrix with non-negative elements. The transformation matrix M is composed of natural moments [32]. $\theta$ is a real parameter, corresponding to the MRT model [30] for $\theta=1$ and a scheme in Ref. [33] for $\theta=0$, respectively. $f_{i}(\textbf{x},t)$ and $f_{i}^{eq}(\textbf{x},t)$ are the distribution function and equilibrium distribution function (EDF) associated with the discrete velocity $\textbf{c}_{i}$ at position x and time $t$ respectively, and $f_{i}(\textbf{x},t)^{+}$ is the distribution function after collision. And to simplify the derivation, only the following linear EDF is considered here,

[TABLE]

where $\omega_{i}$ is the weight coefficient, $c_{s}$ is the so-called lattice sound speed. $R_{i}$ is the discrete source term, and can be defined as

[TABLE]

Firstly, for the D$1$Q$3$ model, the set of discrete velocities are $\textbf{c}=\{-1,0,1\}c$, where $c=\delta_{x}/\delta_{t}$ with $\delta_{x}$ being the lattice spacing. The transformation matrix $\textbf{M}=(\textbf{c}_{ix}^{m})(m=0,1,2)$, which can be expressed as $\textbf{M}=\textbf{C}_{d}\textbf{M}_{0}$ [34],

[TABLE]

[TABLE]

[TABLE]

As for the D2Q9 model, the discrete velocities can be given by

[TABLE]

and the transformation matrix as $\textbf{M}=(\textbf{c}^{m}_{ix}\textbf{c}^{n}_{iy})=\textbf{C}_{d}\textbf{M}_{0}$, $(m,n=0,1,2,m+n\leq 2)$,

[TABLE]

[TABLE]

[TABLE]

In the present MRT model, the macroscopic variable $\phi$ should be computed by

[TABLE]

3 Discrete effects of some boundary schemes

We now analyze the discrete effects of these boundary scheme in the framework of the MRT model for CDE. For simplicity, we conducted an analysis of Dirichlet boundary conditions for the simple steady problems with a parabolic distribution in one direction.

3.1 Equivalent difference equation of the MRT model

Firstly, we consider the D$1$Q$3$ MRT model for one-dimensional steady problems with const $R$, and set the distribution function as $f_{i}^{j}=f_{i}(x_{j})$, with $x_{j}$ being a discrete grid point. To make the derivation easier to understand, we rewrite Eq. (2) as

[TABLE]

where

[TABLE]

After taking some manipulations of the evolution equation, as shown in Fig. 1(see Appendix A for details), we can obtain the following equivalent difference equation of the MRT model,

[TABLE]

where $D=({1}/{s_{1}}-{1}/{2})c^{2}_{s}\delta t$, $c^{2}_{s}=2\omega_{1}c^{2}$. Then we consider the D$2$Q$9$ MRT model for x-direction steady problems with constant $R$, and set the distribution function as $f_{i}^{j}=f_{i}(x_{k},y_{j})$, with $y_{j}$ being a discrete grid point,and $i$ being the direction of distribution function. Eq. (2) can be rewritten as

[TABLE]

where ${f}_{i}^{j,+}=f_{i}(x_{k},y_{j},t)^{+}$ is the distribution function after collision. Then we can take a combination of distribution function as

[TABLE]

where $f^{k}_{ijm}=f^{k}_{i}+f^{k}_{j}+f^{k}_{m}$, $f^{k,eq}_{ijm}=f^{k,eq}_{i}+f^{k,eq}_{j}+f^{k,eq}_{m}$. According to Eq. (12), we can obtain

[TABLE]

Substituting Eq. (18) into Eq. (17a), one can obtain

[TABLE]

where $a_{0}=\omega_{0}+2\omega_{1}$, $a_{1}=\omega_{1}+2\omega_{5}$, $A=-(a_{0}+a_{1}s_{2}\theta)/s_{2}$, with $\omega_{1}=\omega_{2}=\omega_{3}=\omega_{4}$, $\omega_{5}=\omega_{6}=\omega_{7}=\omega_{8}$. With the help of Eq. (19), we can rewritten Eqs. (17b) and (17c) as

[TABLE]

where $B=a_{1}(1-{\theta s_{2}}/{2})-({s_{2}}-{s_{1}})A/2$. Then we can get the following equation according to Eqs. (20a) and (20b)

[TABLE]

With the help of Eq. (19), Eqs. (21a) and (21b) can be written as

[TABLE]

Taking a sum of Eqs. (20a), (20b), (22a) and (22b), one can obtain

[TABLE]

where Eq. (19) has been adopted. Then we can obtain the following equivalent difference equation of the MRT model,

[TABLE]

where $D=({1}/{s_{1}}-{1}/{2})c^{2}_{s}\delta t$, $c^{2}_{s}=2a_{1}c^{2}$, $a_{1}=\omega_{1}+2\omega_{5}$. Here we would like to point out that if we adopt different transform matrix M which is constructed by orthogonal vectors, one can obtain the same equivalent difference equation [35].

Actually, for higher dimensions lattice velocity models (e.g., D3Q27), one can obtain the similar difference scheme as Eq. (24) (see Appendix A for details). Then we will get a useful equation, in the following derivation. When $k=1$, Eq. (21b) can be written as

[TABLE]

Substituting Eq. (19) into Eq. (25), one can obtain

[TABLE]

In addition, substituting Eq. (22a) into Eq. (26) with the help of Eq. (3) gives rise to

[TABLE]

We can rewrite the Eq. (27) as

[TABLE]

3.2 Discrete effect of the ABB boundary scheme

To simplify the analysis on the discrete effect of the ABB boundary scheme, a unidirectional and time-independent diffusion problem is adopted, and it can be described by the following simplified equation and boundary conditions for one dimensional problem

[TABLE]

[TABLE]

where $\phi_{0}$ and $\phi_{L}$ are constant, $L$ is the width and $D$ is the diffusion coefficient. $R$ is a constant source term, and is defined by

[TABLE]

The analytical solution of the problem is given by

[TABLE]

Based on Eq. (15), equivalent difference equation for the MRT model for Eq. (29),

[TABLE]

Then we can obtain the solution of Eq. (33),

[TABLE]

where $a$, $b$ are parameters to be determined. If we consider ABB scheme, the value of $\phi$ at bottom and top boundaries can be given by

[TABLE]

where $\phi_{s}^{0.5}$, $\phi_{s}^{N+0.5}$ are numerical slip caused by ABB scheme, $N$ representing grid number. Substituting Eq. (35) into Eq. (34), we obtain the numerical solution

[TABLE]

In the following, we will focus on how to determine $\phi_{s}^{0.5}$ and $\phi_{s}^{N+0.5}$ from the ABB scheme. As Fig. 2 shown, the unknown distribution functions at the layers $k=1$, $k=N$ can be determined by the following equations [29],

[TABLE]

[TABLE]

where $f^{1,+}_{-1}$, $-f^{N,+}_{1}$ represent the distribution function after collision at the layers $k=1$ and $k=N$ respectively. Following the process in Appendix B, we can get the numerical slip,

[TABLE]

[TABLE]

As we can see, $\phi_{s}^{0.5}$ and $\phi_{s}^{N+0.5}$ have the same expression, thus we denote them by $\phi_{s}$ in the following discussion. If the free parameter $s_{2}$ is chosen to satisfy the relation,

[TABLE]

the discrete effect of the ABB scheme can be eliminated.

Furthermore, when we use the BGK model ($s_{1}=s_{2}$) to deal with the problem, and take the weight coefficients $\omega_{0}$ and $\omega_{1}$ to satisfy Eq. (41), the discrete effect on the ABB boundary scheme can also be eliminated. However, this selection of the weight coefficients in the BGK model is limited due to the fact that the weight coefficients should be greater than 0 and less than 1.

Similarly, for the two-dimensional unidirectional steady problem with a parabolic distribution in one direction, we analyze the discrete effect in D2Q9 MRT model. For the ABB scheme,

[TABLE]

where the $f^{1,+}_{i}=f_{i}(x_{k},y_{1},t)^{+}$ represent the distribution function after the collision. Taking a sum of Eqs. (42a), (42b), and (42c), we obtain

[TABLE]

which can be written as

[TABLE]

with the help of Eqs. (19) and (20b). Substituting

[TABLE]

and Eq. (28) into Eq. (44), we can obtain

[TABLE]

where $a_{0}=\omega_{0}+2\omega_{1}$, $a_{1}=\omega_{1}+2\omega_{5}$ in D2Q9 model.

Similarly, for the three-dimensional unidirectional steady problem with a parabolic distribution in one direction, one can obtain the following results with a similar derivation process,

[TABLE]

where $c_{s}^{2}=2a_{1}c^{2}$, $a_{0}=\omega_{0}+4\omega_{1}+4\omega_{7}$, $a_{1}=\omega_{1}+4\omega_{7}+4\omega_{19}$ in D3Q19 model. Taking the following equation

[TABLE]

in Eq. (47), one can eliminate the discrete effect. The parameters $a_{0}$ and $a_{1}$ in the different lattice model are listed in Table 1, the velocities of D2Q9 and D3Q27 models are presented in Fig. 3, and the relaxation factors $s_{1}$ and $s_{2}$ are associated with first and second moments. We note that when $\theta=1$, $\omega_{i}=1/4(i=1-4)$ in D2Q4 model, $\omega_{i}=1/5(i=0-4)$ in D2Q5 model, $\omega_{0}=4/9,\omega_{1-4}=1/9,\omega_{5-8}=1/36$ in D2Q9 model, Eq. (47) contains the previous works [30]. And Eq. (48) is consist with the recent results [31] when $\theta=1$ in the frame of TRT model. It should be noted that for a specified lattice model, we can determine the explicit expression of $\phi_{s}$ from Eq. (47), but the numerical slip $\phi_{s}$ could not be eliminated since $w_{i}$ is not flexible enough to satisfy Eq. (48). For example, in the D1Q2 model, ($\omega_{0}=0$, $\omega_{1}=1/2$), Eq. (48) can not be satisfied under the condition of $0<s_{1}<2$ and $0<s_{2}<2$.

3.3 Discrete effect of the BB boundary scheme

In this section, we analyze the BB boundary scheme under the same assumptions. For D2Q9 model with the BB boundary scheme [36],

[TABLE]

Summing Eqs. (49a), (49b), and (49c), one can obtain

[TABLE]

Then Eq. (17b) can be written as

[TABLE]

where $a_{1}=\omega_{1}+2\omega_{5}$. One can obtain

[TABLE]

with the help of Eqs. (19) and (20b). Substituting

[TABLE]

and Eq. (28) into Eq. (52), we can obtain

[TABLE]

3.4 Discrete effect of the NEE boundary scheme

For the NEE scheme [37],

[TABLE]

Based on Eq. (55), we obtain

[TABLE]

According to Eq. (21a), we have

[TABLE]

which can be written as

[TABLE]

with the help of Eq. (56). Substituting Eq. (19) into Eq. (58), one can obtain

[TABLE]

Substituting

[TABLE]

and Eq. (28) into Eq. (59), one can obtain

[TABLE]

4 NUMERICAL RESULTS

In this section, some simulations of CDEs are performed to test above analysis, and ABB scheme is employed to treat the Dirichlet boundary conditions. In our simulations, the global relative error (GRE) and maximum error($E_{max}$) are used to measure accuracy, and are defined as

[TABLE]

where $\phi$ and $\phi^{*}$ are the numerical and analytical solutions, respectively. In addition, the following convergent criterion for the steady problems is used,

[TABLE]

In our simulations, $f_{i}^{eq}$ is applied to approximate the initial distribution function $f_{i}$.

4.1 Some unidirectional time-independent real-valued CDEs

4.1.1 A linear time-independent diffusion equation

We first consider a two-dimensional linear time-independent diffusion equation with a constant source term,

[TABLE]

where $\phi_{0}$ and $\phi_{L}$ are two constants, $L$ is the width between the top and bottom boundaries, and $R$ is the source term and is defined by

[TABLE]

The analytical solution of this problem is given by

[TABLE]

Here we consider the popular D2Q9 MRT model with $\theta=1$, the physical parameter $L=1.0$, $u_{x}=0.1$, $u_{y}=0.0$, the diffusion coefficient $D=0.1$, the boundary conditions $\phi_{0}=0$, $\phi_{L}=1$, $\delta_{x}=L/N$ with the grid number $N$ varying from 5 to 17.

First, we would like to verify that the parameters except $s_{1}$ and $s_{2}$ have little effect on numerical results. In our simulations, the value of $s_{1}$ is determined by the diffusion coefficient, while $s_{2}$ is given by Eq. (48). We measured the GREs of the problem under different values of $s_{3}$, and present the results in Table 2 and Table 3. As shown in these table, for the fixed $s_{1}$ and $N$, the relaxation parameter $s_{3}$ has little influence on GREs. For this reason, except $s_{1}$ and $s_{2}$, the other parameters in S are set to be 1.0 in the following simulations. In general, the GRE decreases with the increase of grid number $N$, and as we shown in Table 2 the GRE increases for the accumulation of mechanical errors when the grid number $N$ increases.

After that, we test different weight coefficients in the D2Q9 BGK model when $s_{1}=0.1$ and $0.5$. In Fig. (4), the case 1 is $\omega_{i}=1/9,(i=0-8)$, the case 2 is $\omega_{0}=4/9$, $\omega_{1}=1/9$, $\omega_{5}=1/36$, the case 3 is a set of weight coefficients satisfied Eq. (48). In our simulation, case 3 is $\omega_{0}=1/1083$, $\omega_{1}=1/4332$, $\omega_{5}=1081/4332$ when $s_{1}=0.1$, $\omega_{0}=1/27$, $\omega_{1}=1/108$, $\omega_{5}=25/108$ when $s_{1}=0.5$. We can see that case 3 has more accurate results than case 1 and case 2. As we known, the weight coefficients in the D2Q9 model are given as $\omega_{0}=4/9$, $\omega_{1}=1/9$, $\omega_{5}=1/36$ for Navier-Stokes equations. Actually, weight coefficients in the LB model for CDEs are more flexible and they could be adjusted to give more accurate results. This adjustment has certain limitations because the weight coefficients must be greater than zero. For BGK model, taking $s_{1}=s_{2}$, $0<a_{0}<1$, and $(1/s_{1}-1/2)^{2}=1/(8a_{0})$, then we can get the limitation $1/s_{1}>(1+\sqrt{2})/2\sqrt{2}$. When $1/s_{1}>(1+\sqrt{2})/2\sqrt{2}$, $\phi_{s}$ on ABB boundary scheme can be eliminated with the adjustment of the weight coefficients in BGK model. As for BB and NEE boundary schemes, when we consider BGK model(that is $s_{1}=s_{2}$), $\phi_{s}$ can be eliminated only if $s_{1}=1$. For this reason, the adjustment of ABB scheme is more flexible.

$\phi_{s}$ of ABB boundary scheme depends on $s_{1}$, $s_{2}$, $\theta$ and $c_{s}^{2}$, $\phi_{s}$ of BB boundary scheme depends on $s_{1}$ and $s_{2}$, and $\phi_{s}$ of NEE scheme only depends on $s_{1}$. We test the same problem with NEE scheme, taking $s_{1}=0.6,1.2,1.9$, with different $s_{2}$. As we shown in Fig. (5), $s_{2}$ has little effect on numerical results. And when we change the value of $s_{1}$, we can see that the GRE has a minimum when $s_{1}=1$, which agree with Eq. (61).

Then we consider $\phi_{s}$ of BB scheme as Eq. (54). When $s_{2}=2-s_{1}$, the discrete effect can be eliminated for the unidirectional steady problem with a parabolic distribution in one direction. We take a simulation of the same problem as Eq. (64) with BB scheme, taking $s_{1}=0.1,0.6,1.0,1.9$ respectively, and shown the result in Fig. (7). Under the same lattice size to eliminate the numerical slip in MRT model, we can adjust the parameter $s_{2}$ to satisfy $s_{1}+s_{2}=2$ for BB boundary scheme while in the BGK model $s_{2}$ is determined by diffusion coefficient, and can not be adjusted. As the figures shown, we can adjust $s_{2}$ to get more accurate results.

Then, we consider a three-dimensional linear time-independent diffusion equation with a constant source term,

[TABLE]

The analytical solution of this problem is given by

[TABLE]

Here we consider the popular D3Q19 BGK and MRT model, the physical parameters $L=1.0$, $u_{x}=0.1$, $u_{y}=0.0$, $u_{z}=0.0$, the diffusion coefficient $D=0.1$, the boundary conditions $\phi_{0}=0$, $\phi_{L}=1$, $\delta_{x}=L/N$ with the grid number $N$ varying from 5 to 17. For the ABB boundary scheme, we can adjust the parameter $s_{2}$ to satisfy Eq. (48) to get more accurate results. We perform some simulations with both BGK and MRT models, and present the results in Figs. 8, 9, 10 and 11. In these figures, the values of $s_{1}$ are taken to be $0.1$, $0.6$, $1.9$, and a particular value satisfying Eq. (48) under the condition of $s_{1}=s_{2}$. From the results in Figs. 8, 9, 10 and 11, one can see that when $s_{2}$ satisfies Eq. (48), the numerical results are in good agreement with analytical solutions.

Here we give some comparisons of the GRE and $E_{max}$ among D2Q5 and D2D9, D3Q7 and D3Q19 models in Tables. 4, 5, 6 and 7, and find that there are no apparent differences among D2Q5 and D2D9, D3Q7 and D3Q19 models when we adjust $s_{2}$ to satisfy Eq. (48) for ABB boundary scheme. However, the D2Q5 and D3Q7 models are more efficient since less discrete velocities are included.

4.1.2 Helmholtz equation

We also concidered the following linear Helmholtz equation, as

[TABLE]

with the boundary conditions

[TABLE]

The physical domain is $\Omega=[0,H]\times[0,H]$, $\lambda$ and $\mu$ are two constants. Under above conditions, steady analytical solution of Eq. (69) can be obtained

[TABLE]

which is more complicated than Eq. (64). We conducted some simulations with $\lambda=0$ and $\mu=1.0$, and present the results of D2Q9 MRT model under different values of $s_{1}$ in Figs. 12, 13, 14, where different weight coefficients are used. As we can see, the analytical solution Eq. (71) is time-independent and only depends on y when $\lambda=0$ and $\mu=1.0$. As shown in these figures, the relaxation parameter $s_{2}$ has a significant effect on numerical results, what is more, we can obtain the most accurate results when the value of $s_{2}$ determined by Eq. (48) is adopted.

4.2 A unidirectional time-independent complex-valued CDEs

In this part, we further considered a simple two-dimension complex-valued problem governed by Eq. (64) to verify Eqs. (48) and (54) where $D=1+i$, $R=4i$, $L=1.0$, $u_{x}=0.1$, $u_{y}=0.0$, and the boundary conditions $\phi_{0}=0$, $\phi_{L}=1+i$. In our simulations, $\delta_{x}={L}/{N}$ with the grid number $N$ varying from 5 to 17, the D2Q5 MRT model ($\theta=0$) is used.

The $\tau_{r}$, $\tau_{i}$ are the relaxation times of the real and the imaginary parts respectively, and $S_{r}=diag(s_{0},s_{r1},s_{r1},s_{r2},s_{r2})$ and $S_{i}=diag(s_{0},s_{i1},s_{i1},s_{i2},s_{i2})$ are the diagonal relaxation matrix. Then we have [38]

[TABLE]

where $D=D_{r}+iD_{i}$. In our simulations, we take $s_{0}$ = 0.0, $s_{r1}$ = 1.0, 10.0, 0.501, and $s_{i1}$ is determined by Eq. (72). Substituting $s_{1}=s_{r1}+is_{i1}$ and $s_{2}=s_{r2}+is_{i2}$ into Eq. (48), we have

[TABLE]

[TABLE]

where $a_{0}=\omega_{0}+2\omega_{1}$, $a_{1}=\omega_{1}$ in the D2Q5 model. The $s_{r2}$ and $s_{i2}$ are choose to satisfy Eqs. (73) and (74), and it shows a good accuracy in Tables. 8 and 9.

Then we take the same simulation with BB boundary scheme. $s_{2}$ is satisfied $s1+s_{2}=2$, that is $s_{r2}=2-s_{r1}$, $s_{i2}=-s_{r1}$ and shows the results in Tables. 10 and 11 which have good agreement with analytical solutions.

5 CONCLUSIONS

In this work, we performed a detailed analysis on the discrete effects of ABB, BB and NEE schemes of the popular one- to three- dimensional D$n$Q$q$ MRT LB model for real- and complex-valued CDEs. Firstly, through the analysis with ABB boundary scheme, we obtain a relation with four adjustable parameters the weight coefficient, the relaxation factors $s_{1}$ and $s_{2}$ associated with first and second moments and a model parameter $\theta$, which can be used to eliminate the discrete effect. We would also like to point out that taking $\theta=1$ under some assumption, the relation in [31] in the framework of TRT model would be the special case of Eq. (48). The weight coefficient $\omega$ can be considered as an adjustable parameter makes the general relation Eq. (48) more flexible. Then we analyse the discrete effects of BB and NEE boundary schemes and indicate that the discrete effect of BB scheme can be eliminated when $s_{1}+s_{2}=2$, and the discrete effect of NEE scheme can not be eliminated except $s_{1}=1$. The adjustment of ABB boundary scheme is more flexible than BB and NEE boundary schemes. We also carried out some numerical simulations of several special equations, including the real-valued linear time-independent diffusion equations in two- and three-dimensional space, the real-valued two-dimensional Helmholtz equation, and the complex-valued linear time-independent diffusion equation. The results also show that when the relation Eq. (48) for ABB boundary scheme and $s_{1}+s_{2}=2$ for BB boundary scheme is satisfied, the discrete effect (or numerical slip) can be eliminated.

Acknowledgements

This work is supported by the National Natural Science Foundation of China (Grants No. 51576079 and No. 51836003), and the National Key Research and Development Program of China (Grant No. 2017YFE0100100)

APPENDIX

5.1 Equivalent difference equation of the MRT model

In this Appendix, we show how to derive the equivalent difference equation. Firstly, for the D1Q3 MRT model, from Eq. (13), we can obtain the expressions of the distribution functions,

[TABLE]

where $f^{k}_{i}$, $f^{k,eq}_{i}$ are the distribution function and its equilibrium part at $x=k\delta_{x}$. According to Eqs. (12) and (3), we have

[TABLE]

[TABLE]

Substituting Eq. (76) into Eq. (75b), one can obtain

[TABLE]

Based on Eq. (78), we can get

[TABLE]

Substituting Eq. (79a) into Eq. (75c), and with the help of Eq. (78), we have

[TABLE]

Similarly, substituting Eq. (79b) into Eq. (75a), and with the aid of Eq. (78), one can obtain

[TABLE]

In addition, from Eqs. (80) and (81), we also have

[TABLE]

Summing Eqs. (82a) and (82b), one can derive the following equation,

[TABLE]

where Eqs. (80) and (81) have been used. Substituting Eq. (78) into Eq. (83) yields

[TABLE]

where Eq. (77) has been adopted. From Eq. (84), we can obtain the equivalent difference equation of the MRT model, i.e., Eq. (15).

For the D3Q27 model, we have

[TABLE]

where $f^{k}_{i,j,m}=f^{k}_{i}+f^{k}_{j}+f^{k}_{m}$, $f^{k,eq}_{i,j,m}=f^{k,eq}_{i}+f^{k,eq}_{j}+f^{k,eq}_{m}$. If the parts of $f^{k}_{0,1,3,5,6,15,16,17,18}$, $f^{k}_{2,8,9,13,14,21,22,25,26}$, and $f^{k}_{4,7,10,11,12,19,20,23,24}$ in the D3Q27 model are viewed as $f^{k}_{0}$, $f^{k}_{1}$, and $f^{k}_{-1}$ in the D1Q3 model, $w_{0}+4w_{1}+4w_{7}$ and $w_{1}+4w_{7}+4w_{19}$ in the D3Q27 model are considered as $w_{0}$ and $w_{1}$ in D1Q3 model, we can derive the equivalent different Eq. (15) through the similar process.

5.2 Discrete effect of the ABB boundary condition

In the D1Q3 model, when $k=1$, $k=N$, Eqs. (82b) and (82a) can be written as

[TABLE]

Substituting Eq. (79a) into Eq. (86a), substituting Eq. (79b) into Eq. (86b), we can obtain

[TABLE]

In addition, substituting Eqs. (77) and (82a) into Eqs. (87a) and (87b) gives rise to

[TABLE]

On the other hand, the ABB scheme can be given by

[TABLE]

Substituting Eq. (82b) into Eq. (89a), and substitute Eq. (82a) into Eq. (89b), one can obtain

[TABLE]

Substituting Eqs. (89a) and (89b) into Eqs. (90a) and (90b), we can obtain

[TABLE]

which can also be written as

[TABLE]

From Eq. (36), we have

[TABLE]

Substituting Eqs. (93a) and (93b) into Eq. (92a), and Eqs. (93c) and (93d) into Eq. (92b), we can obtain Eqs. (39) and (40).

6 References