Discrete unified gas kinetic scheme for nonlinear convection-diffusion equations

TL;DR

This paper introduces a discrete unified gas kinetic scheme (DUGKS) for nonlinear convection-diffusion equations, demonstrating its accuracy, stability, and flexibility with non-uniform meshes through analytical and numerical comparisons.

Contribution

The paper develops a novel DUGKS for NCDEs, showing its second-order convergence, compatibility with non-uniform meshes, and superior stability and accuracy over existing lattice Boltzmann models.

Findings

DUGKS achieves second-order convergence.

DUGKS is more accurate than FV-LBM.

DUGKS is more stable than SLBM.

Abstract

In this paper, we develop a discrete unified gas kinetic scheme (DUGKS) for general nonlinear convection-diffusion equation (NCDE), and show that the NCDE can be recovered correctly from the present model through the Chapman-Enskog analysis. We then test the present DUGKS through some classic convection-diffusion equations, and find that the numerical results are in good agreement with analytical solutions and the DUGKS model has a second-order convergence rate. Finally, as a finite-volume method, DUGKS can also adopt the non-uniform mesh. Besides, we performed some comparisons among the DUGKS, finite-volume lattice Boltzmann model (FV-LBM), single-relaxation-time lattice Boltzmann model (SLBM) and multiple-relaxation-time lattice Boltzmann model (MRT-LBM). The results show that the DUGKS model is more accurate than FV-LBM, more stable than SLBM, and almost has the same accuracy as the MRT-LBM. Besides, the using of non-uniform mesh may make DUGKS model more flexible.

Tables

Table 1. Table 1: A comparison of DUGKS, FV-LBM and MRT-LBM

	$Pe=100$	$Pe=1000$	$Pe={10}^7$	$Pe={10}^9$
DUGKS	$3.641\times {10}^{-4}$	$4.109\times {10}^{-4}$	$7.390\times {10}^{-5}$	$7.383\times {10}^{-5}$
FV-LBM	$1.543\times {10}^{-3}$	$1.431\times {10}^{-3}$	$1.436\times {10}^{-3}$	$1.436\times {10}^{-3}$
MRT-LBM	$3.265\times {10}^{-4}$	$1.709\times {10}^{-4}$	$1.453\times {10}^{-4}$	$1.452\times {10}^{-4}$

Table 2. Table 2: G ​ R ​ E 𝐺 𝑅 𝐸 GRE s of the DUGKS model with uniform grid (denoted as DUGKS a ) and rectangular grid (denoted as DUGKS b )

		$b=0.05$	$b=0.1$	$b=0.5$	$b=1.0$
$c=10$	DUGKSa	$1.0826\times {10}^{-3}$	$1.8261\times {10}^{-3}$	$3.5657\times {10}^{-3}$	$2.2113\times {10}^{-3}$
	DUGKSb	$1.1406\times {10}^{-3}$	$1.8489\times {10}^{-3}$	$3.6188\times {10}^{-3}$	$2.1893\times {10}^{-3}$
$c=20$	DUGKSa	$4.0612\times {10}^{-4}$	$5.3645\times {10}^{-4}$	$1.0287\times {10}^{-3}$	$6.9190\times {10}^{-4}$
	DUGKSb	$5.9501\times {10}^{-4}$	$5.8983\times {10}^{-4}$	$1.0694\times {10}^{-3}$	$7.3166\times {10}^{-4}$

Table 3. Table 3: The global relative errors with different values of α 𝛼 \alpha and c 𝑐 c at t = 1.0 𝑡 1.0 t=1.0 and δ = 1.2 𝛿 1.2 \delta=1.2

	$\alpha =0.01$	$\alpha =0.05$	$\alpha =0.1$	$\alpha =1.0$
$c=10$	$1.4858\times {10}^{-3}$	$4.5385\times {10}^{-4}$	$4.3775\times {10}^{-4}$	$3.9429\times {10}^{-4}$
$c=20$	$1.0398\times {10}^{-3}$	$2.5046\times {10}^{-4}$	$1.3507\times {10}^{-4}$	$1.0327\times {10}^{-4}$
$c=100$	$9.9290\times {10}^{-5}$	$1.5214\times {10}^{-4}$	$1.0280\times {10}^{-4}$	$2.5874\times {10}^{-5}$

Table 4. Table 4: The global relative errors at t = 1.0 𝑡 1.0 t=1.0 , δ = 1.2 𝛿 1.2 \delta=1.2 , c = 10.0 𝑐 10.0 c=10.0 and different values of α 𝛼 \alpha and meshes (mesh a denotes the uniform mesh with the lattice size 100 × 100 100 100 100\times 100 , mesh b denotes the non-uniform mesh with the lattice size 100 × 100 100 100 100\times 100 , mesh c denotes the non-uniform mesh with the lattice size 50 × 50 50 50 50\times 50 )

	$\alpha =0.01$	$\alpha =0.05$	$\alpha =0.1$	$\alpha =1.0$
mesha	$1.4858\times {10}^{-3}$	$4.5385\times {10}^{-4}$	$4.3775\times {10}^{-4}$	$3.9429\times {10}^{-4}$
meshb	$3.3793\times {10}^{-4}$	$3.2038\times {10}^{-4}$	$3.3250\times {10}^{-4}$	$3.6957\times {10}^{-4}$
meshc	$1.0070\times {10}^{-3}$	$4.1587\times {10}^{-4}$	$3.5491\times {10}^{-4}$	$3.7804\times {10}^{-5}$

Equations

$${\partial }_t\varphi +\nabla \cdot \mathbf{B}=\nabla \cdot (\alpha \nabla \cdot \mathbf{D})+F,$$

$$\frac{\partial f_i}{\partial t}+{\mathbf{c}}_i\cdot \nabla f_i={\Omega }_i+R_i+F_i,$$

$$\begin{array}{rl}& f_i^{eq}={\omega }_i\left[\varphi +\frac{{\mathbf{c}}_i\cdot \mathbf{B}}{c_s^2}+\frac{(c_s^2\mathbf{D}+\mathbf{C}-c_s^2\varphi \mathbf{I}):({\mathbf{c}}_i{\mathbf{c}}_i-c_s^2\mathbf{I})}{2c_s^4}\right],\\ & R_i={\omega }_i\frac{{\mathbf{c}}_i\cdot ({\partial }_t\mathbf{B}+\nabla \cdot \mathbf{C})}{c_s^2},\\ & F_i={\omega }_{i}F,\end{array}$$

$$\begin{array}{rl}& {\mathbf{c}}_i=(0,1,-1)c\\ & {\omega }_0=\frac{2}{3},{\omega }_1={\omega }_2=\frac{1}{6},c_s=\frac{c}{\sqrt{3}},\end{array}$$

$$\begin{array}{rl}& {\mathbf{c}}_i=\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\\ & {\omega }_0=\frac{4}{9},{\omega }_{1-4}=\frac{1}{9},{\omega }_{1-4}=\frac{1}{36},c_s=\frac{c}{\sqrt{3}},\end{array}$$

$$\text{centering}\begin{array}{rl}& {\mathbf{c}}_i=\left(\text{setcounter}MaxMatrixCols15\begin{array}{ccccccccccccccc}0& 1& -1& 0& 0& 0& 0& 1& 1& 1& -1& -1& -1& -1& 1\\ 0& 0& 0& 1& -1& 0& 0& 1& 1& -1& 1& -1& -1& 1& -1\\ 0& 0& 0& 0& 0& 1& -1& 1& -1& 1& 1& -1& 1& -1& -1\end{array}\right)c\\ & {\omega }_0=\frac{2}{9},{\omega }_{1-6}=\frac{1}{9},{\omega }_{7-14}=\frac{1}{72},c_s=\frac{c}{\sqrt{3}}.\end{array}\text{@add@centering}$$

$$\text{centering}\begin{array}{rl}& {\mathbf{c}}_i=\left(\text{setcounter}MaxMatrixCols19\begin{array}{ccccccccccccccccccc}0& 1& -1& 0& 0& 0& 0& 1& -1& 1& -1& 1& -1& -1& 1& 0& 0& 0& 0\\ 0& 0& 0& 1& -1& 0& 0& 1& -1& -1& 1& 0& 0& 0& 0& 1& -1& 1& -1\\ 0& 0& 0& 0& 0& 1& -1& 0& 0& 0& 0& 1& -1& 1& -1& 1& -1& -1& 1\end{array}\right)c\\ & {\omega }_0=\frac{1}{3},{\omega }_{1-6}=\frac{1}{18},{\omega }_{7-18}=\frac{1}{36},c_s=\frac{c}{\sqrt{3}}.\end{array}\text{@add@centering}$$

$$\begin{array}{rl}& \sum _i{f}_i=\sum _i{f}_i^{eq}=\varphi ,\sum _i{R}_i=0,\sum _i{F}_i=F,\\ & \sum _i{\mathbf{c}}_i{f}_i^{eq}=\mathbf{B},\sum _i{\mathbf{c}}_i{R}_i={\partial }_t\mathbf{B}+\nabla \cdot \mathbf{C},\sum _i{\mathbf{c}}_i{F}_i=0,\\ & \sum _i{\mathbf{c}}_i{\mathbf{c}}_i{f}_i^{eq}=c_s^2\mathbf{D}+\mathbf{C}.\end{array}$$

$$f_i^{n+1}-f_i^n+\frac{\Delta t}{|V_j|}{J}^{n+1/2}=\frac{\Delta t}{2}({\Omega }_i^{n+1}+{\Omega }_i^n)+\Delta t(R_i^n+\frac{\Delta t}{2}{\partial }_t{R}_i^n)+\Delta t(F_i^n+\frac{\Delta t}{2}{\partial }_t{F}_i^n),$$

$$J^{n+1/2}={\int }_{\partial V_j}({\mathbf{c}}_i\cdot \mathbf{n})f_i(\mathbf{x},{\mathbf{c}}_i,t_{n+1/2})d\mathbf{S}$$

$$A_i^n=\frac{1}{|V_j|}{\int }_{V_j}{A}_i({\mathbf{x}}_j,{\mathbf{c}}_i,t_n)d\mathbf{x},A\in \left\{f,R,F,\Omega \right\}$$

$${\tilde{f}}_i=f_i-\frac{\Delta t}{2}{\Omega }_i=\frac{2\lambda +\Delta t}{2\lambda }{f}_i-\frac{\Delta t}{2\lambda }{f}_i^{eq}.$$

$${\tilde{f}}_i^{n+1}={\tilde{f}}_i^{+,n}-\frac{\Delta t}{|V_j|}{J}^{n+1/2}+\Delta t(R_i^n+\frac{\Delta t}{2}{\partial }_t{R}_i^n)+\Delta t(F_i^n+\frac{\Delta t}{2}{\partial }_t{F}_i^n),$$

$${\tilde{f}}_i^{+}=\frac{2\lambda -\Delta t}{2\lambda +\Delta t}{\tilde{f}}_i+\frac{2\Delta t}{2\lambda +\Delta t}{f}_i^{eq}.$$

$${\tilde{f}}_i^{n+1}={\tilde{f}}_i^{+,n}-\frac{\Delta t}{|V_j|}{J}^{n+1/2},$$

$$\begin{array}{rl}& {\tilde{f}}_i=f_i-\frac{\Delta t}{2}{\Omega }_i-\frac{\Delta t}{2}{R}_i-\frac{\Delta t}{2}{F}_i,\\ & {\tilde{f}}_i^{+}=\frac{2\lambda -\Delta t}{2\lambda +\Delta t}{\tilde{f}}_i+\frac{2\Delta t}{2\lambda +\Delta t}{f}_i^{eq}+\frac{2\lambda \Delta t}{2\lambda +\Delta t}{R}_i+\frac{2\lambda \Delta t}{2\lambda +\Delta t}{F}_i.\end{array}$$

$$\begin{array}{rl}& f_i({\mathbf{x}}_b,{\mathbf{c}}_i,t_n+h)-f_i({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)=\frac{h}{2}[{\Omega }_i({\mathbf{x}}_b,{\mathbf{c}}_i,t_n+h)\\ & +{\Omega }_i({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)]+h[R_i({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)+\frac{h}{2}{D}_i{R}_i({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)]\\ & +h[F_i({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)+\frac{h}{2}{D}_i{F}_i({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)],\end{array}$$

$${\bar{f}}_i=f_i-\frac{h}{2}{\Omega }_i=\frac{2\lambda +h}{2\lambda }-\frac{h}{2\lambda }{f}_i^{eq}.$$

$$\begin{array}{rl}& {\bar{f}}_i({\mathbf{x}}_b,{\mathbf{c}}_i,t_n+h)={\bar{f}}_i^{+}({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)+h[R_i({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)\\ & +\frac{h}{2}{D}_i{R}_i({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)]+h[F_i({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)+\frac{h}{2}{D}_i{F}_i({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)],\end{array}$$

$${\bar{f}}_i^{+}=\frac{2\lambda -h}{2\lambda +h}{\bar{f}}_i+\frac{2h}{2\lambda +h}{f}_i^{eq}.$$

$${\bar{f}}_i({\mathbf{x}}_b,{\mathbf{c}}_i,t_n+h)={\bar{f}}_i^{+}({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n),$$

$$\begin{array}{rl}& {\bar{f}}_i=f_i-\frac{h}{2}{\Omega }_i-\frac{h}{2}{R}_i-\frac{h}{2}{F}_i,\\ & {\bar{f}}_i^{+}=\frac{2\lambda -h}{2\lambda +h}{\bar{f}}_i+\frac{2h}{2\lambda +h}{f}_i^{eq}+\frac{2\lambda h}{2\lambda +h}{R}_i+\frac{2\lambda h}{2\lambda +h}{F}_i.\end{array}$$

$$\begin{array}{rl}& {\bar{f}}_i({\mathbf{x}}_b,{\mathbf{c}}_i,t_n+h)={\bar{f}}_i^{+}({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)+\frac{2\lambda h}{2\lambda +h}[R_i({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)\\ & +\frac{h}{2}{\partial }_t{R}_i({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)]+h[F_i({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)+\frac{h}{2}{\partial }_t{F}_i({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)].\end{array}$$

$$\begin{array}{rl}& {\bar{f}}_i^{+}({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)={\bar{f}}_i^{+}({\mathbf{x}}_b,{\mathbf{c}}_i,t_n)-h{\mathbf{c}}_i\cdot \nabla {\bar{f}}_i^{+}({\mathbf{x}}_b,{\mathbf{c}}_i,t_n),\\ & R_i({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)=R_i({\mathbf{x}}_b,{\mathbf{c}}_i,t_n)-h{\mathbf{c}}_i\cdot \nabla R_i({\mathbf{x}}_b,{\mathbf{c}}_i,t_n),\\ & F_i({\mathbf{x}}_b-{\mathbf{c}}_{i}h,{\mathbf{c}}_i,t_n)=F_i({\mathbf{x}}_b,{\mathbf{c}}_i,t_n)-h{\mathbf{c}}_i\cdot \nabla F_i({\mathbf{x}}_b,{\mathbf{c}}_i,t_n),\end{array}$$

$$\begin{array}{rl}& \nabla {\bar{f}}_i^{+}(x_{j+1/2},{\mathbf{c}}_i,t_n)=\frac{{\bar{f}}_i^{+}(x_{j+1},{\mathbf{c}}_i,t_n)-{\bar{f}}_i^{+}(x_j,{\mathbf{c}}_i,t_n)}{x_{j+1}-x_j},\\ & {\bar{f}}_i^{+}(x_{j+1/2},{\mathbf{c}}_i,t_n)={\bar{f}}_i^{+}(x_j,{\mathbf{c}}_i,t_n)+(x_{j+1/2}-x_j)\nabla {\bar{f}}_i^{+}(x_{j+1/2},{\mathbf{c}}_i,t_n).\end{array}$$

$${\bar{f}}_i^{+}=\frac{2\lambda -h}{2\lambda +\Delta t}{\tilde{f}}_i+\frac{3h}{2\lambda +\Delta t}{f}_i^{eq},$$

$$f_i=\frac{2\lambda }{2\lambda +h}\overline{f_i}+\frac{h}{2\lambda +h}{f}_i^{eq},$$

$${\tilde{f}}_i^{+}=\frac{4}{3}{\bar{f}}_i^{+}-\frac{1}{3}{\tilde{f}}_i.$$

$$\begin{array}{rl}\varphi ({\mathbf{x}}_j,t),F({\mathbf{x}}_j,t)& \stackrel{(\text{ref}eq:2_3)}{\to }{\tilde{f}}_i({\mathbf{x}}_j,t),R_i({\mathbf{x}}_j,t),F_i({\mathbf{x}}_j,t)\\ & \stackrel{(\text{ref}eq:2_{1}3),(\text{ref}eq:2_{2}1)}{\to }{\tilde{f}}_i^{+}({\mathbf{x}}_j,t),{\bar{f}}_i^{+}({\mathbf{x}}_j,t)\\ & \stackrel{(\text{ref}eq:2_{2}0)}{\to }{\bar{f}}_i^{+}({\mathbf{x}}_b,t),R_i({\mathbf{x}}_b,t),F_i({\mathbf{x}}_b,t)\\ & \stackrel{(\text{ref}eq:2_{1}8),(\text{ref}eq:2_{1}9)}{\to }{\bar{f}}_i({\mathbf{x}}_b,t+h)\to \varphi ({\mathbf{x}}_b,t+h)\\ & \stackrel{(\text{ref}eq:2_{2}2)}{\to }{f}_i({\mathbf{x}}_b,t+h)\stackrel{(\text{ref}eq:2_9)}{\to }J({\mathbf{x}}_b,t+h)\\ & \stackrel{(\text{ref}eq:2_{1}2)}{\to }{\tilde{f}}_i({\mathbf{x}}_j,t+\Delta t)\to \varphi ({\mathbf{x}}_j,t+\Delta t),F({\mathbf{x}}_j,t+\Delta t)\end{array}$$

$$\begin{array}{rl}& f_i=f_i^{(0)}+ϵf_i^{(1)}+{ϵ}^2{f}_i^{(2)},R_i=ϵR_i^{(1)}+{ϵ}^2{R}_i^{(2)},F_i=ϵF_i^{(1)}+{ϵ}^2{F}_i^{(2)},\\ & {\partial }_t=ϵ{\partial }_{t_1}+{ϵ}^2{\partial }_{t_2},\nabla =ϵ{\nabla }_1,\end{array}$$

Full text

1 Introduction

The convection-diffusion equation (CDE) is usually used to describe the physical phenomena where particles, energy or other physical quantities are transferred inside a physical system, and in particular, plays an important role in the field of heat and mass transfer [1]. However, as a kind of partial differential equation (PDE), CDE is usually so complicated that it is diffucult to get the analytical solution most of time. With the development of computing power, some numerical methods have been developed to solve CDEs, such as finite-element method [2], finite-difference method [3] and finite-volume method [4].

In the past decades, the lattice Boltzmann method (LBM), as a mesoscopic numerical approach, has achieved great success in the simulation of hydrodynamic problems [5, 6, 7, 8, 9, 10, 11]. On the other hand, the LBM has also been extended to solve the CDEs. Dawson et al. [12] first proposed a LB model for CDE, but the model cannot give correct CDE. Shi and Guo [13] developed a lattice Bhatnagar-Gross-Krook (LBGK) model to solve the general nonlinear convection-diffusion equations (NCDEs), where an auxiliary moment $\mathbf{C}$ is used to correctly recover the NCDE. However, in their work, the convection term $\mathbf{B}$ should be a function of $\phi$. Chopard [14] developed a new LB model where a source term related to temporal derivative or spatial derivative is adopted to give correct CDE. We noted that LB models are limited to the isotropic CDEs. To solve the nonlinear anisotropic convection-diffusion equations (NACDEs). The two-relaxation-time (TRT) and multiple-relaxation-time (MRT) LB models are considered by Ginzburg [15, 16, 17, 18, 19], while in these models, some assumptions on the convection and diffusion terms, and the assumptions may not be satisfied for some special NACDEs. Yoshida and Nagaoka [20] also developed a MRT LB model, and did some analysis on different boundary conditions, however, the assumptions were also adopted to recover the CDE. Recently, Chai et al. [21] presented a MRT LB model for general NACDEs without any assumptions on the convection and diffusion terms. Although everything looks perfect, some restrictions still exist in all the above LB models. The first is that the temporal and spatial steps are coupled, causing the selection of parameters to be very limited. The second is that all above LB models must be implemented on uniform grid.

Recently, Guo et al. [22] proposed the discrete unified gas kinetic scheme (DUGKS) for all Knudsen number flows. The DUGKS combines the advantages of LBM and unified gas kinetic scheme (UGKS). Firstly, as a finite volume scheme, DUGKS can adopt the flexible mesh. Secondly, the DUGKS is more accurate than finite-volume LBM, this is because the evaluation of the flux at cell interface is simplified by employing a transformation of distribution function with collision effect, which has also been used in LBM. Finally, the asymptotic preserving (AP) property still exists in the DUGKS. It should be noted that at the beginning, the DUGKS in Ref. [22] is developed based on the Bhatnagar-Gross-Krook (BGK) collision model [23], and the source term is mot included. Then, Wu et al. [24] developed a DUGKS with a force term for incompressible fluid flows and also presented the non-equilibrium extrapolation (NEE) scheme for DUGKS. Recently, Zhang et al. [25] and Yang et al. [26] developed the phase-field based DUGKS for two-phase flows, the difference between their two works is that the Chan-Hilliard (CH) equation [27, 28] is considered in Ref. [25] while the Allen-Cahn (AC) equation [29] is adopted in Ref. [26]. In these works, the DUGKS was used to solve the phase field equations. Huo and Rao [30] uesd the DUGKS to study the solid-liquid phase change problem, in which the energy equation was solved by DUGKS. From above discussion, the DUGKS has been widely used to study single and two-phase flows, and also the phase-field and energy equations. However, it is unclear whether the phase-field and energy equations as some special types of CDEs can be recovered from the DUGKS. Through the DUGKS, we found that the above restrictions of LBM solving CDEs are avoided, perfectly. So whether we can use the present DUGKS model to solve the more general partial differential equations (PDEs)? In this work, we will develop a DUGKS for general NCDEs, and also perform a detailed Chapman-Enskog analysis.

The rest of the paper is organized as follow. In Sec. 2, the DUGKS for the general NCDE is proposed. In Sec. 3, through the Chapman-Enskog analysis, the NCDE is recovered correctly from the present DUGKS. In addition, some special cases and distinct characteristics are also discussed. In Sec. 4, the accuracy and convergence rate of the DUGKS model are tested through some classic CDEs, and some comparisons among the present DUGKS, finite-volume LB model, LBGK model and MRT LB model are conducted. Finally, some conclusions are given in Sec. 5.

2 The DUGKS model for general NCDEs

In this section, we will present a DUGKS for n-dimensional NCDE with variable coefficients

[TABLE]

where $\phi$ is a scalar function of position $\mathbf{x}$ and time t, $\nabla$ is the gradient operator with respect to the position $\mathbf{x}$ in n dimensions. $\mathbf{B}$ and $\mathbf{D}$ are the known convection and diffusion terms, and usually they are related to position $\mathbf{x}$, $\phi$, time t. $\alpha$ and $F$ are the diffusion coefficient and source term, respectively.

Following the idea in the previous work [22], the DUGKS with DnQq lattice (q is the number of discrete directions) for the NCDE is considered here. First, the discrete velocity Boltzmann equation (DBE) can be written as

[TABLE]

where $f_{i}=f_{i}(\mathbf{x},\mathbf{c}_{i},t)$ is the particle distribution function with discrete velocity $\mathbf{c}_{i}$ at time t and position $\mathbf{x}$. $\Omega_{i}=-(f_{i}-f_{i}^{eq})/\lambda$ is the Bhatnagar-Gross-Krook (BGK) collision model [13], $\lambda$ is the relaxation time. $f_{i}^{eq}$ is the equilibrium distribution function, $R_{i}$ and $F_{i}$ are the distribution functions of source term. To derive correctly NCDE (1) from present DUGKS, the distribution functions $f_{i}^{eq}$, $R_{i}$ and $F_{i}$ are given by

[TABLE]

where $\mathbf{I}$ is the unit matrix, $\mathbf{C}$ is a tensor function which can be set to be [math] or $\int\mathbf{B}^{\prime}(\phi)\mathbf{B}^{\prime}(\phi)d\phi$ [13]. $c_{s}$ is the so called sound speed related to discrete velocity. $\omega_{i}$ and $\mathbf{c}_{i}$ are weight coefficient and discrete velocity, and in different discrete velocity models, they can be defined as

D1Q3:

[TABLE]

D2Q9:

[TABLE]

D3Q15:

[TABLE]

D3Q19:

[TABLE]

Based on the conservation law and Eq. (3), we have

[TABLE]

In the DUGKS, we divide the computational domain into a set of control volumes (cells), and $\mathbf{x}_{j}$ is used to denote the cell $j$. Then, integrating Eq. (2) over the volume $V_{j}$ from $t_{n}$ to $t_{n+1}$, and using the midpoint rule, trapezoidal rule and Taylor expansion for the integration of the flux term at cell interface, collision term and source terms inside the cell, we can obtain

[TABLE]

where

[TABLE]

is the flux of cell $j$, $\partial V_{j}$ and $|V_{j}|$ are the surface area and volume of cell $j$, $\mathbf{n}$ is the outward unit vector to the surface, $\Delta t=\alpha\frac{\Delta x}{c}$ is the time step and it is only determined by the Courant-Friedrichs-Lewy (CFL) condition ($\alpha$ is the CFL number and lies between 0 and 1). It should be noted that $f_{i}^{n}$, $R_{i}^{n}$, $F_{i}^{n}$ and $\Omega_{i}^{n}$ in Eq. (9) are the cell-averaged values of the distribution functions and collision term, respectively, i.e.,

[TABLE]

Because the collision term $\Omega_{i}^{n+1}$ involves the unknown variables at $t_{n+1}$, thus the evolution equation Eq. (9) is implicit scheme. In order to remove the implicity, a new distribution function is adopted,

[TABLE]

Then Eq. (9) can be rewritten as

[TABLE]

where

[TABLE]

Based on Eqs. (8) and (12), the conserved variable $\phi$ can be computed by $\phi=\sum_{i}\tilde{f}_{i}$. With this fact, we only need to track the distribution function $\tilde{f}_{i}$ instead of $f_{i}$ in practical computation. Besides, from the computational point of view, if we use $(R_{i}^{n}-R_{i}^{n-1})/\Delta t$ and $(F_{i}^{n}-F_{i}^{n-1})/\Delta t$ to evaluate the values of $\partial_{t}R_{i}^{n}$ and $\partial_{t}F_{i}^{n}$, Eq. (13) would become an explicit format. However, if $(R_{i}^{n+1}-R_{i}^{n})/\Delta t$ and $(F_{i}^{n+1}-F_{i}^{n})/\Delta t$ are used to estimate $\partial_{t}R_{i}^{n}$ and $\partial_{t}F_{i}^{n}$, we can rewrite Eq. (13) as

[TABLE]

with

[TABLE]

This method is the same as Ref. [24]. If necessary, we can also use $(F_{i}^{n}-F_{i}^{n-1})/\Delta t$ to estimate $\partial_{t}F_{i}^{n}$ and use $(R_{i}^{n+1}-R_{i}^{n})/\Delta t$ to estimate $\partial_{t}R_{i}^{n}$.

Now, the key ingredient in updating $\tilde{f}_{i}$ according to Eq. (13) is to evaluate the flux $J^{n+1/2}$. From Eq. (10), we can see that the flux $J^{n+1/2}$ is only determined by the original distribution $f_{i}$ at time $t+\frac{\Delta t}{2}$. In order to compute $f_{i}(\mathbf{x},\mathbf{c}_{i},t_{n+1/2})$, we integrate the discrete velocity Boltzmann equation within a half time step $h=\Delta t/2$ along the characteristic line with the point ($\mathbf{x}_{b}$) located at the cell interface ($\mathbf{x}_{b}=\mathbf{x}_{j+1/2}$ in the one-dimensional case, see Fig. 1),

[TABLE]

where the trapezoidal rule and Taylor expansion are applied to evaluate the collision term and source terms, $D_{i}=\partial_{t}+\mathbf{c}_{i}\cdot\nabla$. Then, similar to the treatment in Eq. (12), another distribution function $\bar{f}$ is introduced to remove the implicity of Eq. (17).

[TABLE]

As a result, the implicit formulation Eq. (17) can be rewritten by

[TABLE]

where

[TABLE]

Here, we can adopt $[R_{i}(\mathbf{x}_{b},\mathbf{c}_{i},t_{n}+h)-R_{i}(\mathbf{x}_{b}-\mathbf{c}_{i}h,\mathbf{c}_{i},t_{n})]/h$ and $[F_{i}(\mathbf{x}_{b},\mathbf{c}_{i},t_{n}+h)-F_{i}(\mathbf{x}_{b}-\mathbf{c}_{i}h,\mathbf{c}_{i},t_{n})]/h$ to evaluate the values of $D_{i}R_{i}(\mathbf{x}_{b}-\mathbf{c}_{i}h,\mathbf{c}_{i},t_{n})$ and $D_{i}F_{i}(\mathbf{x}_{b}-\mathbf{c}_{i}h,\mathbf{c}_{i},t_{n})$, and Eq. (19) can be rewritten as

[TABLE]

with

[TABLE]

The method is the same as Ref. [24]. If the source term is a nonlinear function of the variable $\phi$, we can use the explicit difference method to avoid solving the nonlinear equations, but, the finite-difference scheme for gradient term would destroy the locality of the DUGKS.

Actually, based on the previous works [13, 21], we can rewrite the evolution equation Eq. (19) as

[TABLE]

This evolution equation can avoid calculating the gradient term, and we will explain why we can do this in next section. Now, we focus on the computation of the distribution function $\bar{f}_{i}^{+}(\mathbf{x}_{b}-\mathbf{c}_{i}h,\mathbf{c}_{i},t_{n})$, $R_{i}(\mathbf{x}_{b}-\mathbf{c}_{i}h,\mathbf{c}_{i},t_{n})$ and $F_{i}(\mathbf{x}_{b}-\mathbf{c}_{i}h,\mathbf{c}_{i},t_{n})$.

With the Taylor expansion, the cell interface $\mathbf{x}_{b}$, the distribution function $\bar{f}_{i}^{+}(\mathbf{x}_{b}-\mathbf{c}_{i}h,\mathbf{c}_{i},t_{n})$, $R_{i}(\mathbf{x}_{b}-\mathbf{c}_{i}h,\mathbf{c}_{i},t_{n})$ and $F_{i}(\mathbf{x}_{b}-\mathbf{c}_{i}h,\mathbf{c}_{i},t_{n})$ can be approximated as

[TABLE]

where the distribution functions at $\mathbf{x}_{b}$ and the gradient terms can be approximated by linear interpolations, respectively. For example, as shown in Fig. 1, in one-dimensional case, the reconstructions become

[TABLE]

The estimations of $R_{i}(\mathbf{x}_{b},\mathbf{c}_{i},t_{n})$, $F_{i}(\mathbf{x}_{b},\mathbf{c}_{i},t_{n})$, $\nabla R_{i}(\mathbf{x}_{b},\mathbf{c}_{i},t_{n})$ and $\nabla F_{i}(\mathbf{x}_{b},\mathbf{c}_{i},t_{n})$ are similar to Eq. (25). From Eqs. (18) and (23), we can obtain the conserved variable at the cell interface $\mathbf{x}_{b}$, $\phi(\mathbf{x}_{b},t_{n}+h)=\sum_{i}\bar{f}_{i}(\mathbf{x}_{b},\mathbf{c}_{i},t_{n}+h)$.

Besides, we can get the relationship of distribution functions $f_{i}$, $f_{i}^{eq}$, $\bar{f}_{i}$, $\tilde{f}_{i}$, $\bar{f}_{i}^{+}$ and $\tilde{f}_{i}^{+}$ from Eqs. (12), (14), (18) and (20),

[TABLE]

[TABLE]

[TABLE]

The update of conserved variable $\phi$ in one time step of the present DUGKS can be summarized as follows:

[TABLE]

3 The Chapman-Enskog Analysis

In this part, the present DUGKS for NCDE is analyzed through the Chapman-Enskog (CE) analysis. In general, the using of Chapman-Enskog analysis in LBM is to recover the macroscopic equations from evolution equations [13, 21]. However, we found that the using of CE analysis in DUGKS is to recover the macroscopic equations from DBE [25]. The results of the two ways are the same, such as equilibrium distribution function, moment conditions and so on. In the following, the NCDE will be exactly recovered from Eq. (2) and (23), respectively.

Firstly, we expand the distribution functions $f_{i}$, $R_{i}$, $F_{i}$, the derivatives of time and space as

[TABLE]

where $\epsilon$ is a small parameter and keeps the same order of the Knudsen number.

Substituting Eq. (29) into Eq. (2), some equations at different orders of $\epsilon$ are obtained,

[TABLE]

Summing Eq. (30) over $i$, we can get

[TABLE]

From Eq. (30), we can obtain $\sum_{i}\mathbf{c}_{i}f_{i}^{(1)}$,

[TABLE]

Substituting Eq. (32) into Eq. (31) yields

[TABLE]

Using the conditions in Eq. (8), we have

[TABLE]

Combining above equations at the orders of $O(\epsilon^{1})$ and $O(\epsilon^{2})$, and taking $\alpha=c_{s}^{2}\lambda$, we can get the Eq. (1).

In addition, we redesign Eq. (19) to Eq. (23), and show that Eq. (1) can be exactly recovered from Eq. (23). We expand the distribution functions $\bar{f}_{i}$, $R_{i}$, $F_{i}$, the derivatives of time and space, the same as Eq. (29),

[TABLE]

By applying Taylor expansion to Eq. (23), we have

[TABLE]

where $D_{i}=\epsilon D_{1i}+\epsilon^{2}\partial_{t_{2}}$ and $D_{1i}=\partial_{t_{1}}+\mathbf{c}_{i}\cdot\nabla_{1}$.

Substituting Eq. (35) into Eq. (36), we can derive the following equations at different orders of $\epsilon$,

[TABLE]

Summing Eq. (37) over $i$ and using Eq. (8), we can get

[TABLE]

From Eq. (30), we have

[TABLE]

Then, substituting Eq. (39) into Eq. (38), we can obtain Eq. (34). Similarly, combining equations at the orders $O(\epsilon^{1})$ and $O(\epsilon^{2})$, and taking $\alpha=c_{s}^{2}\lambda$, the NCDE (1) is correctly recovered from Eq. (23).

*Remark 1.*As a model only include one relaxation time, DUGKS is different from single-relaxation-time lattice Boltzmann model (SLBM). The spatial and temporal step in DUGKS are decoupled, the temporal step in DUGKS is determined by Courant-Friedrichs-Lewy(CFL) condition($\Delta t=\alpha\frac{\Delta x}{c}$). In SLBM, however, the spatial and temporal steps are coupled through $\Delta\mathbf{x}=\mathbf{c}_{i}\Delta t$. Therefore, the restriction in SLBM does not exist in DUGKS.

*Remark 2.*As a finite volume scheme, the present DUGKS is different from finite-volume lattice Boltzmann Method (FV-LBM). In DUGKS, the flux $J$ is appropriated by $f_{i}(t_{n+1/2})$ instead of $f_{i}(t_{n})$ in FV-LBM. The analysis in Refs. [31, 32] show that the FV-LBM may suffer from severe numerical dissipation.

*Remark 3.*The tensor function $\mathbf{C}$ in Eq. (3) is an auxiliary-moment. If $\mathbf{B}$ is the function of $\phi$, $u$, $\mathbf{x}$ and t, we can define $\mathbf{C}=0$ so that $R_{i}=\omega_{i}\frac{\mathbf{c}_{i}\cdot\partial_{t}\mathbf{B}}{c_{s}^{2}}$. If $\mathbf{B}$ is only a function of $\phi$, we can define $\mathbf{C}=\int\mathbf{B}^{\prime}(\phi)\mathbf{B}^{\prime}(\phi)d\phi$ so that $R_{i}=\omega_{i}\frac{\mathbf{c}_{i}\cdot\mathbf{B}^{\prime}(\phi)F}{c_{s}^{2}}$. In the second case, we do not have to calculate the temporal derivative, and additionally, the equilibrium distribution function $f_{i}^{eq}$ can also be simplified by $f_{i}^{eq}=\omega_{i}(\phi+\frac{\mathbf{c}_{i}\cdot\mathbf{B}}{c_{s}^{2}})$ if $\mathbf{D}$ is just a function of $\phi$, too. Furthermore, when the linear equilibrium distribution function $f_{i}^{eq}=\omega_{i}(\phi+\frac{\mathbf{c}_{i}\cdot\mathbf{B}}{c_{s}^{2}})$ is considered, we can adopt the DdQ2d+1 discrete velocity model, for instance, D1Q3, D2Q5 and D3Q7.

4 Numerical results and discussion

In this part, some examples, including isotropic CDE with a constant velocity, Burgers-Fisher equation, the nonlinear heat conduction equation (NHCE), Gaussian hill problem and CDE with nonlinear convection and diffusion terms, are adopted to test the accuracy and stability of the present DUGKS. In our simulations, the distribution function $\tilde{f}_{i}$ is initialized by the equilibrium distribution function $f_{i}^{eq}$, i.e., $\tilde{f}_{i}(\mathbf{x},t_{0})=f_{i}^{eq}(\mathbf{x},t_{0})$. Unless otherwise stated, the non-equilibrium extrapolation scheme [24] is used to treat the boundary conditions. The following global relative error ($GRE$) is used to measure the accuracy of the present DUGKS,

[TABLE]

where $\phi_{a}$ and $\phi_{n}$ are the analytical and numerical solutions. In addition, to obtain stable results with present DUGKS, the CFL condition number should be less than 1.

Example 4.1 Two-dimensional isotropic CDE with a constant velocity can be expressed as

[TABLE]

where $u_{x}$ and $u_{y}$ are constants, and set to be 0.1, $\alpha$ is the diffusion coefficient. $F$ is the source term, and is given by

[TABLE]

Under the periodic boundary and following initial conditions,

[TABLE]

the solution of the problem can be expressed as

[TABLE]

When the present DUGKS is used to study this problem, the functions $\mathbf{B}$, $\mathbf{C}$ and $\mathbf{D}$ are given by $\mathbf{B}=\phi\mathbf{u}$ with $\mathbf{u}=(u_{x},u_{y})^{T}$, $\mathbf{C}=\phi\mathbf{u}\mathbf{u}$ and $\mathbf{D}=\phi\mathbf{I}$.

Now, we performed some simulations under different Péclet numbers and different time, where $Pe=Lu_{x}/\alpha$, $L$ is the characteristic length (here $L=2.0$), and the CFL condition number is equal to 0.5. The results are presented in Fig. 2 where $c=1.0$, the uniform grid is $200\times 200$, $\alpha$ can be determined by the specified $Pe$ (100 or 1000). As seen from the figure, the numerical solutions are in good agreement with analytical solutions. Besides, we also measured the values of $GRE$ at time $t=3.0$, and they are $3.641\times 10^{-4}$ for $Pe=100$ and $4.109\times 10^{-4}$ for $Pe=1000$. In addition, to test the capacity of present DUGKS for this problem with a larger $Pe$, some simulations were performed with $Pe=10^{7}$ and $10^{9}$, and the results are presented in Fig. 3. From the Figure, we can find that the numerical solutions still agree well with the analytical solutions, and the values of $GRE$ at time $t=3.0$ are $7.390\times 10^{-5}$ for $Pe=10^{7}$ and $7.383\times 10^{-5}$ for $Pe=10^{9}$. It is clearly that the deviations are small enough. We also performed a comparison among DUGKS, FV-LBM and MRT-LBM under the same conditions, and listed the results in Table 1. As we can see from this table, the performance of FV-LBM is worst, which is mainly caused by the severe numerical dissipation. Besides, the accuracies of DUGKS and MRT-LBM are almost the same.

Finally, the problem is applied to test the convergence rate of the present DUGKS. Since it is a periodic problem, the effect of the boundary conditions can be excluded. To this end, we carried out some simulations with different lattice sizes ($25\times 25$ $\sim$ $200\times 200$), and the time step is fixed at $1.0\times 10^{-5}$. As seen from Fig. 4, that the present DUGKS has a second-order convergence rate.

Example 4.2 The Burgers-Fisher equation in two dimensions [21] can be written as

[TABLE]

The analytical solution of Eq. (45) can be given by [21]

[TABLE]

where $A=-\frac{a\delta}{4b(\delta+1)}$, $m=\frac{a^{2}+2bk(\delta+1)^{2}}{a(\delta+1)}$, $a$, $b$, $k$ and $\delta$ are constants. Different from the first problem, this problem is nonlinear, and boundary conditions are nonperiodic.

For this problem, $\mathbf{B}=(\frac{a}{\delta+1},0)^{T}\phi^{\delta+1}$, and the simulations are performed on [-1,2] $\times$ [-1,2] with a 300 $\times$ 300 uniform grid size. In our simulations, parameters are set as $\delta=2.0$, $k=1.0$, $c=10$, $a=6.0$, $b=0.05$, and the CFL condition number is equal to 0.1. We presented the result in Fig. 5, and found that the numerical solutions are in good agreement with the corresponding analytical solutions.

To test the convergence rate of DUGKS for this problem, some simulations were carried out at time $t=1.0$, the lattice sizes are varied from $25\times 25$ to $100\times 100$, and time step $\Delta t=1.0\times 10^{-5}$. From the results in Fig. 6, it is clearly that the present DUGKS has a second-order convergence rate in space.

As a finite-volume scheme, the DUGKS has the distinct advantage in adopting the non-uniform mesh. To show the advantage more clearly, we also performed some simulations on rectangular grid (lattice size is $300$ $\times$ $150$), and the other parameters are the same as above. We presented a comparison between uniform and non-uniform grids in Table 2 where time $t=0.5$. As seen from this table, the errors of DUGKS with rectangular grid and those with uniform grid are of the same accuracy. While, computational cost of DUGKS with uniform grid ($300\times 300$) is about twice as that of DUGKS rectangular grid ($300\times 150$), the model on rectangular grid is more efficient than that on uniform grid.

Example 4.3 The generalized two dimensions NHCE in [13]

[TABLE]

has the following analytical solution,

[TABLE]

where $\alpha$ and $\delta$ are constants.

For this problem, we take $\mathbf{B}=0$ and $\mathbf{D}=\phi^{\delta}\mathbf{I}$, which leads to the following equilibrium distribution function

[TABLE]

We carried out some simulations on $[0,1]\times[0,1]$ with the lattice size $100\times 100$. As seen from Fig. 7 and 8, the numerical solutions are close to the analytical solutions at different values of $\alpha$, and the gradient term $\nabla\phi$ increases very fast with the decrease of $\alpha$. To see the difference between analytical and numerical solutions, we also measured the global relative errors and present them in Table 3 with CFL condition number equaling to 0.1. To test the convergence rate of DUGKS for this problem, We plotted the global relative errors at different lattice size in Fig. 9 where $dt=1.0\times 10^{-6}$, $\delta=1.2$ and $\alpha=0.01$. From this figure, it is also found that the DUGKS model for the NHCE is of second-order accuracy in space.

For this problem, some simulations were also performed with the non-uniform mesh. The non-uniform mesh is generated by the following transformation,

[TABLE]

[TABLE]

where $k=1.5$, which is used to control the distribution of non-uniform mesh. The grid point in ($\xi,\eta$) plane are defined by $\xi_{i}=i/Nx$ and $\eta_{j}=j/Ny$ for $i=0,1,...,N_{x}$ and $j=0,1,...,N_{y}$. The distributions of the uniform and the non-uniform meshes used in our simulations are shown in Fig. 10. With the same parameters shown in Table 3, we carried out some simulations with the non-uniform mesh, and the results are listed in Table 4.

As seen from Table 4, the $GRE$s with non-uniform meshes are smaller than those of mesha, and the difference becomes more obvious with the decrease of $\alpha$. This illustrates that the appropriate non-uniform mesh can improve the accuracy of the present DUGKS.

Example 4.4 The Gaussian hill problem is described by the following anisotropic convection diffusion equation [21]

[TABLE]

where $\mathbf{u}=(u_{x},u_{y})^{T}$ is a constant velocity, $\mathbf{K}$ is the constant diffusion tensor, and can be defined as

[TABLE]

The analytical solution to this Gaussian hill problem can be expressed as

[TABLE]

where $\mathbf{x}=(x,y)^{T}$, $\sigma=\sigma_{0}^{2}\mathbf{I}+2\mathbf{K}t$, $\sigma^{-1}$ is the inverse matrix of $\sigma$, $\det(\sigma)$ is the determinant of $\sigma$.

To study the Gaussian hill problem, we first write Eq. (52) in an isotropic form,

[TABLE]

where $\mathbf{B}=\phi\mathbf{u}$ and the tensor $\mathbf{D}$ is given by $\mathbf{D}=\mathbf{K}\phi/\kappa$ with $\kappa$ being a positive constant. The physical domain of the problem $[-1,1]\times[-1,1]$ and the periodic boundary conditions are applied to all directions. In our simulations, $\sigma=0.01$, $\mathbf{u}=(0.01,0.01)^{T}$, $\phi_{0}=2\pi\sigma_{0}^{2}$, CFL = 0.5 and the lattice size is $400\times 400$. To test the capacity of the present DUGKS for the Gaussian hill problem, the following three types of diffusion tensor are considered,

[TABLE]

which are corresponding to the isotropic, diagonally anisotropic and fully anisotropic diffusion problems.

We conducted several simulations and presented the numerical solutions at time $t=10$ in Figs. 11, 12 and 13 where $\kappa=0.001$ and $c=1.0$. As shown in these figures, the numerical solutions are consistent with the analytical solutions. In addition, to see the deviation between the numerical and analytical solutions, the $GRE$s of isotropic, diagonally anisotropic and fully anisotropic diffusion problems are also calculated, and they are $1.0829\times 10^{-3}$, $7.8158\times 10^{-4}$, and $1.7746\times 10^{-3}$, which also illustrate that the present DUGKS is accurate for Gaussian hill problem. Fig. 14 shows the accuracy of the present DUGKS for this problem, the lattice size is varied from $200\times 200$ to $500\times 500$ with time step $\Delta t=1.0\times 10^{-4}$. From this figure, we can find that the present DUGKS has a second-order convergence rate in space.

In the early work of Chai et al. [21], we have known that if $\kappa=10^{-4}$ and $\mathbf{u}=(0.1,0.1)^{\intercal}$, the SLBM could not give a stability solution. To test the stability of the present DUGKS, we performed some simulations with $\kappa=10^{-4}$ and $\mathbf{u}=(0.1,0.1)^{\intercal}$. As seen from Fig. 15, the DUGKS can give a stable numerical and accurate solution.

Example 4.5 We continue to consider the following convection-diffusion equation with nonlinear convection and diffusion terms as

[TABLE]

where $m$ and $\alpha$ are two constants, $\mathbf{D}$ is the tensor function of $\phi$, and is given by

[TABLE]

$F$ is the source term, and is defined as

[TABLE]

where $n_{x}$, $n_{y}$ and $A$ are the constants. Under the proper initial and periodic conditions, the analytical solution of this problem can be obtained,

[TABLE]

where $\kappa$ is a constant. For this problem, the functions $\mathbf{B}$ and $\mathbf{C}$ are given by

[TABLE]

We performed the simulations on [0,1] $\times$ [0,1] with the uniform grid 400 $\times$ 400, and the physical parameters are set as $\kappa=1.1$, $A=1.0$, $m=2.0$, $n_{x}=2.0$, $n_{y}=3.0$ and $c=1.0$. Besides, the CFL condition number is equal to 0.5. As seen from Figs. 16 and 17, the numerical solutions at $t=3.0$ and different Péclet numbers ($Pe=Lu_{x}/\alpha$) are in good agreement with analytical solutions, and the $GRE$s are about $3.348\times 10^{-5}$ for $Pe=100$ and $3.039\times 10^{-5}$ for $Pe=1000$. In addition, we also find that the values of $GRE$ are much smaller than $2.865\times 10^{-3}$ and $7.162\times 10^{-4}$ in Ref. [21].

Then the convergence rate of the present DUGKS for this problem is also considered, and the lattice size is varied from 100 $\times$ 100 to 500 $\times$ 500 with a fixed time step fixed $\Delta t=1.0\times 10^{-5}$. As shown in Fig. 18, the present DUGKS also has a second-order convergence rate for this nonlinear convection-diffusion equation.

Example 4.6 We now considered the nonlinear Fokker-Planck equation [33]

[TABLE]

with the initial condition

[TABLE]

and analytical solution

[TABLE]

where $\left\langle{x(t)}\right\rangle=\exp(t+\frac{t^{2}}{2})$, $\eta(t)=1-\exp(-t^{2})$, $\int_{-\infty}^{+\infty}\delta(x-x_{0})=1$ and

[TABLE]

We noted that in the above tests, the convection term $\mathbf{B}$ is only a function of $\phi$. However, in this problem, the convection term $\mathbf{B}$ is the function of $\phi$, $x$ and $t$, thus we have to define the auxiliary moment $\mathbf{C}=0$. For this example, the initial condition of $\phi(x,0)$ is taken as

[TABLE]

where $\eta$ is a small constant and $\Delta x$ is the lattice spacing.

In our simulations, the physical domain is fixed on $[-2,8]$, the uniform grid 400 $\times$ 400 is adopted, $\text{CFL}=0.5$. We presented the results at different time in Fig. 19. From this figure, we can see that the numerical solutions agree well with the analytical solutions.

To test the convergence rate of the DUGKS for this problem, some simulations were carried out at different lattice size ($\Delta x=1/4\sim 1/32$), and the time step is fixed at $\Delta t=1.0\times 10^{-5}$. As shown in Fig. 20, the present DUGKS indeed has a second-order convergence rate in space. Besides, theoretically, the DUGKS should also have a second-order convergence rate in time, to confirm this statement, we also carried out some simulations with a fixed $\Delta x=1/100$, the time step is varied from $8.0\times 10^{-5}$ to $1.0\times 10^{-5}$. As seen from Fig. 21, the present DUGKS does have a second-order convergence rate in time.

5 Conclusion

In this work, the discrete unified gas kinetic scheme is developed to solve general nonlinear convection-diffusion equation. Through Chapman-Enskog analysis, the NCDE can be recovered exactly from the present DUGKS. Through a lot of numerical simulations, we find that the numerical solutions are in good agreement with analytical solutions, and the present DUGKS has a second-order convergence rate in both space and time. In Example 4.1, a comparison was made between DUGKS, FV-LBM and MRT-LBM, the results show that the present DUGKS is more accurate than FV-LBM, and has almost same accuracy with the MRT-LBM. In Example 4.2 and Example 4.3, one can see that the present DUGKS is efficient and can be implemented on the non-uniform meshes. In Example 4.4, the results show that the DUGKS model is more stable than SLBM. In Example 4.5, a more nonlinear equation is considered to test our model, and finally in Example 4.6, we tested the present DUGKS, and found that the DUGKS also has a second-order convergence rate in time.

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)