claimClaim remarkRemark hypothesisHypothesis

    Unsteady GSIS: fast convergence and asymptotic preservingLei Wu

    # General synthetic iterative scheme for unsteady rarefied gas flows
    ††thanks: Submitted to the editors DATE.

    Lei Wu
    Department of Mechanics and Aerospace Engineering, Southern University of Science and Technology, Shenzhen 518055, China ()
<EMAIL_ADDRESS>

    ###### Abstract

    In rarefied gas flows, the spatial grid size could vary by several orders of magnitude in a single flow configuration (e.g., inside the Knudsen layer it is at the order of mean free path of gas molecules, while in the bulk region it is at a much larger hydrodynamic scale). Therefore, efficient implicit numerical method is urgently needed for time-dependent problems.
    However, the integro-differential nature of gas kinetic equations poses a grand challenge, as the gain part of the collision operator is non-invertible. Hence an iterative solver is required in each time step, which usually takes a lot of iterations in the (near) continuum flow regime where the Knudsen number is small; worse still, the solution does not asymptotically preserve the fluid dynamic limit when the spatial cell size is not refined enough.
    Inspired by our general synthetic iteration scheme for steady-state solution of the Boltzmann equation, we propose two numerical schemes to push the multiscale simulation of unsteady rarefied gas flows to a new boundary, that is, the numerical solution not only converges within dozens of iterations in each time step, but also asymptotic preserves the Navier-Stokes-Fourier limit at coarse spatial grid, even when the time step is large (e.g., in the extreme slow decay of two-dimensional Taylor vortex, the time step is at the order of vortex decay time). The properties of fast convergence and asymptotic preserving of the proposed schemes are not only rigorously proven by the Fourier stability analysis, but also demonstrated by solid numerical examples.

    ###### keywords:

    gas kinetic equation, fast convergence, asymptotic preserving

    76P05, 65L04, 65M12

    ##
    1 Introduction

    Rarefied gas flows have attracted significant research interest in the past decades due their wide range of engineering applications. These flows are characterized by the Knudsen number Kn, which is defined as the ratio of mean free path $\lambda$ (or mean collision time $t_{c}$) of gas molecules to the characteristic flow length $L$ (or time/period $T$). Only when the Knudsen number is small can the rarefied gas dynamics be well described by macroscopic equations in the bulk flow region [15], such as the Euler, Navier-Stokes-Fourier (NSF), Burnett [3], and Grad 13 moments equations [7]; see the short review on the performance of macroscopic equations [22]. For general values of Knudsen number, however, the Boltzmann equation or simplified gas kinetic equations should be used.

    Since the velocity distribution function is defined in a six-dimensional phase space, the computational cost of memory and time in the solving of gas kinetic equations is huge. Thus, many numerical methods are proposed to solve the kinetic equations under a numerical scale larger than the kinetic one [12, 5, 6, 4, 10, 23, 9], that is, the spatial grid size $\Delta{x}\gg\lambda$, and/or the time step $\Delta{t}\gg\tau_{c}$. Some schemes asymptotically preserve the Euler limit, as they become a consistent discretization of the Euler equations when $\text{Kn}\rightarrow 0$ [5, 6].
    Nevertheless, from a practical point of view, the Euler equations cannot be applied for some flows, even when the Knudsen number is very small. For instances, in the Poiseuille flow [21] and Rayleigh-Brillouin scattering [19], the flow velocity and density amplitude scales as $1/\text{Kn}$. If the Euler equations are used, they become divergent, which is not physical.
    Therefore, some numerical schemes are designed to asymptotically preserve the NSF limit when $\Delta{t}\gg\tau_{c}$ [4, 10], under the assumption that the spatial derivatives are handled exactly. Recently, it is found that the NSF limit can be captured by the (discrete) unified gas-kinetic scheme, when both the time step and spatial cell size are much larger than the corresponding kinetic scales [23, 9, 8, 26], e.g., $\Delta{x}\sim\sqrt{\text{Kn}}L\gg\lambda$.
    It should be noted that, however, these large numerical scales are used in bulk flow region only, without resolving the Knudsen layer which occupies a spatial region within a few mean free path away from the solid walls.

    In reality, rarefied gas flows are intrinsically multiscale, say, in the two-dimensional thermal edge flow in the (near) continuum flow regime where the Knudsen number is small [19], the spatial grid size varies by several orders of magnitude: inside the Knudsen layer $\Delta{x}\sim\lambda\sim 0.001$, while in the bulk region it is at a much larger hydrodynamic scale: $\Delta{x}\sim{}L\sim 1$. Such a disparate grid cell size is necessary, as the fine grid in the Knudsen layer helps to capture the ghost effect that arises from the rarefaction effects inside the Knudsen layer and determines the velocity field in the whole computational domain [24], while the coarse grid help to save the computational cost; our numerical tests show that, if the Knudsen layer is under-resolved, the vortexes rotate in wrong directions. Therefore, for time-dependent problems, it is highly desired to use the implicit numerical method, otherwise the time step in explicit numeral methods will be restricted by the Courant–Friedrichs–Lewy condition, namely, $\Delta{t}\sim\tau_{c}$, which renders any practical numerical simulation impossible when the Knudsen number is small.

    However, a grand challenge arises for implicit numerical method, since the integro-differential nature of gas kinetic equations makes the gain part of the collision operator non-invertible. Therefore, an iterative solver is required in each time step. A simple way is to use the conventional iterative scheme (CIS), where the time derivative, streaming operator, and loss part of the collision operator are evaluated at the current iteration step, while the gain part of the collision operator is obtained from the previous iteration step. The CIS is efficient when the Knudsen number is large, but it takes a lot of iterations in the (near) continuum flow regime; worse still, the CIS does not asymptotically preserve the NSF limit, so a huge number of spatial cells should be used to capture the gas dynamics [21]. To improve the efficiency of implicit iteration, a fast convergence and asymptotic NSF preserving scheme is urgently needed.

    Thanks to the pioneering work of Larsen in neutron transport [11], the acceleration is possible if the kinetic equation and its moment equations are coupled: the kinetic equation provides high-order moments to moment equations (diffusion equations in the content of neutron transport), while the moment equations provide macroscopic quantities appearing in the gain part of the collision operator.
    The essential idea was extended to gas kinetic system, such as the high-order/low-order (HOLO) method [20] for the Bhatnagar-Gross-Krook (BGK) model [2] and the general synthetic iterative scheme (GSIS) for the Boltzmann equation [18]. In the Sod shock tube problem, it was found that in each time step HOLO finds the converged solution within dozens of iterations [20]. Albeit this promising property, the stability and asymptotic NSF preserving of HOLO have not been rigorously analyzed. On the other hand, while the fast convergence and asymptotic NSF preserving of GSIS has been proven [19, 17, 25], so far it has been limited to steady-state problems.

    In this paper we extend the GSIS to time-dependent problems, retaining its unique properties of fast convergence and asymptotic NSF preserving.
    In GSIS, the kinetic equation and macroscopic synthetic equations can be solved by different numerical methods with different order of accuracy. This brings tremendous numerical convenience because the kinetic equation, which requires discretization in the high-dimensional phase space, are usually time-consuming and hence should be handled by as simple algorithm as possible. The macroscopic synthetic equation, on the contrary, are well studied in computational fluid dynamics and can be handled by sophisticated high-order numerical methods.
    From the practical application point of view, we hope the developed kinetic schemes could use as large time step and spatial cell size as possible. That is, $\Delta{t}$ and $\Delta{x}$ can be chosen as the maximum ones used in the NSF equations to capture continuum flow dynamics accurately.

    In section 2 the gas kinetic equation and its asymptotic behavior at small Knudsen number is introduced. In section 3 the CIS to solve the time-dependent gas kinetic equation is introduced and the convergence rate of iterations is calculated by the Fourier stability analysis, at the whole range of Knudsen number and time step. In sections 4 and 5 two general synthetic iterative schemes are proposed and their convergence rates are calculated. That of HOLO is also calculated and compared. In section 6 the conditions for asymptotic NSF preserving are analyzed. Several challenging numerical examples are used to assess the accuracy and efficiency of the proposed GSIS for time-dependent problems in section 7.
    Finally, a summary is given in section 8.

    ##
    2 Kinetic equations and moment equations

    We take the linearized gas kinetic models to illustrate how the GSIS can be extended to unsteady problems; the extension to nonlinear Boltzmann equation is straightforward [18, 25]. In the absence of external acceleration, the kinetic model reads

    (1)

    $\frac{\partial h}{\partial t}+\bm{v}\cdot\frac{\partial h}{\partial\bm{x}}=\mathcal{L}_{s}(h)\equiv\mathcal{L}^{+}-\delta_{rp}h.$

    Here, $h(t,\bm{x},\bm{v})$ is the perturbation velocity distribution function that depends on the time $t\in\mathbb{R}^{+}$, molecular velocity $\bm{v}=(v_{1},v_{2},v_{3})\in\mathbb{R}^{3}$, and spatial coordinates $\bm{x}=(x_{1},x_{2},x_{3})\in\mathbb{R}^{3}$.
    The rarefaction parameter $\delta_{rp}$ (its inverse is known as the Knudsen number Kn) is a constant in linearized problems. When $\delta_{rp}$ is large, the kinetic model becomes stiff, which requires special care to make the solution converge fast and asymptotically preserve the NSF limit even with coarse spatial grid and large temporal step. The loss part of the collision operator is $\delta_{rp}h$, while the gain part is [13]

    (2)

    $\mathcal{L}^{+}=\delta_{rp}\left[\varrho+2\bm{u}\cdot\bm{v}+\tau\left(v^{2}-\frac{3}{2}\right)+\frac{4(1-\text{Pr})}{5}\bm{q}\cdot{\bm{v}}\left(v^{2}-\frac{5}{2}\right)\right]f_{eq},$

    where the global equilibrium velocity distribution function is

    (3)

    $f_{eq}=\frac{1}{\pi^{3/2}}\exp(-v^{2}),$

    and
    the macroscopic quantities, such as the perturbation density $\rho$, flow velocity $\bm{v}$, perturbation temperature $\tau$, stress $\sigma_{ij}$, and heat flux $\bm{q}$ are defined as the velocity moments of the velocity distribution function:

    (4)

    $\displaystyle\varrho=\int{h}d\bm{v},\quad\bm{u}=\int{\bm{v}h}d\bm{v},\quad\tau=\frac{2}{3}\int\left({v^{2}}-\frac{3}{2}\right)hd\bm{v},$

    $\displaystyle\sigma_{ij}=2\int{v_{\langle{i}}v_{j\rangle}h}d\bm{v}\equiv{}2\int\left(v_{i}v_{j}-\frac{v^{2}}{3}\delta_{ij}\right){h}d\bm{v},\quad\bm{q}=\int\left({v^{2}}-\frac{5}{2}\right)\bm{v}hd\bm{v}.$

    The kinetic model reduces to the linearized BGK [2] and Shakhov [14] models when the Prandtl number is $\text{Pr}=1$ and $2/3$, respectively. For real monatomic gas, the Prandtl number is approximately 2/3. Note that in the BGK model, the heat flux does not appear in the collision operator (2).

    Note that in the kinetic model (1) the spatial coordinates and molecular velocity have been normalized by the characteristic flow length $L$ and most probable speed $v_{m}=\sqrt{2RT_{0}}$, where $R$ is the gas constant and $T_{0}$ is the reference temperature. Also, the time has been normalized by $L/v_{m}$.

    On multiply Eq. (1) by 1, 2$\bm{v}$, and $v^{2}-\frac{3}{2}$, respectively, and integrate the resultant equations with respect to the molecular velocity $\bm{v}$, we obtain the following equations for the evolution of density, velocity, and temperature:

    (5)

    $\displaystyle\frac{\partial{\varrho}}{\partial{t}}+\frac{\partial{u_{i}}}{\partial{x_{i}}}=0,$

    $\displaystyle 2\frac{\partial{u_{i}}}{\partial{t}}+\frac{\partial{\varrho}}{\partial{x_{i}}}+\frac{\partial{\tau}}{\partial{x_{i}}}+\frac{\partial{{\sigma_{ij}}}}{\partial{x_{j}}}=0,$

    $\displaystyle\frac{3}{2}\frac{\partial{\tau}}{\partial{t}}+\frac{\partial{{q_{j}}}}{\partial{x_{j}}}+\frac{\partial{u_{j}}}{\partial{x_{j}}}=0,$

    which are not closed since the expressions for stress $\sigma_{ij}$ and heat flux $\bm{q}$ are not known. To close these moment systems, one either uses the Chapman-Enskog expansion [3] or the Grad moment method [7]. However, the resultant equations are usually limited to finite value of Knudsen number. For instance, to the first order of Knudsen number, the constitutive relation is exactly the Newton’s law for shear stress and Fourier’s law for heat conduction:

    (6)

    $\displaystyle\sigma_{ij}$
    $\displaystyle=-2\delta_{rp}^{-1}\frac{\partial u_{<i}}{\partial{x_{j>}}}\equiv-\delta_{rp}^{-1}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}-\frac{2}{3}\frac{\partial u_{k}}{\partial x_{k}}\delta_{ij}\right),\ q_{i}=-\frac{5}{4\text{Pr}}\delta_{rp}^{-1}\frac{\partial\tau}{\partial x_{i}}.$

    To describe the rarefied gas dynamics over the entire region of Knudsen numbers, it is necessary to solve the gas kinetic equation numerically.

    ##
    3 Conventional iterative scheme

    A direct method to solve the kinetic equation, which is a complicated integro-differential equation, is the use of conventional iteration scheme. In this paper, we consider the following typical temporal discretization of Eq. (1) for unsteady problems:

    (7)

    $\frac{h_{n+1}-h_{n}}{\Delta{t}}+\frac{\bm{v}}{2}\cdot\left(\frac{\partial{h_{n+1}}}{\partial\bm{x}}+\frac{\partial{h_{n}}}{\partial\bm{x}}\right)=r\mathcal{L}_{n+1}+(1-r)\mathcal{L}_{n},$

    where the quantities with the subscript $n$ are evaluated at the time $t_{n}$, $\Delta{t}=t_{n+1}-t_{n}$ is the time step, and the parameter $r$ varies between 0.5 and 1 [20, 26]. When $r=1/2$, the scheme is second-order accuracy in time, while when $r=1$ it is a backward Euler scheme with first-order temporal accuracy.

    Since $\mathcal{L}_{n+1}$ is a function of $h_{n+1}$, Eq. (7) must be solved iteratively. In CIS, given the value of velocity distribution function $h_{n+1}^{k}$ at the $k$-th iteration step (this is often called inner iteration in time-dependent implicit schemes [26]), its value at the next iteration step is calculated by:

    (8)

    $\frac{h^{k+1}_{n+1}-h_{n}}{\Delta{t}}+\frac{\bm{v}}{2}\cdot\left(\frac{\partial{h^{k+1}_{n+1}}}{\partial\bm{x}}+\frac{\partial{h_{n}}}{\partial\bm{x}}\right)=r\left({\mathcal{L}_{n+1}^{+,k}-\delta_{rp}h_{n+1}^{k+1}}\right)+(1-r){\mathcal{L}_{n}},$

    and this process repeats until the relative difference in macroscopic quantities between two consecutive inner iterations are less than a fixed value.

    We adopt the Fourier stability analysis to investigate the efficiency of this inner iteration, that is, to see how fast the error decays when $k$ increases. We define the error functions between velocity distribution functions at two consecutive inner iterations as

    (9)

    $Y^{k+1}(\bm{x},\bm{v})=h_{n+1}^{k+1}(\bm{x},\bm{v})-h_{n+1}^{k}(\bm{x},\bm{v}),$

    and the error functions for macroscopic quantities $M=[\varrho,\bm{u},\tau,\bm{q}]$ between two consecutive inner iteration steps as

    (10)

    $\displaystyle\Phi^{k+1}(\bm{x})\equiv$
    $\displaystyle\left[\Phi^{k+1}_{\varrho},\Phi^{k+1}_{\bm{u}},\Phi^{k+1}_{\tau},\Phi^{k+1}_{\bm{q}}\right]$

    $\displaystyle=$
    $\displaystyle M_{n+1}^{k+1}(\bm{x}\,)-M_{n+1}^{k}(\bm{x}\,)=\int{Y^{k+1}(\bm{x},\bm{v}\,)\phi(\bm{v}\,)}d\bm{v},$

    where

    (11)

    $\displaystyle\phi(\bm{v}\,)=$
    $\displaystyle\left[1,v_{1},v_{2},v_{3},\frac{2}{3}v^{2}-1,v_{1}\left(v^{2}-\frac{5}{2}\right),v_{2}\left(v^{2}-\frac{5}{2}\right),v_{3}\left(v^{2}-\frac{5}{2}\right)\right].$

    From Eq. (8), it can be easily found that $Y^{k+1}(\bm{x},\bm{v})$ satisfies

    (12)

    $\displaystyle\left(1+\frac{1}{r\Delta{t}\delta_{rp}}+\frac{1}{2r\delta_{rp}}\bm{v}\cdot\bm{\nabla}\right)Y^{k+1}=$

    $\displaystyle f_{eq}\left[\Phi^{k}_{\varrho}+2\Phi^{k}_{\bm{u}}\cdot\bm{v}+\Phi^{k}_{\tau}\left(v^{2}-\frac{3}{2}\right)+\frac{4(1-\text{Pr})}{5}\Phi^{k}_{\bm{q}}\cdot{\bm{v}}\left(v^{2}-\frac{5}{2}\right)\right],$

    because in the $(n+1)$-th time step, all variables in the $n$-th time step are fixed and hence eliminated.

    To determine the error decay rate $e$ we perform the Fourier stability analysis by seeking the eigenfunctions $\bar{Y}(\bm{v}\,)$ and $\alpha=[\alpha_{\varrho},\bm{\alpha}_{u},\alpha_{\tau},\bm{\alpha}_{q}]$ of the following forms:

    (13)

    $\displaystyle Y^{k+1}(\bm{x},\bm{v}\,)=e^{k}\bar{Y}(\bm{v}\,)\exp(i\bm{\theta}\cdot{\bm{x}}\,),\quad\Phi^{k+1}(\bm{x}\,)=e^{k+1}\alpha\exp(i\bm{\theta}\cdot{\bm{x}}\,),$

    where $\bm{\theta}=(\theta_{1},\theta_{2},\theta_{3})$ is the wavevector of perturbance and $i$ is the imaginary unit. The iteration is unstable when the error decay rate is larger than unity, while slow (fast) convergence occurs when the error decay rate $|e|$ approaches one (zero). Note that the two exponents in the right-hand-side of Eq. (13) are different, due to the fact that in CIS we first need macroscopic quantities to start the iteration.

    Figure 1:
    The error decay rate as a function of the Knudsen number in CIS, GSIS, and HOLO [20], when $r=1/2$. Note that the iteration is unstable when the error decay rate is larger than one.

    Obviously, from Eqs. (10) and (13) we have

    (14)

    $e\alpha=\int\bar{Y}(\bm{v}\,)\phi(\bm{v}\,)d\bm{v},$

    and from Eqs. (12) and (13), we have

    (15)

    $\displaystyle\bar{Y}(\bm{v})=$
    $\displaystyle\left[\alpha_{\varrho}+2\bm{\alpha}_{u}\cdot\bm{v}+\alpha_{\tau}\left(v^{2}-\frac{3}{2}\right)+\frac{4(1-\text{Pr})}{5}\bm{\alpha}_{q}\cdot\bm{v}\left({v}^{2}-\frac{5}{2}\right)\right]y_{0}(\bm{v}),$

    where

    (16)

    $y_{0}(\bm{v})=\frac{{f_{\text{eq}}}}{1+(r\Delta{t}\delta_{rp})^{-1}+i(2r\delta_{rp})^{-1}\bm{\theta}\cdot\bm{v}}.$

    On multiplying Eq. (15) with $\phi(\bm{v})$ and integrating the resultant equations with respect to the molecular velocity $\bm{v}$, we obtain eight linear algebraic equations for eight unknown elements in $\alpha$ with the help of Eq. (14). These algebraic equations can be written in the matrix form as
    $C_{8}\alpha^{\top}=e\alpha^{\top}$,
    where the superscript $\top$ is the transpose operator, and
    $8\times 8$ matrix is

    (17)

    $C_{8}=\int\left[1,2\bm{v},v^{2}-\frac{3}{2},\frac{4(1-\text{Pr})}{5}\bm{v}\left(v^{2}-\frac{5}{2}\right)\right]^{\top}\phi(\bm{v})y_{0}(\bm{v})d\bm{v}.$

    The error decay rate can be obtained by numerically computing the eigenvalues of the matrix $C_{8}$ and taking the maximum absolute value of $e$; results as a function of the Knudsen number, for both the BGK and Shakhov models, are shown in Fig. 1. It is clear that when the Knudsen number is large, $e$ goes to zero so that the error decays quickly. This means that CIS is very efficient for highly rarefied gas flows, i.e., the converged solution can be found within dozens of iterations [18]. On the contrary, $e\rightarrow 1$ when $\text{Kn}\rightarrow 0$, which means that the CIS is extremely slow in the (near) continuum flow regime.

    ##
    4 The scheme II GSIS

    To expedite the convergence of inner iteration, that is, to reduce the number of $k$ in Eq. (8), macroscopic synthetic equations are needed. There are several ways of constructing these equations, and we first develop the GSIS-II for time-dependent problems due to its relative simplicity [25].

    First, in GSIS, given the value of velocity distribution function $h_{n+1}^{k}$ at the $k$-th iteration step, its value at the intermediate $(k+1/2)$-th step is obtained in a similar way to Eq. (8):

    (18)

    $\frac{h^{k+1/2}_{n+1}-h_{n}}{\Delta{t}}+\frac{\bm{v}}{2}\cdot\left(\frac{\partial{h^{k+1/2}_{n+1}}}{\partial\bm{x}}+\frac{\partial{h_{n}}}{\partial\bm{x}}\right)=r\left({\mathcal{L}_{n+1}^{+,k}-\delta_{rp}h_{n+1}^{k+1/2}}\right)+(1-r){\mathcal{L}_{n}}.$

    This velocity distribution function $h^{k+1/2}_{n+1}$ will be used to construct high-order constitutive relations in macroscopic synthetic equations; and when the synthetic equations are solved to obtain macroscopic quantities, say, $M^{k+1}_{n+1}=[\varrho,\bm{u},\tau,\bm{q}]$, they will be used in the gain term (2) for the next inner iteration, until convergence criterion is met.

    Certainly, the macroscopic synthetic equations should be derived exactly from the gas kinetic model. In GSIS-II, the constitutive relations are constructed, with a free parameter $\delta$, in the following manner [25]:

    (19)

    $\displaystyle\sigma^{k+1}_{ij}=-2\delta^{-1}\frac{\partial u^{k+1}_{<i}}{\partial{x_{j>}}}+\left(\sigma^{k+1/2}_{ij}+2\delta^{-1}\frac{\partial u^{k+1/2}_{<i}}{\partial{x_{j>}}}\right),$

    (20)

    $\displaystyle q^{k+1}_{i}=-\frac{5}{4\mathrm{Pr}}\delta^{-1}\frac{\partial\tau^{k+1}}{\partial x_{i}}+\left(q_{i}^{k+1/2}+\frac{5}{4\mathrm{Pr}}\delta^{-1}\frac{\partial\tau^{k+1/2}}{\partial x_{i}}\right).$

    At first sight it seems that the GSIS-II is very similar to the HOLO [20]. However, in HOLO the effective rarefaction parameter $\delta$ is chosen to be infinity in Eqs. (19) and (20), while in GSIS the NSF constitution relations are explicitly included in the macroscopic synthetic equations by setting $\delta=\delta_{rp}$, which allows fast convergence and unconditional stability when using large time step, as will be shown later.

    The macroscopic synthetic equations are solved by the following Crank-Nicolson scheme (although other schemes can also be used):

    (21)

    $\displaystyle\frac{\varrho^{k+1}}{\Delta{t}}+\frac{1}{2}\frac{\partial{u^{k+1}_{i}}}{\partial{x_{i}}}=$
    $\displaystyle\frac{\varrho_{n}}{\Delta{}t_{n}}-\frac{1}{2}\frac{\partial{u_{i,n}}}{\partial{x_{i}}},$

    $\displaystyle 2\frac{u_{i}^{k+1}}{\Delta{t}}+\frac{1}{2}\frac{\partial}{\partial{x_{i}}}\left(\varrho^{k+1}+{\tau^{k+1}}+\sigma^{k+1}_{ij}\right)=$
    $\displaystyle 2\frac{u_{i,n}}{\Delta{t}}-\frac{1}{2}\frac{\partial}{\partial{x_{i}}}\left({\varrho_{n}}+\tau_{n}+\sigma_{ij,n}\right),$

    $\displaystyle\frac{3}{2}\frac{\tau^{k+1}}{\Delta{t}}+\frac{1}{2}\frac{\partial{{q^{k+1}_{i}}}}{\partial{x_{i}}}+\frac{1}{2}\frac{\partial{{u^{k+1}_{i}}}}{\partial{x_{i}}}=$
    $\displaystyle\frac{3}{2}\frac{\tau_{n}}{\Delta{t}}-\frac{1}{2}\frac{\partial{q_{i,n}}}{\partial{x_{i}}}-\frac{1}{2}\frac{\partial{u_{i,n}}}{\partial{x_{i}}},$

    where terms in the left-hand-side are evaluated at the $(n+1)$-th time step (for clarity the subscript is ignored), while these on the right-hand-side are evaluated at the $n$-th time step.

    Therefore,
    to calculate the convergence rate of the time-dependent GSIS, the error functions are defined as

    (22)

    $\displaystyle Y^{k+1/2}(\bm{x},\bm{{v}}\,)=h_{n+1}^{k+1/2}(\bm{x},\bm{{v}}\,)-h_{n+1}^{k}(\bm{x},\bm{{v}}\,)=e^{k}\bar{Y}(\bm{{v}}\,)\exp(i\bm{\theta}\cdot{\bm{x}}\,),$

    $\displaystyle\Phi^{k+1}(\bm{x}\,)=M_{n+1}^{k+1}(\bm{x}\,)-M_{n+1}^{k}(\bm{x}\,)=e^{k+1}\alpha\exp(i\bm{\theta}\cdot{\bm{x}}).$

    Note that here the macroscopic quantities $M^{k+1}$ are calculated from the synthetic equations (21), rather than directly from the velocity distribution function $h^{k+1/2}$. After some algebra, the error decay rate can be obtained by solving the following linear systems:

    (23)

    $\displaystyle e\left[\frac{2}{\Delta{t}}\alpha_{\varrho}+i\theta_{1}\alpha_{u_{1}}+i\theta_{2}\alpha_{u_{2}}+i\theta_{3}\alpha_{u_{3}}\right]=0,$

    $\displaystyle e\left[\left(\frac{4}{\Delta{t}}+\frac{{\theta}^{2}}{\delta}\right)\alpha_{u_{j}}+i\theta_{j}(\alpha_{\varrho}+\alpha_{\tau})\right]=S_{j+1},$

    $\displaystyle e\left[\frac{3}{\Delta{t}}\alpha_{\tau}+i\theta_{1}(\alpha_{q_{1}}+\alpha_{u_{1}})+i\theta_{2}(\alpha_{q_{2}}+\alpha_{u_{2}})+i\theta_{3}(\alpha_{q_{3}}+\alpha_{u_{3}})\right]=0,$

    $\displaystyle e\left(\frac{5i}{4\text{Pr}}\theta_{j}{\delta^{-1}}\alpha_{\tau}+\alpha_{q_{1}}\right)=S_{j+5},$

    where $j=1,2,3$, and the source terms are

    (24)

    $\displaystyle S_{2}=\int\left[{\delta^{-1}}{\theta}^{2}{}{v}_{1}-2i\theta_{1}v_{\langle{1}}v_{1\rangle}-2i\theta_{2}{v}_{1}{v}_{2}-2i\theta_{3}{v}_{1}{v}_{3}\right]\bar{Y}(\bm{{v}}\,)\mathrm{d}^{3}\bm{{v}},$

    $\displaystyle S_{3}=\int\left[{\delta^{-1}}{\theta}^{2}{}{v}_{2}-2i\theta_{2}v_{\langle{2}}v_{2\rangle}-2i\theta_{1}{v}_{1}{v}_{2}-2i\theta_{3}{v}_{2}{v}_{3}\right]\bar{Y}(\bm{{v}}\,)\mathrm{d}^{3}\bm{{v}},$

    $\displaystyle S_{4}=\int\left[{\delta^{-1}}{\theta}^{2}{}{v}_{3}-2i\theta_{3}v_{\langle{3}}v_{3\rangle}-2i\theta_{1}{v}_{1}{v}_{3}-2i\theta_{2}{v}_{2}{v}_{3}\right]\bar{Y}(\bm{{v}}\,)\mathrm{d}^{3}\bm{{v}},$

    $\displaystyle S_{j+5}=\int\left[\frac{5i}{4\text{Pr}}\theta_{j}{\delta^{-1}}{}\left(\frac{2}{3}{v}^{2}-1\right)+{v}_{j}\left({v}^{2}-\frac{5}{2}\right)\right]\bar{Y}(\bm{{v}}\,)\mathrm{d}^{3}\bm{{v}}.$

    Equations (23) and (24) can be rearranged as
    $L_{8}e\alpha_{M}^{\top}=R_{8}\alpha_{M}^{\top}$,
    and the error decay rate can be obtained from the eigenvalues of the matrix $L_{8}^{-1}R_{8}$. A comparison of GSIS-II and HOLO is shown in Fig. 1, which clearly show that HOLO is unstable at large time step, say, when $\Delta{t}=15$ and $|\bm{\theta}|=1$. When the perturbation wavevector $\theta$ in Eq. (13) is increased, the time step for stable iteration in HOLO is much reduced, while the GSIS is always stable, see the last subfigure.

    ##
    5 The scheme I GSIS

    This time-dependent scheme is modified from Ref. [18] which is initially developed to find steady-state solutions of the Boltzmann equation and kinetic model equations. In addition the the five macroscopic equations from the mass, momentum and energy consecration (5), evolution equations for the stress and heat flux are also included, like that in the Grad 13 moment equations [7]. That is, we multiply Eq. (1) by $v_{\langle{i}}v_{j\rangle}$ and integrate the resultant equation with respect to the molecular velocity $\bm{v}$, and obtain

    (25)

    $\frac{\partial\sigma_{ij}}{\partial{t}}+2\int{v_{\langle{i}}v_{j\rangle}}\bm{v}\cdot\frac{\partial h}{\partial\bm{x}}d\bm{v}=-\delta_{rp}\sigma_{ij}.$

    However, unlike the Grad 13 moment system where the high order term (i.e. the second term in the left-hand-side) is closed by expanding the velocity distribution function into Hermite polynomial of molecular velocity, in GSIS-I Eq. (25) is rearranged in the following form to reflect Newton’s law for shear viscosity, which shall allow numerical stability at large time steps:

    (26)

    $\displaystyle\frac{\partial\sigma_{ij}}{\partial{t}}$
    $\displaystyle+\underbrace{2\int{v_{\langle{i}}v_{j\rangle}}\bm{v}\cdot\frac{\partial h}{\partial\bm{x}}d\bm{v}-2\frac{\partial{u_{<i}}}{\partial{x_{j>}}}}_{\text{HoT}_{\sigma_{ij}}}+\underbrace{2\frac{\partial{u_{<i}}}{\partial{x_{j>}}}=-\delta_{rp}\sigma_{ij}}_{\text{Newton's law}}.$

    This equation is solved, again, by the Crank-Nicolson scheme:

    (27)

    $\displaystyle\left(\frac{1}{\Delta{t}}+\frac{\delta_{rp}}{2}\right)$
    $\displaystyle\sigma^{k+1}_{ij}+\frac{\partial{u^{k+1}_{<i}}}{\partial{x_{j>}}}$

    $\displaystyle=\left[\left(\frac{1}{\Delta{t}}-\frac{\delta_{rp}}{2}\right)\sigma_{ij,n}-\frac{\partial{u_{<i,n}}}{\partial{x_{j>}}}-\frac{\text{HoT}_{\sigma_{ij},n}}{2}\right]-\frac{\text{HoT}^{k+1/2}_{\sigma_{ij}}}{2}.$

    It is noted that, when the inner iteration converges, $\text{HoT}^{k+1/2}_{\sigma_{ij}}$ will be the same as $\text{HoT}^{k+1}_{\sigma_{ij}}$, and hence this numerical scheme is a Crank-Nicolson scheme with second-order accuracy in time.

    Likewise, we multiply Eq. (1) by $v_{i}(v^{2}-5/2)$ and integrate the resultant equation with respect to $\bm{v}$; we obtain

    (28)

    $\displaystyle\frac{\partial q_{i}}{\partial{t}}$
    $\displaystyle+\underbrace{\int{\left(v^{2}-\frac{5}{2}\right)}v_{i}\bm{v}\cdot\frac{\partial h}{\partial\bm{x}}d\bm{v}-\frac{5}{4}\frac{\partial{\tau}}{\partial{x_{i}}}}_{\text{HoT}_{q_{i}}}+\underbrace{\frac{5}{4}\frac{\partial{\tau}}{\partial{x_{i}}}=-\delta_{rp}\text{Pr}q_{i}}_{\text{ Fourier's law}},$

    which can also be solved by the Crank-Nicolson scheme.

    Although the macroscopic synthetic equations (5), (26) and (28) resemble the Grad 13 moment equations [7, 16], no approximations are introduced here since the higher-order terms are computed directly from the velocity distribution function.

    According to the Fourier stability analysis, the error decay rate can be obtained by solving the following equations:

    (29)

    $\displaystyle e\left[\frac{2}{\Delta{t}}\alpha_{\varrho}+i\theta_{1}\alpha_{u_{1}}+i\theta_{2}\alpha_{u_{2}}+i\theta_{3}\alpha_{u_{3}}\right]=0,$

    $\displaystyle e\left[\left(\frac{4}{\Delta{t}}+\frac{{\theta}^{2}}{\delta_{rp}+\frac{2}{\Delta{t}}}\right)\alpha_{u_{j}}+i\theta_{j}(\alpha_{\varrho}+\alpha_{\tau})\right]=S_{j+1},$

    $\displaystyle e\left[\frac{3}{\Delta{t}}\alpha_{\tau}+i\theta_{1}(\alpha_{q_{1}}+\alpha_{u_{1}})+i\theta_{2}(\alpha_{q_{2}}+\alpha_{u_{2}})+i\theta_{3}(\alpha_{q_{3}}+\alpha_{u_{3}})\right]=0,$

    $\displaystyle e\left(\frac{5i\theta_{j}\alpha_{\tau}}{4(\delta_{rp}\text{Pr}+\frac{2}{\Delta{}t})}+\alpha_{q_{j}}\right)=S_{j+5}$

    where $j=1,2,3$, and the source terms are

    (30)

    $\displaystyle S_{2}=\frac{1}{\delta_{rp}+\frac{2}{\Delta{t}}}\int\left[{\theta}^{2}{v}_{1}-{2\Theta}(\theta_{1}v_{\langle{1}}v_{1\rangle}+\theta_{2}{v}_{1}{v}_{2}+\theta_{3}{v}_{1}{v}_{3})\right]\bar{Y}(\bm{{v}}\,)\mathrm{d}^{3}\bm{{v}},$

    $\displaystyle S_{3}=\frac{1}{\delta_{rp}+\frac{2}{\Delta{t}}}\int\left[{\theta}^{2}{v}_{2}-{2\Theta}(\theta_{1}{v}_{1}{v}_{2}+\theta_{2}v_{\langle{2}}v_{2\rangle}+\theta_{3}{v}_{2}{v}_{3})\right]\bar{Y}(\bm{{v}}\,)\mathrm{d}^{3}\bm{{v}},$

    $\displaystyle S_{4}=\frac{1}{\delta_{rp}+\frac{2}{\Delta{t}}}\int\left[{\theta}^{2}{v}_{3}-{2\Theta}(\theta_{1}{v}_{1}{v}_{3}+\theta_{2}{v}_{2}{v}_{3}+\theta_{3}v_{\langle{3}}v_{3\rangle})\right]\bar{Y}(\bm{{v}}\,)\mathrm{d}^{3}\bm{{v}},$

    $\displaystyle S_{j+5}=\frac{i}{\delta_{rp}\text{Pr}+\frac{2}{\Delta{t}}}\int\left[\frac{5}{4}\theta_{j}\left(\frac{2}{3}{v}^{2}-1\right)-{\Theta}{v}_{j}\left({v}^{2}-\frac{5}{2}\right)\right]\bar{Y}(\bm{{v}}\,)\mathrm{d}^{3}\bm{{v}}.$

    and $\Theta=\theta_{1}{v}_{1}+\theta_{2}v_{2}+\theta_{3}v_{3}$.

    Numerical results for the BGK and Shakhov models when $\Delta{t}$=1 and 15 are shown in Fig. 1. For the Shakhov kinetic model, GSIS-I has a smaller error decay rate than GSIS-II, especially when $\text{Kn}\rightarrow 0$ the error decay rate of GSIS-I goes to zero, while that of GSIS-II goes to $1/3$. This is because the heat flux appears in the gain term (2) of the Shakhov model, and the GSIS-I has the synthetic equation (28) to guild the evolution of heat flux. However, GSIS-II does not have this capability, which results in a slower convergence (larger error decay rate $e$) than GSIS-I. For the BGK model, both GSIS schemes have the error decay rate approaching zero when $\text{Kn}\rightarrow 0$, because only the density, velocity and temperature appears in the collision term, and both GSIS schemes have the evolution equations for these macroscopic quantities.

    ##
    6 The property of asymptotic preserving

    From a practical point of view, the property of asymptotic NSF preserving should be investigated based on the numerical scale solving the kinetic equation and macroscopic synthetic equations [8]. Since in GSIS the kinetic equation and synthetic equations can be solved by different numerical methods with different order of accuracy, in this section, we consider the influence of spatial and temporal discretizations in the gas kinetic solver on the accuracy of GSIS, based on the assumptions that the spatial grid size $\Delta{x}$ and the time step $\Delta{t}$ are refined enough to capture the physical solution of NSF equations. Namely, we investigate at what values of $\alpha$ and $\beta$, can the macroscopic synthetic equations be exactly reduced to the NSF equations when Kn is small, through the Chapman-Enskog expansion [3] of the discretized gas kinetic equation, with the following scaling:

    (31)

    $\displaystyle{\Delta{x}}\sim{\text{Kn}^{1/\alpha}},\quad{\Delta{t}}\sim{\text{Kn}^{1/\beta}}.$

    Note that the spatial grid size $\Delta{x}$ and time step $\Delta{t}$ have been normalized by the characteristic flow length $L$ and time ($L/v_{m}$), respectively.
    Here $\alpha$ and $\beta$ denote the order of accuracy in the asymptotic preserving of NSF equations. Clearly, the larger the values of $\alpha$ and $\beta$, the better the numerical scheme. If $\alpha=\infty$, the scheme will capture the hydrodynamical behavior when $\Delta{x}$ is approximately the system size (no matter what the value of Kn is), as long as this size is adequate to capture the flow physics. On the other hand, if $\beta=\infty$, the scheme will capture the hydrodynamical behavior when $\Delta{t}$ is approximately the characteristic time of the system (e.g., the oscillation period of a sound wave).

    Since when the Knudsen number is small, the error decay rate of GSIS is much smaller than unity (i.e., zero in GSIS-I for both BGK and Shakhov models, and 1/3 in GSIS-II scheme when the Shakhov model is used), the converged solution can be found within a few inner iterations. Thus, we have $h_{eq}^{k+1}=h_{eq}^{k}$ and $M^{k+1}=M^{k}$. When the inner iteration is converged, the iterative scheme (8) can be expressed as

    (32)

    $\frac{\partial h}{\partial t}+\bm{v}\cdot\frac{\partial{h}}{\partial\bm{x}}+O(\Delta{t}^{m})\delta_{t}(h)+O(\Delta{x}^{n})\delta_{x}(h)=\mathcal{L}_{s},$

    where $n$ is the order of approximation for the spatial derivative in the kinetic equation, while $\delta_{x}(h)$ is the $(n+1)$-th order derivative of $h$ with respect to spatial coordinates. For instance, if the second-order upwind finite difference scheme, we have $n=2$.
    Similar applies to $m$ and $\delta_{t}(h)$.

    ###
    6.1 Chapman-Enskog expansion

    In the Chapman-Enskog expansion the velocity distribution function is approximated by the Taylor expansion $h=h^{(0)}+\text{Kn}h^{(1)}+\text{Kn}^{2}h^{(2)}+\cdots$,
    so are the stress $\sigma_{ij}=\sum_{\ell=0}^{\infty}\text{Kn}^{\ell}\sigma_{ij}^{(\ell)}$ and heat flux $\bm{q}=\sum_{\ell=0}^{\infty}\text{Kn}^{\ell}\bm{q}^{(\ell)}$,
    where $\sigma^{(\ell)}_{ij}=2\int{v_{\langle{i}}v_{j\rangle}h^{(\ell)}}d\bm{v}$ and
    $\bm{q}^{(\ell)}=\int\left({v^{2}}-\frac{5}{2}\right)\bm{v}h^{(\ell)}d\bm{v}$. However, the five conservative variables $C_{M}=\{\rho,\bm{u},\tau\}$ are calculated only according to the zeroth-order expansion. That is,

    (33)

    $\displaystyle\rho=\int{}h^{(0)}d\bm{v},\quad\bm{u}=\int{}\bm{v}h^{(0)}d\bm{v},\quad\tau=\frac{2}{3}\int{}\left(v^{2}-\frac{3}{2}\right)h^{(0)}d\bm{v},$

    with the compatibility condition $\int{}h^{(\ell)}d\bm{v}=\int{}\bm{v}h^{(\ell)}d\bm{v}=\int{}v^{2}h^{(\ell)}d\bm{v}=0$ for $\ell\geq 1$.
    From Eq. (33) and the compatibility condition, one finds that the time derivatives in Eq. (5) can be formally written as a series in Kn [16]:
    $\frac{\partial}{\partial{}t}=\sum_{\ell=0}^{\infty}\text{Kn}^{\ell}\frac{\partial}{\partial{}t_{\ell}}$.

    ###
    6.2 GSIS-I

    By substituting the Taylor expansion of velocity distribution function into Eq. (32) and collecting terms with the order of $\text{Kn}^{-1}$, we have

    (34)

    $h^{(0)}=\left[\varrho+2\bm{u}\cdot\bm{v}+\tau\left(v^{2}-\frac{3}{2}\right)\right]f_{eq},$

    and $\sigma^{(0)}_{ij}=\bm{q}^{(0)}=0$, with the following largest scaling:

    (35)

    $\displaystyle\Delta{x}\sim{\text{Kn}^{1/\infty}}=O(1),\quad\Delta{t}\sim{}{\text{Kn}^{1/\infty}}=O(1).$

    Under this circumstance, by collecting terms with the order $\text{Kn}^{0}$, we have

    (36)

    $h^{(1)}=-\frac{\partial h^{(0)}}{\partial t_{0}}-\bm{v}\cdot\frac{\partial{h^{(0)}}}{\partial\bm{x}}\underbrace{-O(\Delta{t}^{m})\delta_{t}(h)-O(\Delta{x}^{n})\delta_{x}(h)}_{\text{discretization error terms}}.$

    Note that in the standard Chapman-Enskog expansion, this $h^{(1)}$ term is used to produce the NSF constitutive relations (6). Although there are some error terms in $h^{(1)}$, this does not affect the exact derivation of NSF constitutive relations in GSIS-I. This is because we only need $h^{(0)}$ to evaluate the stress and heat flux according to Eqs. (26) and (28), while the $h^{(1)}$ term will result constitutive relations at the order of $\text{Kn}^{2}$. To be specific, let us take the stress for an example. When the inner iteration is converged, Eq. (26)
    becomes

    (37)

    $\displaystyle\frac{\partial\sigma_{ij}}{\partial{t}}$
    $\displaystyle+2\int{v_{\langle{i}}v_{j\rangle}}\bm{v}\cdot\frac{\partial h}{\partial\bm{x}}d\bm{v}=-\delta_{rp}\sigma_{ij}.$

    Since $\sigma_{ij}\propto\text{Kn}$, the leading order solution is

    (38)

    $\displaystyle\sigma^{(1)}_{ij}$
    $\displaystyle=-\frac{2}{\delta_{rp}}\int{v_{\langle{i}}v_{j\rangle}}\bm{v}\cdot\frac{\partial h^{(0)}}{\partial\bm{x}}d\bm{v}=-\frac{2}{\delta_{rp}}\frac{\partial u_{<i}}{\partial{x_{j>}}},$

    which is exactly Newton’s law of stress.

    Therefore, GSIS-I asymptotically preserves the NSF equations with the largest scalings (35). That is to say, as long as the spatial resolution $\Delta{x}=O(1)$ and temporal resolution $\Delta{t}=O(1)$ are able to capture the physical solution of NSF equations, GSIS-I is able to recover the linearized NSF equations when $\text{Kn}\rightarrow 0$. This means that the overall order of accuracy of GSIS-I depends only on the order of accuracy in solving macroscopic synthetic equations.
    In reality, however, such a large spatial grid size $\Delta{x}=O(1)$ cannot be used in regions with Knudsen layer or shock structure, where the physical solutions require a spatial resolution of $O(\text{Kn})$. Fortunately, these kinetic layers only take up a small fraction of the computational domain, say, in the vicinity of solid walls, which can be captured by implicit schemes with non-uniform spatial discretization. This will be demonstrated in section 7.3 below. Likewise, such a large temporal step $\Delta{t}=O(1)$ will be much reduced to the maximum time step in solving the NSF equations accurately; this will be demonstrated in section 7.1.

    ###
    6.3 GSIS-II

    We now follow the standard procedure to see at what conditions can the NSF equations be recovered from the macroscopic synthetic equations in GSIS-II. At the zeroth-order approximation, we have Eq. (34), and the following Euler equation:

    (39)

    $\displaystyle\frac{\partial{\varrho}}{\partial{t_{0}}}+\frac{\partial{u_{i}}}{\partial{x_{i}}}=0,\quad 2\frac{\partial{u_{i}}}{\partial{t_{0}}}+\frac{\partial{\varrho}}{\partial{x_{i}}}+\frac{\partial{\tau}}{\partial{x_{i}}}=0,\quad\frac{3}{2}\frac{\partial{\tau}}{\partial{t_{0}}}+\frac{\partial{u_{j}}}{\partial{x_{j}}}=0,$

    that relates the time derivative to spatial derivatives.

    With the largest scaling (35) $h^{(1)}$ is given by Eq. (36). If one calculates the stress and heat flux according to Eq. (4), then the constitutive relations are not exactly the NSF ones due to the presence of error terms. That is to say, in general, the NSF equations cannot be exactly recovered at such large temporal step and spatial size. In practice, if the error terns are very small, GSIS-II can also approach the NSF limit in the continuum flow regime. This will be demonstrated in section 7.2.

    In GSIS-II, in order to recover the NSF constitutive relations (6) exactly in general cases,
    the velocity distribution function in the Taylor expansion to the first-order of Knudsen number must be exactly recovered as

    (40)

    $\displaystyle h^{(1)}=$
    $\displaystyle-\frac{\partial h^{(0)}}{\partial t_{0}}-\bm{v}\cdot\frac{\partial{h^{(0)}}}{\partial\bm{x}}=-\left[2v_{\langle{i}}v_{j\rangle}\frac{\partial u_{\langle{i}}}{\partial x_{j\rangle}}+\left(v^{2}-\frac{5}{2}\right)v_{i}\frac{\partial\ln\tau}{\partial x_{i}}\right]f_{eq}.$

    where $h^{(0)}$ is given in Eq. (34), and the time derivative is changed to spatial derivatives with the help of Eq. (39). This requires the following scaling

    (41)

    $\displaystyle\Delta{x}\sim{\text{Kn}^{1/n}},\quad\Delta{t}\sim{\text{Kn}^{1/m}}.$

    Figure 2:
    Schematic of the asymptotic path to the limiting hydrodynamic flow regimes. The NSF equations are valid in the dark gray region. The region below the line $\Delta=O(\text{Kn})$ suggests resolved kinetic scale. The
    red solid line stands for the maximum spatial grid size/time step to solve the NSF equations accurately using some discretization method (certainly $\Delta$ depends on the numerical method for NSF equations and specific flow problems; $\Delta$ could be larger than $O(1)$ because the time step can be proportional to $\text{Kn}^{-1}$ in some flows), below which GSIS-I can capture the hydrodynamic behavior. GSIS-II works at a smaller value of $\Delta$, say, beneath the dash-dotted line $\Delta=O(\text{Kn}^{n})$. Due to the instability at large time step, the validation range of HOLO is smaller than GSIS-II.

    Normally the kinetic equation is solved with second-order accuracy both in temporal and spatial directions, that is, $n=m=2$. Therefore, the NSF equations are recovered in GSIS-II with $\Delta{x}\sim\sqrt{\text{Kn}}$ and $\Delta{t}\sim\sqrt{\text{Kn}}$, like the (discrete) unified gas-kinetic scheme [23, 9, 8].

    The asymptotic path to the limiting hydrodynamic flow regimes are summarized in Figure 2.

    ##
    7 Numerical results

    In this section we propose several solid numerical examples to assess the accuracy and efficiency of both GSIS schemes, as well as HOLO.

    ###
    7.1 Rayleigh-Brillouin scattering

    The coherent Rayleigh-Brillouin scattering is a promising technique to probe the property of gas, in which the wavelike density perturbation in gas is created by a moving optical lattice. We choose this problem because it is a zero-dimensional problem so any one can quickly test/compare our and their methods.

    Applying the Fourier transform in the scattering (say, $x_{1}$) direction, the governing equation for the velocity distribution function $h$ can be written as [22]

    (42)

    $\frac{\partial{h}}{\partial{t}}+2i\pi{}v_{1}h=\mathcal{L}(h)+2v_{1}\cos(2\pi{f_{s}}t)f_{eq},$

    where $f_{s}$ is the frequency of the moving optical lattice; we choose $f_{s}=\sqrt{5/6}$ (i.e. the sound speed normalized by the most probable speed of gas molecules) so that the amplitude of perturbed density will scale as $1/\text{Kn}$ when the Knudsen number is small. If the numerical scheme cannot preserve the NSF limit when $\text{Kn}\rightarrow 0$, then a very small time step is needed to keep a low numerical dissipation; otherwise the amplitude will be much smaller than the converged solution.

    Since the spatial derivative $\partial/\partial{x_{1}}$ is replaced by $2i\pi$ in the Fourier transform, the conservation equations (5) becomes

    (43)

    $\displaystyle\frac{\partial{\varrho}}{\partial t}+2i\pi{{u}_{1}}=0,\quad\frac{3}{2}\frac{\partial{\tau}}{\partial t}+2i\pi(u_{1}+q_{1})=0,$

    (44)

    $\displaystyle 2\frac{\partial{u}_{1}}{\partial t}+2i\pi(\varrho+\tau+\sigma_{11})=2\cos(2\pi{f_{s}}t),$

    and the same transform is applied to the synthetic equations for stress and heat flux in Eqs. (19), (20), (26), and (28).

    Figure 3:
    The imaginary part of the density perturbation $\varrho$ in the coherent Rayleigh-Brillouin scattering, when the rarefaction parameter is $\delta_{rp}=1000$. Inset shows the iteration number, where the iteration is terminated when the relative error in macroscopic quantities between two consecutive iterations is less than $10^{-10}$. The kinetic equations are solved with a time step of $\Delta{t}=0.01$, by the backward Euler method and Crank-Nicolson scheme, respectively. The macroscopic synthetic equations are always solved by the Crank-Nicolson scheme with the same time step.

    In the numerical simulation, the molecular velocity space $\bm{v}$ is discretized by $6\times 6\times 6$ Gauss-Hermite quadrature, which is accurate to resolve the velocity distribution function when the Knudsen number is small: in the continuum flow regime the velocity distribution function contains $v^{3}_{i}$, therefore, the integrand in the high order term in Eq. (28) has the highest polynomial of $v^{7}_{i}$, while the Gauss-Hermite quadrature of order 6 is accurate for polynomial up to the order of 11. Starting from the zero initial values for the velocity distribution function and macroscopic quantities, the synthetic equations are solved by the Crank-Nicolson scheme, while the kinetic equation is solved by backward Euler and Crank-Nicolson schemes, respectively. The inner iteration terminates when the relative error in density between two consecutive steps are less than $10^{-10}$.

    It should be noted that in HOLO, in order to make the mesoscopic and macroscopic equations consistent, consistency terms are introduced to enslave the solution of synthetic equations to that of the kinetic equation [20]. In this problem, since the spatial derivative is handled exactly by the Fourier transform, and the molecular velocity space is discretized by adequate quadrature, these consistency terms vanish. In GSIS, however, we enslave the solution of kinetic equation to that of the synthetic equations, since the kinetic equation converges so slowly that false convergence might happen [1, 18]. That is, when the inner iteration is converged (judged by the relative error in macroscopic quantities), the velocity distribution function is updated to reflect the converged macroscopic quantities:

    (45)

    $\displaystyle h^{k+1}=h^{k+1/2}$
    $\displaystyle+\left({\varrho}^{k+1}-\varrho^{k+1/2}\right)+2\left({\bm{u}}^{k+1}-\bm{u}^{k+1/2}\right)\cdot\bm{v}$

    $\displaystyle+\left({\tau}^{k+1}-\tau^{k+1/2}\right)\left(v^{2}-\frac{3}{2}\right).$

    Numerical results of the GSIS and HOLO for coherent Rayleigh-Brillouin scattering is shown in Fig. 3. From the figure we see that, when the kinetic equation is solved with second-order accuracy, GSIS-I, GSIS-II, and HOLO yield identical solutions. For the BGK model, in each time step, CIS needs around 200 inner iterations to find converged solutions (not shown), while GSIS-I, GSIS-II, and HOLO needs about 5 iterations. This is consistent with the Fourier stability analysis presented in Fig. 1, since the wavevector in this problem is always unity. For the Shakhov model, GSIS-I still needs 5 iterations, while GSIS-II and HOLO need about 20 iterations at each time step. This is because the former has synthetic equation for the evolution of heat flux that appears in the gain term of kinetic model so that the error decay rate approaches zero when $\text{Kn}\rightarrow 0$, while the latter has the error decay rate approaching $1/3$.

    We also find that, when the kinetic equation is solved with first-order temporal accuracy, both GSIS schemes yield accurate solutions. However, HOLO exhibits visible dissipation as the density amplitude is smaller. This is because the time step should be no larger than $O(\text{Kn})$ for first-order scheme, otherwise the artificial viscosity is comparable or larger than the physical viscosity, which leads to inaccurate solutions. It is very interesting to note that, if we turn off the correction of velocity distribution function in Eq. (45), both GSIS schemes produce the same inaccurate result as HOLO. This clearly demonstrates that we should enslave the solution of kinetic equation to that of the synthetic equation, while the HOLO does the reverse way and get inaccurate results.

    Finally, it should be noted that although GSIS-I allows a time step of $O(1)$ to asymptotic preserve the NSF limit, here it is chosen as $\Delta{t}=0.01$ because the NSF equations cannot be accurately solved by a second-order scheme when the time step is larger than 0.01 in this specific problem.

    Figure 4:
    Velocity profiles in the two-dimensional Taylor vortex flow at the time $t=200$, when the BGK model with $\delta_{rp}=10000$ is solved by GSIS. The profiles with time step $\Delta{t}=50$ are shifted upward for clarity. DVM 1st (2nd) order means that the kinetic equation is solved with first (second) order temporal accuracy. Inset: velocity contour at $t=200$.

    ###
    7.2 Decay of two-dimensional Taylor vortex

    To test the property of fast converging and asymptotic NSF preserving of GSIS and HOLO at the time step $\Delta{t}\sim{O(1)}$, we consider the decay of two-dimensional incompressible Taylor vortex within a periodic domain $0\leq{x_{1}},x_{2}\leq 1$. In the continuum flow regime, the flow is governed by the incompressible Navier-Stokes equations and has the following analytical solution:

    (46)

    $\displaystyle u_{1}(x_{1},x_{2},t)=$
    $\displaystyle-\cos(2\pi{x_{1}})\sin(2\pi{x_{2}})\exp\left(-{4\pi^{2}}t/{\delta_{rp}}\right),$

    $\displaystyle u_{2}(x_{1},x_{2},t)=$
    $\displaystyle~{}~{}~{}\sin(2\pi{x_{1}})\cos(2\pi{x_{2}})\exp\left(-{4\pi^{2}}t/{\delta_{rp}}\right),$

    where the density and temperature are always zero. Therefore, the BGK model is used.

    This is an ideal test case to assess the asymptotic NSF preserving property, since when $\delta_{rp}$ is large (Kn is small), the Taylor vortex takes a long time to decay, so the time step can be made very large if the numerical scheme for the kinetic system asymptotic preserves the NSF equations.

    Without loss of generality we choose $\delta_{rp}=10^{4}$, and in order to focus only on the error in temporal discretization, we solve the kinetic equation and macroscopic synthetic equation by the Fourier spectral method in the spatial directions, with the spatial grid size of $\Delta{x}=1/16$. The kinetic equation is solved by the backward Euler and Crank-Nicolson scheme for the first- and second-order temporal accuracy, respectively, with the initial condition
    $h(\bm{v},t=0)=2[v_{1}u_{1}(x_{1},x_{2},0)+v_{2}u_{2}(x_{1},x_{2},0)]f_{eq}$,
    while the synthetic equations are solved by the second-order Crank-Nicolson scheme. The molecular velocity space $\bm{v}$ is discretized by $6\times 6\times 6$ Gauss-Hermite quadrature, which is sufficient accurate to evaluate the high order term in Eq. (26).

    Numerical results are summarized in Fig. 4. When $\Delta{t}=1$, both GSIS and HOLO yield accurate results (not shown for clarity). However, when $\Delta{t}$ is increased to 10, HOLO becomes unstable, which is consistent with the result from the Fourier stability analysis in Fig. 1 (note that in this case the wavevector is $\theta=8\pi^{2}$). Both GSIS schemes produce stable results, and converged solution in each time step is found within 4 iterations even when the time step is as large as $\Delta{t}=50$. From the figure one can also find that when the kinetic equation is solved with second-order temporal accuracy, accurate results are obtained. However, when the kinetic equation is solved with first-order temporal accuracy, GSIS-II generates inaccurate solutions when the time step is large, say, $\Delta{t}=50$. This is because GSIS-II do not have the property of asymptotic NSF preserving at large temporal step.

    To be specific, we find that in GSIS-II, when the streaming operator is treated exactly by the Fourier spectral method, the first-order term in the Taylor expansion of velocity distribution function is given by Eq. (36).
    Suppose $\Delta{t}\sim{O(1)}$, it is estimated that the last error term has a contribution to the shear viscosity at the order of

    (47)

    $\Delta{t}\left(\frac{4\pi^{2}}{\delta_{rp}}\right)^{m+1}=\underbrace{\left(\Delta{t}\frac{16\pi^{4}}{\delta_{rp}}\right)}_{O(1)}\left(\frac{4\pi^{2}}{\delta_{rp}}\right)^{m-1}\frac{1}{\delta_{rp}},$

    where the underbraced term can be made $O(1)$ since this is largest time step to get the Navier-Stokes equation correct when solved by the second-order temporal scheme.
    Therefore, when the kinetic equation is solved by the backward Euler scheme, we have $m=1$ and the error term is at the same order with the viscosity coefficient, hence the GSIS-II is not accurate at large time step. However, when the kinetic equation is solved by the second-order temporal accuracy, we have $m=2$. Therefore, the error term is much smaller than the viscosity coefficient when the Knudsen number is small. In this case, GSIS-II can yield accurate NSF solutions even at large time steps, as observed in Fig. 4.

    ###
    7.3 Oscillatory Couette flow

    Consider the oscillatory Couette flow between two parallel plates and investigate the combined effect of spatial and temporal discretizations. The plate located at $x_{1}=0$ oscillates in the $x_{2}$ direction with the velocity
    $u_{wall}=\sin(2\pi{f_{s}}t)$,
    while the other plate at $x_{1}=1$ is stationary. The rarefaction parameter is chosen to be $\delta_{rp}=1000$, while the oscillation frequency is $2\pi{f_{s}}=0.1$. This problem can be greatly simplified. In fact, we have $\rho=\tau=0$, and the only the evolution equations for the velocity $u_{2}$ and stress $\sigma_{12}$ are needed.

    Figure 5:
    Velocity profiles in the oscillating Couette flow at the time $t=60$, when the BGK model with $\delta_{rp}=1000$ is solved by the two GSIS schemes with different spatial resolution, with the time step $\Delta{t}=1$. The reference solution is obtained from GSIS-I with $N=121$ and $\Delta{t}=0.1$. Both kinetic and synthetic equations are solved by the Crank-Nicolson scheme. Inset: the inner iteration number at each time step.

    The spatial coordinate $x_{1}\in[0,1]$ is discretized non-uniformly with $N$ points, in the following manner

    (48)

    $x_{1}=\frac{1}{2}+\frac{\tanh(8m)}{2\tanh(4)},\quad j=\frac{0,1,\cdots,N-1}{N-1}-\frac{1}{2},$

    so that the Knudsen layer can be well resolved. Later it will be shown this affects the efficiency of inner iterations. The molecular velocity space $v_{2}\times{v_{3}}$ is discretized by $6\times 6$ Gauss-Hermite quadrature, while $v_{1}$ is discretized non-uniformly:

    (49)

    $v_{1}=\frac{6}{(N_{v}-1)^{3}}(-N_{v}+1,-N_{v}+3,\cdots,{N_{v}-1})^{3},$

    where $N_{v}=32$ is the total number of discretized velocity in the $v_{1}$ direction.

    The kinetic equation is solved by the Crank-Nicolson scheme with second-order spatial and temporal accuracy, with the initial condition $h(x_{1},\bm{v},t=0)=0$, and the boundary conditions are $h(x_{1}=0,v)=2u_{wall}v_{2}f_{eq}$ for $v_{1}>0$ and $h(x_{1}=1,v)=0$ for $v_{1}<0$.
    The synthetic equations are solved by the second-order temporal accuracy, while the spatial derivative is approximated by central finite-difference with 5 stencils; hence in the solving the synthetic equations $j$ in Eq. (48) is taken from 2 to $N-3$, while the 4 left points, obtained from the kinetic equation, provide boundary conditions to synthetic equations.

    Figure 6:
    The error decay rate as a function of the rarefaction parameter in GSIS and HOLO. For wall bounded problems such as the Poiseuille and Couette flows, the perturbation from the wall bears a wavevector of $|\bm{\theta}|=\delta_{rp}$ if the Knudsen layer is resolved. Therefore, HOLO is only stable when the time step is small, while the robustness of two GSISs is manifested at any time step.

    Numerical solutions of the oscillatory flow at $\delta_{rp}=1000$ are shown in Fig. 5 for different spatial resolutions. It is seen that when the number of spatial grid is decreased from $N=121$ to 61 and then to 21, the accuracy of GSIS-I barely reduces. However, in GSIS-II when $N=21$, large difference to the reference solution is observed. This is for the GSIS-II the NSF equations can only be derived exactly when the spatial resolution is about $\Delta{x}={O(\sqrt{\text{Kn}})}\approx 0.032$. When $N=121,61$, and 21, the maximum grid size according to Eq. (48) are 0.033, 0.066, and 0.190, respectively. The case of $\max(\Delta{x})=0.190$ is certainly too large.

    On the other hand, we find that, when the same spatial resolution is used, generally speaking, GSIS-I needs less iteration numbers than GSIS-II. When $N$ is decreased, the iteration number increases. This is because this flow is driven by the oscillation of left plate; smaller value of $N$ means coarser spatial grid size and hence the information from the plate cannot be effectively passed to the bulk flow regime. Therefore, the use of non-uniform spatial grid (48) not only allows the capture of flow dynamics around the Knudsen layer, but also facilitates fast convergence.

    Note that the result of HOLO is not shown, since we find that it is not stable when $\Delta{t}>0.004$. This can be understood as follows. When the Knudsen layer is resolved, the oscillating plate will generate a perturbation with a wavelength at the order of mean free path, that is, the wavevector is about $\theta\approx 2\pi\delta$. With this information, we re-calculate the error decay rate of GSIS-I, GSIS-II and HOLO in Fig. 6 with $\theta=\delta_{rp}$. It is seen that GSIS is stable, while HOLO is unstable when $\Delta{t}>0.01$ at $\delta_{rp}=1000$. This is consistent with our numerical observations. If a coarse spatial grid is used, the stability region of HOLO increases, however, the convergence speed is much reduced.

    Finally, we assess the temporal accuracy of GSIS-I. It is seen from Fig. 7 that even when the time step is much large than the mean collision time of gas molecules, say, one sixth of the oscillation frequency, the phase of the velocity is preserved after 10 oscillation periods.

    Figure 7:
    Velocity profiles in the oscillating Couette flow at $x_{1}=0.15$ and 0.5, when the BGK model with $\delta_{rp}=1000$ is solved by GSIS-I with $\Delta{t}=5$ and 10, respectively. The reference solution is obtained from GSIS-I with $\Delta{t}=0.1$. Both kinetic and synthetic equations are solved by the Crank-Nicolson scheme, when the spatial region is discretized by Eq. (48) with $N=121$.

    ##
    8 Conclusions

    In summary, we have proposed two general synthetic iterative schemes to solve the gas kinetic equations efficiently and accurately. First, our rigorous Fourier stability analysis has shown that the GSIS enables fast convergence in the entire range of Knudsen number and time step. Second, the Chapman-Enskog expansion has been used to prove that, in the continuum flow regime, GSIS-I asymptotically preserves the Navier-Stokes-Fourier equations when the spatial and temporal step is $\Delta{t},\Delta{x}\sim{O(1)}$, provided that the NSF equations can capture the hydrodynamics at this spatial grid size and temporal step. For GSIS-II, the spatial and temporal step should be made smaller to recover the NSF equations exactly, e.g., $\Delta{t},\Delta{x}\sim{O(\sqrt{\text{Kn}})}$ when the kinetic equation is solved by a second-order accuracy temporal and spatial scheme.

    From the analytical and numerical results we conclude that, in order to construct a fast converging and asymptotic preserving scheme for rarefied gas flows, it is necessary to couple the gas kinetic equation with the macroscopic synthetic equations. And in order to construct the GSIS-I scheme, it is necessary to reach at least the Grad 13 moment equations, where the highest order velocity moments are calculated directly from the numerical solution of kinetic equation. For other kinetic systems, say, for radiative heat transfer, phonon dynamics, we need the next-level moment system beyond the equations derived from the conservation laws. We plan to test this conjecture in other kinetic systems in the near future.

    ## References

      *
    [1]

    M. L. Adams and E. W. Larsen, Fast iterative methods for
    discrete-ordinates particle transport calculations, Progr. Nucl. Energy, 40
    (2002), pp. 3–159.

      *
    [2]

    P. L. Bhatnagar, E. P. Gross, and M. Krook, A model for collision
    processes in gases. I. Small amplitude processes in charged and neutral
    one-component systems, Phys. Rev., 94 (1954), pp. 511–525.

      *
    [3]

    S. Chapman and T. Cowling, The Mathematical Theory of Non-uniform
    Gases, Cambridge University Press, 1970.

      *
    [4]

    G. Dimarco and L. Pareschi, Implicit-explicit linear multistep
    methods for stiff kinetic equations, Siam J. Numer. Anal., 55 (2017),
    pp. 664–690.

      *
    [5]

    F. Filbet and S. Jin, A class of asymptotic-preserving schemes for
    kinetic equations and related problems with stiff sources, J. Comput. Phys.,
    229 (2010), p. 7625.

      *
    [6]

    F. Filbet and S. Jin, An asymptotic preserving scheme for the
    ES-BGK model of the Boltzmann equation, J. Sci. Comput., 46 (2011),
    pp. 204–224.

      *
    [7]

    H. Grad, On the kinetic theory of rarefied gases, Comm. Pure Appl.
    Math., 2 (1949), pp. 331–407.

      *
    [8]

    Z. L. Guo, J. Q. Li, and K. Xu, On unified preserving properties of
    kinetic schemes, arXiv:1909.04923, (2019).

      *
    [9]

    Z. L. Guo, K. Xu, and R. J. Wang, Discrete unified gas kinetic
    scheme for all Knudsen number flows: Low-speed isothermal case, Physical
    Review E, 88 (2013), p. 033305.

      *
    [10]

    J. W. Hu and X. X. Zhang, On a class of implicit-explicit
    Runge-Kutta schemes for stiff kinetic equations preserving the Navier-Stokes
    limit, J. Sci. Comput., 73 (2017), pp. 797–818.

      *
    [11]

    E. W. Larsen, On numerical solutions of transport problems in the
    diffusion limit, Nuc. Sci. Eng., 83 (1983), pp. 90–99.

      *
    [12]

    L. Mieussens,  Discrete velocity models and numerical schemes for
    the Boltzmann-BGK equation in plane and axisymmetric geometries, J. Comput.
    Phys., 162 (2000), pp. 429–466.

      *
    [13]

    E. M. Shakhov, Approximate kinetic equations in rarefied gas
    theory, Fluid Dynamics, 3 (1968), pp. 112–115.

      *
    [14]

    E. M. Shakhov, Generalization of the Krook kinetic relaxation
    equation, Fluid Dyn., 3 (1968), pp. 95–96.

      *
    [15]

    Y. Sone, Kinetic theory and fluid dynamics, Birkhauser Boston,
    2002.

      *
    [16]

    H. Struchtrup, Macroscopic Transport Equations for Rarefied Gas
    Fows: Approximation Methods in Kinetic Theory, Heidelberg, Germany:
    Springer, 2005.

      *
    [17]

    W. Su, Y. H. Zhang, and L. Wu, Multiscale simulation of molecular
    gas flows by the general synthetic iterative scheme, Comput. Methods Appl.
    Mech. Engrg., 373 (2021), p. 113548.

      *
    [18]

    W. Su, L. H. Zhu, P. Wang, Y. H. Zhang, and L. Wu, Can we find
    steady-state solutions to multiscale rarefied gas flows within dozens of
    iterations?, J. Comput. Phys., 407 (2020), p. 109245.

      *
    [19]

    W. Su, L. H. Zhu, and L. Wu, Fast convergence and asymptotic
    preserving of the general synthetic iterative scheme, SIAM J. Sci. Comput.,
    42 (2020), pp. B1517–B1540.

      *
    [20]

    W. T. Taitano, D. A. Knoll, L. Chacón, J. M. Reisner, and A. K.
    Prinja, Moment-based acceleration for neutral gas kinetics with BGK
    collision operator, Journal of Computational and Theoretical Transport, 43
    (2014), pp. 83–108.

      *
    [21]

    P. Wang, M. T. Ho, L. Wu, Z. Guo, and Y. Zhang, A comparative study
    of discrete velocity methods for low-speed rarefied gas flows, Computers &
    Fluids, 161 (2018), pp. 33 – 46.

      *
    [22]

    L. Wu and X. J. Gu, On the accuracy of macroscopic equations for
    linearized rarefied gas flows, Advances in Aerodynamics, 2 (2020), p. 2.

      *
    [23]

    K. Xu and J. C. Huang, A unified gas-kinetic scheme for continuum
    and rarefied flows, J. Comput. Phys., 229 (2010), pp. 7747–7764.

      *
    [24]

    Y. Yoshio, Flows induced by temperature fields in a rarefied gas and
    their ghost effect on the behavior of a gas in the continuum limit, Annu.
    Rev. Fluid Mech., 32 (2000), pp. 779–811.

      *
    [25]

    L. H. Zhu, X. C. Pi, W. Su, Z. H. Li, Y. H. Zhang, and L. Wu, General synthetic iterative scheme for nonlinear gas kinetic simulation of
    multi-scale rarefied gas flows, J. Comput. Phys., (2021).

      *
    [26]

    Y. J. Zhu, C. W. Zhong, and K. Xu, An implicit unified gas-kinetic
    scheme for unsteady state solutions in all flow regimes, J. Comput. Phys.,
    386 (2019), pp. 190–217.