A multi-prediction implicit scheme for steady state solutions of gas flow in all flow regimes

TL;DR

This paper introduces a multiscale implicit scheme with multiple macroscopic predictions for steady gas flow solutions across all regimes, significantly accelerating convergence especially in continuum flows.

Contribution

The paper develops a novel multiscale implicit method with multiple macroscopic predictions within the DVM framework, enhancing convergence speed in continuum flow regimes.

Findings

Method is one order of magnitude faster than previous schemes in continuum flows.

Effective across all flow regimes with improved convergence.

Utilizes a multiscale flux similar to DUGKS.

Abstract

An implicit multiscale method with multiple macroscopic prediction for steady state solutions of gas flow in all flow regimes is presented. The method is based on the finite volume discrete velocity method (DVM) framework. At the cell interface a multiscale flux with a construction similar to discrete unified gas-kinetic scheme (DUGKS) is adopted. The idea of the macroscopic variable prediction is further developed and a multiple prediction structure is formed. A prediction scheme is constructed to give a predicted macroscopic variable based on the macroscopic residual, and the convergence is accelerated greatly in the continuum flow regime. Test cases show the present method is one order of magnitude faster than the previous implicit multiscale scheme in the continuum flow regime.

Figures

Tables

Table 1. Table 1: Comparison of the efficiency between the implicit method in reference [ 17 ] and the present method for cavity flow simulations.

Case	Velocity space	Implicit method [17]		Present		Speedup
		Steps	Time (s)	Steps	Time (s)
Re=1000	1192	865	1102	23	74.4	14.8
Kn=0.075	1192	117	150	28	93.3	1.6
Kn=10	1192	173	225	33	114.3	2.0

Equations

$$\frac{\partial f}{\partial t}+\vec{u}\cdot \frac{\partial f}{\partial \vec{x}}=\frac{g-f}{\tau },$$

$$g=\rho {\left(\frac{\lambda }{\pi }\right)}^{\frac{3}{2}}{e}^{-\lambda {\vec{c}}^2},$$

$$g^{\ast }=\rho {\left(\frac{\lambda }{\pi }\right)}^{\frac{3}{2}}{e}^{-\lambda {\vec{c}}^2}\left[1+\frac{4(1-\mathrm{Pr}){\lambda }^2\vec{q}\cdot \vec{c}}{5\rho }(2\lambda {\vec{c}}^2-5)\right],$$

$$\vec{W}=\int \vec{\psi }fd\Xi ,$$

$$\mathbf{P}=\int \vec{c}\vec{c}fd\Xi ,$$

$$\vec{q}=\int \frac{1}{2}\vec{c}{\vec{c}}^{2}fd\Xi .$$

$$\int \vec{\psi }(g-f)d\Xi =\vec{0}.$$

$$\underset{\partial \Omega }{\int }\vec{u}\cdot \vec{n}fdA=\underset{\Omega }{\int }\frac{g-f}{\tau }dV,$$

$$\underset{\partial \Omega }{\int }\vec{F}dA=\vec{0},$$

$$\vec{F}=\int \vec{u}\cdot \vec{n}\vec{\psi }fd\Xi .$$

$$\sum _{j\in N\left(i\right)}{A}_{ij}{\vec{u}}_k\cdot {\vec{n}}_{ij}{f}_{ij,k}=V_i\frac{g_{i,k}-f_{i,k}}{{\tau }_i},$$

$$\sum _{j\in N\left(i\right)}{A}_{ij}{\vec{u}}_k\cdot {\vec{n}}_{ij}{f}_{ij,k}^{n+1}=V_i\frac{{\tilde{g}}_{i,k}^{n+1}-f_{i,k}^{n+1}}{{\tilde{\tau }}_i^{n+1}},$$

$$f_{ij,k}^{n+1}=\frac{{\tilde{\tau }}_{ij}^{n+1}}{{\tilde{\tau }}_{ij}^{n+1}+h_{ij}}f\left({\vec{x}}_{ij}-{\vec{u}}_k{h}_{ij},0,{\vec{u}}_k\right)+\frac{h_{ij}}{{\tilde{\tau }}_{ij}^{n+1}+h_{ij}}\tilde{g}\left({\vec{x}}_{ij},0,{\vec{u}}_k\right),$$

$$f({\vec{x}_{ij}}-{\vec{u}_{k}}{h_{ij}},0,{\vec{u}_{k}})=\left\{{\begin{array}[]{*{20}{l}}{f_{i,k}^{n+1}+({{\vec{x}}_{ij}}-{{\vec{x}}_{i}}-{{\vec{u}}_{k}}{h_{ij}})\nabla f_{i,k}^{n+1}{\mkern 1.0mu}{\kern 1.0pt}{\mkern 1.0mu}{\mkern 1.0mu}{\kern 1.0pt},\;{\kern 1.0pt}\;{\kern 1.0pt}\;{\kern 1.0pt}{{\vec{u}}_{k}}\cdot{{\vec{n}}_{ij}}\geq 0,}\\ {f_{j,k}^{n+1}+({{\vec{x}}_{ij}}-{{\vec{x}}_{j}}-{{\vec{u}}_{k}}{h_{ij}})\nabla f_{j,k}^{n+1}{\mkern 1.0mu}{\kern 1.0pt}{\mkern 1.0mu}{\kern 1.0pt},\;{\kern 1.0pt}\;{\kern 1.0pt}\;{\kern 1.0pt}{{\vec{u}}_{k}}\cdot{{\vec{n}}_{ij}}<0.}\end{array}}\right.$$

$${\widetilde{\stackrel{⃗}{W}}}_{ij}^{n+1}={\int }_{\vec{u}\cdot {\vec{n}}_{ij}\ge 0}\vec{\psi }{\tilde{g}}_{ij}^{\mathrm{l},\mathrm{n}+1}d\Xi +{\int }_{\vec{u}\cdot {\vec{n}}_{ij}<0}\vec{\psi }{\tilde{g}}_{ij}^{\mathrm{r},\mathrm{n}+1}d\Xi ,$$

$${\tilde{\tau }}_{ij}^{n+1}=\frac{\mu ({\widetilde{\stackrel{⃗}{W}}}_{ij}^{n+1})}{p({\widetilde{\stackrel{⃗}{W}}}_{ij}^{n+1})}+\frac{\left|p_{ij}^{\mathrm{l},n+1}-p_{ij}^{\mathrm{r},n+1}\right|}{\left|p_{ij}^{\mathrm{l},n+1}+p_{ij}^{\mathrm{r},n+1}\right|}{h}_{ij},$$

$$h_{ij}=\min (h_i,h_j).$$

$$h_i=\frac{V_i}{\underset{k}{\max }\left(\sum _{j\in N(i)}\left({\vec{u}}_k\cdot {\vec{n}}_{ij}{A}_{ij}\mathrm{H}[{\vec{u}}_k\cdot {\vec{n}}_{ij}]\right)\right)}\mathrm{CFL},$$

$${\rm{H}}[x]=\left\{\begin{array}[]{l}0,\quad x<0,\\ 1,\quad x\geq 0.\end{array}\right.$$

$$r_{i,k}^{n+1,(l)}=\frac{{\tilde{g}}_{i,k}^{n+1}-f_{i,k}^{n+1,(l)}}{{\tilde{\tau }}_i^{n+1}}-\frac{1}{V_i}\sum _{j\in N\left(i\right)}{A}_{ij}{\vec{u}}_k\cdot {\vec{n}}_{ij}{f}_{ij,k}^{n+1,(l)}.$$

$$r_{i,k}^{n+1,(l)}+\Delta r_{i,k}^{n+1,(l+1)}=\frac{1}{\Delta {\xi }_{i,k}^{n+1,(l+1)}}\Delta f_{i,k}^{n+1,(l+1)},$$

$$\left(\frac{1}{\Delta {\xi }_{i,k}^{n+1,(l+1)}}+\frac{1}{{\tilde{\tau }}_i^{n+1}}\right)\Delta f_{i,k}^{n+1,(l+1)}=r_{i,k}^{n+1,(l)}-\frac{1}{V_i}\sum _{j\in N(i)}{A}_{ij}{\vec{u}}_k\cdot {\vec{n}}_{ij}\Delta f_{ij,k}^{n+1,(l+1)}.$$

$$\Delta f_{ij,k}^{n+1,(l+1)}=\left\{\begin{array}[]{l}\frac{{\tilde{\tau}_{ij}^{n+1}}}{{\tilde{\tau}_{ij}^{n+1}+{h_{ij}}}}\Delta f_{i,k}^{n+1,(l+1)},\quad{{\vec{u}}_{k}}\cdot{{\vec{n}}_{ij}}\geq 0\\ \frac{{\tilde{\tau}_{ij}^{n+1}}}{{\tilde{\tau}_{ij}^{n+1}+{h_{ij}}}}\Delta f_{j,k}^{n+1,(l+1)},\quad{{\vec{u}}_{k}}\cdot{{\vec{n}}_{ij}}<0\end{array}\right.,$$

$$\Delta f_{ij,k}^{n+1,(1)}=\left\{{\begin{array}[]{*{20}{l}}{\Delta f_{i,k}^{n+1,(1)},\quad{{\vec{u}}_{k}}\cdot{{\vec{n}}_{ij}}\geq 0}\\ {\Delta f_{j,k}^{n+1,(1)},\quad{{\vec{u}}_{k}}\cdot{{\vec{n}}_{ij}}<0}\end{array}}\right..$$

$$\left(\frac{1}{\Delta {\xi }_{i,k}^{n+1,(l+1)}}+\frac{1}{{\tilde{\tau }}_i^{n+1}}+\frac{1}{V_i}\sum _{j\in N_k^{+}(i)}\frac{{\tilde{\tau }}_{ij}^{n+1}}{{\tilde{\tau }}_{ij}^{n+1}+h_{ij}}{A}_{ij}{\vec{u}}_k\cdot {\vec{n}}_{ij}\right)\Delta f_{i,k}^{n+1,(l+1)}$$

$$=$$

$${\vec{F}}_{ij}^{n+1}=\sum _k{\vec{\psi }}_k{\vec{u}}_k\cdot {\vec{n}}_{ij}{f}_{ij,k}^{n+1}\Delta {\Xi }_k,$$

$${\vec{R}}_i^{n+1}=-\frac{1}{V_i}\sum _{j\in N(i)}{A}_{ij}{\vec{F}}_{ij}^{n+1}.$$

$${\vec{W}}_i^{n+1}=\sum _k{\vec{\psi }}_k{f}_{i,k}^{n+1}\Delta {\Xi }_k+{\widetilde{\stackrel{⃗}{W}}}_i^{n+1}-\sum _k{\vec{\psi }}_k{\tilde{g}}_{i,k}^{n+1}\Delta {\Xi }_k,$$

$$\frac{1}{\Delta t_i^{n+1}}\left({\widetilde{\stackrel{⃗}{W}}}_i^{n+1}-{\vec{W}}_i^n\right)={\vec{R}}_i^n+\Delta {\widetilde{\stackrel{⃗}{R}}}_i^{n+1}.$$

$$\Delta {\widetilde{\stackrel{⃗}{R}}}_i^{n+1}=-\frac{1}{V_i}\sum _{j\in N(i)}{A}_{ij}{\widetilde{\stackrel{⃗}{\mathcal{F}}}}_{ij}^{n+1}+\frac{1}{V_i}\sum _{j\in N(i)}{A}_{ij}{\vec{\mathcal{F}}}_{ij}^n,$$

Full text

Keywords: implicit scheme, rarefied flow, kinetic scheme, multiscale scheme

1 Introduction

Rarefied gas flow simulation is always a research hotspot of computational fluid dynamics (CFD). Recent years, multiscale gas-kinetic methods based on the discrete velocity method (DVM, [1, 2, 3, 4, 5]) framework for nonequilibrium rarefied flow simulation have been developed, like the unified gas-kinetic scheme (UGKS) [6] by Xu and Huang, the discrete unified gas-kinetic scheme (DUGKS) [7, 8] by Guo et al. These multiscale methods overcome the time step and cell size restrictions of the original DVM method which requires time step and cell size of the order of mean collision time and mean free path, and thus have attracted more and more researchers’ attention. It is worth pointing out that although UGKS and DUGKS can adopt time step and cell size comparable to the traditional macroscopic Navier-Stokes (NS) method, they still involve large amount of computation due to the curse of dimensionality. Hence, many researches on the acceleration of these multiscale methods have been carried out, including Mao et al.’s implicit UGKS [9], Zhu et al.’s prediction based implicit UGKS [10, 11], Zhu et al.’s implicit multigrid UGKS algorithm [12], Yang et al.’s memory saving implicit multiscale scheme [13], Pan et al.’s implicit DUGKS [14], etc. Following these previous works, it is quite valuable to further develop the fast algorithm for the multiscale method.

In this paper, a multiple prediction implicit multiscale method for steady state calculation of gas flow in all flow regimes is proposed. The idea of macroscopic prediction presented by Zhu et al. [10] is further developed. A prediction solver is used to predict the macroscopic variable based on the macroscopic residual, and a multiple prediction procedure is constructed. The prediction solver is designed to ensure the accuracy of the predicted macroscopic variable in the continuum flow regime and the stability of the numerical system in all flow regimes, which makes the method very efficient in the continuum flow regime and stable in all flow regimes. Our test cases show that the present method is one order of magnitude faster than the previous implicit multiscale method in the continuum flow regime.

2 Numerical method

In this paper, the monatomic gas is considered and the governing equation is BGK-type equation [15],

[TABLE]

where $f$ is the gas particle velocity distribution function, $\vec{u}$ is the particle velocity, $\tau$ is the relaxation time calculated as $\tau=\mu/p$ ($\mu$ and $p$ are the viscosity and pressure). $g$ is the equilibrium state which has a form of Maxwellian distribution,

[TABLE]

or if the Shakhov model [16] is used

[TABLE]

where $\vec{c}$ is the peculiar velocity $\vec{c}=\vec{u}-\vec{U}$ and $\vec{U}$ is the macroscopic gas velocity, $\vec{q}$ is the heat flux, $\lambda$ is a variable related to the temperature $T$ by $\lambda=1/(2RT)$. Pr is the Prandtl number and has a value of $2/3$ for monatomic gas. $f$ is related to the macroscopic variables by

[TABLE]

where $\vec{W}=(\rho,\rho\vec{U},\rho E)^{T}$ is the vector of the macroscopic conservative variables, $\vec{\psi}$ is the vector of moments $\vec{\psi}={\left({1,\vec{u},\frac{1}{2}{{\vec{u}}^{2}}}\right)^{T}}$, $d\Xi=du_{x}du_{y}du_{z}$ is the velocity space element. The stress tensor $\boldsymbol{P}$ and the heat flux $\vec{q}$ can also be calculated by $f$ as

[TABLE]

[TABLE]

Moreover, $f$ and $g$ obey the conservation law,

[TABLE]

Adopting the integral form, the steady state of the governing equation Eq. 1 is

[TABLE]

where $\Omega$ is the control volume, $dV$ is the volume element, $dA$ is the surface area element and $\vec{n}$ is the outward normal unit vector. Take the moment of Eq. 8 for $\vec{\psi}={\left({1,\vec{u},\frac{1}{2}{{\vec{u}}^{2}}}\right)^{T}}$, the corresponding macroscopic governing equation can be written as

[TABLE]

where the flux $\vec{F}$ has the relation with the distribution function $f$ by

[TABLE]

This paper is about the numerical method of determining the steady state defined by Eq. 8. It is time-consuming to directly solve Eq. 8 through a microscopic scheme involving discretization in both physical space and velocity space. The main idea of the present method is summarized as that, using the accurate but expensive scheme to calculate the residual of the system deviating from the steady state, then utilizing this residual and using the less accurate but efficient scheme to do the evolution. More distinctly, an accurate multiscale microscopic scheme based on the DVM framework is used to handle the microscopic numerical system with Eq. 8, and a fast prediction scheme is used to do the evolution of the macroscopic variables. The prediction scheme can be some kind of macroscopic scheme based on macroscopic variables or even a scheme based on the DVM framework but with less velocity points. The schematic of the general algorithmic framework for the present method is shown in Fig. 1. The method consists of several loops in different layers. The outermost loop is denoted by $n$. One iteration of the $n$ loop includes a loop denoted by $m$ and a loop denoted by $l$, where the macroscopic variable $\vec{W}^{n}_{i}$ and the residual $\vec{R}^{n}_{i}$ are given as the input, the new $\vec{W}^{n+1}_{i}$ and $\vec{R}^{n+1}_{i}$ are the output. In the $m$ loop, the predicted macroscopic variable $\tilde{\vec{W}}^{n+1}_{i}$ is determined by the prediction scheme and by the numerical smoothing process. In the $l$ loop, the microscopic variable $f^{n+1}_{i,k}$ is calculated and the new $\vec{W}^{n+1}_{i}$ and $\vec{R}^{n+1}_{i}$ are obtained. The present method is a development of the prediction method of Zhu et al. [10], and has a structure similar to the multigrid method of Zhu et al. [12], therefore we call it as “multiple prediction method”. The method is detailed in following paragraphs.

2.1 Construction of the $l$ loop

In the $l$ loop, residuals of the numerical system are evaluated through the microscopic scheme and the microscopic variables (the discrete distribution function) are updated through an implicit method (the numerical smoothing process). The microscopic scheme is very important because it determines the final steady state of the whole numerical system and thus determines the nature of the present numerical method.

The microscopic scheme is based on Eq. 8. Discretizing the physical space by finite volume method and discretizing the velocity space into discrete velocity points, the microscopic governing equation Eq. 8 can be expressed as

[TABLE]

where the signs $i,k$ correspond to the discretizations in physical space and velocity space respectively. $j$ denotes the neighboring cell of cell $i$ and $N\left(i\right)$ is the set of all of the neighbors of $i$. Subscript $ij$ denotes the variable at the interface between cell $i$ and $j$. $A_{ij}$ is the interface area, ${\vec{n}_{ij}}$ is the outward normal unit vector of interface $ij$ relative to cell $i$, and $V_{i}$ is the volume of cell $i$. The $l$ loop aims to find the solution of Eq. 11 with the input predicted variable $\tilde{\vec{W}}^{n+1}_{i}$, therefore Eq. 11 can be written more exactly as

[TABLE]

where the symbol $\sim$ denotes the predicted variables at the $(n+1)$th step. $\tilde{g}_{i,k}^{n+1}$ and $\tilde{\tau}_{i}^{n+1}$ can be directly calculated from the input variable $\tilde{\vec{W}}^{n+1}_{i}$. The distribution function $f_{ij,k}^{n+1}$ at the interface $ij$ is very important to ensure the multiscale property of the scheme. In this paper, following the idea of DUGKS [7, 8], the construction of $f_{ij,k}^{n+1}$ in reference [17] is adopted, i.e.

[TABLE]

where

[TABLE]

In above equations, $\nabla f_{i,k}^{n+1}$ and $\nabla f_{j,k}^{n+1}$ can be obtained through the reconstruction of the distribution function data. $\tilde{g}\left({{{\vec{x}}_{ij}},0,{{\vec{u}}_{k}}}\right)$ and ${\tilde{\tau}_{ij}^{n+1}}$ are calculated by the same way as the method of GKS [18] and they can be both calculated from the predicted macroscopic variable $\tilde{\vec{W}}^{n+1}_{i}$. For $\tilde{g}({\vec{x}_{ij}},0,{\vec{u}_{k}})$, it is determined by the interface macroscopic variables $\tilde{\vec{W}}^{n+1}_{ij}$, which can be calculated as

[TABLE]

where the superscripts $\rm{l}$ and $\rm{r}$ denote variables at the left and right sides of the interface, ${\tilde{g}_{ij}^{\rm{l},n+1}}$ and ${\tilde{g}_{ij}^{\rm{r},n+1}}$ can be determined after the spacial reconstruction of $\tilde{\vec{W}}^{n+1}_{i}$. For ${\tilde{\tau}_{ij}^{n+1}}$, it is calculated as

[TABLE]

where the pressure ${p_{ij}^{{\rm{l}},n+1}}$, ${p_{ij}^{{\rm{r}},n+1}}$ at two sides of the interface can be obtained from the reconstruction and the second term on the right is for artificial viscosity. $h_{ij}$ in above equations is calculated from the physical local time step

[TABLE]

The physical local time step $h_{i}$ for the cell $i$ is determined by the local CFL condition as

[TABLE]

where ${\rm{H}}[x]$ is the Heaviside function defined as

[TABLE]

For more details about the construction of the interface distribution function $f_{ij,k}^{n+1}$ please refer to reference [17].

Eq. 12 is solved by iterations. The microscopic residual $r_{i,k}^{n+1,(l)}$ at the $l$th iteration can be defined as

[TABLE]

According to the previous descriptions, $r_{i,k}^{n+1,(l)}$ can be calculated from $f_{i,k}^{n+1,(l)}$ and $\tilde{\vec{W}}^{n+1}_{i}$ through the spatial data reconstruction. The increment equation to get the microscopic variable $f_{i,k}^{n+1,(l+1)}$ at the iteration $l+1$ is constructed by backward Euler method,

[TABLE]

where ${\Delta\xi_{i,k}^{n+1,(l+1)}}$ is the pseudo time step and ${\Delta\xi_{i,k}^{n+1,(l+1)}}$ is always set to be $\infty$ in the present study. Combined with the residual expression Eq. 20, Eq. 21 can be written as

[TABLE]

For the increment of the interface distribution function ${\Delta f_{ij,k}^{n+1,(l+1)}}$, it is simply handled by a modified upwind scheme

[TABLE]

where the coefficient ${\frac{{\tilde{\tau}_{ij}^{n+1}}}{{\tilde{\tau}_{ij}^{n+1}+{h_{ij}}}}}$ is the corresponding coefficient multiplied by $f\left({{{\vec{x}}_{ij}}-{{\vec{u}}_{k}}{h_{ij}},0,{{\vec{u}}_{k}}}\right)$ in Eq. 13. This coefficient is multiplied because during the whole $l$ loop the term $\tilde{g}\left({{{\vec{x}}_{ij}},0,{{\vec{u}}_{k}}}\right)$ in Eq. 13 is calculated by the predicted macroscopic variable $\tilde{\vec{W}}^{n+1}_{i}$ and therefore is an invariant, so the variation of the microscopic variable $f_{i,k}^{n+1,(l+1)}$ only influences the term $f\left({{{\vec{x}}_{ij}}-{{\vec{u}}_{k}}{h_{ij}},0,{{\vec{u}}_{k}}}\right)$, which is multiplied by the coefficient ${\frac{{\tilde{\tau}_{ij}^{n+1}}}{{\tilde{\tau}_{ij}^{n+1}+{h_{ij}}}}}$. It is worth noting that in an actual implementation of the method, the interface distribution function $f_{ij,k}^{n}$ at the $n$th step may be taken as the initial value $f_{ij,k}^{n+1,(0)}$ at $l=0$ for the step $n+1$ to reduce computation cost, in this situation the variation $\Delta f_{i,k}^{n+1,(1)}$ should also account for the variation of $\tilde{g}\left({{{\vec{x}}_{ij}},0,{{\vec{u}}_{k}}}\right)$, and the coefficient ${\frac{{\tilde{\tau}_{ij}^{n+1}}}{{\tilde{\tau}_{ij}^{n+1}+{h_{ij}}}}}$ shouldn’t be multiplied at the first iteration of the $l$ loop, i.e.

[TABLE]

In this situation, after the first iteration of the $l$ loop, the interface distribution function $f_{ij,k}^{n+1,(l>0)}$ will be calculated with the newly predicted $\tilde{\vec{W}}^{n+1}_{i}$ and Eq. 23 is used to handle ${\Delta f_{ij,k}^{n+1,(l+1)}}$ again. Without loss of generality, substituting Eq. 23 into Eq. 22 will yield

[TABLE]

where $N_{k}^{+}(i)$ is the set of $i$’s neighboring cells satisfying ${\vec{u}_{k}}\cdot{\vec{n}_{ij}}\geq 0$ while for $N_{k}^{-}(i)$ it satisfies ${\vec{u}_{k}}\cdot{\vec{n}_{ij}}<0$. For simplicity, Eq. 25 is solved by the Symmetric Gauss-Seidel (SGS) method, or also known as the Point Relaxation Symmetric Gauss-Seidel (PRSGS) method [19, 20]. In each time of the SGS iteration, a forward sweep from the first to the last cell and a backward sweep from the last to the first cell are implemented, during which the data of a cell is always updated by the latest data of its adjacent cells through Eq. 25. Such a SGS iteration procedure is totally matrix-free and easy to implement.

After several times of SGS iterations for solving Eq. 25, an evaluation of $f_{i,k}^{n+1,(l+1)}$ with a certain precision can be obtained. Then the residual $r_{i,k}^{n+1,(l+1)}$ at the $(l+1)$th iteration of the $l$ loop can be computed from $f_{i,k}^{n+1,(l+1)}$ and $\tilde{\vec{W}}^{n+1}_{i}$, and a new turn of the $l$ loop will be performed. After several iterations of the $l$ loop, an evaluation of $f_{i,k}^{n+1}$ with a certain precision can be obtained, and the interface distribution function $f_{ij,k}^{n+1}$ can be calculated by Eq. 13. Then the macroscopic numerical flux ${\vec{F}_{ij}^{n+1}}$ at the interface can be got by numerical integral in the discrete velocity space

[TABLE]

and the macroscopic residual $\vec{R}_{i}^{n+1}$ defined by the macroscopic governing equation Eq. 9 at the $(n+1)$th step can be calculated from the flux by

[TABLE]

Note that in the $l$ loop we solve the microscopic system Eq. 12 which is under the condition of the predicted variable $\tilde{\vec{W}}^{n+1}_{i}$, so $\vec{R}_{i}^{n+1}$ is not zero even if the microscopic system is solved sufficiently accurately. Finally, the macroscopic variable ${\vec{W}}^{n+1}_{i}$ is calculated by numerical integral as

[TABLE]

where the term $\tilde{\vec{W}}_{i}^{n+1}-\sum\limits_{k}{{{\vec{\psi}}_{k}}\tilde{g}_{i,k}^{n+1}\Delta{\Xi_{k}}}$ is the integral error compensation term to make the scheme conservative, more details about this term please refer to reference [17].

The iteration of the $l$ loop is similar to the numerical smoothing process in multigrid method [12]. The computation procedure of the $l$ loop is listed as follows:

Step 1.

Set the initial value $f_{i,k}^{n+1,(0)}=f_{i,k}^{n}$.

Step 2.

Calculate the interface distribution function $f_{ij,k}^{n+1,(l)}$ by Eq. 13 from $f_{i,k}^{n+1,(l)}$ and $\tilde{\vec{W}}_{i}^{n+1}$ through data spatial reconstruction. Calculate the microscopic residual $r_{i,k}^{n+1,(l)}$ by Eq. 20.

Step 3.

Make judgement: if the residual $r_{i,k}^{n+1,(l)}$ meets the convergence criterion, or if the iteration number of the $l$ loop meets the maximum limit, break out of the $l$ loop and go to Step 5.

Step 4.

Solve Eq. 25 by several times of SGS iterations, obtain $f_{i,k}^{n+1,(l+1)}$, and go to Step 2.

Step 5.

By Eq. 26, Eq. 27 and Eq. 28, do numerical integral in the velocity space to get ${\vec{W}}_{i}^{n+1}$ and $\vec{R}_{i}^{n+1}$ for the step $n+1$.

2.2 Construction of the $m$ loop

In the $m$ loop, based on the macroscopic variable ${\vec{W}}_{i}^{n}$ and the macroscopic residual $\vec{R}_{i}^{n}$ at the $n$th step, a reasonable estimation for the macroscopic variable $\tilde{\vec{W}}_{i}^{n+1}$ is obtained through a fast prediction scheme to accelerate convergence. Theoretically speaking, the prediction scheme can be either a macroscopic scheme based on macroscopic variables or a microscopic scheme based on the DVM framework but with less velocity points. In this paper, a macroscopic scheme is designed to do the prediction. The process of the $m$ loop has certain similarity to the coarse grid correction in the multigrid method [12].

2.2.1 Framework

The macroscopic residual has the form of Eq. 27. To reduce the residual, a prediction equation is constructed by backward Euler formula

[TABLE]

$\Delta{t_{i}^{n+1}}$ is the local prediction time step, the purpose of this time step is to constrain the marching time depth of the prediction process to make the scheme stable in the extreme case. The predicted residual increment $\Delta\tilde{\vec{R}}_{i}^{n+1}$ is calculated by

[TABLE]

where $\vec{\mathcal{F}}_{ij}^{n}$ and $\tilde{\vec{\mathcal{F}}}_{ij}^{n+1}$ are fluxes calculated by the prediction solver from $\vec{W}_{i}^{n}$ and the predicted $\tilde{\vec{W}}_{i}^{n+1}$ with data reconstruction. This prediction solver is well-designed to balance between accuracy and stability, and will be presented later in the next section.

The aim of the $m$ loop is to solve Eq. 29 and give an estimation for $\tilde{\vec{W}}_{i}^{n+1}$ with a certain precision. Like what we do in the $l$ loop, Eq. 29 is also solved by iterations. The residual $\vec{\mathcal{R}}_{i}^{n+1,(m)}$ at the $m$th iteration can be defined by Eq. 29 and expressed as

[TABLE]

and the corresponding increment equation for $\tilde{\vec{W}}_{i}^{n+1,(m+1)}$ is

[TABLE]

where ${\Delta\eta_{i}^{n+1,(m+1)}}$ is the pseudo time step. Considering Eq. 31, the increment of the residual $\Delta\vec{\mathcal{R}}_{i}^{n+1,(m+1)}$ can be expressed as

[TABLE]

the variation of the flux $\Delta\tilde{\vec{\mathcal{F}}}_{ij}^{n+1,(m+1)}$ is further approximated by

[TABLE]

where $\vec{\mathsf{F}}_{ij}$ has the form [21] of the well-known Roe’s flux function

[TABLE]

Here ${\vec{\mathbb{F}}_{ij}}(\vec{W})$ is the Euler flux

[TABLE]

and ${\mathfrak{r}_{ij}}$ is

[TABLE]

where $a_{ij}$ is the acoustic speed at the interface and $\Delta{x_{ij}}$ is the distance between cell center $i$ and $j$. Substitute Eq. 33, Eq. 34 and Eq. 35 into Eq. 32, approximate ${{\mathfrak{r}}_{ij}^{n+1,(m+1)}}$ by ${\mathfrak{r}_{ij}^{n+1,(m)}}$, and note that $\sum\limits_{j\in N(i)}{{A_{ij}}{\vec{\mathbb{F}}_{ij}}({{\vec{W}}_{i}})}=\vec{0}$ holds, we can get

[TABLE]

Eq. 38 is solved by several times’ SGS iterations. An estimation of $\tilde{\vec{W}}_{i}^{n+1,(m+1)}$ with a certain precision can be obtained from Eq. 38, then ${\mathfrak{r}}_{ij}^{n+1,(m+1)}$ and the residual $\vec{\mathcal{R}}_{i}^{n+1,(m+1)}$ at the $(m+1)$th iteration of the $m$ loop can be calculated. After several turns of the $m$ loop, the predicted macroscopic variable $\tilde{\vec{W}}_{i}^{n+1}$ with a certain precision can be determined.

In fact, utilizing ${\vec{W}}_{i}^{n}$ and $\tilde{\vec{W}}_{i}^{n+1}$, a prediction for the microscopic variable ${\tilde{f}}_{i,k}^{n+1}$ can also be obtained to accelerate the convergence of the microscopic numerical system (i.e. the $l$ loop). The increment of the distribution function $\Delta{\tilde{f}}_{i,k}^{n+1}$ can be calculated from the Chapman-Enskog expansions [22] based on macroscopic variables ${\vec{W}}_{i}^{n}$ and $\tilde{\vec{W}}_{i}^{n+1}$. This strategy will increase the complexity of the algorithm and thus is not adopted in the present method.

Likewise, as one can see, the process of the $m$ loop is similar to the numerical smoothing process in multigrid method [12]. The computation procedure of the $m$ loop is listed as follows:

Step 1.

Set the initial value $\tilde{\vec{W}}_{i}^{n+1,(0)}={\vec{W}}_{i}^{n}$.

Step 2.

Calculate the residual $\vec{\mathcal{R}}_{i}^{n+1,(m)}$ by Eq. 31 from $\vec{R}_{i}^{n}$, $\vec{W}_{i}^{n}$ and $\tilde{\vec{W}}_{i}^{n+1,(m)}$ (data reconstruction is implemented).

Step 3.

Make judgement: if the residual $\vec{\mathcal{R}}_{i}^{n+1,(m)}$ meets the convergence criterion, or if the iteration number of the $m$ loop meets the maximum limit, break out of the $m$ loop and the predicted macroscopic variable $\tilde{\vec{W}}_{i}^{n+1}$ is determined.

Step 4.

Solve Eq. 38 by several times of SGS iterations, obtain $\tilde{\vec{W}}_{i}^{n+1,(m+1)}$, and go to Step 2.

2.2.2 Prediction solver

The prediction solver used to calculate the fluxes $\vec{\mathcal{F}}_{ij}^{n}$ and $\tilde{\vec{\mathcal{F}}}_{ij}^{n+1}$ in Eq. 30 requires careful design. For the continuum flow, the prediction solver should be as accurate as a traditional NS solver. For the rarefied flow, it’s unrealistic for a fast solver based on macroscopic variables to provide a very precise flux, but the solver should be stable so that the present method can be applied to all flow regimes. Thus, there are two principles for the prediction solver: accurate in the continuum flow regime, stable in all flow regimes.

We start constructing the solver from the view of gas kinetic theory. Based on the famous Chapman-Enskog expansion [22], the distribution function $f$ obtained from the BGK equation Eq. 1 to the first order of $\tau$ is

[TABLE]

Suppose there is an interface in $x$ direction. If the interface distribution function has the form of Eq. 39, take moments of $u_{x}\vec{\psi}$ to Eq. 39 and ignore second (and higher) order terms of $\tau$, we can get the NS flux [18, 23], where the term $g$ corresponds to the Euler flux and terms with $\tau$ (i.e. terms except $g$) correspond to viscous terms in the NS flux. Flux directly calculated from Eq. 39 will lead to divergence in many cases, and we introduce some modifications below.

The Euler flux often causes stability issue. Inspiring by gas-kinetic scheme (GKS) or also known as BGK-NS scheme [18], we replace it by a weighting of Euler flux and the flux of kinetic flux vector splitting (KFVS) [24]. That is, we replace the term $g$ in Eq. 39 and the interface distribution function is expressed as

[TABLE]

where $g^{\rm{lr}}$ is

[TABLE]

which is determined by the reconstructed macroscopic variables on the two side of the interface. The interface macroscopic variable $\vec{W}$ is calculated as

[TABLE]

and $g$ is obtained from $\vec{W}$. The weight factors $\tau^{\prime}/(\tau^{\prime}+h)$ and $h/(\tau^{\prime}+h)$ share the same forms as those in Eq. 13 (for how these weight factors are constructed please refer to [17]), and $\tau^{\prime}$ is calculated by

[TABLE]

where ${\tau_{{\rm{artificial}}}}$ is for artificial viscosity. $h$ is the local CFL time step and is equal to $h_{ij}$ in Eq. 13. Eq. 40 has a form similar to the interface distribution function of GKS [18], except that the viscous term is not upwind split and the weight factor is constructed following the thought of DUGKS [7, 8]. Because the KFVS scheme is very robust, the flux obtained from Eq. 40 makes the numerical system more stable than directly using Eq. 39. In the continuum flow regime, $h\gg\tau$, if the flow is continuous the term ${\tau_{{\rm{artificial}}}}$ for artificial viscosity will be negligible and Eq. 40 will recover the NS flux, while if the flow is discontinuous the term ${\tau_{{\rm{artificial}}}}$ will be activated and Eq. 40 will work as a stable KFVS solver. In the rarefied flow simulation, $\tau>h$ and the inviscid part of Eq. 40 generally provides a KFVS flux, which can increase the stability of the scheme.

The flux obtained from Eq. 40 works well in the continuum flow regime. However, in the case of large Kn number, the numerical system based on Eq. 40 is very stiff due to the large NS-type linear viscous term and the scheme is easy to blow up. Therefore, we multiply the viscous term by a limiting factor ${\mathfrak{q}}(\kappa)$ and Eq. 40 is transformed into

[TABLE]

Here we emphasize that the limiting factor ${\mathfrak{q}}(\kappa)$ aims not to accurately calculate the flux, but to increase the stability in the case of large Kn number. One can view it as an empirical parameter. The limiting factor ${\mathfrak{q}}(\kappa)$ is constructed considering the form of nonlinear coupled constitutive relations (NCCR) [25, 26], and is expressed as

[TABLE]

which has $\mathop{\lim}\limits_{\kappa\to 0}{\mathfrak{q}}(\kappa)=1$ and $\mathop{\lim}\limits_{\kappa\to+\infty}{\mathfrak{q}}(\kappa)=0$. $\kappa$ is related to the viscous term and calculated as

[TABLE]

where $-k\nabla T$ and $2\mu{S_{ij}}$ correspond to the heat flux and stress in NS equation, $C_{p}$ is the specific heat at constant pressure. $\beta$ is a molecular model coefficient [26] involved in the variable soft sphere (VSS) model [27, 28] and is calculated as

[TABLE]

where the molecular scattering factor $\alpha$ and the heat index $\omega$ depend on the type of gas molecule. The limiting factor ${\mathfrak{q}}(\kappa)$ is constructed to weaken the viscous term in large Kn number case to make the scheme stable. It can be seen from Eq. 45 and Eq. 46 that, when the stress and heat flux are small, ${\mathfrak{q}}(\kappa)$ is approaching to $1$ and we can get the NS viscous term in Eq. 44, when the stress and heat flux are large, ${\mathfrak{q}}(\kappa)$ is approaching to [math] and the viscous term is weakened. Here we further reveal the mechanism of ${\mathfrak{q}}(\kappa)$ through a simple one-dimensional case where there is no stress but only heat flux, i.e. $k\partial T/\partial x\neq 0$ and $\partial{U}/\partial x=0$. In this case $\kappa$ is

[TABLE]

and the heat flux from the viscous term of Eq. 44 is

[TABLE]

If the magnitude of the NS heat flux $\left|{q_{{\rm{NS}}}}\right|$ approaches [math], ${\mathfrak{q}}(\kappa)$ will approach $1$ and $q$ approaching $q_{{\rm{NS}}}$ holds true for Eq. 49. If $\left|{q_{{\rm{NS}}}}\right|$ approaches $+\infty$, in this case the heat flux $q$ from Eq. 49 goes to

[TABLE]

where $M$ is

[TABLE]

On the other hand, in the NCCR relation [26], for one-dimensional case, if $\partial{U}/\partial x=0$, the heat flux is calculated as

[TABLE]

where ${\kappa_{{\rm{NCCR}}}}$ is

[TABLE]

If the magnitude of $\left|{q_{{\rm{NCCR}}}}\right|$ approaches [math], similarly $q_{{\rm{NCCR}}}$ approaching $q_{{\rm{NS}}}$ holds true, i.e. the NS heat flux is recovered. If the magnitude of $\left|{q_{{\rm{NCCR}}}}\right|$ approaches $+\infty$, in this limiting case the magnitude of the heat flux can be deduced from Eq. 52 and Eq. 53 as

[TABLE]

where $M$ has the same definition as Eq. 51. Comparing Eq. 50 and Eq. 54, one can find that $q$ and $q_{{\rm{NCCR}}}$ are identical in the limiting case. The above derivation implies that $q$ and $q_{{\rm{NCCR}}}$ are very similar when their magnitudes are very small or very large. Of course, instead of the above special case, for more general multidimensional case, $\vec{q}$ from the viscous term of Eq. 44 and $\vec{q}_{{\rm{NCCR}}}$ based on the NCCR relation [26] are not exactly same when their magnitudes approach $+\infty$, but they are generally of the same order of magnitude when they are large. All in all, the viscous term of Eq. 44 recovers the NS viscous term when the stress and heat flux are small, and this viscous term will be reduced to the same order of magnitude as the NCCR viscous term when the stress and heat flux are large. Thus, in small Kn number case the flux obtained from Eq. 44 is accurate as the NS flux, while in large Kn number case the viscous term of Eq. 44 is suppressed to make the numerical system more stable.

Finally, take moments of $u_{x}\vec{\psi}$ to Eq. 44 (ignore second and higher order terms of $\tau$, i.e. $\int{\vec{\psi}\left({\frac{{\partial g}}{{\partial t}}+\vec{u}\cdot\frac{{\partial g}}{{\partial\vec{x}}}}\right)d\Xi}=\vec{0}$ is used to transform time derivatives into spatial derivatives), the prediction flux is

[TABLE]

For the calculation about the moments of the Maxwellian distribution function, one can refer to reference [18] for some instruction.

The present prediction solver based on Eq. 55 is efficient compared to the solver of GKS [18]. It is accurate as an NS solver in the continuum flow regime and has enhanced stability in large Kn number case. It is not accurate for rarefied flow calculation, but as mentioned before, it’s unrealistic for a solver based on macroscopic variables to provide a very precise flux in large Kn number case. As a prediction solver, stability is the most important. The accuracy of the final solution obtained from the present method only depends on the microscopic scheme described in Section 2.1.

3 Numerical results and discussions

More test cases will be added during the preparation of the final paper.

3.1 Lid-driven cavity flow

The test case of lid-driven cavity flow is performed to test the efficiency of the present method, and to test if the viscous effect can be correctly simulated by the present method. Three cases Re=1000 and Kn=0.075, 10 are considered, involving gas flows from continuum regime to free molecular regime. The Mach number, which is defined by the upper wall velocity $U_{\rm{wall}}$ and the acoustic velocity, is 0.16. The Shakhov model is used and the Prandtl number Pr=$2/3$. The hard sphere (HS) model is used, with heat index $\omega$=0.5 and molecular scattering factor $\alpha$=1. On the wall of the cavity, the diffuse reflection boundary condition with full thermal accommodation [23] is implemented. For the physical space discretization, as shown in Fig. 2, a nonuniform 61$\times$61 mesh with a mesh size $0.004L$ ($L$ is the width of the cavity) near the wall is used for the case Re=1000 while a uniform 61$\times$61 mesh is used for the cases Kn=0.075, 10. For the velocity space discretization, as shown in Fig. 3, a 1192 cells’ unstructured mesh is used, where the central area is refined to reduce the ray effect. For the iteration strategy, in each step $n$, 60 turns of $m$ loop and 3 turns of $l$ loop are performed, while 40 times and 6 times of SGS iterations are executed for each turn of the $m$ loop and the $l$ loop respectively. The prediction step $\Delta t_{i}^{n+1}$ in Eq. 29 is set as $+\infty$ in this set of test cases. The convergence criterion is that the global root-mean-square of the infinite norm about the macroscopic residual vector defined by Eq. 27 is less than $10^{-9}$. Computations are run on a single core of a computer with Intel(R) Xeon(R) CPU E5-2678 v3 @ 2.50GHz. The computational efficiency compared with the implicit multiscale method in reference [17] is shown in Tab. 1. It can be seen that in the continuum flow regime (case of Re=1000), the present method is one order of magnitude faster than the implicit method of reference [17]. Considering that the implicit method of reference [17] is two orders of magnitude faster than explicit UGKS (discussed in [17]) in the continuum flow regime, the present method should be thousands of times faster than explicit UGKS in the continuum flow regime. For the cases Kn=0.075, 10, the present method is only one to two times faster than the method of reference [17]. This efficiency is reasonable because in the continuum flow regime the prediction scheme (the $m$ loop) gives very accurate predicted macroscopic variables and the numerical system converges rapidly, while for the rarefied flow the prediction scheme fails to be so precise and the slight efficiency increase compared with the method of reference [17] is due to the improved iteration strategy (namely, the $l$ loop) of the present method. The results of the present method for this set of test cases are shown in Fig. 4, Fig. 5 and Fig. 6. The present results agree very well with the results of GKS and UGKS.