# A flux reconstruction kinetic scheme for the Boltzmann equation

Tianbai Xiao<EMAIL_ADDRESS>Karlsruhe Institute of Technology,
Karlsruhe, Germany

###### Abstract

It is challenging to solve the Boltzmann equation accurately due to the
extremely high dimensionality and nonlinearity. This paper addresses the idea
and implementation of the first flux reconstruction method for high-order
Boltzmann solutions. Based on the Lagrange interpolation and reconstruction,
the kinetic upwind flux functions are solved simultaneously within physical
and particle velocity space. The fast spectral method is incorporated to solve
the full Boltzmann collision integral with a general collision kernel. The
explicit singly diagonally implicit Runge-Kutta (ESDIRK) method is employed as
time integrator and the stiffness of the collision term is smoothly overcome.
Besides, we ensure the shock capturing property by introducing a self-adaptive
artificial dissipation, which is derived naturally from the effective cell
Knudsen number at the kinetic scale. As a result, the current flux
reconstruction kinetic scheme can be universally applied in all flow regimes.
Numerical experiments including wave propagation, normal shock structure, one-
dimensional Riemann problem, Couette flow and lid-driven cavity will be
presented to validate the scheme. The order of convergence of the current
scheme is clearly identified. The capability for simulating cross-scale and
non-equilibrium flow dynamics is demonstrated.

###### keywords:

Boltzmann equation, computational fluid dynamics, high-order methods, flux
reconstruction, discontinuous Galerkin

##  1 Introduction

The computational fluid dynamics (CFD) has been in voracious self-evolution in
recent decades. A highly visible direction is the development of high-order
numerical methods. In spite of the benefits from being intuitive, robust and
easy for implementation, the traditional second-order methods have proven to
be insufficiently accurate under a comparable computational cost [1]. The high
resolution and low dissipation inherited by high-order methods enable high-
fidelity simulation of intricate flows in turbulence, acoustics, plasma
physics, etc. It is more complex to implement high-order methods and they are
basically less robust than first- and second-order schemes.

High-order methods have been developed in the context of the finite difference
(FD), finite volume (FV) and finite element (FE) formulations. By extending
the difference stencils, higher-order finite difference methods can be
constructed and it is feasible to construct compact stencils [2]. However, the
straightforward extensions are restricted to problem domains with regular
geometry only [3]. The finite volume methods can handle complex geometries in
design, and a series of high-order extensions have been developed with regular
and irregular geometries [4, 5, 6, 7]. The reconstructions in FV methods are
mostly based on cell-averaged values, resulting in non-compact stencils.

The thriving finite element methods provide an alternative to design high-
order methods. The discontinuous Galerkin (DG) method is arguably one of the
most studied high-order FE algorithms [8, 9, 10, 11], which originates from
the work on neutron transport problem by Reed and Hill [12]. The basic idea of
the DG methods lies in the unified consideration of spatial discretization and
spectral decomposition. Within each element, the solutions are represented via
polynomial basis functions and are allowed to be discontinuous across cell
boundaries, which encourages the method to capture sharp wave structures that
arise in fluid mechanics. Thanks to the in-cell polynomials, it is
straightforward to extend the DG methods to arbitrarily order of accuracy for
smooth solutions. As a special case of DG methods, the nodal DG scheme employs
Lagrange polynomials as basis functions to interpolate solutions between
distinct nodal points [13]. Such idea is implemented similarly in another
class of algorithms named the spectral difference (SD) methods [14, 15], but
based on the differential form of governing equations.

Huynh’s work on the flux reconstruction (FR) approach provides profound
insight into constructing high-order methods for any advection-diffusion type
equation [16]. It establishes a general framework, where many existing
approaches such as the nodal DG and spectral difference methods can be
understood as its particular cases. Jameson used the FR formulation to prove
that the SD method is uniformly stable in a norm of Sobolev type provided that
the flux collocation points are placed at the zeros of the corresponding
Legendre polynomial [17]. The essential connections between FR and DG methods
have been analyzed in [18, 19]. A series of flux reconstruction methods have
been developed correspondingly [20, 21, 22, 23, 24]. Specifically, Vincent et
al. proposed a new class of energy stable flux reconstruction methods based on
Huynh’s approach, which is often referred as Vincent-Castonguay-Jameson-Huynh
(VCJH) schemes [20]. In what follows, we refer the terminology Flux
Reconstruction corresponding to Vincent’s formulation if unspecified.

Another hot topic in computational fluid dynamics research might go into the
study of multi-scale and non-equilibrium flow dynamics. As an example, the
Boltzmann equation provides a statistical description of particle transports
and collisions at the mesoscopic scale, i.e. the molecular mean free path and
collision time. The evolution of one-particle probability density function is
followed within the phase space. Compared to macroscopic fluid equations, the
Boltzmann equation provides many more degrees of freedom and thus can be used
to describe both equilibrium and non-equilibrium systems. Hilbert’s 6th
problem [25] served as an intriguing beginning of trying to link the behaviors
of an interacting many-particle system across different scales. It has been
shown since then that hydrodynamic equations can be recovered from the
asymptotic limits of the Boltzmann solutions [26, 27].

Continuous Efforts efforts have been devoted to the construction of numerical
solvers for the Boltzmann and its related equations [28, 29, 30, 31, 32, 33,
34]. It is challenging to solve the Boltzmann equation accurately and
efficiently due to the extremely high dimensionality and nonlinearity, and
most of the solvers above are no more than second order accuracy. The existing
attempts on constructing high-order Boltzmann solvers are very limited.
Boscheri and Dimarco [35] developed a class of central WENO implicit-explicit
Runge Kutta schemes for the simplified BGK model of the Boltzmann equation.
Jaiswal et al. [36] and Su el al. [37] developed the discontinuous Galerkin
methods for the Boltzmann equation. To the best of the author’s knowledge, no
preliminary work has been done on developing the flux reconstruction method
for the Boltzmann equation.

In this paper, a novel flux reconstruction kinetic scheme (FRKS) is presented
for the Boltzmann equation. Based on the Lagrange interpolation and
reconstruction, the kinetic upwind flux functions are solved simultaneously
within physical and particle velocity space. The fast spectral method is
incorporated into the FR framework to solve the full Boltzmann collision
integral. The explicit singly diagonally implicit Runge-Kutta (ESDIRK) method
[38] is incorporated as numerical integrator and thus the stiffness of the
collision operator in the continuum flow regime can be overcome. We ensure the
shock capturing property by introducing a self-adaptive artificial
dissipation, which is derived from the effective cell Knudsen number at the
kinetic scale. As a result, the FRKS is able to capture the cross-scale flow
dynamics where resolved and unresolved regions coexist inside a flow field.

The rest of this paper is organized as follows. Section 2 is a brief
introduction of the kinetic theory of gases. Section 3 presents the
formulation of the solution algorithm and its detailed implementation. Section
4 includes numerical experiments to demonstrate the performance of the flux
reconstruction kinetic scheme. The last section is the conclusion.

##  2 Kinetic Theory

The gas kinetic theory describes the time-space evolution of particle
distribution function $f(t,\mathbf{x},\mathbf{v})$. With a separate modeling
of particle transport and collision processes, the Boltzmann equation of
dilute monatomic gas in the absence of external force is

$\frac{\partial f}{\partial
t}+\mathbf{v}\cdot\nabla_{\mathbf{x}}f=Q(f,f)=\int_{\mathbb{R}^{3}}\int_{\mathbb{S}^{2}}\left[f\left(\mathbf{v}^{\prime}\right)f\left(\mathbf{v}_{*}^{\prime}\right)-f(\mathbf{v})f\left(\mathbf{v}_{*}\right)\right]\mathcal{B}(\cos\theta,g)d\mathbf{\Omega}d\mathbf{v}_{*},$
(1)

where $\\{\mathbf{v},\mathbf{v}_{*}\\}$ are the pre-collision velocities of
two classes of colliding particles, and
$\\{\mathbf{v}^{\prime},\mathbf{v}_{*}^{\prime}\\}$ are the corresponding
post-collision velocities. The collision kernel $\mathcal{B}(\cos\theta,g)$
measures the probability of collisions in different directions, where $\theta$
is the deflection angle and $g=|\mathbf{g}|=|\mathbf{v}-\mathbf{v}_{*}|$ is
the magnitude of relative pre-collision velocity. The solid angle
$\mathbf{\Omega}$ is the unit vector along the relative post-collision
velocity $\mathbf{v}^{\prime}-\mathbf{v}_{*}^{\prime}$, and the deflection
angle satisfies the relation $\theta=\mathbf{\Omega}\cdot\mathbf{g}/g$. With
the collision frequency defined as

$\nu(\mathbf{v})=\int_{\mathbb{R}^{3}}\int_{\mathbb{S}^{2}}f\left(\mathbf{v}_{*}\right)\mathcal{B}\left(\cos\theta,g\right)d\mathbf{\Omega}d\mathbf{v}_{*},$
(2)

The Boltzmann collision integral can be written as a combination of gain and
loss, i.e.,

$Q(f,f)=Q^{+}+Q^{-}=\int_{\mathbb{R}^{3}}\int_{\mathbb{S}^{2}}f\left(\mathbf{v}^{\prime}\right)f\left(\mathbf{v}_{*}^{\prime}\right)\mathcal{B}(\cos\theta,g)d\mathbf{\Omega}d\mathbf{v}_{*}-\nu(\mathbf{v})f(\mathbf{v}).$
(3)

A particle distribution function is related to unique macroscopic state. The
conservative flow variables can be obtained from the velocity moments of
distribution function, i.e.

$\mathbf{W}(t,\mathbf{x})=\left(\begin{array}[]{c}\rho\\\ \rho\mathbf{V}\\\
\rho E\end{array}\right):=\int_{\mathbb{R}^{3}}f\psi d\mathbf{v},$ (4)

where $\psi=[1,\mathbf{v},\mathbf{v}^{2}/2]^{T}$ is a vector of collision
invariants satisfying $\int_{\mathbb{R}^{3}}Q(f,f)\psi d\mathbf{v}=0$, and
temperature is defined as

$\frac{3}{2}kT=\frac{1}{2n}\int(\mathbf{v}-\mathbf{V})^{2}fd\mathbf{v},$ (5)

where $k$ is the Boltzmann constant and $n$ is the number density of gas.

Substituting the function $H=\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}f\log
fd\mathbf{v}d\mathbf{x}$ into the Boltzmann equation, we have

$\displaystyle\frac{dH}{dt}$
$\displaystyle=-\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}(\log
f+1)\mathbf{v}\cdot\nabla_{\mathbf{x}}fd\mathbf{v}d\mathbf{x}+\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}(\log
f+1)Q(f,f)d\mathbf{v}d\mathbf{x}$ (6)
$\displaystyle=\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\log
fQ(f,f)d\mathbf{v}d\mathbf{x}.$

From the H-theorem [26], we learn that the $H$ function is minimal only if $f$
is a Maxwellian,

$f=\mathcal{M}:=\rho\left(\frac{m}{2\pi
kT}\right)^{3/2}\exp(-\frac{m}{2kT}\left(\mathbf{v}-\mathbf{V})^{2}\right),$
(7)

where $m$ is the molecular mass.

##  3 Solution Algorithm

###  3.1 Formulation

Considering the domain $\bm{\Omega}$ with $N$ non-overlapping cells

$\bm{\Omega}=\bigcup_{i=1}^{N}\bm{\Omega}_{i},\quad\bigcap_{i=1}^{N}\bm{\Omega}_{i}=\emptyset,$
(8)

we represent the solution of the Boltzmann equation with piecewise
polynomials. Within each element $\bm{\Omega}_{i}$, the particle distribution
function is approximated by a polynomial of degree $m$ denoted $f_{i}\simeq
f_{i}^{\delta}(t,\mathbf{x},\mathbf{v})$, and the corresponding flux function
is approximated of degree $m+1$, i.e. $F_{i}\simeq
F_{i}^{\delta}(t,\mathbf{x},\mathbf{v})$. Therefore, the total approximate
solutions are

$f^{\delta}=\bigoplus_{i=1}^{N}f_{i}^{\delta}\approx f,\quad
F^{\delta}=\bigoplus_{i=1}^{N}F_{i}^{\delta}\approx F.$ (9)

For convenience, a standard coordinate can be introduce locally as
$\bm{\Omega}_{S}=\\{\mathbf{r}|\mathbf{r}\in[-1,1]^{3}\\}$. The transformation
of coordinates is made by the mapping

$\mathbf{r}=\Gamma_{i}(\mathbf{x})=\left[\begin{array}[]{c}2\left(\frac{x-x_{i-1/2}}{x_{i+1/2}-x_{i-1/2}}\right)-1\\\
2\left(\frac{y-y_{j-1/2}}{y_{j+1/2}-y_{j-1/2}}\right)-1\\\
2\left(\frac{z-z_{k-1/2}}{z_{k+1/2}-z_{k-1/2}}\right)-1\\\
\end{array}\right].$ (10)

Here we take structured mesh for illustration, while the mapping in
unstructured mesh can be found in [21, 39]. And thus the Boltzmann equation in
the local coordinate system becomes

$\frac{\partial\hat{f}^{\delta}}{\partial
t}=-\nabla_{\mathbf{r}}\cdot\hat{\mathbf{F}}^{\delta}+\hat{Q}^{\delta},$ (11)

where $\hat{\mathbf{F}}^{\delta}$ and $\hat{Q}^{\delta}$ are the numerical
flux and collision term respectively.

###  3.2 Flux

####  3.2.1 Discontinuous flux

In the flux reconstruction method, the solution is approximated by piecewise
polynomials. For brevity, let us take one-dimensional geometry as example,
while the extension to multi-dimensional case is straightforward via
tensorization. We define the following Lagrange polynomials of degree $m$

$l_{p}=\prod_{q=0,q\neq p}^{m}\left(\frac{r-r_{q}}{r_{p}-r_{q}}\right),$ (12)

and the particle distribution function can be represented on the basis of
$m+1$ solution points

$\hat{f}^{\delta}=\sum_{p=0}^{m}\hat{f}_{p}^{\delta}l_{p}.$ (13)

For the Boltzmann equation, the flux function is defined as

$F(t,x,v)=vf(t,x,v),$ (14)

and thus the local flux function is

$\hat{F}(t,r,v)=\frac{F\left(t,\Gamma_{i}^{-1}(r),v\right)}{J_{i}},$ (15)

where $J_{i}=(x_{i+1/2}-x_{i-1/2})/2$ is the Jacobian. Therefore, the flux
polynomials can be constructed as

$\hat{F}^{\delta D}=\sum_{p=0}^{m}\hat{F}_{p}^{\delta D}l_{p},$ (16)

where $\hat{F}_{p}^{\delta D}$ is the flux calculated by Eq.(15) at solution
point $r_{p}$. The notation $\delta D$ implies that such a flux is basically
discontinuous since it is derived directly from piecewise discontinuous
solutions of $\hat{f}^{\delta}$.

####  3.2.2 Interactive flux

The discontinuou flux polynomials in Eq.(15) is of the same degree of freedom
$m$ as solutions, which doesn’t meet the accuracy requirement. Besides, it
doesn’t take the information from adjacent cells into consideration and can by
no means deal with gas-surface interactions. A natural idea is to introduce a
degree $m+1$ correction flux to the approximate transformed discontinuous
flux, i.e.

$\hat{F}^{\delta}=\hat{F}^{\delta D}+\hat{F}^{\delta C}.$ (17)

The total flux is expected to equal the correct interactive fluxes at cell
boundaries, and to preserve a similar in-cell profile of discontinuous flux. A
feasible approach, as proposed in [16], is to introduce two symmetric
auxiliary functions $\\{g_{L},g_{R}\\}$, which satisfy

$\displaystyle g_{L}(r)=g_{R}(-r),$ (18) $\displaystyle g_{L}(-1)=1,\
g_{R}(-1)=0,$ $\displaystyle g_{L}(1)=0,\ g_{R}(1)=1.$

The corresponding correction flux can be reconstructed as

$\hat{F}^{\delta C}=(\hat{F}^{\delta I}_{L}-\hat{F}^{\delta
D}_{L})g_{L}+(\hat{F}^{\delta I}_{R}-\hat{F}^{\delta D}_{R})g_{R}.$ (19)

Here $\\{\hat{F}^{\delta D}_{L},\hat{F}^{\delta D}_{R}\\}$ are the
reconstructed discontinuous fluxes from the polynomial representation at the
left and right boundary of element, and $\\{\hat{F}^{\delta
I}_{L},\hat{F}^{\delta I}_{R}\\}$ are the interactive fluxes at the
boundaries. In the Boltzmann equation, we can clearly identify the flight
directions of particle transports, and the corresponding upwind flux can be
evaluated as

$\displaystyle\hat{F}^{\delta I}_{i,L}=\hat{F}^{\delta
D}_{i-1,R}H(v)+\hat{F}^{\delta D}_{i,L}H(1-v),$ (20)
$\displaystyle\hat{F}^{\delta I}_{i,R}=\hat{F}^{\delta
D}_{i,R}H(v)+\hat{F}^{\delta D}_{i+1,L}H(1-v),$

where $H(x)$ is the heaviside step function.

####  3.2.3 Total flux

Given the total flux $\hat{F}^{\delta}$, the its derivatives can be expressed
as

$\displaystyle\frac{\partial\hat{F}^{\delta}}{\partial
r}=\frac{\partial\hat{F}^{\delta D}}{\partial r}+\frac{\partial\hat{F}^{\delta
C}}{\partial r}.$ (21)

It can be evaluated by calculating the divergences of the Lagrange polynomials
and the correction functions at each solution point $r_{p}$, i.e.

$\frac{\partial\hat{F}^{\delta}}{\partial
r}\left(r_{p}\right)=\sum_{q=0}^{m}\hat{F}_{q}^{\delta
D}\frac{\mathrm{d}l_{q}}{\mathrm{~{}d}r}\left(r_{p}\right)+\left(\hat{f}_{L}^{\delta
I}-\hat{f}_{L}^{\delta
D}\right)\frac{\mathrm{d}g_{L}}{\mathrm{~{}d}r}\left(r_{p}\right)+\left(\hat{f}_{R}^{\delta
I}-\hat{f}_{R}^{\delta
D}\right)\frac{\mathrm{d}g_{R}}{\mathrm{~{}d}r}\left(r_{p}\right).$ (22)

###  3.3 Collision

It is challenging to solve the Boltzmann collision integral due to the
extremely high dimensionality and nonlinearity. The numerical Boltzmann
solvers are pioneered by Goldstein et al. [28]. The early approaches solves
the fivefold integral with Eulerian grid points and interpolations [40, 41].
Given the two-body collision model, the computational cost of these methods
are of $O(N^{6})$, where $N$ is the number of discrete velocity points in one
direction, and only half order of convergence is realized.

Another class of methods solves the Boltzmann equation with the Fourier
transform. Bobylev made a preliminary attempt of such method for the Maxwell
molecules in a homogeneous flow field [42]. In 2006, Mouhot and Pareschi
proposed a fast spectral method based on the Carleman-type Boltzmann equation
[43], with the spectral accuracy and the computational cost of
$O(M^{2}N^{3}\log N)$ [44]. Here $M$ is the number of grid points for
discretizing polar angles, which is much smaller than the number of velocity
grids $N$ in each direction. The advantageous efficiency and efficiency
guarantee its dominance in wide applications [45, 46, 36, 37, 47]. Here, we
briefly introduce the basic idea of this method.

The Carleman-type Boltzmann equation can be obtained with the following
transformations,

$\displaystyle Q(f,f)$
$\displaystyle=\int_{\mathbb{R}^{3}}\int_{\mathbb{S}^{2}}\Theta
g\left[f\left(\mathbf{v}^{\prime}\right)f\left(\mathbf{v}_{*}^{\prime}\right)-f(\mathbf{v})f\left(\mathbf{v}_{*}\right)\right]d\mathbf{\Omega}d\mathbf{v}_{*}$
(23) $\displaystyle=\int_{\mathbb{R}^{3}}\int_{\mathbb{S}^{2}}\Theta
g\left[f\left(\mathbf{v}+\frac{g\Omega-\mathbf{g}}{2}\right)f\left(\mathbf{v}_{*}-\frac{g\Omega-\mathbf{g}}{2}\right)-f(\mathbf{v})f\left(\mathbf{v}_{*}\right)\right]d\Omega
d\mathbf{v}_{*}$
$\displaystyle=2\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\Theta\delta\left(2\mathbf{y}\cdot\mathbf{g}+\mathbf{y}^{2}\right)\left[f\left(\mathbf{v}+\frac{\mathbf{y}}{2}\right)f\left(\mathbf{v}_{1}-\frac{\mathbf{y}}{2}\right)-f(\mathbf{v})f\left(\mathbf{v}_{*}\right)\right]d\mathbf{y}d\mathbf{v}_{*}$
$\displaystyle=4\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\Theta\delta\left(\mathbf{y}\cdot\mathbf{g}+\mathbf{y}^{2}\right)\left[f(\mathbf{v}+\mathbf{y})f\left(\mathbf{v}_{*}-\mathbf{y}\right)-f(\mathbf{v})f\left(\mathbf{v}_{*}\right)\right]d\mathbf{y}d\mathbf{v}_{*}$
$\displaystyle=4\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\Theta\delta(\mathbf{y}\cdot\mathbf{z})[f(\mathbf{v}+\mathbf{y})f(\mathbf{v}+\mathbf{z})-f(\mathbf{v})f(\mathbf{v}+\mathbf{y}+\mathbf{z})]d\mathbf{y}d\mathbf{z},$

where $\Theta=\mathcal{B}/g$ is the differential cross section,
$\mathbf{y}=(g\mathbf{\Omega}-\mathbf{g})/2$ and
$\mathbf{z}=\mathbf{v}_{*}-\mathbf{v}-\mathbf{y}=-\mathbf{g}-\mathbf{y}$.

As proposed in [45], a general collision kernel can be represented as

$\mathcal{B}=C_{\alpha,\gamma}^{\prime\prime}\sin^{\alpha+\gamma-1}\left(\frac{\theta}{2}\right)\cos^{-\gamma}\left(\frac{\theta}{2}\right)g^{\alpha},$
(24)

with

$C_{\alpha,\gamma}^{\prime\prime}=\frac{\Gamma[(7+\alpha)/2]}{6\Gamma[(3+\alpha+\gamma)/2]\Gamma(2-\gamma/2)}C_{\alpha},$
(25)

where $C_{\alpha}$ is the factor given in [48]. For commonly-used molecular
models, the above equation can be simplified. For example, it reduces to

$\mathcal{B}=\frac{1}{4}gd^{2}$ (26)

for the hard-sphere (HS) model.

In the fast spectral method, the particle distribution function is discretized
with $N=[N_{u},N_{v},N_{w}]^{T}$ quadrature points and periodized in a
truncated domain $\mathcal{D}=[-L,L]^{3}$. For the particle distribution
functions in a standard element $\bm{\Omega}_{S}$ of the flux reconstruction
scheme, the Fourier series can be constructed as,

$\displaystyle\hat{f}^{\delta}(t,\mathbf{r},\mathbf{v})=\sum_{k=-N/2}^{N/2-1}\hat{f}^{\delta}_{k}(t,\mathbf{r})\exp\left(i\xi_{k}\cdot\mathbf{v}\right),$
(27)
$\displaystyle\hat{f}^{\delta}_{k}=\frac{1}{(2L)^{3}}\int_{D_{L}}\hat{f}^{\delta}(t,\mathbf{r},\mathbf{v})\exp\left(-i\xi_{k}\cdot\mathbf{v}\right)d\mathbf{v},$

where $i$ is the imaginary unit and $\xi_{k}=k\pi/L$ is the frequency
component. Similarly, the collision operator can be expanded as

$\hat{Q}^{\delta}_{k}=\sum_{l,m=-N/2,(l+m=k)}^{N/2-1}\hat{f}_{l}^{\delta}\hat{f}_{m}^{\delta}[\beta(l,m)-\beta(m,m)],$
(28)

where $l=[l_{x},l_{y},l_{z}]^{T}$ and $m=[m_{x},m_{y},m_{z}]^{T}$. The kernel
mode is

$\displaystyle\beta(l,m)$
$\displaystyle=4\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\delta(\mathbf{y}\cdot\mathbf{z})\Theta\exp\left(i\xi_{l}\cdot\mathbf{y}+i\xi_{m}\cdot\mathbf{z}\right)d\mathbf{y}d\mathbf{z}$
(29)
$\displaystyle=\iint\delta\left(\mathbf{e}\cdot\mathbf{e}^{\prime}\right)\Theta\left[\int_{-R}^{R}\rho\exp\left(i\rho\xi_{l}\cdot\mathbf{e}\right)d\rho\right]\left[\int_{-R}^{R}\rho^{\prime}\exp\left(i\rho^{\prime}\xi_{m}\cdot\mathbf{e}^{\prime}\right)d\rho^{\prime}\right]d\mathbf{e}d\mathbf{e}^{\prime}$
$\displaystyle=4\int\left[\int_{0}^{R}\rho\cos\left(\rho\xi_{l}\cdot\mathbf{e}\right)d\rho\right]\left[\int\delta\left(\mathbf{e}\cdot\mathbf{e}^{\prime}\right)\int_{0}^{R}\rho^{\prime}\Theta\cos\left(\rho^{\prime}\xi_{m}\cdot\mathbf{e}^{\prime}\right)d\rho^{\prime}d\mathbf{e}^{\prime}\right]d\mathbf{e}$

with $\mathbf{e}$ and $\mathbf{e}^{\prime}$ being the vectors on the unit
sphere $\mathbb{S}^{2}$. The splitting of $l$ and $m$ kernel modes can be
achieved by numerical quadrature. Afterwards the convolution from the
Boltzmann collision integral can be evaluated efficiently with the fast
Fourier transform in the frequency domain. To avoid tedious repetition, we
refer [44] for more details of the fast spectral method.

###  3.4 Integrator

After finishing the evaluations of fluxes and collision terms, we get the time
derivatives of particle distribution function at the solution points
$\\{\mathbf{r}_{i},\mathbf{v}_{j}\\}$ from Eq.(11), i.e.

$\frac{\partial\hat{f}_{i,j}^{\delta}}{\partial
t}=-\nabla_{\mathbf{r}}\cdot\hat{\mathbf{F}}_{i,j}^{\delta}+\hat{Q}_{i,j}^{\delta}=\hat{\mathcal{C}}^{\delta}_{i,j},$
(30)

where
$\hat{\mathcal{C}}_{i,j}^{\delta}={\mathcal{C}}_{j}(\hat{f}^{\delta}_{i})$
denotes a combination of flux and collision operators. Note that
$\hat{Q}_{i,j}^{\delta}$ can become stiff in the continuum limit, when the
particle distribution function is close to the Maxwellian [31]. To circumvent
the CFL restriction, an appropriate time integrator needs to be chosen in hope
that it is efficient and A- or L-stable for stiff and oscillatory problems.

A prevailing family of integration methods for stiff differential equations is
the backward differentiation formula (BDF) thanks to its ease of
implementation [49]. As linear multi-step methods, the BDF methods with an
order greater than two cannot be A-stable. In spite of the attempts on
constructing higher-order A-stable methods by introducing additional stages
[50], these methods haven’t been proven to be universally effective and thus
the most commonly used method is BDF-2.

An alternative integrator is the multi-stage implicit Runge–Kutta (IRK)
methods [51]. In the original IRK methods, a fully coupled nonlinear system
needs to be solved at each step or each stage. To reduce the computational
complexity, the diagonally implicit Runge–Kutta (DIRK) and singly diagonally
implicit Runge–Kutta (SDIRK) methods have been proposed [38]. As a further
simplification, the explicit singly diagonally implicit Runge-Kutta (ESDIRK)
method employs an explicit first step and thus reduces the degree of the
nonlinear systems from SDIRK by one. The comparisons from compressible Navier-
Stokes equations indicated that the ESDIRK methods are more efficient than the
BDF methods [52, 53].

In this paper, the A-L stable ESDIRK method is employed to construct the flux
reconstruction scheme. We provide a brief introduction of this method, while
the comprehensive numerical implementation can be found in [54]. The SDIRK
method with $s$ stages can be written into the following form,

$\displaystyle t^{p}=t^{n}+c_{p}\Delta t^{n},$ (31)
$\displaystyle\hat{f}^{p}_{i,j}=\hat{f}^{n}_{i,j}+\Delta
t^{n}\sum_{q=1}^{p}a_{pq}\mathcal{C}_{j}\left(\hat{f}_{i}^{q}\right),\
p=1,\ldots,s,$ $\displaystyle\hat{f}^{n+1}_{i,j}=\hat{f}^{n}+\Delta
t^{n}\sum_{p=1}^{s}b_{p}\mathcal{C}_{j}\left(\hat{f}^{p}_{i}\right).$

For the ESDIRK method, the first step is explicit and thus $a_{11}=0$. The
Butcher tableau of SDIRK and ESDIRK methods is presented in Table 1, where the
$s$ or $s-1$ nonlinear equation systems to be solved are clearly identified.

Table 1: Butcher tableau of SDIRK (left) and ESDIRK (right) methods.

$c_{1}$ | $a_{11}$ | 0 | $\cdots$ | 0
---|---|---|---|---
$c_{2}$ | $a_{21}$ | $a_{22}$ | $\cdots$ | 0
$\vdots$ | $\vdots$ | $\vdots$ | $\ddots$ | $\vdots$
$c_{s}$ | $a_{s1}$ | $a_{s2}$ | $\cdots$ | $a_{ss}$
| $b_{1}$ | $b_{2}$ | $\cdots$ | $b_{s}$

$0$ | $0$ | 0 | $\cdots$ | 0
---|---|---|---|---
$c_{2}$ | $a_{21}$ | $a_{22}$ | $\cdots$ | 0
$\vdots$ | $\vdots$ | $\vdots$ | $\ddots$ | $\vdots$
$c_{s}$ | $a_{s1}$ | $a_{s2}$ | $\cdots$ | $a_{ss}$
| $b_{1}$ | $b_{2}$ | $\cdots$ | $b_{s}$

The advantage of DIRK-type integrators is that the computation of the stage
vectors is decoupled. With $N_{\mathbf{r}}$ physical solution points and
$N_{\mathbf{v}}$ velocity grid points, instead of solving one nonlinear system
with $s\times N_{\mathbf{r}}\times N_{\mathbf{v}}$ unknowns, $s$ nonlinear
systems with $N_{\mathbf{r}}\times N_{\mathbf{v}}$ unknowns are solved in
practice. The solution algorithm at each implicit stage can be written as

$\frac{\hat{f}^{p}_{i,j}-\hat{f}^{n}_{i,j}}{a_{pp}\Delta
t^{n}}=\mathcal{C}_{j}\left(\hat{f}_{i}^{p}\right)+\frac{1}{a_{pp}}\sum_{q=1}^{p-1}a_{pq}\mathcal{C}_{j}\left(\hat{f}^{q}_{i}\right),$
(32)

where the stage vectors and derivatives can be obtained via

$\displaystyle{s}^{p}=\hat{f}^{n}_{i,j}+\Delta
t^{n}\sum_{q=1}^{p-1}a_{pq}\mathcal{C}_{j}\left(\hat{f}_{i}^{q}\right),$ (33)
$\displaystyle\mathcal{C}_{j}\left(\hat{f}_{i}^{p}\right)=\frac{1}{a_{pp}\Delta
t^{n}}\left(\hat{f}^{p}_{i,j}-{s}^{p}\right).$

###  3.5 Artificial dissipation

Robust shock capturing is the critical factor for evaluating high-order
methods in hyperbolic conservation laws. In the vicinity of discontinuities in
a self-evolving flow field, oscillations tend to appear due to the Gibbs
phenomenon and lead spurious or unstable solutions. Given the less dissipation
by nature, such effects are usually more severe for higher-order methods. For
the Boltzmann equation, the shock capturing is not as strongly desired since
the shock structures can be resolved at particle mean free path level. Here we
still consider the handling of this issue as we expect to design a universal
approach that can be applied in both resolved and unresolved regions. We show
that the high dimensional information from the Boltzmann equation can be
extracted to inject more physically consistent artificial dissipation.

The issue of introducing artificial dissipation into high-order methods has
been around for a long time. The basic ideas can be categorized as follows.

Limiting: The idea is to co-opt the slope or flux limiters from finite volume
methods based on certain rules, e.g. the total variation diminishing (TVD) or
total variation bounded (TVB) principle. Preliminary work has been done in the
context of discontinuous Galerkin methods [55, 56, 57]. In principle, such
methods smear the discontinuous or sharp solutions across several adjacent
cells, which significantly diminishes the significance of introducing solution
points inside elements. Besides, a naive usage of limiters can easily lead to
descending order of accuracy around local extrema.

Artificial viscosity: An alternative way is to introduce artificial
dissipative term around discontinuous regions. The pioneer work was done in
the Jameson-Schmidt-Turkel (JST) schemes [58]. The idea has been implemented
in the discontinuous Galerkin [59] and spectral difference methods [60]. The
artifacts are expected to vanish in smooth regions and therefore the
artificial viscosity coefficients in front of the even-order derivatives need
to be solution or grid dependent [61]. It is extremely hard to set up a
viscosity that can be universally applied to different equations or
geometries. Also, the design of boundary conditions is ambiguous for the
artificial viscosity.

Filtering: Filters are commonly used in Galerkin-type methods [62, 63]. The
idea is to damp the high-order coefficients of the polynomials to eliminate
high-frequency oscillations. For the finite element type methods, the solution
inside each element is a polynomial in essence, and can be expressed
equivalently with an orthogonal polynomial basis of the same degree [64].
Different filter functions can be constructed and applied to the orthogonal
polynomials, e.g. the $L^{1}$, $L^{2}$ and the exponential filters [65]. The
filtering plays basically as a separate step in the solution algorithm and is
easy for implementation. However, if the filter is applied everywhere as
limiter in the domain, its strength needs to be very carefully chosen so that
it doesn’t destroy the key solution structure while mitigating the Gibbs
phenomenon. As a result, it is more often used locally with a detector of
trouble cells [66].

As we stand on top of the Boltzmann equation, it provides us a different point
of view to construct the artificial dissipation from the underlying kinetic
physics. Let us introduce the following dimensionless variables

$\tilde{\mathbf{x}}=\frac{\mathbf{x}}{L_{0}},\
\tilde{t}=\frac{t}{L_{0}/V_{0}},\
\tilde{\mathbf{v}}=\frac{\mathbf{v}}{V_{0}},\
\tilde{f}=\frac{f}{n_{0}V_{0}^{3}}$ (34)

where $V_{0}=\sqrt{2kT_{0}/m}$ is the most probable molecular speed, and the
Boltzmann equation becomes

$\frac{\partial\tilde{f}}{\partial\tilde{t}}+\tilde{\mathbf{v}}\cdot\nabla_{\tilde{\mathbf{x}}}{\tilde{f}}=\frac{1}{\mathrm{Kn}}Q(\tilde{f},\tilde{f}).$
(35)

The Knudsen number is defined as

$\mathrm{Kn}=\frac{V_{0}}{L_{0}\nu_{0}}=\frac{\ell_{0}}{L_{0}},$ (36)

where $\ell_{0}$ and $\nu_{0}$ are the molecular mean free path and mean
collision frequency in the reference state. For brevity, we drop the tilde
notation to denote dimensionless variables henceforth.

As calculated in [67], the thickness of a weak shock wave is around 10
molecular mean free paths, and is therefore of $O(10L_{0}\mathrm{Kn})$. When
the shock is resolved by the cell resolution, the Boltzmann equation is able
to recover the physical solution profile. Instead, if the fluid dynamics is
solved at a coarser level, the shock becomes under resolved and thus performs
as a discontinuity. The sampling theorem indicates the best numerical
solutions that can be captured under certain numerical resolution. In this
case, the physical shock thickness will be replaced by the numerical one
anyway, where the finest discontinuity length equals the distance between two
solution points. An effective numerical dissipation can be introduce following
this principle.

We modify the dimensionless Boltzmann equation as follows

$\frac{\partial f}{\partial
t}+{\mathbf{v}}\cdot\nabla_{{\mathbf{x}}}{f}=\frac{1}{\mathrm{Kn}_{c}}Q(f,f),$
(37)

where a cell Knudsen number is introduce instead of the original one. The
definition is given by

$\mathrm{Kn}_{c}=\mathrm{Kn}+\frac{\Delta x_{m}}{L_{c}}.$ (38)

Here $\Delta_{m}$ denotes the minimal distance between two adjacent solution
points with polynomials of degree $m$. A characteristic length scale of local
cell is introduced based on the gradient,

$L_{c}=\frac{\phi}{\nabla_{\mathbf{x}}\phi},$ (39)

where $\phi$ is a physical quantity of interest. Here we choose pressure as
criterion of gradient,

$p=\frac{1}{3}\int_{\mathcal{R}^{3}}(\mathbf{v}-\mathbf{V})fd\mathbf{v},$ (40)

and the evaluation of derivatives is conducted the same way as section 3.2.1.

The modified Knudsen number is related to an augmented viscosity. Let us
define the symmetric linearized operator first,

$\mathcal{L}_{g}(f):=Q(g,f)+Q(f,g),$ (41)

where $g$ is another class of particle distribution functions. Considering a
small Knudsen number $\mathrm{Kn}_{c}=\varepsilon$, we can apply the Chapman-
Enskog expansion to approximate the particle distribution function [26],

$f\simeq f_{\varepsilon}=\mathcal{M}+\sum_{n=1}^{\infty}\varepsilon^{n}g_{n},$
(42)

where $g_{n}\in\mathcal{R}(\mathcal{L}_{\mathcal{M}})$. As proved in [68],
$\mathcal{L}_{\mathcal{M}}$ is self-adjoint with respect to
$(f,g)_{\mathcal{M}}=\int_{\mathbb{R}^{3}}f(v)g(v)/\mathcal{M}(v)dv$.

Let $\mathcal{P}$ denote the projection
$\mathcal{N}(\mathcal{L}_{\mathcal{M}})$. Then applying $\mathcal{P}$ to the
Boltzmann equation leads

$\mathcal{P}\left(\frac{\partial f}{\partial
t}+\mathbf{v}\cdot\nabla_{\mathbf{x}}f\right)=\frac{1}{\varepsilon}\mathcal{P}\left(Q(f,f)\right).$
(43)

It is noticeable that
$(Q(f,f),\psi)_{\mathcal{M}}=(Q(f,f),\psi\mathcal{M})_{\mathcal{M}}=0$ holds,
so
$Q(f,f)\in\mathcal{N}\left(\mathcal{L}_{\mathcal{M}}\right)^{\perp}=\mathcal{R}\left(\mathcal{L}_{\mathcal{M}}\right)$
implies

$\mathcal{P}\left({\partial_{t}f}+\mathbf{v}\cdot\nabla_{\mathbf{x}}f\right)=0.$
(44)

Inserting the Chapman-Enskog expansion into the equation above yields

$\mathcal{P}\left(\partial_{t}+\mathbf{v}\cdot\nabla_{\mathbf{x}}\right)\mathcal{M}=-\mathcal{P}\left(\partial_{t}+\mathbf{v}\cdot\nabla_{\mathbf{x}}\right)\left(\varepsilon
g_{1}+\varepsilon^{2}g_{2}+\ldots\right).$ (45)

By matching the coefficients on the terms of order $\varepsilon$, we come to

$\left(\partial_{t}+\mathbf{v}\cdot\nabla_{\mathbf{x}}\right)\mathcal{M}=\mathcal{L}_{\mathcal{M}}\left(g_{1}\right).$
(46)

This equation has solution only if
$\mathcal{L}_{\mathcal{M}}\left(g_{1}\right)\in\mathcal{R}(\mathcal{L}_{\mathcal{M}})$.
Enforcing the projection onto $\mathcal{R}(\mathcal{L}_{\mathcal{M}})$, and
making use of the invertibility of $\mathcal{L}_{\mathcal{M}}$, we get

$\displaystyle\mathcal{P}\left(\partial_{t}+\mathbf{v}\cdot\nabla_{\mathbf{x}}\right)\mathcal{M}$
$\displaystyle=-\varepsilon\mathcal{P}\left(\partial_{t}+\mathbf{v}\cdot\nabla_{\mathbf{x}}\right)\mathcal{L}_{\mathcal{M}}^{-1}(I-\mathcal{P})\left(\partial_{t}+\mathbf{v}\cdot\nabla_{\mathbf{x}}\right)\mathcal{M}$
(47)
$\displaystyle=-\varepsilon\mathcal{P}\left(\mathbf{v}\cdot\nabla_{\mathbf{x}}\right)L_{\mathcal{M}}^{-1}(I-\mathcal{P})\left(\mathbf{v}\cdot\nabla_{\mathbf{x}}\right)\mathcal{M},$

which is the compact form of the Navier-Stokes equations. As is shown, the
augmented Knudsen number plays an equivalent role as artificial viscosity at
the Navier-Stokes level.

###  3.6 Summary

The flowchart of the current solution algorithm is summarized in Fig. 1.

Figure 1: Flowchart of solution algorithm.

##  4 Numerical Experiments

In this section, we will conduct numerical experiments to validate the current
scheme. In order to demonstrate the cross-scale computing capability of the
algorithm, the results at different degrees of gas rarefaction are presented.
As introduced in section 3.5, the dimensionless variables are used in all the
numerical simulations.

###  4.1 Wave propagation

First we study the order of accuracy of the flux reconstruction kinetic
scheme. The propagation of an one-dimensional traveling wave is used as the
test case. The initial particle distribution function is set as Maxwellian in
correspondence with the following macroscopic variables

$\left[\begin{array}[]{c}\rho\\\ U\\\ V\\\ W\\\ p\\\
\end{array}\right]_{{t=0}}=\left[\begin{array}[]{c}1+\alpha\sin(2\pi x)\\\
1\\\ 0\\\ 0\\\ 0.5\\\ \end{array}\right],$

and the detailed computational setup is presented in Table 2.

Table 2: Computational setup of wave propagation problem. $t$ | $x$ | $N_{x}$ | Polynomial | Degree | Points | Correction
---|---|---|---|---|---|---
$(0,1]$ | $[0,1]$ | $[5,40]$ | Lagrange | $[2,3]$ | Legendre | Radau
$\mathbf{v}$ | $N_{u}$ | $N_{v}$ | $N_{w}$ | Quadrature | Kn | $\alpha$
$[-8,8]^{3}$ | 80 | 28 | 28 | Rectangular | $[0.0001,0.1]$ | 0.1
Integrator | Boundary | CFL |  |  |  |
ESDIRK | Periodic | 0.1 |  |  |  |

As listed, the Lagrange polynomials of degree $2$ and $3$ are used in the
computation, resulting in third and fourth order of accuracy by design. Fig. 2
shows the traveling wave solutions with $N_{x}=20$ and polynomial degree $3$
at different reference Knudsen numbers. The reference solutions are produced
by the fast spectral method [45] with 1000 cells. With the increasing
molecular mean free path, the enhanced particle transports result in stronger
viscous dissipation and the smeared wave structure. Table 3 to 10 provide the
convergence orders of the flux reconstruction kinetic schemes in design of the
third (FRKS3) and fourth order of accuracy (FRKS4). It is clear that the
current method preserves the desired accuracy in all Knudsen regimes.

Table 3: Errors and convergences of FRKS3 in the wave propagation problem at $\mathrm{Kn}=0.0001$. $\Delta x$ |  $L^{1}$ error | Order |  $L^{2}$ error | Order |  $L^{\infty}$ error | Order
---|---|---|---|---|---|---
0.2 | 1.688797E-3 |  | 5.120427E-4 |  | 2.215626E-4 |
0.1 | 2.678810E-4 | 2.66 | 5.438995E-5 | 3.23 | 1.593844E-5 | 3.79
0.05 | 3.263045E-5 | 3.04 | 4.729391E-6 | 3.52 | 1.014085E-6 | 3.97
0.025 | 4.045788E-6 | 3.01 | 4.135537E-7 | 3.51 | 6.213239E-8 | 4.03
Table 4: Errors and convergences of FRKS3 in the wave propagation problem at $\mathrm{Kn}=0.001$. $\Delta x$ |  $L^{1}$ error | Order |  $L^{2}$ error | Order |  $L^{\infty}$ error | Order
---|---|---|---|---|---|---
0.2 | 1.657782E-3 |  | 4.984405E-4 |  | 2.105969E-4 |
0.1 | 2.105969E-4 | 2.98 | 2.105969-5 | 4.56 | 1.499155E-5 | 3.81
0.05 | 2.621523E-5 | 3.01 | 4.315510E-6 | 2.29 | 9.015856E-7 | 4.06
0.025 | 3.276879-6 | 3.00 | 3.657462E-7 | 3.56 | 5.422131E-8 | 4.06
Table 5: Errors and convergences of FRKS3 in the wave propagation problem at $\mathrm{Kn}=0.01$. $\Delta x$ |  $L^{1}$ error | Order |  $L^{2}$ error | Order |  $L^{\infty}$ error | Order
---|---|---|---|---|---|---
0.2 | 1.419118E-3 |  | 4.085279E-4 |  | 1.671561E-4 |
0.1 | 1.785668E-4 | 2.99 | 3.648516E-5 | 3.49 | 1.069315E-5 | 3.97
0.05 | 2.141419E-5 | 3.06 | 3.119967E-6 | 3.55 | 6.628941E-7 | 4.01
0.025 | 2.579123E-6 | 3.05 | 2.655507E-7 | 3.55 | 4.047194E-8 | 4.03
Table 6: Errors and convergences of FRKS3 in the wave propagation problem at $\mathrm{Kn}=0.1$. $\Delta x$ |  $L^{1}$ error | Order |  $L^{2}$ error | Order |  $L^{\infty}$ error | Order
---|---|---|---|---|---|---
0.2 | 3.939578E-4 |  | 1.146447E-4 |  | 4.812587E-5 |
0.1 | 4.804452E-5 | 3.04 | 9.801638E-6 | 3.55 | 2.943710E-6 | 4.03
0.05 | 5.964627E-6 | 3.01 | 8.636624E-7 | 3.50 | 1.827464E-7 | 4.01
0.025 | 7.427253E-7 | 3.01 | 7.592175E-8 | 3.51 | 1.140249E-8 | 4.00
Table 7: Errors and convergences of FRKS4 in the wave propagation problem at $\mathrm{Kn}=0.0001$. $\Delta x$ |  $L^{1}$ error | Order |  $L^{2}$ error | Order |  $L^{\infty}$ error | Order
---|---|---|---|---|---|---
0.2 | 1.677265E-4 |  | 4.413744E-5 |  | 1.877038E-5 |
0.1 | 1.027635E-5 | 4.03 | 1.826210E-6 | 4.60 | 4.743682E-7 | 5.31
0.05 | 6.395205E-7 | 4.01 | 8.050237E-8 | 4.50 | 1.508162E-8 | 4.98
0.025 | 4.122186E-8 | 3.96 | 3.654658E-9 | 4.46 | 4.724505E-10 | 5.00
Table 8: Errors and convergences of FRKS4 in the wave propagation problem at $\mathrm{Kn}=0.001$. $\Delta x$ |  $L^{1}$ error | Order |  $L^{2}$ error | Order |  $L^{\infty}$ error | Order
---|---|---|---|---|---|---
0.2 | 1.647013E-4 |  | 4.236274E-5 |  | 1.668588E-5 |
0.1 | 1.013567E-5 | 4.02 | 1.788745-6 | 4.57 | 4.867548E-7 | 5.10
0.05 | 6.243648E-7 | 4.02 | 7.876072E-8 | 4.51 | 1.501070E-8 | 5.02
0.025 | 3.962609-8 | 3.98 | 3.528171E-9 | 4.48 | 4.704976E-10 | 5.00
Table 9: Errors and convergences of FRKS4 in the wave propagation problem at $\mathrm{Kn}=0.01$. $\Delta x$ |  $L^{1}$ error | Order |  $L^{2}$ error | Order |  $L^{\infty}$ error | Order
---|---|---|---|---|---|---
0.2 | 1.444726E-4 |  | 3.450912E-5 |  | 1.172881E-5 |
0.1 | 8.303901E-6 | 4.12 | 1.491799E-6 | 4.53 | 4.560153E-7 | 4.68
0.05 | 5.107348E-7 | 4.02 | 6.450361E-8 | 4.53 | 1.445434E-8 | 4.98
0.025 | 3.210330E-8 | 3.99 | 2.853924E-9 | 4.50 | 4.355419E-10 | 5.05
Table 10: Errors and convergences of FRKS4 in the wave propagation problem at $\mathrm{Kn}=0.1$. $\Delta x$ |  $L^{1}$ error | Order |  $L^{2}$ error | Order |  $L^{\infty}$ error | Order
---|---|---|---|---|---|---
0.2 | 3.957461E-5 |  | 9.537192E-6 |  | 3.142089E-6 |
0.1 | 2.353530E-6 | 4.07 | 4.209624E-7 | 4.50 | 1.209994E-7 | 4.70
0.05 | 1.470697E-7 | 4.00 | 1.853411E-8 | 4.51 | 3.762843E-9 | 5.01
0.025 | 9.404308E-9 | 3.97 | 8.354799E-10 | 4.47 | 1.189082E-10 | 4.98

###  4.2 Normal shock structure

We continue considering a well-resolved problem, i.e. the normal shock wave
structure. The initial particle distribution function is set as Maxwellian in
correspondence with the following macroscopic variables

$\left[\begin{array}[]{c}\rho\\\ U\\\ V\\\ W\\\ T\\\
\end{array}\right]_{{t=0,L}}=\left[\begin{array}[]{c}\rho_{-}\\\ U_{-}\\\ 0\\\
0\\\ T_{-}\\\ \end{array}\right],\quad\left[\begin{array}[]{c}\rho\\\ U\\\
V\\\ W\\\ T\\\
\end{array}\right]_{{t=0,R}}=\left[\begin{array}[]{c}\rho_{+}\\\ U_{+}\\\ 0\\\
0\\\ T_{+}\\\ \end{array}\right].$

Based on the reference frame of shock wave, the upstream and downstream gases
are related with the well-known Rankine-Hugoniot relation,

$\displaystyle\frac{\rho_{+}}{\rho_{-}}$
$\displaystyle=\frac{(\gamma+1)\mathrm{Ma}^{2}}{(\gamma-1)\mathrm{Ma}^{2}+2},$
(48) $\displaystyle\frac{U_{+}}{U_{-}}$
$\displaystyle=\frac{(\gamma-1)\mathrm{Ma}^{2}+2}{(\gamma+1)\mathrm{Ma}^{2}},$
$\displaystyle\frac{T_{+}}{T_{-}}$
$\displaystyle=\frac{\left((\gamma-1)\mathrm{Ma}^{2}+2\right)\left(2\gamma\mathrm{Ma}^{2}-\gamma+1\right)}{(\gamma+1)^{2}\mathrm{Ma}^{2}},$

where $\rm Ma$ is the upstream Mach number, and $\gamma=5/3$ is the specific
heat ratio of monatomic molecule. The reference state is set with the upstream
flow conditions. As the shock profile is resolved in this case, no stiffness
will be introduced and thus we employ the Bogacki-Shampine integrator, which
is a third-order explicit Runge-Kutta method [69]. The detailed computation
setup is presented in Table 11.

Fig. 3 provides the profiles of density, $U$-velocity and temperature at
different upstream Mach numbers. The reference solutions are produced by the
fast spectral method [45] with 200 cells. As is shown, excellent agreement has
been achieved between the flux reconstruction solutions under a coarse mesh
and the reference results. It demonstrates the capability of the current
scheme to simulate the evolution of non-equilibrium particle distributions.

Table 11: Computational setup of normal shock structure. $x$ | $N_{x}$ | Polynomial | Degree | Points | Correction
---|---|---|---|---|---
$[-25,25]$ | $50$ | Lagrange | $2$ | Legendre | Radau
$\mathbf{v}$ | $N_{u}$ | $N_{v}$ | $N_{w}$ | Quadrature | Kn
$[-14,14]^{3}$ | 64 | 32 | 32 | Rectangular | $1.0$
Ma | Integrator | Boundary | CFL |  |
$[2,3]$ | Bogacki–Shampine | Dirichlet | 0.2 |  |

###  4.3 Riemann problem

Now we shift our gaze to the problem where resolved and unresolved regions
coexist in the solution domain. We employ the Sod shock tube problem, which is
a standard one-dimensional Riemann problem. The particle distribution function
is initialized as Maxwellian, which corresponds to the following macroscopic
variables

$\left[\begin{array}[]{c}\rho\\\ U\\\ V\\\ W\\\ p\\\
\end{array}\right]_{{t=0,L}}=\left[\begin{array}[]{c}1\\\ 0\\\ 0\\\ 0\\\
0.5\\\ \end{array}\right],\quad\left[\begin{array}[]{c}\rho\\\ U\\\ V\\\ W\\\
p\\\ \end{array}\right]_{{t=0,R}}=\left[\begin{array}[]{c}0.125\\\ 0\\\ 0\\\
0\\\ 0.1\\\ \end{array}\right].$

To test the capability of the current scheme to solve resolved/unresolved wave
structures and the corresponding multi-scale performance, simulations are
performed with different reference Knudsen numbers $\mathrm{Kn}=[0.0001,1]$,
with respect to typical continuum, transition, and free molecular flow
regimes. The detailed computation setup is listed in Table 12.

Table 12: Computational setup of Sod shock tube. $t$ | $x$ | $N_{x}$ | Polynomial | Degree | Points | Correction
---|---|---|---|---|---|---
$[0,0.15]$ | $[0,1]$ | $50$ | Lagrange | $2$ | Legendre | Radau
$\mathbf{v}$ | $N_{u}$ | $N_{v}$ | $N_{w}$ | Quadrature | Kn | CFL
$[-8,8]^{3}$ | 64 | 32 | 32 | Rectangular | $[0.0001,1]$ | $0.15$
Integrator | Boundary |  |  |  |  |
ESDIRK-3 | Dirichlet |  |  |  |  |

Fig. 4 presents the profiles of density, $U$-velocity and temperature inside
the shock tube at the output instant $t=0.2$ under different Knudsen numbers.
The reference solutions are derived from the Euler and collisionless Boltzmann
equations. In the continuum regime with $\mathrm{Kn}=0.0001$, the molecular
mean free path is much less than the grid size, and thus the current method
becomes a shock capturing method under limited resolution in space and time.
As shown in Fig. 4a, oscillatory solutions from the original flux
reconstruction method emerge around the shock wave front due to the Gibbs
phenomenon. Conversely, the adaptive artificial dissipation introduced in the
current scheme eliminates the oscillations effectively while preserving the
sharp wave structures. With increasing Knudsen number and molecular mean fee
path, the enhanced transport phenomena widen the waves and reduce the
gradients of characteristic variables. Therefore, the solution profiles
becomes resolvable under the current resolution. From $\mathrm{Kn}=0.0001$ to
$\mathrm{Kn}=1$, the artificial dissipation doesn’t destroy the solutions from
the current method, and a smooth transition is recovered from the Euler
solutions of Riemann problem to the collisionless Boltzmann solutions.

###  4.4 Couette flow

The former cases consider only periodic or Dirichlet boundary conditions. In
this case, we employ the Couette flow as an example to test the gas-surface
interactions in the flux reconstruction kinetic scheme. The initial particle
distribution is set as Maxwellian based on the homogeneous fluids,

$\left[\begin{array}[]{c}\rho\\\ U\\\ V\\\ W\\\ T\\\
\end{array}\right]_{{t=0}}=\left[\begin{array}[]{c}1\\\ 0\\\ 0\\\ 0\\\ 1\\\
\end{array}\right].$

The boundary temperature at both ends of the domain are set as $T_{w}=1$, and
the velocities differ as $\mathbf{V}_{wL}=[0,-1,0]^{T}$,
$\mathbf{V}_{wR}=[0,1,0]^{T}$. Maxwell’s diffusive boundary is adopted to
model the gas-surface interaction. The detailed computational setup is
recorded in Table 13.

Table 13: Computational setup of Couette flow. $x$ | $N_{x}$ | Polynomial | Degree | Points | Correction
---|---|---|---|---|---
$[-1,1]$ | $30$ | Lagrange | $2$ | Legendre | Radau
$\mathbf{v}$ | $N_{u}$ | $N_{v}$ | $N_{w}$ | Quadrature | Kn
$[-8,8]^{3}$ | 72 | 72 | 28 | Rectangular | $[0.2/\sqrt{\pi},20/\sqrt{\pi}]$
Integrator | Boundary | CFL |  |  |
Bogacki–Shampine | Maxwell | 0.15 |  |  |

Fig. 5a shows the macroscopic $V$-velocity profiles in the transition regimes
with three Knudsen numbers
$\mathrm{Kn}=\\{0.2/\sqrt{\pi},2/\sqrt{\pi},20/\sqrt{\pi}\\}$. The current
numerical solutions agree perfectly with the reference solutions, which are
produced by the information-preserving DSMC method [70]. Fig. 5b draws the
relation of surface shear stress versus Knudsen number, where the
collisionless solution is used to determine the normalization factor
$\tau_{0}$. It is clear that the current solutions fall exactly on the
linearized Boltzmann solutions [71] across different Knudsen regimes.

###  4.5 Lid-driven cavity

In the last case, we test the current scheme with multi-dimensional geometry.
The lid-driven cavity is employed as the test problem. The initial particle
distribution is set as Maxwellian with the homogeneous fluids,

$\left[\begin{array}[]{c}\rho\\\ U\\\ V\\\ W\\\ T\\\
\end{array}\right]_{{t=0}}=\left[\begin{array}[]{c}1\\\ 0\\\ 0\\\ 0\\\ 1\\\
\end{array}\right].$

The solution domain is enclosed by four solid walls with $T_{w}=1$. The upper
wall moves in the tangent direction with $\mathbf{V}_{w}=[0.15,0,0]^{T}$, and
the rest three walls are kept still. Maxwell’s diffusive boundary is adopted
to all the walls. The detailed computational setup is provided in Table 14.

Table 14: Computational setup of lid-driven cavity. $x$ | $y$ | $N_{x}$ | $N_{x}$ | Polynomial | Degree
---|---|---|---|---|---
$[0,1]$ | $[0,1]$ | $15$ | 15 | Lagrange | $2$
Points | Correction | $\mathbf{v}$ | $N_{u}$ | $N_{v}$ | $N_{w}$
Legendre | Radau | $[-8,8]^{3}$ | 32 | 32 | 28
Quadrature | Kn | Integrator | Boundary | CFL |
Rectangular | $[0.2/\sqrt{\pi},20/\sqrt{\pi}]$ | Bogacki–Shampine | Maxwell | $0.15$ |

Fig. 6 shows the contours of $U$-velocity with streamlines and temperature
with heat flux vectors inside the cavity. As explained in [72], the anti-
Fourier’s heat flux driven by stress is clearly identified. Fig. 7 the
velocity profiles along the vertical and horizontal central lines of the
cavity. The DSMC solutions with $60\times 60$ physical mesh are plotted for
comparison. The quantitative comparison demonstrates that the current scheme
is able to provide equivalent DSMC solutions in the transition regime with
much coarser mesh.

##  5 Conclusion

Non-equilibrium statistical mechanics is profoundly built upon the Boltzmann
equation. For the first time, a high-order kinetic scheme based on flux
reconstruction is proposed for solving the Boltzmann equation in this paper.
The upwind flux solver is integrated with flux reconstruction formulation
seamlessly throughout the phase space. The fast spectral method is constructed
to solve the exact Boltzmann collision integral with an arbitrary collision
kernel. Besides, the explicit singly diagonally implicit Runge-Kutta method
ensures the compatible accuracy in time direction and overcome the stiffness
of collision term in the continuum flow regime. The current method provides an
accurate and efficient tool for the study of cross-scale and non-equilibrium
flow phenomena. It shows the potential to be extended to other complex
systems, e.g. astrophysics [73], plasma physics [74], uncertainty
quantification [75], etc.

## References

  * [1] Zhi Jian Wang, Krzysztof Fidkowski, Rémi Abgrall, Francesco Bassi, Doru Caraeni, Andrew Cary, Herman Deconinck, Ralf Hartmann, Koen Hillewaert, Hung T Huynh, et al.  High-order CFD methods: current status and perspective.  International Journal for Numerical Methods in Fluids, 72(8):811–845, 2013.
  * [2] Sanjiva K Lele.  Compact finite difference schemes with spectral-like resolution.  Journal of computational physics, 103(1):16–42, 1992.
  * [3] Tadeusz Liszka and Janusz Orkisz.  The finite difference method at arbitrary irregular grids and its application in applied mechanics.  Computers & Structures, 11(1-2):83–95, 1980.
  * [4] Ami Harten, Bjorn Engquist, Stanley Osher, and Sukumar R Chakravarthy.  Uniformly high order accurate essentially non-oscillatory schemes, III.  In Upwind and high-resolution schemes, pages 218–290. Springer, 1987.
  * [5] Rémi Abgrall.  On essentially non-oscillatory schemes on unstructured meshes: analysis and implementation.  Journal of Computational Physics, 114(1):45–58, 1994.
  * [6] Xu-Dong Liu, Stanley Osher, and Tony Chan.  Weighted essentially non-oscillatory schemes.  Journal of computational physics, 115(1):200–212, 1994.
  * [7] Guang-Shan Jiang and Chi-Wang Shu.  Efficient implementation of weighted ENO schemes.  Journal of computational physics, 126(1):202–228, 1996.
  * [8] Bernardo Cockburn and Chi-Wang Shu.  Runge–Kutta discontinuous Galerkin methods for convection-dominated problems.  Journal of scientific computing, 16(3):173–261, 2001.
  * [9] Bernardo Cockburn and Chi-Wang Shu.  The local discontinuous Galerkin method for time-dependent convection-diffusion systems.  SIAM Journal on Numerical Analysis, 35(6):2440–2463, 1998.
  * [10] Jaime Peraire and P-O Persson.  The compact discontinuous Galerkin (CDG) method for elliptic problems.  SIAM Journal on Scientific Computing, 30(4):1806–1824, 2008.
  * [11] Hong Luo, Luqing Luo, Robert Nourgaliev, Vincent A Mousseau, and Nam Dinh.  A reconstructed discontinuous galerkin method for the compressible navier–stokes equations on arbitrary grids.  Journal of Computational Physics, 229(19):6961–6978, 2010.
  * [12] William H Reed and TR Hill.  Triangular mesh methods for the neutron transport equation.  Technical report, Los Alamos Scientific Lab., N. Mex.(USA), 1973.
  * [13] Jan S Hesthaven and Tim Warburton.  Nodal discontinuous Galerkin methods: algorithms, analysis, and applications.  Springer Science & Business Media, 2007.
  * [14] David A Kopriva and John H Kolias.  A conservative staggered-grid Chebyshev multidomain method for compressible flows.  Journal of computational physics, 125(1):244–261, 1996.
  * [15] Yen Liu, Marcel Vinokur, and Zhi Jian Wang.  Spectral difference method for unstructured grids I: Basic formulation.  Journal of Computational Physics, 216(2):780–801, 2006.
  * [16] Hung T Huynh.  A flux reconstruction approach to high-order schemes including discontinuous Galerkin methods.  In 18th AIAA Computational Fluid Dynamics Conference, page 4079, 2007.
  * [17] Antony Jameson.  A proof of the stability of the spectral difference method for all orders of accuracy.  Journal of Scientific Computing, 45(1):348–358, 2010.
  * [18] Y Allaneau and Antony Jameson.  Connections between the filtered discontinuous Galerkin method and the flux reconstruction approach to high order discretizations.  Computer Methods in Applied Mechanics and Engineering, 200(49-52):3628–3636, 2011.
  * [19] D De Grazia, G Mengaldo, D Moxey, PE Vincent, and SJ Sherwin.  Connections between the discontinuous Galerkin method and high-order flux reconstruction schemes.  International journal for numerical methods in fluids, 75(12):860–877, 2014.
  * [20] Peter E Vincent, Patrice Castonguay, and Antony Jameson.  A new class of high-order energy stable flux reconstruction schemes.  Journal of Scientific Computing, 47(1):50–72, 2011.
  * [21] Patrice Castonguay, Peter E Vincent, and Antony Jameson.  A new class of high-order energy stable flux reconstruction schemes for triangular elements.  Journal of Scientific Computing, 51(1):224–256, 2012.
  * [22] Kartikey Asthana and Antony Jameson.  High-order flux reconstruction schemes with minimal dispersion and dissipation.  Journal of Scientific Computing, 62(3):913–944, 2015.
  * [23] Ray Vandenhoeck and Andrea Lani.  Implicit high-order flux reconstruction solver for high-speed compressible flows.  Computer Physics Communications, 242:1–24, 2019.
  * [24] Ji Li, Chengwen Zhong, and Sha Liu.  High-order kinetic flow solver based on the flux reconstruction framework.  Physical Review E, 102(4):043306, 2020.
  * [25] David Hilbert.  Mathematical problems.  Bulletin of the American Mathematical Society, 8(10):437–479, 1902\.
  * [26] Sydney Chapman, Thomas George Cowling, and David Burnett.  The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases.  Cambridge university press, 1990.
  * [27] Harold Grad.  On the kinetic theory of rarefied gases.  Communications on pure and applied mathematics, 2(4):331–407, 1949\.
  * [28] D Goldstein, B Sturtevant, and JE Broadwell.  Investigations of the motion of discrete-velocity gases.  Progress in Astronautics and Aeronautics, 117:100–117, 1989.
  * [29] A Bobylev and S Rjasanow.  Difference scheme for the boltzmann equation based on the fast fourier transform.  European Journal of Mechanics-B/Fluids, 16(2):293–306, 1997.
  * [30] Luc Mieussens.  Discrete-velocity models and numerical schemes for the Boltzmann-BGK equation in plane and axisymmetric geometries.  Journal of Computational Physics, 162(2):429–466, 2000.
  * [31] Francis Filbet and Shi Jin.  A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources.  Journal of Computational Physics, 229(20):7625–7648, 2010.
  * [32] Kun Xu and Juan-Chen Huang.  A unified gas-kinetic scheme for continuum and rarefied flows.  Journal of Computational Physics, 229(20):7747–7764, 2010.
  * [33] Tianbai Xiao, Qingdong Cai, and Kun Xu.  A well-balanced unified gas-kinetic scheme for multiscale flow transport under gravitational field.  Journal of Computational Physics, 332:475–491, 2017.
  * [34] Tianbai Xiao, Kun Xu, and Qingdong Cai.  A unified gas-kinetic scheme for multiscale and multicomponent flow transport.  Applied Mathematics and Mechanics, 40(3):355–372, 2019.
  * [35] Walter Boscheri and Giacomo Dimarco.  High order central WENO-Implicit-Explicit Runge Kutta schemes for the BGK model on general polygonal meshes.  Journal of Computational Physics, 422:109766, 2020.
  * [36] Shashank Jaiswal, Alina A Alexeenko, and Jingwei Hu.  A discontinuous Galerkin fast spectral method for the full Boltzmann equation with general collision kernels.  Journal of Computational Physics, 378:178–208, 2019.
  * [37] Wei Su, Peng Wang, Yonghao Zhang, and Lei Wu.  Implicit discontinuous Galerkin method for the Boltzmann equation.  Journal of Scientific Computing, 82(2):1–35, 2020.
  * [38] Christopher Alan Kennedy and Mark Huitt Carpenter.  Diagonally Implicit Runge-Kutta Methods for Ordinary Differential Equations, a Review.  National Aeronautics and Space Administration, Langley Research Center, 2016.
  * [39] Freddie D Witherden, Antony M Farrington, and Peter E Vincent.  PyFR: An open source framework for solving advection–diffusion type problems on streaming architectures using the flux reconstruction approach.  Computer Physics Communications, 185(11):3028–3040, 2014.
  * [40] Yoshio Sone, Taku Ohwada, and Kazuo Aoki.  Temperature jump and Knudsen layer in a rarefied gas over a plane wall: Numerical analysis of the linearized Boltzmann equation for hard-sphere molecules.  Physics of Fluids A: Fluid Dynamics, 1(2):363–370, 1989.
  * [41] Taku Ohwada.  Structure of normal shock waves: Direct numerical analysis of the Boltzmann equation for hard-sphere molecules.  Physics of Fluids A: Fluid Dynamics, 5(1):217–234, 1993.
  * [42] AV Bobylev.  The theory of the nonlinear spatially uniform Boltzmann equation for Maxwell molecules.  Mathematical physics reviews, 7:111–233, 1988.
  * [43] T Carleman.  L’Intégrale de Fourier et Questions que s’y Rattachent, volume 1 of Publications Scientifiques de l’Institut Mittag-Leffler.  Almqvist & Wiksells, Uppsala, 1944.
  * [44] Clément Mouhot and Lorenzo Pareschi.  Fast algorithms for computing the Boltzmann collision operator.  Mathematics of computation, 75(256):1833–1852, 2006.
  * [45] Lei Wu, Craig White, Thomas J Scanlon, Jason M Reese, and Yonghao Zhang.  Deterministic numerical solutions of the Boltzmann equation using the fast spectral method.  Journal of Computational Physics, 250:27–52, 2013.
  * [46] Irene M Gamba, Jeffrey R Haack, Cory D Hauck, and Jingwei Hu.  A fast spectral method for the Boltzmann collision operator with general collision kernels.  SIAM Journal on Scientific Computing, 39(4):B658–B674, 2017.
  * [47] Tianbai Xiao, Chang Liu, Kun Xu, and Qingdong Cai.  A velocity-space adaptive unified gas kinetic scheme for continuum and rarefied flows.  Journal of Computational Physics, 415:109535, 2020.
  * [48] Graeme Austin Bird.  Molecular gas dynamics and the direct simulation of gas flows.  Clarendon, 1994.
  * [49] Charles Francis Curtiss and Joseph O Hirschfelder.  Integration of stiff equations.  Proceedings of the National Academy of Sciences of the United States of America, 38(3):235, 1952.
  * [50] Mark H Carpenter, Eric J Nielsen, and Sally A Viken.  Efficiency of High Order Temporal Schemes.  AIAA, 2003.
  * [51] Antony Jameson.  Evaluation of fully implicit Runge Kutta schemes for unsteady flow calculations.  Journal of Scientific Computing, 73(2):819–852, 2017.
  * [52] Hester Bijl, Mark H Carpenter, Veer N Vatsa, and Christopher A Kennedy.  Implicit time integration schemes for the unsteady compressible Navier–Stokes equations: laminar flow.  Journal of Computational Physics, 179(1):313–329, 2002.
  * [53] Lai Wang and Meilin Yu.  Comparison of ROW, ESDIRK, and BDF2 for Unsteady Flows with the High-Order Flux Reconstruction Formulation.  Journal of Scientific Computing, 83:1–27, 2020.
  * [54] Christopher A Kennedy and Mark H Carpenter.  Additive Runge–Kutta schemes for convection–diffusion–reaction equations.  Applied numerical mathematics, 44(1-2):139–181, 2003.
  * [55] Bernardo Cockburn and Chi-Wang Shu.  TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. II. General framework.  Mathematics of computation, 52(186):411–435, 1989.
  * [56] Lilia Krivodonova, Jianguo Xin, J-F Remacle, Nicolas Chevaugeon, and Joseph E Flaherty.  Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws.  Applied Numerical Mathematics, 48(3-4):323–338, 2004.
  * [57] Lilia Krivodonova.  Limiters for high-order discontinuous Galerkin methods.  Journal of Computational Physics, 226(1):879–896, 2007.
  * [58] Antony Jameson, Wolfgang Schmidt, and Eli Turkel.  Numerical solution of the Euler equations by finite volume methods using Runge Kutta time stepping schemes.  In 14th fluid and plasma dynamics conference, page 1259, 1981.
  * [59] Per-Olof Persson and Jaime Peraire.  Sub-cell shock capturing for discontinuous Galerkin methods.  In 44th AIAA Aerospace Sciences Meeting and Exhibit, page 112, 2006\.
  * [60] Sachin Premasuthan, Chunlei Liang, and Antony Jameson.  Computation of flows with shocks using the spectral difference method with artificial viscosity, I: basic formulation and application.  Computers & Fluids, 98:111–121, 2014.
  * [61] Andrew W Cook and William H Cabot.  A high-wavenumber viscosity for high-resolution numerical methods.  Journal of Computational Physics, 195(2):594–601, 2004.
  * [62] David Flad, Andrea Beck, and Claus-Dieter Munz.  Simulation of underresolved turbulent flows by adaptive filtering using the high order discontinuous Galerkin spectral element method.  Journal of Computational Physics, 313:1–12, 2016.
  * [63] Martin Frank, Cory Hauck, and Kerstin Kuepper.  Convergence of filtered spherical harmonic equations for radiation transport.  Communications in Mathematical Sciences, 14(5):1443–1465, 2016\.
  * [64] Abhishek Sheshadri.  An analysis of stability of the flux reconstruction formulation with applications to shock capturing.  PhD thesis, Ph. D. thesis, Stanford University, 2016.
  * [65] John P Boyd.  Chebyshev and Fourier spectral methods.  Courier Corporation, 2001.
  * [66] Mathea J Vuik and Jennifer K Ryan.  Multiwavelet troubled-cell indicator for discontinuity detection of discontinuous galerkin schemes.  Journal of Computational Physics, 270:138–160, 2014.
  * [67] Cheng-Shu Wang.  On the theory of the thickness of weak shock waves.  University of Michigan, 1956.
  * [68] Laure Saint-Raymond.  Hydrodynamic limits of the Boltzmann equation.  Springer Science & Business Media, 2009.
  * [69] Przemyslaw Bogacki and Lawrence F Shampine.  A 3(2) pair of Runge-Kutta formulas.  Applied Mathematics Letters, 2(4):321–325, 1989.
  * [70] Jing Fan and Ching Shen.  Statistical simulation of low-speed rarefied gas flows.  Journal of Computational Physics, 167(2):393–412, 2001.
  * [71] Y Sone, S Takata, and T Ohwada.  Numerical analysis of the plane couette flow of a rarefied gas on the basis of the linearized boltzmann equation for hard-sphere molecules.  European Journal of Mechanics B Fluids, 9(3):273–288, 1990.
  * [72] Benzi John, Xiao-Jun Gu, and David R Emerson.  Effects of incomplete surface accommodation on non-equilibrium heat transfer in cavity flow: A parallel dsmc study.  Computers & fluids, 45(1):197–201, 2011.
  * [73] Tianbai Xiao, Kun Xu, Qingdong Cai, and Tiezheng Qian.  An investigation of non-equilibrium heat transport in a gas system under external force field.  International Journal of Heat and Mass Transfer, 126:362–379, 2018\.
  * [74] Tianbai Xiao and Martin Frank.  A stochastic kinetic scheme for multi-scale plasma transport with uncertainty quantification.  Journal of Computational Physics, 432:110139, 2021.
  * [75] Tianbai Xiao and Martin Frank.  A stochastic kinetic scheme for multi-scale flow transport with uncertainty quantification.  arXiv preprint arXiv:2002.00277, 2020.

Figure 2: Travelling wave solutions with $N_{x}=20$ and polynomial degree $3$
at different reference Knudsen numbers.

(a) $\mathrm{Ma}=2$ (b) $\mathrm{Ma}=3$

Figure 3: Profiles of density, $U$-velocity and temperature across normal
shock wave at different upstream Mach numbers.

(a) $\mathrm{Kn}=0.0001$ (b) $\mathrm{Kn}=0.01$ (c) $\mathrm{Kn}=1$

Figure 4: Profiles of density, $U$-velocity and temperature at $t=0.15$ in the
Sod shock tube at different reference Knudsen numbers.

Figure 5: Profiles of $V$-velocity and surface shear stress in the Couette
flow at different reference Knudsen numbers.

Figure 6: Contours of $U$-velocity with streamlines and temperature with heat
flux vectors inside the cavity.

Figure 7: Velocity profiles along vertical and horizontal central lines inside
the cavity.