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

Wei Su11footnotemark: 1 Lianhua Zhu22footnotemark: 2 Peng Wang33footnotemark:
3 Yonghao Zhang44footnotemark: 4 Lei Wu55footnotemark: 5
<EMAIL_ADDRESS>James Weir Fluids Laboratory, Department of Mechanical
and Aerospace Engineering, University of Strathclyde, Glasgow G1 1XJ, UK

###### 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.
The key ingredient of our scheme is that the macroscopic equations, which are
solved together with the Boltzmann equation and help to adjust the velocity
distribution function, not only asymptotically preserves the Navier-Stokes
limit in the framework of Chapman-Enskog expansion, but also contain Newton’s
law for stress and Fourier’s law for heat conduction explicitly. For this
reason, like implicit UGKS, the constraint that the numerical cell size should
be smaller than the mean free path of gas molecules is removed, but we do not
need the complex evaluation of numerical flux at the cell interface. What’s
more, 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.

††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

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

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

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

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:

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

where the linearized Boltzmann collision operator is [37]:

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

and the equilibrium collision frequency is

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

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]

$B(|{\bm{u}}|,\theta)=\frac{|{\bm{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)},$
(6)

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:

$f(t,\bm{x},\bm{v})=\frac{2\int_{v_{n}^{\prime}<0}|v_{n}^{\prime}|f(t,\bm{x},\bm{v}^{\prime})\mathrm{d}\bm{v}^{\prime}}{\pi{}T_{w}^{2}}\exp\left(-\frac{|\bm{v}-\bm{U}_{w}|^{2}}{T_{w}}\right),$
(7)

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

$\displaystyle\rho=\int{h}\mathrm{d}\bm{v},\quad\bm{U}=\int{\bm{v}h}\mathrm{d}\bm{v},\quad
T=\frac{2}{3}\int{|\bm{v}|^{2}}h\mathrm{d}\bm{v}-\rho,\quad p=\rho+T$ (8)
$\displaystyle\sigma_{ij}=2\int{\left(v_{i}v_{j}-\frac{|\bm{v}|^{2}}{3}\delta_{ij}\right)h}\mathrm{d}\bm{v},\quad\bm{q}=\int{\bm{v}|\bm{v}|^{2}h}\mathrm{d}\bm{v}-\frac{5}{2}\bm{U},$
(9)

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]:

$\nu_{eq}^{(k)}h^{(k+1)}+\bm{v}\cdot\frac{\partial{h}^{(k+1)}}{\partial{\bm{x}}}=L^{+}(h^{(k)},f_{eq}),$
(10)

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

$\delta_{rp}=\frac{H}{\lambda},~{}\lambda=\frac{\mu{(T_{0})}v_{m}}{n_{0}k_{B}T_{0}},$
(11)

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

$\mathrm{Kn}=\frac{\sqrt{\pi}}{2\delta_{rp}}.$ (12)

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:

$\displaystyle\frac{\partial{\rho}}{\partial{t}}+\frac{\partial{U_{i}}}{\partial{x_{i}}}=0,$
(13) $\displaystyle
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,$
$\displaystyle\frac{3}{2}\frac{\partial{T}}{\partial{t}}+\frac{\partial{q_{j}}}{\partial{x_{j}}}+\frac{\partial{U_{j}}}{\partial{x_{j}}}=0,$

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

$\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>}}},\quad
q_{i}=-\frac{5}{4\mathrm{Pr}}\delta_{rp}^{-1}\frac{\partial T}{\partial
x_{i}},$ (14)

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:

$\displaystyle\sigma_{ij}=-2\delta_{rp}^{-1}\frac{\partial
U_{<i}}{\partial{x_{j>}}}+\text{HoT}_{\sigma_{ij}},$ (15) $\displaystyle
q_{i}=-\frac{5}{4\mathrm{Pr}}\delta_{rp}^{-1}\frac{\partial T}{\partial
x_{i}}+\text{HoT}_{q_{i}},$ (16)

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

$\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}\bm{v},$
(17)

where

$L_{s}=\delta_{rp}\left\\{\left[\rho+2\bm{U}\cdot\bm{v}+T\left(|\bm{v}|^{2}-\frac{3}{2}\right)+\frac{4\left(1-\mathrm{Pr}\right)}{5}\bm{q}\cdot{\bm{v}}\left(|\bm{v}|^{2}-\frac{5}{2}\right)\right]f_{eq}-h\right\\}$
(18)

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

$\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.$ (19)

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

$\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}|\bm{v}|^{2}}\mathrm{d}\bm{v},$
(20)

where

$\text{HoT}_{q_{i}}=\frac{\partial}{\partial{x_{i}}}\int\left[(v_{i}^{2}-C_{q})\left(|\bm{v}|^{2}-\frac{3}{2}\right)-v_{i}^{2}\right]h\mathrm{d}\bm{v}+\sum_{j\neq{i}}\frac{\partial}{\partial{x_{j}}}\int{}v_{i}v_{j}\left(|\bm{v}|^{2}-\frac{5}{2}\right)h\mathrm{d}\bm{v},$
(21)

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:

  * 1.

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:

${\nu_{eq}^{(k)}}h^{(k+1/2)}+\bm{v}\cdot\frac{\partial{h}^{(k+1/2)}}{\partial{\bm{x}}}=L^{+}(h^{(k)},f_{eq}),$
(22)

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,

$\displaystyle h^{(k+1)}(\bm{x},\bm{v})=$ $\displaystyle
h^{(k+1/2)}(\bm{x},\bm{v})+\left[2\lambda_{\bm{U}}(\bm{x})\cdot{\bm{v}}+\frac{4}{5}{\lambda_{\bm{q}}}(\bm{x})\cdot\bm{v}\left(|\bm{v}|^{2}-\frac{5}{2}\right)\right]f_{eq}$
(23)
$\displaystyle+\left[\lambda_{\rho}(\bm{x})+\lambda_{T}(\bm{x})\left(|\bm{v}|^{2}-\frac{3}{2}\right)+\lambda_{\sigma_{ij}}(\bm{x})\left(v_{i}v_{j}-\frac{|\bm{v}|^{2}}{3}\delta_{ij}\right)\right]f_{eq},$

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

$h(x_{2}=0,\bm{v})=\left[1-\frac{|\bm{v}|^{2}}{2}-2\sqrt{\pi}\int_{v_{2}<0}v_{2}h(x_{2}=0,\bm{v})\mathrm{d}v_{2}\right]f_{eq},\
\ \text{when~{}}{v_{2}>0},$ (24)

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

$h(v_{1},v_{2},v_{3})=-h(v_{1},-v_{2},v_{3}),$ (25)

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

$\bm{U}=0,\quad\sigma_{ij}=0~{}\text{when}~{}i\neq{j},\quad q_{1}=q_{3}=0,$
(26)

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

$\frac{\partial T}{\partial
x_{2}}=-\frac{4\delta_{rp}}{9C_{q}}q_{2}+\underbrace{\frac{2}{3C_{q}}\int{}v_{2}|\bm{v}|^{2}(L-L_{s})\mathrm{d}\bm{v}}_{H_{1}^{(k)}(x_{2})}-\underbrace{\frac{2}{3C_{q}}\frac{\partial}{\partial
x_{2}}\int{}(v_{2}^{2}-C_{q})\left(|\bm{v}|^{2}-\frac{3}{2}\right)h\mathrm{d}\bm{v}}_{H_{2}^{(k+1/2)}(x_{2})},$
(27)

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

$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}),$
(28)

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

$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].$
(29)

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

$\rho+T+\sigma_{22}=\int{2v_{2}^{2}}hd\mathbf{v},$ (30)

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

$\sigma_{22}=-\frac{\frac{\partial}{\partial
x_{2}}\int{}2\left(v_{2}^{2}-\frac{|\bm{v}|^{2}}{3}\right)v_{2}h\mathrm{d}\bm{v}}{\delta_{rp}}+\frac{2}{\delta_{rp}}\int{(L-L_{s})v_{2}^{2}}\mathrm{d}\bm{v}.$
(31)

Figure 1: The density and temperature profiles at different iteration steps
obtained from the CIS (a, b) and GSIS (c, d), when $\delta_{rp}=50$. Circles
show the converged solution obtained from the GSIS. The linearized Shakhov
model is used with the initial condition $h(x_{2},\bm{v})=0$. The spatial
region is discretized by $N_{2}=51$ equidistant points. The iteration stops
when $\epsilon$ in Eq. (33) is less than $10^{-5}$. Data in the legends are
the iteration steps.

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]:

$v_{2}=\frac{6}{(N_{v}-1)^{\imath}}[(-N_{v}+1)^{\imath},(-N_{v}+3)^{\imath},\cdots,(N_{v}-1)^{\imath}],$
(32)

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

$\epsilon=\max\left\\{\int{}\left|\frac{\rho^{(k+1)}}{\rho^{(k)}}-1\right|\mathrm{d}x_{2},\int{}\left|\frac{T^{(k+1)}}{T^{(k)}}-1\right|\mathrm{d}x_{2},\int{}\left|\frac{q_{2}^{(k+1)}}{q_{2}^{(k)}}-1\right|\mathrm{d}x_{2}\right\\}$
(33)

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 2: The decay of the error $\epsilon$ as a function of the iteration
step, for the Fourier flow between two parallel plates described by the
linearized Shakhov model. The spatial region is discretized by $N_{2}=51$
equidistant points.

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.

Figure 3: The influence of the spatial discretization on the accuracy of both
the CIS and GSIS, for the Fourier flow between two parallel plates described
by the linearized Shakhov model with $\delta_{rp}=50$. The iteration
terminates when $\epsilon<10^{-6}$. The reference solutions (i.e. $\rho_{ref}$
and $q_{2,ref}$) are obtained from the GSIS with $N_{2}=251$, that is, the
spatial cell size is about one tenth of the mean free path of gas molecules.

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

$x_{2}=(10-15s+6s^{2})s^{3},\quad s=(0,1,\cdots,N_{s}-1)/2(N_{s}-1)$ (34)

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

Figure 4: The Knudsen layer function $T_{s}$ for the temperature profile in
the Fourier flow between two parallel plates obtained from GSIS.

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:

$T_{s}\left({x_{2}}\delta_{rp}\mathrm{Pr}\right)=\delta_{rp}\mathrm{Pr}\frac{T_{NS}(x_{2})-T(x_{2})}{k_{1}},$
(35)

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

${\zeta}_{T}=\frac{\delta_{rp}}{2}\left(\frac{5}{4\delta_{rp}\mathrm{Pr}|q_{2}|}-1\right).$
(36)

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

$\displaystyle
h\left(x_{1}=0,\bm{v}\right)=-2\sqrt{\pi}f_{eq}\int_{v_{1}<0}v_{1}h\left(x_{1}=0,\bm{v}\right)\mathrm{d}\bm{v},\quad\text{
when }v_{1}>0,$ (37) $\displaystyle
h\left(x_{1}=1,\bm{v}\right)=2\sqrt{\pi}f_{eq}\int_{v_{1}>0}v_{1}h\left(x_{1}=1,\bm{v}\right)\mathrm{d}\bm{v},\quad\text{
when }v_{1}<0,$ $\displaystyle
h\left(x_{2}=0,\bm{v}\right)=-2\sqrt{\pi}f_{eq}\int_{v_{2}<0}v_{2}h\left(x_{2}=0,\bm{v}\right)\mathrm{d}\bm{v},\quad\text{
when }v_{2}>0,$ $\displaystyle
h\left(x_{2}=1,\bm{v}\right)=\left[\sqrt{\pi}+2v_{1}+2\sqrt{\pi}\int_{v_{2}>0}v_{2}h\left(x_{2}=1,\bm{v}\right)\mathrm{d}\bm{v}\right]f_{eq},\quad\text{
when }v_{2}<0.$

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

$x_{1,2}=(10-15s+6s^{2})s^{3},\quad
s=(0,1,2,\ldots,N_{s}-1)/\left(N_{s}-1\right),$ (38)

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

$h^{(k+1)}(\bm{x},\bm{v})=h^{(k+1/2)}(\bm{x},\bm{v})+\frac{\delta_{rp}}{\max(10,\delta_{rp})}\left[\lambda_{\rho}(\bm{x})+2\lambda_{\bm{U}}(\bm{x})\cdot{\bm{v}}+\lambda_{T}(\bm{x})\left(|\bm{v}|^{2}-\frac{3}{2}\right)\right]f_{eq},$

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

$\epsilon=\iint\left|\frac{|\bm{U}^{(k+1)}|}{|\bm{U}^{(k)}|}-1\right|\mathrm{d}x_{1}\mathrm{d}x_{2}<10^{-5}.$
(39)

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.

Figure 5: The decay of error $\epsilon$ as a function of the iteration step in the lid-driven cavity flow described by the linearized Shakhov model. The iteration is assumed to be converged when $\epsilon$ defined in Eq. (39) is less than $10^{-5}$. 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

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.

Figure 6: Accuracy comparisons between the CIS and GSIS for the lid-driven
cavity flow. In each plot, the right half are the reference solution (GSIS
results on the $N_{s}=61$ grid). In the upper rows, the left halves of the
plots are CIS solutions on grids of $N_{s}=21$, 41 and 101 from left to right.
In the lower rows, the left halves are GSIS solution on grids of $N_{s}=21$,
41 and 61 from lest to right. The contour plot is pressure ($\rho+T$), with
contour levels of -0.2, -0.1, -0.05, -0.02, -0.005, 0, 0.005, 0.02, 0.05, 0.1
and 0.2.

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

$h\left(\bm{x},\bm{v}\right)=\left[2\bm{t}_{w}\cdot\bm{v}-2\sqrt{\pi}\int_{\bm{v}^{\prime}\cdot\bm{n}_{w}<0}\bm{v}^{\prime}\cdot\bm{n}_{w}h\left(\bm{x},\bm{v}^{\prime}\right)\mathrm{d}\bm{v}^{\prime}\right]f_{eq},\quad\text{when}\
\bm{v}\cdot\bm{n}_{w}>0,$ (40)

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}$.

Figure 7: Schematic of the geometry and structured triangular mesh for shear-
driven flow between two eccentric cylinders.

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.

Figure 8: Comparisons of the CIS and GSIS results for shear-driven flow
between two eccentric cylinders. (a) Contours of $U_{1}$ and streamlines at
$\delta_{rp}=1000$; (b) Contours of $U_{2}$ and streamlines at
$\delta_{rp}=1000$; (c) Contours of $U_{1}$ and streamlines at
$\delta_{rp}=10$; (d) Contours of $U_{2}$ and streamlines at $\delta_{rp}=10$.
In each sub-figures, the GSIS results are plotted in the left half domain
while the CIS ones are illustrated in the right half domain. In (a) and (b)
the velocity contours obtained by only solving the Navier-Stokes equations
with non-slip velocity boundary are also included, which are indicated by the
white dashed lines.

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]:

$f=f_{eq}(\bm{v})+A\Re\left[\exp({i{\mathrm{St}}t})h(\bm{x},\bm{v})\right],$
(41)

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

$i\mathrm{St}h+\bm{v}\cdot\frac{\partial{h}}{\partial{\bm{x}}}=L(h,f_{eq}).$
(42)

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}$

${\mathrm{St}}=\frac{\varpi{}H}{v_{m}}$ (43)

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

$S_{s}(\delta_{rp},f_{s})=\Re\left(\int\hat{h}\mathrm{d}\bm{v}\right),$ (44)

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)

$2\pi{i}(f_{s}-v_{2})\hat{h}={L}^{+}(\hat{h})-\nu_{eq}\hat{h}+f_{eq}.$ (45)

Figure 9: Comparisons of the SRBS spectrum (a) and iteration numbers (b)
between the CIS and GSIS when the rarefaction parameter is large. The HS
molecular model is used in the LBE. The molecular velocity space $[-6,6]^{3}$
is discretized by $24\times 32\times 24$ uniformly distributed points. The
solutions are believed to be converged when the relative error in $\hat{\rho}$
between two consecutive iteration steps is less than $10^{-7}$.

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:

$\hat{h}^{(k+1)}(\bm{v})=\frac{{L}^{+}(\hat{h}^{(k)})+{f_{eq}(\bm{v})}}{2\pi{}i(f_{s}-v_{2})+\nu_{eq}(\bm{v})},$
(46)

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:

$\left[\begin{array}[]{cccccc}2i\pi{f_{s}}&-2i\pi&0&0&0\\\ \vskip 6.0pt plus
2.0pt minus 2.0pt\cr-2i\pi&2i\pi{f_{s}}&-2i\pi&-2i\pi&0\\\ \vskip 6.0pt plus
2.0pt minus 2.0pt\cr 0&-2i\pi&3i\pi{f_{s}}&0&-2i\pi\\\ \vskip 6.0pt plus 2.0pt
minus 2.0pt\cr 0&-\frac{8}{3}i\pi&0&2i\pi{f_{s}}+\delta_{rp}&0\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt\cr
0&0&-3i\pi{C_{q}}&0&2i\pi{f_{s}}+\frac{2}{3}\delta_{rp}\end{array}\right]\left[\begin{array}[]{cccccc}\hat{n}^{(k+1)}\\\
\hat{U}_{2}^{(k+1)}\\\ \hat{T}^{(k+1)}\\\ \hat{\sigma}_{22}^{(k+1)}\\\
\hat{q}_{2}^{(k+1)}\end{array}\right]=\left[\begin{array}[]{cccccc}1\\\ 0\\\
0\\\ \text{R}_{4}\\\ \text{R}_{5}\end{array}\right],$ (47)

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.

Figure 10: Comparisons of the SRBS spectrum and iteration numbers (inset)
between the CIS and GSIS when the rarefaction parameter is small. The Maxwell
molecular model is used in the LBE. The molecular velocity space $[-6,6]^{3}$
is discretized by $24\times 192\times 24$ uniformly distributed points due to
high rarefaction effects. The solutions are believed to be converged when the
relative error in $\hat{\rho}$ between two consecutive iteration steps is less
than $10^{-7}$.

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

$\displaystyle\text{R}_{4}=$ $\displaystyle
2i\pi{}\text{HoT}_{\sigma_{22}}^{(k+1/2)}+2\int{}(L-L_{s})v_{2}^{2}\mathrm{d}\bm{v}+(\bar{\delta}_{rp}-\delta_{rp})\hat{\sigma}_{22}^{(k+1/2)},$
(48) $\displaystyle\text{R}_{5}=$ $\displaystyle
2i\pi{}\text{HoT}_{q_{2}}^{(k+1/2)}+\int{}(L-L_{s})v_{2}|\bm{v}|^{2}\mathrm{d}\bm{v}+\frac{2}{3}(\bar{\delta}_{rp}-\delta_{rp})\hat{q}_{2}^{(k+1/2)}.$

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.

Figure 11:  (a) The decay of the relative error
$\epsilon=|\hat{\rho}^{(k+1)}/\hat{\rho}^{(k)}-1|$ between two consecutive
iteration steps, and (b) the SRBS spectra obtained at different level of
convergence criterion. The reference solution is obtained from the GSIS when
$\epsilon=10^{-7}$. The linearized Boltzmann equation with HS molecular model
is used, with the rarefaction parameter $\delta_{rp}=50$.

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

$\left|\frac{\hat{\rho}^{(k+1)}}{\hat{\rho}^{(k)}}-1\right|=\epsilon$ (49)

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

$\left|\frac{\hat{\rho}}{\hat{\rho}^{(k+1)}}-1\right|\approx\frac{\gamma}{1-\gamma}\epsilon,$
(50)

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

$U_{w,1}=\Re\left[U_{0}\exp({i\mathrm{St}t})\right].$ (51)

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

$\displaystyle
2i\mathrm{St}U_{1}+\frac{\partial\sigma_{12}}{\partial{x_{2}}}=0,$ (52)
$\displaystyle i\mathrm{St}\sigma_{12}+\text{HoT}_{\sigma_{12}}+\frac{\partial
U_{1}}{\partial{x_{2}}}=-\delta_{rp}\sigma_{12}+2\int{(L-L_{s})}v_{1}v_{2}\mathrm{d}\bm{v},$

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:

$2i\mathrm{St}(i\mathrm{St}+\delta_{rp})U_{1}^{(k+1)}-\frac{\partial^{2}U_{1}^{(k+1)}}{\partial{x_{2}}^{2}}=\frac{\partial}{\partial{x_{2}}}\left[2\int{(L^{(k)}-L_{s}^{(k)})}v_{1}v_{2}\mathrm{d}\bm{v}-\text{HoT}_{\sigma_{12}}^{(k)}\right].$
(53)

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).

Figure 12: Comparisons of the amplitude of the shear stress exerting on the
oscillating plate and iteration numbers (inset) between the CIS and GSIS, for
the Oscillating Couette flow. The Shakhov model is solved, where the solution
is converged when
$\int{}\left|\frac{U_{1}^{(k+1)}}{U_{1}^{(k)}}-1\right|\mathrm{d}x_{2}<10^{-5}$.

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.

$\text{Kn}_{t}=\frac{\varpi}{v_{m}/\lambda}=\frac{\text{St}}{\delta_{rp}},$
(54)

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:

$\displaystyle
2i\mathrm{St}(i\mathrm{St}+\bar{\delta}_{rp})U_{1}^{(k+1)}-\frac{\partial^{2}U_{1}^{(k+1)}}{\partial{x_{2}}^{2}}=$
$\displaystyle\frac{\partial}{\partial{x_{2}}}\left[2\int{(L^{(k)}-L_{s}^{(k)})}v_{1}v_{2}\mathrm{d}\bm{v}-\text{HoT}_{\sigma_{12}}^{(k)}\right]$
(55)
$\displaystyle+2i\mathrm{St}(\bar{\delta}_{rp}-\delta_{rp})U_{1}^{(k+1/2)},$

where

$\bar{\delta}_{rp}=\max(\delta_{rp},\text{St}).$ (56)

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.

Figure 13: Comparisons of the velocity profiles in the oscillating Couette
flow with different spatial discretizations, when (a, b) $\delta_{rp}=50$ and
$\text{St}=1$ and (c, d) $\delta_{rp}=500$ and $\text{St}=0.1$. The reference
solution is obtained from the GSIS, where the spatial domain is discretized by
Eq. (38), with $N_{s}=100$ when $\delta_{rp}$=50 and $N_{s}=500$ when
$\delta_{rp}=500$.

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.

Figure 14: Comparisons of (a) the amplitude of normal pressure exerting on the
oscillating plate and iteration numbers (b) between the CIS and GSIS, for the
sound propagation problem. The Shakhov model is solved, where the solution is
converged when
$\max\left\\{\int{}\left|\frac{\rho^{(k+1)}}{\rho^{(k)}}-1\right|\mathrm{d}x_{2},\int{}\left|\frac{U_{2}^{(k+1)}}{U_{2}^{(k)}}-1\right|\mathrm{d}x_{2},\int{}\left|\frac{T^{(k+1)}}{q^{(k)}}-1\right|\mathrm{d}x_{2}\right\\}<10^{-5}$.

###  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]

$\displaystyle h(x_{2}=0,\bm{v})=$
$\displaystyle\left[\sqrt{\pi}+2v_{2}-2\sqrt{\pi}\int_{v_{2}<0}v_{2}h(x_{2}=0,\bm{v})\mathrm{d}\bm{v}\right]f_{eq},\
\ \text{when~{}}{v_{2}>0},$ (57) $\displaystyle h(x_{2}=1,\bm{v})=$
$\displaystyle
2\sqrt{\pi}f_{eq}\int_{v_{2}<0}v_{2}h(x_{2}=1,\bm{v})\mathrm{d}\bm{v},\ \
\text{when~{}}{v_{2}<0}.$

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

$\displaystyle i\text{St}\rho+\frac{\partial{U_{2}}}{\partial{x_{2}}}=0,$ (58)
$\displaystyle
2i\text{St}U_{2}+\frac{\partial{\rho}}{\partial{x_{2}}}+\frac{\partial{T}}{\partial{x_{2}}}+\frac{\partial{\sigma_{22}}}{\partial{x_{2}}}=0,$
(59)
$\displaystyle\frac{3}{2}i\text{St}T+\frac{\partial{q_{2}}}{\partial{x_{2}}}+\frac{\partial{U_{2}}}{\partial{x_{2}}}=0,$
(60) $\displaystyle
i\text{St}\sigma_{22}+\text{HoT}_{\sigma_{22}}+\frac{4}{3}\frac{\partial
U_{2}}{\partial{x_{2}}}=-\delta_{rp}\sigma_{22}+2\int{(L-L_{s})}\left(v_{2}^{2}-\frac{|\bm{v}|^{2}}{3}\right)\mathrm{d}\bm{v},$
(61) $\displaystyle
i\text{St}q_{2}+\text{HoT}_{q_{2}}+\frac{3C_{q}}{2}\frac{\partial
T}{\partial{x_{2}}}=-\frac{2}{3}\delta_{rp}q_{2}+\int{(L-L_{s})}v_{2}|\bm{v}|^{2}\mathrm{d}\bm{v}.$
(62)

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

$\displaystyle\delta_{rp}\sigma_{22}^{(k+1)}=$
$\displaystyle\bar{\delta}_{rp}\sigma_{22}^{(k+1)}+(\delta_{rp}-\bar{\delta}_{rp})\sigma_{22}^{(k+1/2)},$
(63) $\displaystyle\delta_{rp}q_{2}^{(k+1)}=$
$\displaystyle\bar{\delta}_{rp}q_{2}^{(k+1)}+(\delta_{rp}-\bar{\delta}_{rp})q_{2}^{(k+1/2)},$

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).

Figure 15: Profiles of macroscopic quantities in the sound propagation problem
with different spatial discretizations, when $\delta_{rp}=50$ and
$\text{St}=2.5$. The reference solution is obtained from the GSIS, where the
spatial domain is discretized by $N_{s}=200$ uniform grids; the corresponding
CIS results overlap with these lines and are not shown here. The normal
pressure is defined as $P_{22}=2\int{v_{2}^{2}}h\mathrm{d}\bm{v}$.

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

$\displaystyle
2i\text{St}U_{1}+\frac{\partial\sigma_{12}}{\partial{x_{2}}}+\frac{\partial\sigma_{13}}{\partial{x_{3}}}=0,$
(64) $\displaystyle
i\text{St}\sigma_{12}+\text{HoT}_{\sigma_{12}}+\frac{\partial
U_{1}}{\partial{x_{2}}}=-\delta_{rp}\sigma_{12}+2\int{(L-L_{s})}v_{1}v_{2}\mathrm{d}\bm{v},$
$\displaystyle i\text{St}\sigma_{13}+\text{HoT}_{\sigma_{13}}+\frac{\partial
U_{1}}{\partial{x_{3}}}=-\delta_{rp}\sigma_{13}+2\int{(L-L_{s})}v_{1}v_{3}\mathrm{d}\bm{v},$

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

$\displaystyle
2i\text{St}(i\text{St}+\bar{\delta}_{rp})U_{1}^{(k+1)}-\left(\frac{\partial^{2}}{\partial{x_{2}}^{2}}+\frac{\partial^{2}}{\partial{x_{3}}^{2}}\right)U_{1}^{(k+1)}=\text{Source}+2i\text{St}(\bar{\delta}_{rp}-\delta_{rp})U_{1}^{(k+1/2)},$
(65)

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

$\displaystyle\text{Source}=$
$\displaystyle\frac{\partial}{\partial{x_{2}}}\left[2\int{(L^{(k)}-L_{s}^{(k)})}v_{1}v_{2}\mathrm{d}\bm{v}-\text{HoT}^{(k+1/2)}_{\sigma_{12}}\right]$
(66)
$\displaystyle+\frac{\partial}{\partial{x_{3}}}\left[2\int{(L^{(k)}-L_{s}^{(k)})}v_{1}v_{3}\mathrm{d}\bm{v}-\text{HoT}^{(k+1/2)}_{\sigma_{13}}\right].$

Figure 16:  (a) Schematic of the oscillatory flow in a 3D rectangular cavity,
where ‘O’ is the origin of the coordinate. The top lid oscillates in the
$x_{1}$ direction periodically. (b) The amplitude of shear force exerting on
the oscillating lid that is normalized by the aspect ration
$\text{Asp}=OH/OA$, see Eq. (67); Inset shows the iteration number when the
relative error in $U_{1}$ between two consecutive iteration is less than
$10^{-5}$. (c, d, e) The distribution of shear stress when $\text{St}=0$, 10,
and 50, respectively, and $\text{Asp}=2$. (f, g, h) Same as (c, d, e),
respectively, but with $\text{Asp}=0.05$. The linearized Shakhov model is used
with $\delta_{rp}=50$ in all cases.

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

$\bar{\sigma}_{12}=\frac{2|\int_{0}^{\text{Asp}/2}\sigma_{12}(x_{2}=1)\mathrm{d}x_{3}|}{\text{Asp}},$
(67)

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:

  1. 1.

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

$\displaystyle\nabla\cdot\left[\bm{\mathcal{G}}_{\text{c}}+\bm{\mathcal{G}}_{\text{d}}\right]=0,$
(A.68) $\displaystyle\bm{L}-\nabla\bm{u}-\bm{\Pi}=0,$
$\displaystyle\bm{E}-\nabla T-\bm{\Theta}=0,$

where

$\displaystyle\bm{\mathcal{G}_{\text{c}}}=\left[\begin{array}[]{c}\bm{U}\\\
p\bm{I}\\\
\bm{0}\end{array}\right],\quad\bm{\mathcal{G}_{\text{d}}}=\left[\begin{array}[]{c}\bm{0}\\\
-\frac{1}{\delta_{rp}}\left(\bm{L}+\bm{L}^{\text{T}}-\frac{2}{3}\mathrm{tr}\left(\bm{L}\right)\bm{I}\right)\\\
-\frac{5}{4\delta_{rp}\mathrm{Pr}}\bm{E}\end{array}\right],$ (A.69)
$\displaystyle\bm{\Pi}=\left[\begin{array}[]{cc}\text{HoT}_{\sigma_{11}}+\frac{1}{2}\text{HoT}_{\sigma_{22}}&\frac{1}{2}\text{HoT}_{\sigma_{12}}\\\
\frac{1}{2}\text{HoT}_{\sigma_{12}}&\frac{1}{2}\text{HoT}_{\sigma_{11}}+\text{HoT}_{\sigma_{22}}\end{array}\right],\quad\bm{\Theta}=\left[\begin{array}[]{c}\frac{4}{5}\text{HoT}_{q_{1}}\\\
\frac{4}{5}\text{HoT}_{q_{2}}\end{array}\right]$

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

$\sigma_{ij}=-\frac{1}{\delta_{rp}}\left(L_{ij}+L_{ji}-\frac{2}{3}L_{kk}\delta_{ij}\right),\quad
q_{i}=-\frac{5}{4\delta_{rp}\mathrm{Pr}}E_{i}$ (A.70)

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

$\displaystyle\mathcal{V}=\\{\varphi:\
\varphi|_{\Delta_{i}}\in\mathcal{P}^{k}(\Delta_{i}),\ \forall\
\Delta_{i}\subset\Delta\\},$ (A.71) $\displaystyle\mathcal{W}=\\{\psi:\
\psi|_{\Gamma_{c}}\in\mathcal{P}^{k}(\Gamma_{c}),\ \forall\
\Gamma_{c}\subset\Gamma\\},$

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

$\displaystyle-\left(\bm{\mathcal{G}_{\text{c}}}+\bm{\mathcal{G}_{\text{d}}},\nabla\bm{r}\right)_{\Delta_{i}}+\langle\hat{\bm{\mathcal{F}}}\cdot\bm{n},\bm{r}\rangle_{\partial\Delta_{i}}=0,$
(A.72)
$\displaystyle\left(\bm{L},\bm{w}\right)_{\Delta_{i}}+\left(\bm{U},\nabla\cdot\bm{w}\right)_{\Delta_{i}}-\langle\hat{\bm{U}},\bm{w}\cdot\bm{n}\rangle_{\partial\Delta_{i}}=\left(\bm{w},\bm{\Pi}\right)_{\Delta_{i}}$
$\displaystyle\left(\bm{E},\bm{z}\right)_{\Delta_{i}}+\left(T,\nabla\cdot\bm{z}\right)_{\Delta_{i}}-\langle\hat{T},\bm{z}\cdot\bm{n}\rangle_{\partial\Delta_{i}}=\left(\bm{z},\bm{\Theta}\right)_{\Delta_{i}}$

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]

$\hat{\bm{\mathcal{F}}}\cdot\bm{n}=\left[\begin{array}[]{c}\bm{U}\\\
\hat{p}\bm{I}-\frac{1}{\delta_{rp}}\left(\bm{L}+\bm{L}^{\mathrm{T}}-\frac{2}{3}\mathrm{tr}\left(\bm{L}\right)\bm{I}\right)\\\
-\frac{5}{4\delta_{rp}\mathrm{Pr}}\bm{E}\par\end{array}\right]\cdot\bm{n}+\left[\begin{array}[]{ccc}\tau&&\\\
&\frac{\tau}{\delta_{rp}}&\\\
&&\frac{5\tau}{4\delta_{rp}\mathrm{Pr}}\end{array}\right]\left[\begin{array}[]{c}p-\hat{p}\\\
\bm{U}-\hat{\bm{U}}\\\ T-\hat{T}\end{array}\right].$ (A.73)

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

$\langle\bm{\left(\hat{\mathcal{F}}}\cdot\bm{n}\right)^{+},\bm{\psi}\rangle_{\Gamma_{c}}+\langle\bm{\left(\hat{\mathcal{F}}}\cdot\bm{n}\right)^{-},\bm{\psi}\rangle_{\Gamma_{c}}=0,\quad\text{on}\
\Gamma_{c}\in\Gamma\backslash\partial\Delta,$ (A.74)

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

$\langle\hat{\bm{Q}}-\bm{Q}^{+}_{VDF},\bm{\psi}\rangle_{\Gamma_{c}},\quad\text{on}\
\Gamma_{c}\in\Gamma\cap\partial\Delta,$ (A.75)

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

$\left[\begin{array}[]{cccc}A_{Q}&A_{L}&A_{E}&A_{\hat{Q}}\\\
B_{Q}&B_{L}&B_{E}&B_{\hat{Q}}\\\ C_{Q}&C_{L}&C_{E}&C_{\hat{Q}}\\\
D_{Q}&D_{L}&D_{E}&D_{\hat{Q}}\\\
\end{array}\right]\left[\begin{array}[]{c}\mathbb{Q}\\\ \mathbb{L}\\\
\mathbb{E}\\\
\hat{\mathbb{Q}}\end{array}\right]=\left[\begin{array}[]{c}S_{Q}\\\ S_{L}\\\
S_{E}\\\ S_{\hat{Q}}\end{array}\right]$ (A.76)

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.

## References

## References

  * [1] S. Chapman, T. G. Cowling, The Mathematical Theory of Non-uniform Gases, Cambridge University Press, 1970.
  * [2] G. A. Bird, Direct simulation and the Boltzmann equation, Physics of Fluids 13 (11) (1970) 2676–2681.
  * [3] V. V. Aristov, Direct Methods for Solving the Boltzmann Equation and Study of Nonequilibrium Flows, Springer Netherlands, Dordrecht, 2001.
  * [4] J. M. Burt, I. D. Boyd, A low diffusion particle method for simulating compressible inviscid flows, Journal of Computational Physics 227 (9) (2008) 4653 – 4670.
  * [5] L. Wu, J. Zhang, H. H. Liu, Y. H. Zhang, J. M. Reese, A fast iterative scheme for the linearized Boltzmann equation, J. Comput. Phys. 338 (2017) 431–451.
  * [6] L. Mieussens, Discrete velocity models and numerical schemes for the Boltzmann-BGK equation in plane and axisymmetric geometries, J. Comput. Phys. 162 (2000) 429–466.
  * [7] F. G. Tcheremissine, Direct numerical solution of the Boltzmann equation, Rar. Gas. Dynam (2005) 677–685.
  * [8] C. Mouhot, L. Pareschi, Fast algorithms for computing the Boltzmann collision operator, Math. Comput. 75 (2006) 1833–1852.
  * [9] K. Xu, J. C. Huang, An improved unified gas-kinetic scheme and the study of shock structures, IMA J. Appl. Math. 76 (2011) 698–711.
  * [10] J. C. Huang, K. Xu, P. B. Yu, A unified gas-kinetic scheme for continuum and rarefied flows II: Multi-dimensional cases, Commun. Comput. Phys. 12 (2012) 662–690.
  * [11] Z. Guo, K. Xu, R. Wang, Discrete unified gas kinetic scheme for all Knudsen number flows: Low-speed isothermal case, Phys. Rev. E 88 (2013) 033305.
  * [12] Z. Guo, R. Wang, K. Xu, Discrete unified gas kinetic scheme for all Knudsen number flows. II. Thermal compressible case, Physical Review E 91 (3) (2015) 033313.
  * [13] Y. Zhu, C. Zhong, K. Xu, Implicit unified gas-kinetic scheme for steady state solutions in all flow regimes, J. Comput. Phys. 315 (2016) 16–38.
  * [14] S. Chen, K. Xu, A comparative study of an asymptotic preserving scheme and unified gas-kinetic scheme in continuum flow limit, Journal of Computational Physics 288 (2015) 52 – 65.
  * [15] P. Wang, M. T. Ho, L. Wu, Z. Guo, Y. Zhang, A comparative study of discrete velocity methods for low-speed rarefied gas flows, Computers & Fluids 161 (2018) 33 – 46.
  * [16] R. Roveda, D. B. Goldstein, P. L. Varghese, Hybrid euler/particle approach for continuum/rarefied flows, Journal of Spacecraft and Rockets 35 (3) (1998) 258–265.
  * [17] W.-L. Wang, I. Boyd, Hybrid DSMC-CFD Simulations of Hypersonic Flow over Sharp and Blunted Bodies.
  * [18] Q. H. Sun, I. D. Boyd, G. V. Boyd, A hybrid continuum/particle approach for micro-scale gas flows, Rar. Gas. Dynam (2003) 663.
  * [19] H. S. Wijesinghe, R. D. Hornung, A. L. Garcia, N. G. Hadjiconstantinou, F. M. Hydrodynamics, Three-dimensional hybrid continuum-atomistic simulations for multiscale hydrodynamics, Journal of Fluids Engineering 126 (2004) 768 – 777\.
  * [20] T. Schwartzentruber, L. Scalabrin, I. Boyd, A modular particle–continuum numerical method for hypersonic non-equilibrium gas flows, Journal of Computational Physics 225 (1) (2007) 1159 – 1174.
  * [21] V. Kolobov, R. Arslanbekov, V. Aristov, A. Frolova, S. Zabelok, Unified solver for rarefied and continuum flows with adaptive mesh and algorithm refinement, Journal of Computational Physics 223 (2) (2007) 589 – 608.
  * [22] Z. Tang, B. He, G. Cai, Investigation on a coupled navier–stokes–direct simulation monte carlo method for the simulation of plume flowfield of a conical nozzle, International Journal for Numerical Methods in Fluids 76 (2) (2014) 95–108.
  * [23] M. L. Adams, E. W. Larsen, Fast iterative methods for discrete-ordinates particle transport calculations, Progr. Nucl. Energy 40 (2002) 3–159.
  * [24] D. Valougeorgis, S. Naris, Acceleration schemes of the discrete velocity method: Gaseous flows in rectangular microchannels, SIAM Journal on Scientific Computing 25 (2) (2003) 534–552.
  * [25] J. Lihnaropoulos, S. Naris, D. Valougeorgis, Formulation and stability analysis of rapidly convergent iteration schemes for the 2‐d linearized BGK equation, Transport Theory and Statistical Physics 36 (4-6) (2007) 513–528.
  * [26] L. Szalmás, D. Valougeorgis, A fast iterative model for discrete velocity calculations on triangular grids, Journal of Computational Physics 229 (11) (2010) 4315 – 4326.
  * [27] S. Naris, D. Valougeorgis, F. Sharipov, D. Kalempa, Discrete velocity modelling of gaseous mixture flows in mems, Superlattices and Microstructures 35 (3) (2004) 629 – 643, eurotherm 75 ’Microscale Heat Transfer 2’.
  * [28] S. Naris, D. Valougeorgis, D. Kalempa, F. Sharipov, Gaseous mixture flow between two parallel plates in the whole range of the gas rarefaction, Physica A: Statistical Mechanics and its Applications 336 (3) (2004) 294 – 318\.
  * [29] S. Naris, D. Valougeorgis, D. Kalempa, F. Sharipov, Flow of gaseous mixtures through rectangular microchannels driven by pressure, temperature, and concentration gradients, Phys. Fluids 17 (2005) 100607.
  * [30] L. Szalmás, Accelerated discrete velocity method for axial-symmetric flows of gaseous mixtures as defined by the mccormack kinetic model, Computer Physics Communications 184 (11) (2013) 2430 – 2437.
  * [31] L. Szalmás, An accelerated discrete velocity method for flows of rarefied ternary gas mixtures in long rectangular channels, Computers & Fluids 128 (2016) 91 – 97.
  * [32] W. Su, P. Wang, H. Liu, L. Wu, Accurate and efficient computation of the boltzmann equation for couette flow: Influence of intermolecular potentials on knudsen layer function and viscous slip coefficient, Journal of Computational Physics 378 (2019) 573 – 590.
  * [33] J. Fan, C. Shen., Statistical simulation of low-speed rarefied gas flows., J. Comput. Phys. 167 (2001) 393–412.
  * [34] J. Zhang, J. Fan, J. Z. Jiang, Multiple temperature model for the information preservation method and its application to nonequilibrium gas flows, J. Comput. Phys. 230 (2011) 7250–7265.
  * [35] F. Fei, J. Fan, A diffusive information preservation method for small Knudsen number flows, J. Comput. Phys. 243 (2013) 179–193.
  * [36] P. Degond, G. Dimarco, L. Pareschi, The moment guided Monte Carlo method, Int. J. Numerical Methods in Fluids 67 (2011) 189–213.
  * [37] L. Wu, J. M. Reese, Y. H. Zhang, Solving the Boltzmann equation by the fast spectral method: application to microflows, J. Fluid Mech. 746 (2014) 53–84.
  * [38] L. Wu, C. White, T. J. Scanlon, J. M. Reese, Y. H. Zhang, Deterministic numerical solutions of the Boltzmann equation using the fast spectral method, J. Comput. Phys. 250 (2013) 27–52.
  * [39] L. Wu, H. H. Liu, Y. H. Zhang, J. M. Reese, Influence of intermolecular potentials on rarefied gas flows: Fast spectral solutions of the Boltzmann equation, Phys. Fluids 27 (2015) 082002.
  * [40] T. Ohwada, Y. Sone, K. Aoki, Numerical analysis of the shear and thermal creep flows of a rarefied gas over a plane wall on the basis of the linearized Boltzmann equation for hard-sphere molecules, Phys. Fluids A 1 (1989) 1588–1599.
  * [41] W. Su, P. Wang, Y. H. Zhang, L. Wu, A high-order hybridizable discontinuous Galerkin method with fast convergence to steady-state solutions of the gas kinetic equation, J. Comput. Phys. 376 (2019) 973–991.
  * [42] W. Su, P. Wang, Y. Zhang, L. Wu, Implicit Discontinuous Galerkin Method for the Boltzmann Equation, arXiv: 1901.01865.
  * [43] D. Valougeorgis, S. Naris, Acceleration schemes of the discrete velocity method: Gaseous flows in rectangular microchannels, SIAM J. Sci. Comput. 25 (2003) 534–552.
  * [44] Y. J. Zhu, C. W. Zhong, K. Xu, Implicit unified gas-kinetic scheme for steady state solutions in all flow regimes, J. Comput. Phys. 315 (2016) 16–38.
  * [45] L. M. Yang, C. Shu, W. M. Yang, Z. Chen, H. Dong, An improved discrete velocity method (DVM) for efficient simulation of flows in all flow regimes, Phys. Fluids 30 (2018) 062005.
  * [46] L. M. Yang, Z. Chen, C. Shu, W. M. Yang, J. Wu, L. Q. Zhang, Improved fully implicit discrete-velocity method for efficient simulation of flows in all flow regimes, Phys. Rev. E 98 (2019) 063313.
  * [47] E. M. Shakhov, Approximate kinetic equations in rarefied gas theory, Fluid Dynamics 3(1) (1968) 112–115.
  * [48] H. Grad, On the kinetic theory of rarefied gases, Comm. Pure Appl. Math. 2 (1949) 331–407.
  * [49] H. Struchtrup, Macroscopic Transport Equations for Rarefied Gas Fows: Approximation Methods in Kinetic Theory, Heidelberg, Germany: Springer, 2005.
  * [50] L. Szalmás, Accelerated discrete velocity method for axial-symmetric flows of gaseous mixtures as defined by the McCormack kinetic model, Comput. Phys. Commun. 184 (2013) 2430–2437.
  * [51] L. Szalmás, D. Valougeorgis, A fast iterative model for discrete velocity calculations on triangular grids, J. Comput. Phys. 229 (2010) 4315–4326.
  * [52] W. Su, P. Wang, H. H. Liu, L. Wu, Accurate and efficient computation of the Boltzmann equation for Couette flow: Influence of intermolecular potentials on Knudsen layer function and viscous slip coefficient, J. Comput. Phys. 378 (2019) 573–590.
  * [53] W. Su, S. Lindsay, H. H. Liu, L. Wu, Comparative study of the discrete velocity and lattice Boltzmann methods for rarefied gas flows through irregular channels, Phys. Rev. E 96 (2017) 023309.
  * [54] F. Sharipov, D. Kalempa, Oscillatory Couette flow at arbitrary oscillation frequency over the whole range of the Knudsen number, Microfluidics and Nanofluidics 4 (2008) 363–374.
  * [55] D. Kalempa, F. Sharipov, Sound propagation through a rarefied gas confined between source and receptor at arbitrary Knudsen number and sound frequency, Phys. Fluids 21 (2009) 103601.
  * [56] L. Wu, J. M. Reese, Y. H. Zhang, Oscillatory rarefied gas flow inside rectangular cavities, J. Fluid Mech. 748 (2014) 350–367.
  * [57] A. Sugawara, S. Yip, L. Sirovich, Spectrum of density fluctuations in gases., Phys. Fluids 11 (1968) 925–932.
  * [58] L. Wu, C. White, T. J. Scanlon, J. M. Reese, Y. H. Zhang, A kinetic model of the Boltzmann equation for non-vibrating polyatomic gases, J. Fluid Mech. 763 (2015) 24–50.
  * [59] P. Wang, W. Su, Y. H. Zhang, Oscillatory rarefied gas flow inside a three dimensional rectangular cavity, Phys. Fluids 30 (2018) 102002.
  * [60] G. A. Radtke, N. G. Hadjiconstantinou, Variance-reduced particle simulation of the Boltzmann transport equation in the relaxation-time approximation, Phys. Rev. E 79 (2009) 056711.
  * [61] G. A. Radtke, N. G. Hadjiconstantinou, W. Wagner, Low-noise Monte Carlo simulation of the variable hard sphere gas, Phys. Fluids 23 (3) (2011) 030606\.
  * [62] D. R. Ladiges, J. E. Sader, Frequency-domain Monte Carlo method for linear oscillatory gas flows, J. Comput. Phys. 284 (2015) 351–366.
  * [63] C. Liu, Y. J. Zhu, K. Xu, Unified gas-kinetic wave-particle methods I: Continuum and rarefied gas flow, arXiv:1811.07141v1.
  * [64] Y. J. Zhu, C. Liu, C. W. Zhong, K. Xu, Unified gas-kinetic wave-particle methods II: Multiscale simulation on unstructured mesh, arXiv1903.11861v1.
  * [65] L. Wu, Y. Zhang, J. M. Reese, Fast spectral solution of the generalized Enskog equation for dense gases, J. Comput. Phys. 303 (2015) 66–79.
  * [66] L. Wu, H. H. Liu, J. M. Reese, Y. H. Zhang, Non-equilibrium dynamics of dense gas under tight confinement, J. Fluid Mech. 794 (2016) 252–266.
  * [67] A. Frezzotti, L. Gibelli, S. Lorenzani, Mean field kinetic theory description of evaporation of a fluid into vacuum, Phys. Fluids 17 (2005) 012102.
  * [68] B. Cockburn, N. C. Nguyen, J. Peraire, A comparison of hdg methods for stokes flow, Journal of Scientific Computing 45 (1) (2010) 215–237.
  * [69] J. Peraire, N. Nguyen, B. Cockburn, A Hybridizable Discontinuous Galerkin Method for the Compressible Euler and Navier-Stokes Equations, in: 48th AIAA Aerosp. Sci. Meet. Incl. New Horizons Forum Aerosp. Expo., 2010, pp. 1–11.  doi:10.2514/6.2010-363.