Can we find steady-state solutions to multiscale rarefied gas flows within dozens of iterations?

TL;DR

This paper introduces a general synthetic iteration scheme (GSIS) that rapidly finds steady-state solutions for multiscale rarefied gas flows across all Knudsen numbers within dozens of iterations, surpassing existing methods.

Contribution

The paper presents a novel, model-independent iterative scheme (GSIS) that significantly accelerates convergence for steady-state solutions in multiscale rarefied gas flows.

Findings

GSIS achieves convergence within dozens of iterations across all Knudsen numbers.

The method is extendable to mixture flows and chemically reacting flows.

GSIS can enhance the efficiency of Monte Carlo simulations for near-continuum flows.

Abstract

One of the central problems in the study of rarefied gas dynamics is to find the steady-state solution of the Boltzmann equation quickly. When the Knudsen number is large, i.e. the system is highly rarefied, the conventional iteration scheme can lead to convergence within a few iterations. However, when the Knudsen number is small, i.e. the flow falls in the near-continuum regime, hundreds of thousands iterations are needed, and yet the "converged" solutions are prone to be contaminated by accumulated error and large numerical dissipation. Recently, based on the gas kinetic models, the implicit unified gas kinetic scheme (UGKS) and its variants have significantly reduced the iterations in the near-continuum flow regime, but still much higher than that of the highly rarefied gas flows. In this paper, we put forward a general synthetic iteration scheme (GSIS) to find the steady-state solutions of general rarefied gas flows within dozens of iterations at any Knudsen number. As the GSIS does not rely on the specific kinetic model/collision operator, it can be naturally extended to quickly find converged solutions for mixture flows and even flows involving chemical reactions. These two superior advantages are also expected to accelerate the slow convergence in simulation of near-continuum flows via the direct simulation Monte Carlo method and its low-variance version.

Tables

Table 1. Table 1: Number of iteration steps and CPU time to reach convergence for the lid-driven cavity flow.

$\delta$	$N^2$	$N_{v_1}{N}_{v_2}{N}_{v_3}$	Iteration steps		Total CPU time (s)
			CIS	GSIS	CIS	GSIS
0.1	$20\times 20$	$48\times 48\times 24$	14	13	28.5	32.6
1	$20\times 20$	$48\times 48\times 24$	14	16	28.2	38.4
10	$20\times 20$	$24\times 24\times 12$	99	31	121.4	47.7
100	$40\times 40$	$16\times 16\times 16$	1823	36	3176.1	144.1
1000	$60\times 60$	$8\times 8\times 8$	—	36	—	492.5

Equations

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

$$f(t,\mathbf{x},\mathbf{v})=f_{eq}(\mathbf{v})+\alpha h(t,\mathbf{x},\mathbf{v}),$$

$$\frac{\partial h}{\partial t}+\mathbf{v}\cdot \frac{\partial h}{\partial \mathbf{x}}=L(h,f_{eq}),$$

$$L=\underset{L^{+}}{\underbrace{\iint B(\theta ,|\mathbf{u}|)[f_{eq}({\mathbf{v}}^{\prime })h({\mathbf{v}}_{\ast }^{\prime })+f_{eq}({\mathbf{v}}_{\ast }^{\prime })h({\mathbf{v}}^{\prime })-f_{eq}(\mathbf{v})h({\mathbf{v}}_{\ast })]\mathrm{d}\Omega \mathrm{d}{\mathbf{v}}_{\ast }}}-{\nu }_{eq}(\mathbf{v})h(\mathbf{v}),$$

$${\nu }_{eq}(\mathbf{v})=\iint B(|\mathbf{u}|,\theta )f_{eq}({\mathbf{v}}_{\ast })\mathrm{d}\Omega \mathrm{d}{\mathbf{v}}_{\ast }.$$

$$B(|\mathbf{u}|,\theta )=\frac{|\mathbf{u}{|}^{2(1-\omega )}}{K}{\sin }^{\frac{1}{2}-\omega }\left(\frac{\theta }{2}\right){\cos }^{\frac{1}{2}-\omega }\left(\frac{\theta }{2}\right),$$

$$f(t,\mathbf{x},\mathbf{v})=\frac{2{\int }_{v_n^{\prime }<0}|v_n^{\prime }|f(t,\mathbf{x},{\mathbf{v}}^{\prime })\mathrm{d}{\mathbf{v}}^{\prime }}{\pi T_w^2}\exp \left(-\frac{|\mathbf{v}-{\mathbf{U}}_w{|}^2}{T_w}\right),$$

$$\rho =\int h\mathrm{d}\mathbf{v},\mathbf{U}=\int \mathbf{v}h\mathrm{d}\mathbf{v},T=\frac{2}{3}\int |\mathbf{v}{|}^{2}h\mathrm{d}\mathbf{v}-\rho ,p=\rho +T$$

$${\sigma }_{ij}=2\int \left(v_i{v}_j-\frac{|\mathbf{v}{|}^2}{3}{\delta }_{ij}\right)h\mathrm{d}\mathbf{v},\mathbf{q}=\int \mathbf{v}|\mathbf{v}{|}^{2}h\mathrm{d}\mathbf{v}-\frac{5}{2}\mathbf{U},$$

$${\nu }_{eq}^{(k)}{h}^{(k+1)}+\mathbf{v}\cdot \frac{\partial h^{(k+1)}}{\partial \mathbf{x}}=L^{+}(h^{(k)},f_{eq}),$$

$${\delta }_{rp}=\frac{H}{\lambda },\lambda =\frac{\mu (T_0)v_m}{n_0{k}_B{T}_0},$$

$$\mathrm{Kn}=\frac{\sqrt{\pi }}{2{\delta }_{rp}}.$$

$$\frac{\partial \rho }{\partial t}+\frac{\partial U_i}{\partial x_i}=0,$$

$$2\frac{\partial U_i}{\partial t}+\frac{\partial \rho }{\partial x_i}+\frac{\partial T}{\partial x_i}+\frac{\partial {\sigma }_{ij}}{\partial x_j}=0,$$

$$\frac{3}{2}\frac{\partial T}{\partial t}+\frac{\partial q_j}{\partial x_j}+\frac{\partial U_j}{\partial x_j}=0,$$

$${\sigma }_{ij}=-{\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)\equiv -2{\delta }_{rp}^{-1}\frac{\partial U_{<i}}{\partial x_{j>}},q_i=-\frac{5}{4\mathrm{Pr}}{\delta }_{rp}^{-1}\frac{\partial T}{\partial x_i},$$

$${\sigma }_{ij}=-2{\delta }_{rp}^{-1}\frac{\partial U_{<i}}{\partial x_{j>}}+{\text{HoT}}_{{\sigma }_{ij}},$$

$$q_i=-\frac{5}{4\mathrm{Pr}}{\delta }_{rp}^{-1}\frac{\partial T}{\partial x_i}+{\text{HoT}}_{q_i},$$

$$\frac{\partial {\sigma }_{ij}}{\partial t}+{\text{HoT}}_{{\sigma }_{ij}}+\underline{2\frac{\partial U_{<i}}{\partial x_{j>}}=-{\delta }_{rp}{\sigma }_{ij}}+2\int (L-L_s)v_i{v}_j\mathrm{d}\mathbf{v},$$

$$L_s={\delta }_{rp}\left\{\left[\rho +2\mathbf{U}\cdot \mathbf{v}+T\left(|\mathbf{v}{|}^2-\frac{3}{2}\right)+\frac{4\left(1-\mathrm{Pr}\right)}{5}\mathbf{q}\cdot \mathbf{v}\left(|\mathbf{v}{|}^2-\frac{5}{2}\right)\right]f_{eq}-h\right\}$$

$$\text{HoT}_{\sigma_{ij}}=\left\{\begin{array}[]{lr}&\frac{\partial}{\partial x_{i}}\int{}(2v_{i}^{2}-1)v_{j}h\mathrm{d}\bm{v}+\frac{\partial}{\partial x_{j}}\int{}(2v_{j}^{2}-1)v_{i}h\mathrm{d}\bm{v}+\frac{\partial}{\partial x_{k}}\int{}2v_{1}v_{2}v_{3}h\mathrm{d}\bm{v},\\ &\text{for~{}}i\neq{j},k\neq{i},k\neq{j},\\ &\frac{\partial}{\partial x_{i}}\int{}2(v_{i}^{2}-\frac{|\bm{v}|^{2}}{3}-\frac{2}{3})v_{i}h\mathrm{d}\bm{v}+\sum_{k}\frac{\partial}{\partial x_{k}}\int{}2(v_{i}^{2}-\frac{|\bm{v}|^{2}}{3}+\frac{1}{3})v_{k}h\mathrm{d}\bm{v},\\ &\text{for~{}}i=j,k\neq{i}.\end{array}\right.$$

$$\frac{\partial q_i}{\partial t}+{\text{HoT}}_{q_i}+\underline{\frac{3C_q}{2}\frac{\partial T}{\partial x_i}=-\frac{2}{3}{\delta }_{rp}{q}_i}+\int (L-L_s)v_i|\mathbf{v}{|}^2\mathrm{d}\mathbf{v},$$

$${\text{HoT}}_{q_i}=\frac{\partial }{\partial x_i}\int \left[(v_i^2-C_q)\left(|\mathbf{v}{|}^2-\frac{3}{2}\right)-v_i^2\right]h\mathrm{d}\mathbf{v}+\sum _{j\ne i}\frac{\partial }{\partial x_j}\int v_i{v}_j\left(|\mathbf{v}{|}^2-\frac{5}{2}\right)h\mathrm{d}\mathbf{v},$$

$${\nu }_{eq}^{(k)}{h}^{(k+1/2)}+\mathbf{v}\cdot \frac{\partial h^{(k+1/2)}}{\partial \mathbf{x}}=L^{+}(h^{(k)},f_{eq}),$$

$$h^{(k+1)}(\mathbf{x},\mathbf{v})=$$

$$+\left[{\lambda }_{\rho }(\mathbf{x})+{\lambda }_T(\mathbf{x})\left(|\mathbf{v}{|}^2-\frac{3}{2}\right)+{\lambda }_{{\sigma }_{ij}}(\mathbf{x})\left(v_i{v}_j-\frac{|\mathbf{v}{|}^2}{3}{\delta }_{ij}\right)\right]f_{eq},$$

$$h(x_2=0,\mathbf{v})=\left[1-\frac{|\mathbf{v}{|}^2}{2}-2\sqrt{\pi }{\int }_{v_2<0}{v}_{2}h(x_2=0,\mathbf{v})\mathrm{d}{v}_2\right]f_{eq},\text{when }{v}_2>0,$$

$$h(v_1,v_2,v_3)=-h(v_1,-v_2,v_3),$$

$$\mathbf{U}=0,{\sigma }_{ij}=0\text{when}i\ne j,q_1=q_3=0,$$

$$\frac{\partial T}{\partial x_2}=-\frac{4{\delta }_{rp}}{9C_q}{q}_2+\underset{H_1^{(k)}(x_2)}{\underbrace{\frac{2}{3C_q}\int v_2|\mathbf{v}{|}^2(L-L_s)\mathrm{d}\mathbf{v}}}-\underset{H_2^{(k+1/2)}(x_2)}{\underbrace{\frac{2}{3C_q}\frac{\partial }{\partial x_2}\int (v_2^2-C_q)\left(|\mathbf{v}{|}^2-\frac{3}{2}\right)h\mathrm{d}\mathbf{v}}},$$

$$T^{(k+1)}(x_2)=-\frac{4{\delta }_{rp}{q}_2}{9C_q}\left(x_2-\frac{1}{2}\right)+{\int }_{1/2}^{x_2}{H}_1^{(k)}(x_2)\mathrm{d}{x}_2-H_2^{(k+1/2)}(x_2),$$

$$q_2=\frac{9C_q}{2{\delta }_{rp}}\left[T^{(k+1/2)}(x_2=0)+H_2^{(k+1/2)}(x_2=0)-H_1^{(k)}(x_2=0)\right].$$

$$\rho +T+{\sigma }_{22}=\int 2v_2^{2}hd\mathbf{v},$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Multiscale rarefied gas flows involving a wide range of Knudsen number have been encountered in massive engineering problems, e.g. high-altitude aerothermodynamics of space vehicles, microelectromechanical systems, and gas transportation in ultra-tight shale strata. A gas flow can be described by either the macroscopic or the microscopic model. At the macroscopic level, the gas is assumed as a continuous medium and the evolution of gas system is described in terms of the spatial and temporal variations of the familiar flow properties such as density, velocity, pressure and temperature. The mathematical description of any macroscopic model is grounded in two primary aspects: 1) the conservation laws that describe how the mass, momentum and energy must be conserved during transport processes; 2) the constitutive equations that describe how the fluxes of mass dissipation, momentum diffusion and heat conduction response to various stimuli such as pressure difference, gradients of temperature and velocity, and external force. The Navier-Stokes equations provide the conventional mathematical model for a gas as a continuum, in which the conservation laws are closed by the famous constitutive equations of the Newton’s law of viscosity and Fourier’s law of heat conduction. Since the transport terms are expressed in terms of the first-order macroscopic quantities, the Navier-Stokes equations are only valid when the length scale of the gradients of the macroscopic variables is much larger than the mean free path of gas molecules, i.e. the Knudsen number is far smaller than one [1].

The microscopic model postulates that the gas is not continuous but is composed of a finite number of molecules. The molecules rush hither and thither at large speed, and strike with boundary and collide with each other. Actually, the macroscopic transport phenomena stem no other than the random motions of the gas molecules. The mathematical model at the microscopic level is the Boltzmann equation, which governs the evolution of the one-particle velocity distribution function providing information on the state of every molecule at all times [2]. Then, the macroscopic flow properties are identified with average values of the molecular quantities. Note that the Boltzmann equation is applicable for the entire range of Knudsen number.

The Boltzmann equation can be numerically solved either in discretized molecular velocity space via the discrete velocity method (DVM)[3], or by applying the direct simulation Monte Carlo (DSMC) method that uses a collection of particles to represent random points in the molecular velocity space [2]. Compared to the traditional computational fluid dynamic (CFD) techniques for solving macroscopic equations, the Boltzmann equation (or simplified kinetic model equations) is much more expensive to solve in terms of computation time and memory. This is mainly due to the following facts. First of all, additional dimensions of the molecular velocity space are required to be discretized in DVM and particles are required to generate in DSMC. Second, since the random behaviors of gas molecules are modeled on length and time scales comparable to the cell size and simulation time interval, respectively, in order to suppress numerical diffusion errors it is suggested that the size of grid cell and the time interval should be smaller than the molecule’s mean free path and the mean collision time, respectively [4]. As a consequence, the computational cost dramatically increases as the gas flow approaches the near-continuum flow regime. Finally, in DVM, the conventional iteration scheme (CIS) to find steady-state solution converges extremely slowly for flows at low Knudsen numbers, since the exchange of information (e.g. perturbation in the flow field) through molecular steaming becomes very inefficient when binary collisions dominate [5]. Worse still, the “converged” solutions are prone to be contaminated by numerical errors, e.g. the accumulated error from finite discrete molecular velocities [6] and error stemming from the evaluation of molecular collisions, say, by the projection method [7] and the fast spectral method [8]. In DSMC, the simulation time also increases significantly due to this inefficient information exchange process in the near-continuum flow. Note that unified gas-kinetic scheme (UGKS) [9, 10, 11, 12, 13] can remove the restrictions on cell size and time step by simultaneously handling free streaming and collision of gas molecules during transport processes. However, as information exchanging relays on the evolution of velocity distribution function, UGKS still needs a large number of iterations to obtained steady-state solutions in near-continuum flows [14, 15].

There has been a tremendous growth of researches on multiscale hybrid numerical methods that combine multiple models defined at fundamentally different length and time scales within the same overall spatial and temporal domain. Specifically for the flow of interest, the continuum CFD methods are used in regions where the Navier-Stokes equations are valid, while methods based on gas kinetic theory are applied in regions where the continuum equations fail [16, 17, 18, 19, 20, 21, 22]. However, intrinsic difficulties arises when coupling the two different models. First, the mechanism for continuum breakdown is unclear and the criterion to determine where the continuum model is valid relies on empirical parameters that varies for different flow conditions [20]. Second, the continuum-kinetic coupling is strictly required to lie in the region that can be accurately modeled by Navier-Stokes equations, so that the Boltzmann equation is still employed in low-Knudsen-number regions. Therefore, the CIS still needs lots of iterations to achieve convergence; also, DSMC still needs small cell size and time step and hence large evolution steps to find the converged solutions.

In recent years, the synthetic iterative scheme (SIS), which is initially developed for the radiation transport processes [23], has been extended to achieve high efficiency and accuracy in DVM, in particular with fast convergence property across the whole gas flow regimes [24, 25]. In this scheme, the gas kinetic equations and macroscopic equations are solved simultaneously on the same grids in the entire domain. Since the velocity distribution function is guided by the macroscopic flow quantities solved from diffusion-type equations at each iterative step, information propagates accurately and fast even when Knudsen number is small. When the Knudsen number is small, the synthetic macroscopic equations reduce to the Navier-Stokes equations. However, the macroscopic equation contains high-order terms to take into account rarefaction effects, thus the SIS also preserves accuracy in high Knudsen number regimes. The SIS has been successfully applied to Poiseuille flow in channels of arbitrary shapes using the Bhatnagar-Gross-Krook kinetic model for single-species gases [26], and flows of binary and ternary gas mixtures driven by local pressure, temperature and concentration gradients using the McCormak model [27, 28, 29, 30, 31]. The SIS has also been extended to solve the linearized Boltzmann equation, where the role of realistic intermolecular potentials in Poiseuille, Couette and thermal transpiration flows has been analyzed [5, 32].

It is interesting to note that the similar idea of SIS has also been used in DSMC, that is, in addition to traditional DSMC, macroscopic variables are solved and updated according to macroscopic rules/equations. For instances, in the information preservation (IP)-DSMC, the information velocity is introduced to compute macroscopic velocity and shear stress, with the aim of removing “the statistical fluctuation source inherent in the DSMC method that results from the randomness of the thermal velocity” [33, 34, 35], although the rule of updating the information velocity and/or other macroscopic variables is not exactly derived from the Boltzmann equation. On the other hand, the moment guided DSMC is also proposed to reduce the statistical error, where the density, velocity and temperature are updated by five exact macroscopic equations from the conservation law, but with the pressure tensor and heat flux calculated from the DSMC [36].

In DVM, the SIS can not only asymptotically achieve the Navier-Stokes limit with fast convergence rate, but also preserve accuracy in high Knudsen number regimes. The critical point to develop this scheme is that the macroscopic equations must explicitly contain both the constitutive relations predicting the transport phenomena at the continuum level, and high-order terms taking into account rarefaction effects. To the author’s awareness, the SIS is still limited to simple rarefied gas flows such as the Poiseuille, Couette and thermal transpiration flows, where the flow velocity is perpendicular to the computational domain, we refer to [24] for example. In this paper, we intend to put forward the general SIS (GSIS) with the aim to find the steady-state solutions of general rarefied gas flows within dozens of iterations at any Knudsen number. For simplicity we considered linearized flows but the methodology can be extended to nonlinear flows easily.

The remainder of the paper is organized as follows. In Section 2, the linearized Boltzmann equation (LBE) is introduced. In Section 3, the GSIS for general rarefied gas flow is proposed. Numerical tests to assess the efficiency and accuracy of the proposed scheme are presented for stationary problems in Section 4 and for periodic oscillatory problems in Section 5. The paper closes with some finial comments in Section 6.

2 The linearized Boltzmann equation

In kinetic theory, the state of a gas system is described by the one-particle velocity distribution $f\left(t,\bm{x},\bm{v}\right)$. Evolution of the velocity distribution function to the independent variables, i.e. time $t$, spatial position $\bm{x}=\left(x_{1},x_{2},x_{3}\right)$ and molecular velocity $\bm{v}=\left(v_{1},v_{2},v_{3}\right)$, is governed by the Boltzmann equation [1]. When the system deviates slightly from the global equilibrium described by

[TABLE]

the velocity distribution function of gas molecules can be linearized around $f_{eq}$ as:

[TABLE]

where $\alpha h(t,\bm{x},\bm{v})$ is the small perturbance satisfying $|\alpha h/f_{eq}|\ll 1$ with $\alpha$ being a small constant related to the amplitude of perturbation. The velocity distribution function $h(t,\bm{x},\bm{v})$, however, is not necessary smaller than the equilibrium distribution function $f_{eq}$. The LBE for $h(t,\bm{x},\bm{v})$ is:

[TABLE]

where the linearized Boltzmann collision operator is [37]:

[TABLE]

and the equilibrium collision frequency is

[TABLE]

Note that the relative velocity of the two molecules before binary collision is $\bm{u}=\bm{v}-\bm{v}_{\ast}$, and $\Omega$ is a unit vector along the relative post-collision velocity $\bm{v}^{\prime}-\bm{v}^{\prime}_{\ast}$. The deflection angle $\theta$ between the pre- and post-collision relative velocities satisfies $\cos\theta=\Omega\cdot{\bm{u}}/|\bm{u}|$, $0\leq\theta\leq\pi$. Finally, $B(\theta,|\bm{u}|)=|\bm{u}|\sigma$ is the collision kernel, with $\sigma$ being the differential cross-section that is determined by the intermolecular potential. In the present paper, we consider the inverse power-law potentials, where the collision kernels are modeled as [38, 37]

[TABLE]

with $\omega$ being the viscosity index (i.e. the shear viscosity $\mu$ of the gas is proportional to $T^{\omega}$) and $K$ some normalization constants [37]. HS and Maxwell molecules have $\omega=0.5$ and 1, respectively. The details of implementation of the Lennard-Jones potentials in fast spectral method can be found in Ref. [39], but here only the collision kernel (6) will be used to demonstrated that the GSIS works for the LBE.

Note that we present the governing system in terms of dimensionless variables. The coordinate $\bm{x}$ is normalized by the characteristic flow length $H$, the molecular velocity $\mathbf{v}$ is normalized by the most probable speed ${v_{m}}=\sqrt{2k_{B}T_{0}/m}$, the time $t$ is normalized by $H/v_{m}$, and velocity distribution functions $f_{eq}$ and $h$ are normalized by $n_{0}/v_{m}^{3}$, where $n_{0}$ is the average number density of the gas molecules, $T_{0}$ is the reference temperature, $k_{B}$ is the Boltzmann constant, and $m$ is the mass of gas molecules.

To fully determine the gas dynamics in spatially-inhomogeneous problems, the gas-surface boundary condition should be specified. In this paper, the Maxwell diffuse boundary condition will be used: the velocity distribution function $f(t,\bm{x},\bm{v})$ of the reflected gas molecules at the solid surface satisfies the following equation:

[TABLE]

where $T_{w}$ is the wall temperature normalized by the reference temperature $T_{0}$, $\bm{U_{w}}$ is the wall velocity normalized by the most probable speed $v_{m}$, and $v_{n}$ is the normal component of the peculiar velocity $\bm{v}-\bm{U}_{w}$ redirected into the gas.

The macroscopic quantities of interest including the number density $\rho$, bulk velocity $\bm{U}$, temperature $T$, pressure $p$, stress tensor $\sigma_{ij}$ and heat flux $\bm{q}$, which are further normalized by the dimensionless constant $\alpha$, can be calculated as

[TABLE]

where $\delta$ is the Kronecker delta function, and $i,j=1,2,3$ represent the three orthogonal spatial directions in the Cartesian coordinates.

3 The general synthetic iteration scheme

The steady state solution of the integro-differential system (3) is usually solved by the CIS. Given the value of $h^{(k)}(\bm{x},\bm{v})$ at the $k$-th iteration step, the velocity distribution function at the next iteration step is calculated by solving the following equation [40, 38, 39]:

[TABLE]

where the derivative with respect to x can be approximated by any conventional CFD schemes such as the finite difference, finite volume, or Discontinuous Galerkin (DG) methods [41, 42], and the collision operator in Eq. (4) can be calculated by the fast spectral method [37, 39] based on the velocity distribution function at the $k$-th iteration step. The process is repeated until relative differences between successive estimates of macroscopic quantities are less than a convergence criterion $\epsilon$.

A key parameter in the rarefied gas flow is the rarefaction parameter, which is defined as

[TABLE]

where $\mu(T_{0})$ is the shear viscosity of the gas at the reference temperature, and $\lambda$ is the mean free path of the gas molecules. Alternatively, the Knudsen number is defined as

[TABLE]

The CIS is efficient for highly rarefied gas flows when $\delta_{rp}$ is very small, where converged solutions can be quickly found after several iterations. However, the number of iteration increases significantly with the rarefaction parameter [43, 5]. This is due to the frequent collisions of gas molecules, which quickly smear the perturbance and hinder the fluid information exchange. In order to enhance the information exchange across the whole computational domain, synthetic equations for the evolution of macroscopic flow variables that are asymptotic preserving the Navier-Stokes limit should be developed [5].

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

[TABLE]

which are not closed, since expressions for the shear stress $\sigma_{ij}$ and heat flux $\bm{q}$ are not known. One way to close Eq. (13) is to use the Chapman-Enskog expansion, where the distribution function is expressed in the power series of $\mathrm{Kn}$ [1]: $h=\mathrm{Kn}h^{(1)}+\mathrm{Kn}^{2}h^{(2)}+\cdots$. When $f=f^{(0)}$, we have $\sigma_{ij}=q_{i}=0$, and Euler equations are recovered. When the distribution function is truncated at the first-order of $\mathrm{Kn}$, that is, $h=\mathrm{Kn}h^{(1)}$, we have

[TABLE]

and Eq. (13) reduces to Navier-Stokes equations with $\mathrm{Pr}$ being the Prandtl number. Higher-order macroscopic equations can be obtained successively but they are not stable. On the other hand, even the obtained high-order macroscopic equations are stable, they are only the approximate solutions of the Boltzmann equation, rather than the exact solutions.

It should be noted that in implicit UGKS [44] and other variants [45, 46], both the gas kinetic equation and macroscopic equations (13) are solved, where $\sigma_{ij}$ and $\bm{q}$ are obtained according to Eq. (9). These methods are efficient when the Knudsen number is large, like the CIS. However, in the near-continuum flow regime, the number of iterations are still large, at the order of thousands iterations. The reason for the relative slow convergence is that, if the iteration starts from the global equilibrium state where $\sigma_{ij}$ and $\bm{q}$ are zero, in most of the time the Euler equations, rather than the Navier-Stokes equations that dominates the steady-state flow dynamics, are solved, due to the fact that perturbance from the wall boundary takes a long time to reach the bulk region for near-continuum flows. Even when the shear stress and heat flux are non-zero, solutions of Eq. (13) deviate from that of the Navier-Stokes equations in the near-continuum flow regime unless they nearly converge to the steady-state solutions. As a matter of fact, the authors have checked, in the linearized Poiseuille flow, that Eq. (13) cannot boost convergence [5];

Bearing this in mind, to develop an ultra-fast convergence scheme, the macroscopic equations must reduce to the Navier-Stokes equation in the near-continuum flow regime, and must contain the Newton’s law for stress and Fourier’s law for heat conduction explicitly to recover the macroscopic transport mechanism; that is, the shear stress and heat flux should be expressed as follows:

[TABLE]

where $\text{HoT}_{\sigma_{ij}}$ and $\text{HoT}_{q_{i}}$ are the high-order terms containing contributions of all the orders $O(Kn^{\alpha})$ with $\alpha=2,3,\cdots,\infty$.

To obtain (15), we multiply Eq. (3) by $2(v_{i}v_{j}-\delta_{ij}|\bm{v}|^{2}/3)$ and integrate the resultant equation with respect to $\bm{v}$, and obtain

[TABLE]

where

[TABLE]

is the linearized collision operator of the Shakhov kinetic model equation [47], and

[TABLE]

Note that this derivation is rather simple as we just separate the underlined term in Eq. (17) from high-order moments $\int{}2(v_{i}v_{j}-\delta_{ij}|\bm{v}|^{2}/3)v_{k}hd\mathbf{v}$, and the purpose of introducing $L_{s}$ is only to recover the term $\delta_{rp}\sigma_{ij}$, so that the Newton’s law of stress is recovered explicitly. It should also be noted that, for the linearized Boltzmann collision operator, the term $2\int{(L-L_{s})v_{i}v_{j}}d\mathbf{v}$ is negligible small when compared to $\delta_{rp}\sigma_{ij}$. For instances, for the Maxwell model, this term is zero, while for the HS molecular model, this term is less than 2% of $\delta_{rp}\sigma_{ij}$, see page no. 169 in the third edition of the book [1].

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

[TABLE]

where

[TABLE]

and for the linearized Boltzmann collision operator, the term $\int{(L-L_{s})v_{i}|\bm{v}|^{2}}d\bm{v}$ is negligible small when compared to $\delta_{rp}q_{i}$, i.e. within 3% of $\delta_{rp}q_{i}$ [1]. If we choose $C_{q}=5/9\mathrm{Pr}$, then the under-braced term in Eq. (20) recovers the Fourier’s heat conduction law in Eq. (14). Since for monatomic gas the Prandtl number is very close to $2/3$, in the following paper we choose $C_{q}=5/6$.

Note that the macroscopic equations  (13), (17) and (20) resemble the Grad 13 moment equations [48, 49]. However, since the higher-order terms (19) and (21) are computed directly from the velocity distribution function, no approximations are introduced here. If the velocity distribution function is approximated by the Gauss-Hermite polynomials to the third order, where the coefficients before those polynomials are determined by the first 13 moments of the velocity distribution function, then G13 moment equations will be recovered. Since the first-order Chapman-Enskog expansion to G13 equations leads to Eqs. (13) and (14), that is, only the underlined terms in Eqs. (17) and (20) are retained, the derived synthetic equations (13), (17) and (20) are asymptotic preserving the Navier-Stokes limit. Thus, they should be able to boost the convergence to the steady-state solution of the LBE significantly, as in the bulk region (a few mean free path of gas molecules away from solid surfaces) we are effectively solving the Navier-Stokes equations.

With these macroscopic equations to update the macroscopic quantities and the velocity distribution function, we devise the following iteration scheme to find the steady-state solution of the LBE (3) efficiently:

Step 1. When the velocity distribution function $h^{(k)}$ and the corresponding macroscopic quantities in Eqs. (8) and (9) are known at the $k$-th iteration, we calculate $2\int{(L-L_{s})v_{i}v_{j}}\mathrm{d}\bm{v}$ in Eq. (17) and $\int{(L-L_{s})v_{i}|\bm{v}|^{2}}\mathrm{d}\bm{v}$ in Eq. (20). We also calculate the velocity distribution function $h^{(k+1/2)}$ according to the conventional iteration scheme (10), that is, we solve the following equation:

[TABLE]

by a second-order upwind finite difference in the bulk and a first-order upwind scheme at the solid surface [40] or the DG method [41, 42].

• 2.

Step 2. From $h^{(k+1/2)}$, we calculate the density $\rho^{(k+1/2)}(\bm{x})$, flow velocity $\bm{U}^{(k+1/2)}(\bm{x})$, the temperature $T^{(k+1/2)}(\bm{x})$, the shear stress $\sigma_{ij}^{(k+1/2)}(\bm{x})$, the heat flux $\bm{q}^{(k+1/2)}(\bm{x})$, and the high-order terms $\text{HoT}_{\sigma_{ij}}$ and $\text{HoT}_{q_{i}}$ defined in Eqs. (19) and (21), respectively.

• 3.

Step 3. We obtain the macroscopic quantities $\rho^{(k+1)}(\bm{x})$, $\bm{U}^{(k+1)}(\bm{x})$, $T^{(k+1)}(\bm{x})$, $\sigma_{ij}^{(k+1)}(\bm{x})$, and $\bm{q}^{(k+1)}(\bm{x})$ by solving the synthetic equations (13), (17) and (20), That is, for the steady-state problems the shear stress and heat flux can be solved from Eq. (17) and (20), which will then be substituted to Eq. (13) to form the Navier-Stokes equations with source terms related to the higher-order terms defined in Eqs. (19) and (21). These equation can be solved by the SIMPLE algorithm and/or DG method easily in the bulk region, where the boundary values in the vicinity of the walls for the density, velocity, temperature are obtained from the step 2. The detailed DG algorithm to solve the synthetic equations can be found in the Appendix.

• 4.

Step 4. The velocity distribution function $h$ is modified to incorporate the change of macroscopic quantities. That is,

[TABLE]

where $\lambda_{\bm{U}}(\bm{x})=\bm{U}^{(k+1)}(\bm{x})-\bm{U}^{(k+1/2)}(\bm{x})$, $\lambda_{\bm{q}}(\bm{x})=\bm{q}^{(k+1)}(\bm{x})-\bm{q}^{(k+1/2)}(\bm{x})$, $\lambda_{\rho}(\bm{x})=\rho^{(k+1)}(\bm{x})-\rho^{(k+1/2)}(\bm{x})$, $\lambda_{T}(\bm{x})=T^{(k+1)}(\bm{x})-T^{(k+1/2)}(\bm{x})$, and $\lambda_{\sigma_{ij}}(\bm{x})=B\sigma_{ij}^{(k+1)}(\bm{x})-B\sigma_{ij}^{(k+1/2)}(\bm{x})$, with $B=3/2$ when $i=j$ and $B=2$ otherwise.

• 5.

Step 5. The above steps are repeated until convergence.

Since the gas kinetic equation is solved together with the macroscopic equations (13), (17) and (20) for general rarefied gas flows, the above scheme is called the GSIS. Note that although the SIS has been widely applied to the radiation transport processes [23] and rarefied gas flows driven by local pressure, temperature, and concentration gradients [43, 29, 50, 51, 41] to overcome the slow convergence in the near-continuum flow regime, it is the first time that the GSIS is developed for general rarefied gas flows described by the LBE. Also, it is with no doubt that such a methodology can be directly applied to construct the GSIS for the nonlinear Boltzmann equation.

4 Numerical results for stationary problems

Numerical simulations are carried out to assess the efficiency and accuracy of the GSIS. To this end, we consider the one-dimensional heat transfer between two parallel plates, two-dimensional lid-driven cavity flow and shear-driven flow between two eccentric cylinders. The reason is that in previous cases the special SIS is only applicable for rarefied gas flows [43, 29, 50, 51, 41, 52], where the flow velocity is perpendicular to the computational domain. Here we investigate the performance of the GSIS for typical general rarefied gas flows, where the flow velocity (or other macroscopic variables) varies within the computational domain.

4.1 Heat transfer between two parallel plates

Consider the steady Fourier flow of a gas between two infinite parallel plates with a distance $H$, located at $x_{2}=0$ and $x_{2}=1$. The two plates are stationary, the one at $x_{2}=0$ has a temperature $T_{0}-\Delta{T}/2$, while that at $x_{2}=H$ has a temperature $T_{0}+\Delta{T}/2$. We assume that the temperature difference $\Delta{T}$ is negligible compared to $T_{0}$, so that the problem is symmetrical around $x_{2}=1/2$. Therefore, in numerical simulations only the region $x_{2}\in[0,1/2]$ is considered. The Boltzmann equation is linearized by choosing $\alpha=\Delta{T}/{T_{0}}$ in Eq. (2). The boundary condition at $x_{2}=0$, as according to Eqs. (2) and (7), is

[TABLE]

while that at $x_{2}=0.5$ is

[TABLE]

due to the symmetry of this linearized problem.

From the synthetic equations (13), (17) and (20), as well as the symmetry condition (25), we know

[TABLE]

the heat flux perpendicular to the two plates $q_{2}$ is a constant, and the variation of the perturbed temperature satisfies

[TABLE]

whose solution at the $(k+1)$-th iteration step is given by

[TABLE]

where the constant heat flux $q_{2}$ is

[TABLE]

When the temperature is known, the density variation can be easily obtained by solving the following equation

[TABLE]

where the term at the right-hand-side of Eq. (30) is zero due to the symmetry condition (25), and according to Eq. (17) the stress $\sigma_{22}$ can be calculated as

[TABLE]

We first test the efficiency of the GSIS based on the Shakhov model, that is, in Eq. (3) we let the linearized Boltzmann collision operator equal to that of the linearized Shakhov model (18). We choose the rarefaction parameter $\delta_{rp}=50$ and discretize the half spatial space into $N_{2}$ even-spaced points, where the derivative with respect to $x_{2}$ is approximated by a second-order upwind finite difference. The molecular velocity space in the $v_{1}$ and $v_{3}$ directions is truncated to the region $[-6,6]$ by $24\times 24$ equidistant points, while the molecular velocity $v_{2}$ is truncated to $[-6,6]$ and approximated by the non-uniform points [37, 53]:

[TABLE]

which is useful to capture the discontinuity in the velocity distribution function near $v_{2}\sim 0$. In this test we take $\imath=3$ and $N_{v}=64$. The iterations in both CIS and GSIS are terminated when

[TABLE]

is less than a certain value. Note that since $\rho$ and $T$ at $x_{2}=1/2$ are excluded in the above equation since they are zero.

Figure 1 compares the convergence history of the GSIS and CIS when the rarefaction parameter is $\delta_{rp}=50$, that is, the flow is in the near-continuum regime. Starting from the initial guess $h(x_{2},\bm{v})=0$, the perturbance from the solid surface quickly changes the density and temperature near the solid surface in the CIS (about one molecular mean free path away from the wall). However, due to the frequent collision between gas molecules, those in the bulk region takes a long time (i.e. iteration steps) to feel this change. From example, from Fig. 1(b) we see that it takes about 50 iteration steps for the temperature at $x_{2}=0.5$ to feel this change. Moreover, such a change does not necessary lead to the final converged state monotonically, but it could be deviate further away from the final steady state: from Fig. 1(a) we see that the density perturbance in the bulk region is even negative after 50 iterations, while the final steady state the density is always non-negative in the region of $x_{2}\in[0,0.5]$. This is also evidenced in Fig. 2 that the error does not decay monotonically but oscillates several times. Such a slow convergence is completely changed in the GSIS, where the temperature and density are corrected according to the synthetic equations (27) and (30), which the dominated parts are respectively $\frac{\partial T}{\partial x_{2}}=-\frac{4\delta_{rp}}{9C_{q}}q_{2}$ and $\rho=-T$ when $\delta_{rp}$ is large. This means that the temperature and density in the bulk region are corrected to be nearly linear immediately. As we can see from Fig. 1(d), after the first iteration, the temperature from the GSIS at $x_{2}=0$ is the same as that from the CIS, but the temperature from the GSIS in the bulk region varies linearly, while that from the CIS is still zero. From Fig. 1(c) we see that the density also varies linearly in the bulk, while at the solid surface it is more close to the final state than that obtained from the CIS. Since the diffusion-type macroscopic equation (27) allows the efficient exchange of information, fast convergence is realized in the whole computational domain, see Fig. 1(c) and (d).

Figure 2 demonstrates how fast the solution is converged at different values of rarefaction parameter. When $\delta_{rp}$ is small, the errors in both the CIS and GSIS decays at the same rate, which means that the two schemes are as efficient as each other. As $\delta_{rp}$ increases so that the flow enters the transition and near-continuum regimes, the error in the CIS oscillates several times before it decays monotonically. As a consequence, the iteration number of CIS increases rapidly with the rarefaction parameter, which nearly scales as $\delta_{rp}^{2}$. For the GSIS, however, the error is monotonically decreasing, and the rarefaction parameter does not influence the error decay rate, where the converged solutions are obtained within the same number of iterations (here 20 iterations) for each rarefaction parameter from the free molecular to continuum flow regimes. At $\delta=50$, the GSIS is about 100 times more efficient than the CIS, and it can be expected that the gain of using GSIS becomes larger and larger as $\delta_{rp}$ further increases.

Another important property of the GSIS is that the numerical error caused by the spatial discretization is much reduced when compared to that of the CIS. From Fig. 3 we see that when $N_{2}$ is decreased from 251 to 6, that is, when the spatial cell size is respectively about $1/10$ and 5 times of the mean free path of gas molecules, the relative error in the density profile increases from 0.3% to 9%, while that in the heat flux increases from 0.3% to 16% in the CIS. However, the relative error in the GSIS always remain within 1%, even when the cell size is about 5 times larger than the gas mean free path. Note that even when $\delta_{rp}=500$, the heat flux obtained from the GSIS only changes from $3.721\times 10^{-3}$ when $N_{2}=551$ to $3.726\times 10^{-3}$ when $N_{2}=6$. The reason for this excellent performance is that the GSIS is asymptotically preserving the Navier-Stokes limit, while in the CIS the “numerical” thermal conductivity may be different to the physical one. Besides, in the CIS, the false convergence, e.g. the non-uniform distribution of heat flux in in Fig. 3(b), may be reached when the spatial resolution is not enough. The superior GSIS, however, does not suffer this problem.

It should be noted that the implicit UGKS [44] and other variants [45, 46] can also produce accurate results when the cell size is much larger than the molecular mean free path. This is achieved through a complex evaluation of the numerical flux at the cell interface to spontaneously treat the molecular streaming and collision. The GSIS, however, does not need complex flux evaluation.

Using the accurate and efficient GSIS, the LBE is solved for different molecular collision models (6) and the corresponding Knudsen layer functions are obtained. In the numerical simulation, we set the rarefaction parameter to be $\delta_{rp}=60$, so that the distance between two plates is about 60 times as large as the mean free path of gas molecules; thus, the interference between the Knudsen layers near each plate is avoided. In the fast spectral approximation of the linearized Boltzmann collision operator (4), the integral with respect to the solid angle $\Omega$ is calculated by the Gauss-Legendre quadrature with $M=6$, see Eq. (39) in Ref. [38]. In the spatial discretization we let

[TABLE]

with $N_{s}=200$. The iterations terminate when $\epsilon<10^{-6}$.

When the steady-state solution is obtained, the temperature profile in the bulk region (i.e. $0.4\leq{}x_{2}\leq 0.5$) is linearly fitted by $T_{NS}=k_{1}(x_{2}-1/2)$ in the dimensionless form, where $k_{1}$ is the coefficient from the least square fitting. Then the KLF is calculated according to the following equation:

[TABLE]

and the temperature jump coefficient, according to Sharipov’s review paper, is calculated as

[TABLE]

The GSIS results for the LBE with the HS and Maxwell molecules reaches the steady-state after 22 and 27 iterations, respectively, and the temperature jump coefficients are respectively 1.892 and 1.954, which do not vary a lot to the collision model. However, the Knudsen layer functions shown in Fig. 4 has larger difference. It is amazing that the small terms $2\int{(L-L_{s})v_{i}v_{j}}\mathrm{d}\bm{v}$ in Eq. (17) and $\int{(L-L_{s})v_{i}|\bm{v}|^{2}}\mathrm{d}\bm{v}$ in Eq. (20) significantly affect the Knudsen layer function.

4.2 Two-dimensional lid-driven cavity flow

The two-dimensional lid-driven cavity flow is a canonical test for the algorithms of both Navier-Stokes equations and gas-kinetic equations. The flow domain is a square with size of $1\times 1$, with the left and right walls locate at $x_{1}$ = 0 and $x_{1}=1$, bottom and top walls at $x_{2}=0$, $x_{2}=1$. The top wall, i.e., the lid moves in the $x-$direction with a constant velocity of $U_{w}$, while the other sides are static walls. All of the walls are kept at uniform temperature of $T_{0}$. To demonstrate the accuracy and efficient of the GSIS, the Shakhov kinetic equation is linearized by choosing $\alpha=U_{w}/v_{m}$ in Eq. (2). The boundary conditions are

[TABLE]

The problem is solved on non-uniform Cartesian grids, which are discretized by

[TABLE]

where $N_{s}$ is the number of grid nodes in both $x_{1}$ and $x_{2}$ axis. The linearized Shakhov equation is solved by DVM with the 2nd-order upwind finite-difference scheme, where the distribution functions stored at the centers of grid cells. In the synthetic acceleration step, the continuity equation and the momentum equations in Eq. (13) are solved using a finite-difference version of the well known Semi-Implicit Method for Pressure Linked Equations (SIMPLE). In each SIMPLE iteration, we solve four discrete diffusive equations (the two velocity components, pressure correction, and temperature) using the Jacobi iteration methods.

When the macroscopic flow variables are solved by SIMPLE algorithm, the velocity distribution function is updated as

[TABLE]

because (i) the update of the shear stress and heat flux does not affect the accuracy and efficiency of the GSIS, and (ii) for highly rarefied gas flows, the high-order terms are very large and the macroscopic synthetic equations become stiff near the solid corners due to the small value of $\delta_{rp}$, hence the limiter ${\delta_{rp}}/{\max(10,\delta_{rp})}$ is introduced to increase the numerical stability.

We first test the converging speeds of the CIS and GSIS for the cases of $\delta_{rp}=0.1$, 1, 10, 100 and 1000. The corresponding spatial grids are non-uniform with $N_{s}=21$, 21, 21, 41, 61 respectively. For the cases of $\delta_{rp}=0.1$, 1 and 10, the molecular velocity in both $v_{1}$ and $v_{2}$ are discretized by Eq. (32), with $\imath=3$, and $N_{v}=48$, 48 and 24, respectively. For $v_{3}$, 24, 24 and 12 uniform points in the range of $[-6,6]$ are used. While for $\delta=100$ and 1000, the 6- and 8-point Gauss-Hermite quadrature nodes are used in all three velocity components. The iterations in both CIS and GSIS are assumed to be converged when

[TABLE]

Figure 5 compares the decay of error $\epsilon$ as a function of the number of iteration steps in CIS and GSIS for solutions of flows at different values of rarefaction parameter, while Table 1 summarizes the number of iteration steps and the total CPU time of the calculations with a single threaded Matlab 2018 code on Intel Xeon-E5-2680 v4 CPU. Similar to the test case of the Fourier heat transfer, at small $\delta_{rp}$ (0.1 and 1), the errors in both GSIS and CIS decay with the same rate and converged in less than 20 steps. In the cases of larger $\delta_{rp}$, the iteration step in GSIS slightly increases, but it is less than 40 steps even for the case of $\delta_{rp}=1000$. In contrast, the convergence of the CIS iteration deteriorates severely as $\delta_{rp}$ increase. The iteration step reaches 1823 for the case of $\delta_{rp}=100$. Due to the slow convergence of CIS for near continuum flows, the case of $\delta=1000$ is not simulated.

With significantly faster convergence rate, the GSIS takes much less CPU time than the CIS for cases of large $\delta_{rp}$ as shown in Table 1. Note that although the iteration number is reduced in GSIS, the time for each iteration increases as the cost to solve the synthetic equation is non-negligible. see the last column of Table 1. This is because the segregated approach of the SIMPLE algorithm can take up to several hundreds of iterations to converge, depending on the value of $\delta_{rp}$. We note that using a coupled algorithm to solve the discretized pressure and velocity components in a single linear equation system would be much faster than the segregated approach, especially for high $\delta_{rp}$ cases, as have been studied in the incompressible CFD theories. For example, in the following section we find that if the kinetic synthetic equations are solved by the DG, the cost of DG for synthetic equations is negligible since pressure, velocity, and temperature are solved simultaneously.

To compare the accuracy of the GSIS with the CIS, we simulated the case of $\delta_{rp}=100$ with different non-uniform physical grids, including $N_{s}=21$, 41, 61 and 101. Figure 6 presents the comparisons of the pressure fields and streamlines predicted by both the CIS and GSIS on various grids, in which the reference solutions are taken as the GSIS results on grid of $N_{s}=61$. The results show that the GSIS solution on the coarsest grid ($N_{s}=21$) is much more accurate than the CIS counterpart, especially in terms of the pressure field. From Fig. 6(d) to (f), we can observe that the short contour lines near the bottom wall are accurately captured by the GSIS even on the coarsest mesh, while the CIS can capture them only with the finest mesh ($N_{s}=100$, see Fig. 6 (c).)

4.3 Shear-driven flow between two eccentric cylinders

In this section, we consider a shear-driven gas flow between two noncoaxial cylinders. This test case is used to show that the proposed synthetic iterative scheme can be efficiently implemented through other CFD method rather than the finite difference algorithm to deal with more complicated geometries. As shown in Figure 7, the outer cylinder with a radius of 2 rotates clockwise at a constant speed of $U_{w}$, while the inner cylinder with a radius of 1 keeps static. The centers of the outer cylinder and inner cylinder are at $\bm{x}=(0,0.5)$ and the origin, respectively. The cylinders are of a constant temperature $T_{0}$. It is assumed that $U_{w}$ is much smaller than the most probable speed $v_{m}$, thus the gas system can be linearized with $\alpha=U_{w}/v_{m}$. The distribution function for the reflected molecules at the outer cylinder is given by

[TABLE]

where $\bm{n}_{w}$ and $\bm{t}_{w}$ denote the outward unit normal vector and tangential vector of the solid surface. The boundary condition at the inner cylinder is similar but without the term $\bm{t}_{w}\cdot\bm{v}$.

Using both the GSIS and CIS, the shear-driven flow is resolved on structured triangular mesh, in which the grid nodes along the radial direction is described by Eq. (38). The high-order DG methods are employed to seek solutions of the linearized Shakhov model equation and the synthetic macroscopic equations in pecewise polynomial spaces of degree of 3. The detailed DG scheme for the gas kinetic equation can be found in [42], while the hybridizable DG algorithm to solve the synthetic macroscopic equations is listed in the Appendix.

The resultant velocity contours and streamlines are illustrated in Fig. 8 for two selected rarefaction parameters $\delta_{rp}=1000$ and 10, in which the GSIS solutions are plotted in the left half domain and the CIS ones are plotted in the right half domain. The results at $\delta_{rp}=1000$ are obtained on 2400 triangles with cell size (characterized by the height of triangle) varying from 3 to 260 times the mean free path of gas molecules. The molecule velocity space is discretized by 8-point Gauss-Hermite quadrature nodes in $v_{1}$ and $v_{2}$ and 12 equidistant nodes in the range of $[-4,4]$ in $v_{3}$. The results at $\delta_{rp}=10$ are obtained on 1600 triangles with cell size varying from 0.1 to 3 times the mean free path of gas molecules. The molecule velocity space is discretized in the domain of $[-4,4]^{3}$ by 32 non-uniform nodes in $v_{1}$ and $v_{2}$ and 24 equidistant nodes in $v_{3}$. The solutions are believed to be converged when the relative error in velocity magnitude $|\bm{U}|$ between two consecutive iteration steps is less than $10^{-5}$. The streamlines show that, as the gas rotates clockwise from the top to the bottom, due to the shrink of the flow pass, part of the gas near the outer surface is squeezed into the bottom narrow space while the other part of the gas flows back along the surface of the inner cylinder; as a consequence, a vortex appears above the inner cylinder.

Large discrepancies in the velocity contours are observed between the GSIS and CIS results at $\delta_{rp}=1000$. To test the accuracy of both schemes, we also include the results obtained by only solving the Navier-Stokes equations with the non-slip velocity boundary condition, which are illustrated by the white dashed lines in Fig. 8(a) and (b). The GSIS results coincide with the ones from the Navier-Stokes equations, thus the GSIS can asymptotically preserve the Navier-Stokes limit. However, the CIS cannot predict accurate solutions due to the large numerical dissipation on such a coarse mesh, i.e. the maximum cell size is about 260 times of the molecular mean free path. As the rarefaction parameter decreases to 10, the GSIS and CIS can produce close solutions on the same mesh.

Consider the rate of convergence to the steady-state solution, the GSIS cost only 26 iterative steps to reach the convergence criterion for both the cases of $\delta_{rp}=1000$ and 10, while the CIS consumes 49454 and 296 steps, respectively. Since compared to that of solving the kinetic equation, the computational consumption for DG to solve the macroscopic equations is negligible, since the number of degrees of freedom for the latter one is much smaller. Therefore, the GSIS can be nearly 2000 and 10 times faster than the CIS when $\delta_{rp}=1000$ and 10, respectively.

5 Numerical results for periodic oscillation problems

For linearized problems, if the external force that drives the flow changes periodically in time, then the velocity distribution function can be expressed as [54, 55, 56]:

[TABLE]

where $\Re$ is a real part of a variable and $h$ satisfies the following linearized Boltzmann equation:

[TABLE]

Note that here $i$ is the imaginary unit and $h$ is a complex function, so are the macroscopic quantities defined in Eqs. (8) and (9). These complex values will introduce phase shifts relative to that of the external disturbance. The Strouhal number $\mathrm{St}$

[TABLE]

is the oscillation frequency $\varpi$ normalized by $v_{m}/H$.

The solutions to these oscillating problems can also be accelerated by the GSIS; the corresponding macroscopic synthetic equations can be derived if we replace $\partial/\partial{t}$ in Sec. 3 by $i\text{St}$. Several numerical examples are given to demonstrate the accuracy and efficient of the present method.

5.1 Spontaneous Rayleigh-Brillouin scattering

In the spontaneous Rayleigh-Brillouin scattering (SRBS), light propagating through the gas is scattered by the thermal motion of gas molecules, where the spectrum of the scattered light contains the information of gas such as temperature, speed, and viscosity. Thus, SRBS provides a non-intrusive way to probe the gas properties in a remote way. Theoretically, the SRBS spectrum can be obtained by solving the LBE (3) with the initial condition $h(t=0,{x_{2}},\bm{v})\propto{}\delta({x_{2}})f_{eq}(\bm{v})$, which represents a density impulse [57, 58]. To be more specific, the SRBS spectrum is calculated as

[TABLE]

where $\hat{h}(\bm{v})$, the Laplace and Fourier transforms of $h$ in the temporal and spatial directions, respectively, satisfies (suppose the scattered light propagates in the $x_{2}$ direction)

[TABLE]

Note that in Eq. (44) and (45), the rarefaction parameter $\delta_{rp}$ is defined when the characteristic flow length $H$ is $\lambda_{L}/2\sin(\theta_{s}/2)$, with $\lambda_{L}$ being the wavelength of laser and $\theta_{s}$ the angle of light scattering, and $f_{s}(=\text{St}/2\pi)$ is the frequency shift in the scattering process normalized by the characteristic frequency $v_{m}/H$, and the hat denotes the Laplace-Fourier transform of the corresponding quantity. Also note that terms in the left-hand-side of Eq. (45) appear because operators $\partial/\partial{t}$ and $\partial/\partial{x_{2}}$ in Eq. (3) are replaced by $2i\pi{f_{s}}$ and $-2i\pi$, respectively. Finally, the source term $f_{eq}$ in Eq. (45) is from the Laplace transform of the initial density impulse. This term will change the first equation in Eq. (13) to $\frac{\partial{\rho}}{\partial{t}}+\frac{\partial{U_{i}}}{\partial{x_{i}}}=1$, while other synthetic equations remain unchanged.

In CIS, the velocity distribution function is obtained by solving the following equation iteratively:

[TABLE]

which converges fast when $\delta_{rp}$ is small, but extremely slow when $\delta_{rp}$ is large as the flow enters the near-continuum regimes.

In the GSIS, the synthetic equations can be obtained by solving the following matrix at the $(k+1)$-th iteration step:

[TABLE]

where $\text{R}_{4}=2i\pi{}\text{HoT}_{\sigma_{22}}^{(k+1/2)}+2\int{}(L-L_{s})v_{2}^{2}\mathrm{d}\bm{v}$ and $\text{R}_{5}=2i\pi{}\text{HoT}_{q_{2}}^{(k+1/2)}+\int{}(L-L_{s})v_{2}|\bm{v}|^{2}\mathrm{d}\bm{v}$.

In the numerical simulation, started from the zero perturbance at each frequency difference, solutions are believed to be converged when the relative error in $\hat{\rho}$ between two consecutive iteration steps is less than $10^{-7}$. Results in Fig. 9(a) show that GSIS and CIS generate almost the same SRBS spectra, except at $\delta_{rp}=50$ the CIS has a false converged solution (i.e. the discontinuous spectrum) when the frequency difference is around 0.68. As usual, the iteration number in CIS increases significantly with the rarefaction parameter $\delta_{rp}$, while in GSIS this remains nearly unchanged and is far less than that of the CIS. For example, the iteration number of the GSIS is about 10 and 100 times less than that of the CIS when $\delta_{rp}=10$ and 50, respectively. We have also tested that, even when $\delta_{rp}=500$, converged solutions are obtained within 20 steps in the GSIS for every frequency difference.

However, when $\delta_{rp}$ is small and $\mathrm{St}$ is large, the GSIS does not converge or even blows up. This is because the eigenvalue of the matrix in Eq. (47) has large complex values so that any inappropriate initial guess lead to large oscillations that decay rather slow or even blow up. Whereas, physically speaking, the solution should decay fast due to the large rarefaction effect. To remedy this, the small value $\delta_{rp}$ in the left-hand side of Eq. (47) is replaced by a relative large value $\bar{\delta}_{rp}=\max(\delta_{rp},10)$, while the right-hand side terms are modified correspondingly as

[TABLE]

This simple treatment helps to decay non-physical solutions at initial few iteration steps. When the solution of the new system converges, it can be proven that it satisfies Eq. (47). Therefore, no approximation is introduced to the converged solution. This point is proven in Fig. 10, where the GSIS and CIS solutions agree perfectly with each other, and from the inset we see that the GSIS needs slightly less iteration steps than CIS in most of frequency differences.

Another remarkable property of the GSIS is that, at the same level of convergence criterion, the GSIS provides more accurate numerical solutions. One example is given in Fig. 11, where one can see that the relative error between two consecutive iteration steps $\epsilon=|\hat{\rho}^{(k+1)}/\hat{\rho}^{(k)}-1|$ decays rather fast in the GSIS, while in the CIS it decreases slowly with many oscillations. As a consequence, the GSIS finds the correct spectrum profile even when the relative error in density is less $10^{-2}$, while the CIS can only reach the correct solution when the error is less than $10^{-6}$. This can be explained below. According to the analysis of Adam and Larsen for radiation transfer problem [23], if one stops at the $(k+1)$-th step with

[TABLE]

in the CIS, then the relative difference from the true solution $\hat{\rho}$ is

[TABLE]

where $\gamma$ is the spectral radius of the iteration operator. For problem with slow convergence, $\gamma$ is very close to one (see Figure 1 in Ref. [5] for the kinetic BGK model equation), which could make the difference from true solution magnified by thousands of times.

5.2 Oscillatory Couette flow between two parallel plates

Consider the rarefied gas dynamics between two infinite parallel plates with a distance $H$, located at $x_{2}=0$ and $x_{2}=1$. Both plates have a temperature $T_{0}$, the one at $x_{2}=1$ is stationary, while that at $x_{2}=0$ oscillating in the $x_{1}$ direction with the velocity

[TABLE]

The Boltzmann equation is linearized by choosing $\alpha=U_{0}/v_{m}$ in Eq. (42). If we consider the diffuse boundary condition, then we have $h(x_{2}=0,\bm{v})=2v_{1}f_{eq}$ when $v_{2}>0$, and $h(x_{2}=1,\bm{v})=0$ when $v_{2}<0$ [54]. The synthetic equations (13), (17), and (20) can be simplified to

[TABLE]

where the moments involving even order of $v_{1}$ are all zero, and we do not consider the heat flux $q_{1}$ in this problem as it does not affect the rate of convergence. It is noted that the above equations reduce to the synthetic equation developed in Ref. [52] when $\text{St}=0$.

The two equations in Eq. (52) can be combined to produce the following diffusion equation for the flow velocity $U_{1}$ in the $(k+1)$-th iteration step:

[TABLE]

In the numerical simulation, the spatial space is discretized by Eq. (38) with $N_{s}=100$. The kinetic equation (22) is solved by the second-order upwind scheme, while the derivative in Eq. (53) is approximated by the central finite difference scheme with 5 stencils, and the resulting linear algebraic system for $U_{1}$ is solved exactly in the bulk region (i.e. at least three spatial points away from the boundary).

The comparison in the accuracy and efficiency between the CIS and GSIS is summarized in Fig. 12, where the molecular velocity space is discretized in the same way as that in Sec. 4, but with $N_{v}=96$ in Eq. (32). The relative difference in the amplitude of the shear stress $\sigma_{12}$ is within 1%. When the rarefaction parameter is $\delta_{rp}=50$, we see that the number of iterations in the CIS decreases from 30,000 to 100 when Strouhal number increases from 0 to 50. The reason for this reduction can be understood in the following way. The temporal Knudsen number $\text{Kn}_{t}$, which is defined as the ratio of characteristic oscillation frequency to the mean collision frequency of gas molecules, i.e.

[TABLE]

increases with St. Therefore, even when $\delta_{rp}$ is large, that is, when the spatial Knudsen number is small, the large temporal Knudsen number can also make the flow rarefied, and the more rarefied the gas is, the fast the iteration to the steady-state. Even with this effect, the GSIS is still faster than the CIS: only about 20 iterations are needed in the GSIS for each Strouhal number considered.

However, for the GSIS in oscillating problems, there is a problem, like the one encountered in Sec.5.1. From Eq. (53) we see that the eigenvalue of this second-order differential equation is imaginary, which means that when $\delta_{rp}$ is small and St is large, the solution will change quasi-periodically in the spatial direction with large frequency, whereas physically the solution should decay fast from the oscillating sources as the dissipation is huge due to the large values of both spatial and temporal Knudsen numbers. Mathematically speaking, for highly oscillating solutions, any slight inaccurate boundary conditions will lead to completely different solutions. Therefore, in the numerical simulation, when we solve Eq. (53) directly, the solution is either wrong or blows up. To fix this problem, again we introduce a relative large value of $\bar{\delta}_{rp}$ to decay the fast oscillation. That is, instead of solving Eq. (53), we solve the following diffuse-type equation:

[TABLE]

where

[TABLE]

It can be proven that, when the solution of Eq. (55) converges, Eqs. (55) and (53) are equivalent. This treatment does not affect the accuracy and efficiency of the GSIS when $\delta_{rp}$ is small, while when $\delta_{rp}$ is large, the solution from the synthetic equations are always stable, and we see in Fig. 12 that in most cases the GSIS needs slightly less iterations than CIS.

In addition to the significant reduction of iteration number, the GSIS needs less spatial grids than that of the CIS. Two examples are given in Fig. 13, where one can see that the GSIS can yield accurate results even when the cell sizes are respectively about 6.6 and 50 times of the molecular mean free path, while the CIS has large error due to the strong numerical dissipation.

5.3 Sound propagation between two parallel plates

Consider the sound propagation through a gas between two infinite parallel plates with a distance $H$, located at $x_{2}=0$ and $x_{2}=1$. The two plates have a temperature $T_{0}$, the one at $x_{2}=1$ is stationary, while that $x_{2}=0$ oscillating in the $x_{2}$ direction with the speed $U_{w,2}=\Re\left[U_{0}\exp({i\mathrm{St}t})\right]$. The Boltzmann equation is linearized by choosing $\alpha=U_{0}/v_{m}$ in Eq. (42). The boundary conditions are [55]

[TABLE]

The synthetic equations (13), (17), and (20) can be simplified to

[TABLE]

These synthetic equations can be combined to form two diffusion equations for the flow velocity $U_{2}$ and temperature $T$. To quickly decay the non-physical oscillations when $\delta_{rp}$ is small and St is large, in numerical iterations we set

[TABLE]

where $\bar{\delta}_{rp}$ is given in Eq. (56). When $U_{2}$ and $T$ are solved, the perturbed density, shear stress and heat flux can be solved from Eqs. (58), (61), and (62).

Typical numerical results are shown in Fig. 14 when the spatial region $x_{2}\in[0,1]$ is discretized by 200 uniformly-distributed points, while the velocity grids are the same as that used in Sec. 5.2. For the CIS, it is very hard to find the converged solution when the Strouhal number St is small, where the iteration number scales roughly as $\text{St}^{-1.5}$. However, this problem does not exist in the GSIS, as the Strouhal number has little effect on the number of iterations. The effect of spatial resolution on the fidelity of the solution is demonstrated in Fig. 15 when $\text{St}=2.5$, where the sound waves between two plates resonance. It is seen that the GSIS needs less spatial grids than CIS. Again, this example proves the accuracy and efficiency of the GSIS.

5.4 Two-dimensional oscillatory Couette flow

Finally we consider the oscillatory flow in a three-dimensional cavity shown in Fig. 16(a). We assume the side length OD is much larger than OH and OA, so that the problem is quasi two-dimensional. The characteristic length $H$ is chosen as the side length OA, and the aspect ratio is defined as $\text{Asp}=OH/OA$. If $\text{Asp}=\infty$, the problem is just the oscillatory Couette flow between two parallel plates studied in Sec. 5.2. This problem is interesting because it displays a counter-intuitive phenomenon that the shear force exerting on the oscillating lid in two-dimensional cavity could be even smaller than that of the one-dimensional Couette flow [56]. The full three-dimensional oscillatory flow was studied in Ref. [59], but not all the parameter region are covered, for example, the case with OA much larger than OA and OH.

The synthetic equations (13), (17), and (20) can be simplified to

[TABLE]

which leads to the following diffusion-type equation for the flow velocity $U_{1}$ that is solved in an stable iterative manner:

[TABLE]

where $\bar{\delta}_{rp}$ is given in Eq. (56), and

[TABLE]

In numerical simulations, the molecular velocities $v_{2}$ and $v_{3}$ are discretized non-uniformly according to Eq. (32) with 48 points in each direction, while $v_{1}$ is truncated into the region of $[-6,6]$ and discretized by 24 uniformly-distributed points. Due to the symmetry $h(x_{2},x_{3},v_{1},v_{2},v_{3})=h(x_{2},\text{Asp}-x_{3},v_{1},v_{2},-v3)$, we only consider the domain where $0\leq{}x_{2}\leq 1$ and $0\leq{}x_{3}\leq\text{Asp}/2$, which are discretized by 50 and 60 points according to Eqs. (38) and (34), respectively. The velocity distribution function entering the domain from the stationary walls is zero, while that from the oscillating wall is $2v_{1}f_{eq}$. The diffusion-type equation (65) is approximated by the central finite difference with 5 stencils, which are solved exactly by rewriting it in the matrix form. From the inset of Fig. 16(b) we see that the GSIS is very efficient as converged solutions are obtained within 40 iterations.

We are interested in how the average shear force exerting on the oscillating lid change with the normalized oscillation frequency St. Hence $\sigma_{13}$ is not considered here as it is anti-symmetric along the line $x_{3}=\text{Asp}/2$ so its overall contribution to the friction is zero. The amplitude of the average shear force on the oscillating lid is defined as

[TABLE]

which is shown in Fig. 16(b) for different aspect ratios of the cavity over a wide range of the oscillation frequency, when $\delta_{rp}=50$. It can be seen when $\text{Asp}=2$, the average shear force is the same as that of $\text{Asp}=\infty$, except that it is slightly larger when St is small. This is seen more clearly in Fig. 16(c) that the two lateral walls, i.e. the left and right walls in Fig. 16(a), increases the shear stress from a nearly small constant to a high rise near the left top corner. When St increases, the shear stress quickly decay from the oscillating lid to the zero value at the bottom surface, and its value at the oscillating lid is nearly uniform, see Fig. 16(d) and (e). As the aspect ratio of the cavity reduces, the average shear force increases when St is small, see Fig. 16(f) for an example; this is easy to understood as the lateral walls increase the total friction according to our daily life experience. However, from Fig. 16(f) and (g) we can see that the shear stress quickly saturates, as the increase of oscillation frequency only slightly increases the shear stress at the lid, such that the average shear stress on the lid remains nearly constant over a wide range of St; and the smaller the aspect ratio is, the wider this region is. This may be useful to design a micro-electro-mechanical system where the shear force remains constant in a certain wide range of oscillation frequency. It is this efficient algorithm we are able to find this new phenomenon which are missed in Ref. [59]. Another counter-intuitive thing is that, when St is large, the average shear force at small values of cavity aspect ratio is slightly smaller than that of the one-dimensional cavity, although the relative difference is within 5%.

6 Conclusions and outlooks

In summary, we have developed a general synthetic iteration scheme to find the steady-state solution of the linearized Boltzmann equation efficiently and accurately. Various numerical results have demonstrated that our scheme is able to find the converged solution within about 20 iterations at any Knudsen number, due to the fact that the synthetic macroscopic equations not only asymptotically preserves the Navier-Stokes limit in the framework of Chapman-Enskog expansion, but also explicitly contains the constitutive laws for the shear stress and heat flux at the first order approximation in the the Knudsen number to the linearized Boltzmann equation. As a consequence, accurate solutions that are not contaminated by large numerical dissipation and accumulated error, can be obtained when the cell size is much larger than the mean free path of gas molecules. Moreover, the numerical error in the synthetic iteration scheme decays very fast and the convergence criterion can be set at a much smaller value than the conventional iteration scheme. These factors enables our synthetic iteration scheme to find the steady-state solution in 10-ish iterations.

This paper provides a framework to solve the general rarefied gas flow problems. The advantages and future works are highlighted below:

Compared to implicit UGKS [44] and it variants [45, 46], we conclude that in order to develop efficient multiscale numerical schemes, macroscopic equations must be solved together with the Boltzmann or kinetic model equations. While in Refs. [44, 45, 46] only five equations from the conservation law are used so that complex flux evaluation across the cell interface must be adopted to asymptotically preserve the Navier-Stokes limit, our scheme needs no complex flux evaluation as the Navier-Stokes equations are recovered explicitly. Thus, the numerical implementation is much easy and the convergence to steady-state solution is much faster. More importantly, our scheme does not depend on the specific form of the collision operator, while that in Refs. [44, 45, 46] relies only on the BGK-type kinetic equations to enable exact evaluation of numerical flux. 2. 2.

Since the limitation on the cell size is removed and fast convergence is enabled, the present synthetic iteration scheme can be directly applied to low-variance [60, 61] and even frequency-domain [62] DSMC that solves the linearized Boltzmann/kinetic model equations to improve the computational efficiency, especially in the near-continuum flow regime. 3. 3.

The present work can be extended to multi-species and compressible flow easily. The key is to construct macroscopic equations which recovers the compressible Navier-Stokes equation to the first order of Knudsen number. As a matter of fact, the Grad 13 moment equations [48, 49] can be directly used if the high-order velocity moments are calculated from the numerical solution of the Boltzmann equation, rather than closed by making assumption on the form of velocity distribution function. Actually the authors have implemented the synthetic iteration scheme for nonlinear Fourier heat transfer, and started from the global equilibrium distribution converged solution at arbitrary Knudsen number is found within 20 iterations. 4. 4.

It is noted that recently the gas-kinetic wave-particle (UGKWP) method, which uses the essential idea of UGKS that the streaming and collision should be treated spontaneously, has been applied in the framework of DSMC to remove the constraint on the cell size when the Knudsen number is small [63, 64]; the BGK kinetic model is solved and the complex and time-consuming particle sorting is used to enable the asymptotically preserving property. We believe that the synthetic iteration scheme can also be applied to DSMC to remove the limitation on cell size and boost convergence, and the advantage is clear: it relies on no specific collision operator so that can be extended naturally to multi-species flows and even flows involving chemical reactions.

With these new development implemented, it is foreseen that in the near future that the problem of numerical simulation of multiscale rarefied gas flows will be solved completely. Also, the same idea can be applied to other kinetic equations such as the Enskog equation for dense gases dynamics with applications to shale gas extraction and non-equilibrium evaporation and condensation [65, 66, 67].

Appendix

Here, some details to solve the synthetic macroscopic equations using the high-order hybridizable discontinuous Galerkin (HDG) method [68] on arbitrary triangular mesh are presented. The steady-state governing equations can be written in the following mixed form as a system of first-order equations

[TABLE]

where

[TABLE]

with $\bm{I}$ being the identity matrix. The auxiliary variables $\bm{L}$ and $\bm{E}$ are introduced to approximate the combination of the velocity gradient $\nabla\bm{U}$, temperature gradient $\nabla T$ and the high-order moments. Then, the stress tensor and heat flux are evaluated as

[TABLE]

Let $\Delta\in\mathbb{R}^{2}$ be an two-dimensional domain with boundary $\partial\Delta$ in the $x_{1}-x_{2}$ plane. Then, $\Delta$ is partitioned in $M$ disjoint regular triangles $\Delta_{i}$: $\Delta=\cup^{M}_{i}\Delta_{i}$. The boundaries $\partial\Delta_{i}$ of the triangles define a group of $N$ faces $\Gamma_{c}$: $\Gamma=\cup^{M}_{i}\{\partial\Delta_{i}\}=\cup^{N}_{c}\{\Gamma_{c}\}$. For HDG discretization, two types of discontinuous finite element approximation space, one for solutions within $\Delta_{i}$ and the other for traces of solution on $\Gamma_{c}$, are defined as

[TABLE]

where $\mathcal{P}^{k}(D)$ denotes the space of $k-$th order polynomials on a domain $D$.

The HDG method solves the system in two steps. First, a global problem is set up to determine the traces of the flow properties $\hat{\bm{Q}}=\left[\hat{p},\hat{\bm{U}},\hat{T}\right]$ on the faces $\Gamma$. Then, a local problem with $\hat{\bm{Q}}$ as the boundary condition on $\partial\Delta_{i}$ is solved element-by-element to obtain the solutions for the flow properties $\bm{Q}=\left[p,\bm{Q},T\right]$, as well as the ones for the auxiliary variables $\bm{L}$ and $\bm{E}$. Generally speaking, when moving from the interior of the triangle element $\Delta_{i}$ to its boundary $\partial\Delta_{i}$, the traces defines what the values of field variables on the boundary should be. In the HDG method, it is assumed that the traces are singled-valued on each face.

We introduce the notations $\left(a,b\right)_{D}=\int_{D\in\mathbb{R}^{2}}\left(a\odot b\right)\mathrm{d}x_{1}\mathrm{d}x_{2}$ and $\langle a,b\rangle_{D}=\int_{D\in\mathbb{R}^{1}}\left(a\odot b\right)\mathrm{d}\Gamma$, where $\odot$ can be either the dot product $\cdot$ or tensor product $\otimes$. The local problem is stated as: find $\left(\bm{Q},\bm{L},\bm{E}\right)\in\left[\mathcal{V}\right]^{4}\times\left[\mathcal{V}\right]^{4}\times\left[\mathcal{V}\right]^{2}$ such that

[TABLE]

for all $\left(\bm{r},\bm{w},\bm{z}\right)\in\left[\mathcal{V}\right]^{4}\times\left[\mathcal{V}\right]^{4}\times\left[\mathcal{V}\right]^{2}$. The numerical flux $\hat{\bm{\mathcal{F}}}\cdot\bm{n}$ is defined as [69]

[TABLE]

Here $\bm{n}$ being the outward unit normal vector of $\partial\Delta_{i}$. $\tau$ is the stabilization parameter that have important effects on the accuracy and convergence of the HDG method. In this work, we chosen $\tau=1/H_{\text{min}}$, with $H_{\text{min}}$ the minimum height of the triangles $\Delta_{i}$.

The global problem is set up by enforcing the continuity of the numerical flux over all the interior faces. It is stated as: find $\hat{\bm{Q}}\in\left[\mathcal{W}\right]^{4}$ such that

[TABLE]

for all $\bm{\psi}\in\left[\mathcal{W}\right]^{4}$. Here the superscripts $\pm$ denote the numerical fluxes obtained from the triangles on both sides of the face. Note that the traces on boundary faces are calculated as

[TABLE]

where $\bm{Q}^{+}_{VDF}$ is the field solutions directly calculated from the approximated velocity distribution function (see Eq. (8)) within the triangle where the boundary face $\Gamma_{c}$ belongs to.

By assembling the local problem (A.72) and global problem (A.74) and (A.75) over all the triangles and faces, we can obtain a matrix system

[TABLE]

where $\mathbb{Q}$, $\mathbb{L}$, $\mathbb{E}$ and $\hat{\mathbb{Q}}$ are the vectors of degrees of freedom of the flow properties $\bm{Q}$, the auxiliary variables $\bm{L}$ and $\bm{E}$, and the trace of the flow properties $\hat{\bm{Q}}$, respectively. Note that the degrees of freedom for $\bm{Q}$, $\bm{L}$ and $\bm{E}$ are grouped together and ordered element-by-element, and the corresponding coefficient matrix $\left[A_{Q},A_{L},A_{E};B_{Q},B_{L},B_{E};C_{Q},C_{L},C_{E}\right]$ has block-diagonal structure. Therefore, we can eliminate $\bm{Q}$, $\bm{L}$ and $\bm{E}$ to obtained a reduced linear system involving only $\hat{\mathbb{Q}}$. Once $\hat{\mathbb{Q}}$ is determined, $\bm{Q}$, $\bm{L}$ and $\bm{E}$ are reconstructed corresponding to the local problem (A.72) in an element-wise fashion, while the stress tensor and heat flux are calculated as Eq. (A.70).

Acknowledgments

This work is supported in the UK by the Engineering and Physical Sciences Research Council under grant EP/R041938/1. L. Zhu acknowledges the financial support of European Union’s Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie grant agreement number 793007.