# Discrete-velocity-direction models of BGK-type with minimum entropy: II.
Weighted models

Yihong Chen  Department of Mathematical Sciences, Tsinghua University
Beijing, 100084, China<EMAIL_ADDRESS>,  Qian Huang
Department of Energy and Power Engineering, Tsinghua University
Beijing, 100084, China<EMAIL_ADDRESS>and  Wen-An Yong  Department
of Mathematical Sciences, Tsinghua University
Beijing, 100084, China<EMAIL_ADDRESS>
###### Abstract.

In this series of works, we develop a discrete-velocity-direction model (DVDM)
with collisions of BGK-type for simulating gas flows, where the molecular
motion is confined to some prescribed directions but the speed is still a
continuous variable in each orientation. In this article, we introduce a
weighted function in each orientation when recovering the macroscopic
parameters. Moreover, the internal molecular degrees of freedom are
considered. With this weighted DVDM, we develop three submodels by
incorporating the discrete velocity method, the Gaussian-extended quadrature
method of moments and the Hermite spectral method in each direction. These
spatial-time submodels are novel multidimensional versions corresponding to
the three approaches. Numerical tests with a series of 1-D and 2-D flow
problems show the efficiency of the weighted DVDM.

###### Key words and phrases:

BGK equation; minimum entropy principle; discrete-velocity model; extended
quadrature method of moments; Hermite spectral method

* Corresponding author

##  1\. Introduction

Gas flows with both continuum and rarefied regimes have found wide
applications in areas including lunar-lander-induced dusty plumes [33],
dynamics of ultrafine particles [37] and vacuum technology [39]. The motion of
rarefied gases can be properly modelled by kinetic equations describing the
time evolution of problem-specific distribution functions [16, 40]. Moreover,
the spirit of kinetic theory has been used to describe many interacting
systems out of equilibrium, including the airborne aerosols [11, 21] and
active matter systems [2].

Numerical solution of the kinetic equations, including the Boltzmann equation,
is rather expensive due to the binary collision terms as well as the high
dimensionality of the distribution functions (of space, velocity and other
intrinsic properties) [7]. To overcome such difficulties, various
approximations of the kinetic equations have been proposed. Among them, the
BGK model is a formal simplification of the Boltzmann equation with the binary
collision term replaced by a relaxation process towards local equilibriums
(Maxwellian) [3, 42]. This model retains key properties of the Boltzmann
equation, including the conservation laws of mass, momentum and energy, the
$H$-theorem, and the correct hydrodynamic limit [16, 36, 7]. Moreover, the BGK
model becomes so influential that its variants, such as the stochastic
particle BGK method [8] and the BGK-type lattice Boltzmann method [19, 25],
have been developed.

Furthermore, there are lasting efforts to seek efficient numerical methods
solving the BGK equation [7]. At this point, it should be mentioned that the
nonlinearity and non-local dependency of the Maxwellian have posed significant
challenges to velocity discretization. The issue is resolved after Mieussens
presented an elegant discrete-velocity BGK (DV-BGK) model based on the minimum
entropy principle [31, 32]. However, the computational cost raised by high
dimensions remains the major obstacle for a direct discretization of the phase
space.

On the other hand, the moment methods may be a cure by prescribing the form of
the distribution functions. This approach leads to moment closure systems with
their hydrodynamic counterpart being the Euler equations. For instance, the
well-known Grad’s 13-moment theory [13] rests on an expansion of the
Maxwellian, but it then becomes less reliable far from equilibrium. In
contrast, the quadrature-based method of moments (QBMM) [5, 10] is free from
such limitations at its very heart of model assumptions. Particularly, the
Gaussian extended quadrature method of moments (Gaussian-EQMOM) reconstructs
the velocity distribution as a summation over several Gaussian functions, and
the resultant nonlinear moment system for BGK satisfies the desired structural
stability conditions [22]. However, despite the well-posedness of Gaussian-
EQMOM in 1-D velocity space, it is difficult to deduce a multidimensional
version of QBMM.

Aiming at a higher-dimensional version of Gaussian-EQMOM, we proposed a
discrete-velocity-direction model for the BGK equation by forcing all
particles to move in a set of prescribed directions (denoted BGK-DVDM or
simply DVDM) in [23]. The discrete equilibrium is determined with the minimum
entropy principle. We have managed to demonstrate the existence, uniqueness
and numerical feasibility of thus defined equilibriums. This novel semi-
continuous model is flexible enough to incorporate various strategies of
removing the 1-D velocity dependence to generate spatial-time models.
Remarkably, the extension of Gaussian-EQMOM by the DVDM yields a hyperbolic
multidimensional moment closure system [23]. Moreover, combining the DVDM with
the discrete velocity method (DVM) in each direction gives a radial pattern of
discrete velocity nodes, rather than the conventional nodes selected on a
uniform cubic lattice [31, 32]. Let us mention that the DVDM seems a common
practice in solving radiative transfer equations (termed as the discrete
ordinates method therein) [9], while the DVDM for the Boltzmann equation was
developed in [44].

In this paper, we improve the BGK-DVDM in two aspects. First, the internal
molecular degrees of freedom is included so that more realistic fluid
properties can be realized [6]. More importantly, we introduce a weighted
function in each orientation when recovering the macroscopic parameters, as
opposed to the previous treatment. For the new weighted DVDM, the established
properties of the well-behaved discrete equilibrium in [23] still hold. Then
the 1-D DVM, Gaussian-EQMOM and a Hermite spectral method [20] are applied to
generate various DVDM submodels. Our numerical tests show the significance of
the weighted function.

In doing the numerical simulation, we specify proper gas-surface boundary
conditions [40] and second-order finite-volume spatial-time schemes [35] for
the new DVDM submodels. The performance of our DVDM submodels is examined with
a series of 1-D and 2-D problems covering a wide range of flow regimes,
including 2-D Riemann problems and lid-driven cavity flows. The numerical
results show that our BGK-DVDM with internal degrees of freedom is a promising
multiscale flow solver. It might be worthwhile to study the DVDM with the
moment methods further.

The rest of the paper is organized as follows. In Section 2, we introduce the
weighted BGK-DVDM with internal degrees of freedom. Section 3 develops three
DVDM submodels, including DVD-DVM in Section 3.1, DVD-EQMOM in Section 3.2 and
DVD-HSM in Section 3.3. Boundary conditions for the DVDM submodels are
presented in Section 4. Sections 5.1 and 5.2 contain the algorithms to compute
the discrete equilibria and the space-time discretization schemes,
respectively. Numertical results are reported in Section 6. Finally, some
conclusions are given in Section 7.

##  2\. Model development

###  2.1. BGK equation with internal degrees of freedom

We start with the BGK equation for the density function
$f=f(t,\bm{x},\bm{\xi},\bm{\zeta})$:

$\partial_{t}f+\bm{\xi}\cdot\nabla_{\bm{x}}f=\frac{1}{\tau}\left(\mathcal{E}[f]-f\right).$
(2.1)

Here, $(\bm{x},\bm{\xi})\in\mathbb{R}^{D}\times\mathbb{R}^{D}$ with $D=2$ or
$3$, $\bm{\zeta}\in\mathbb{R}^{L}$ represents the possible internal molecular
degrees of freedom, and $\tau$ is a characteristic collision time. The
internal molecular degrees of freedom determine the specific heat ratio of
ideal gases [15]. For example, the diatomic gases have two degrees of freedom:
rotation and vibration. For planar flows ($D=2$), $\bm{\xi}$ is just a part of
the $3$-dimensional molecule velocity and the rest part is included in
$\bm{\zeta}$.

In the right-hand side (RHS) of Eq. (2.1), the local equilibrium state
$\mathcal{E}[f]=\mathcal{E}[f](t,\bm{x},\bm{\xi},\bm{\zeta})$ is modeled as

$\mathcal{E}[f]=\frac{\rho}{(2\pi\theta)^{(D+L)/2}}\exp\left(-\frac{|\bm{\xi}-\bm{U}|^{2}+|\bm{\zeta}|^{2}}{2\theta}\right).$
(2.2)

Here $|\bm{U}|$ denotes the Euclidean length of the vector $\bm{U}$. The
classical fluid quantities including density $\rho$, velocity $\bm{U}$, energy
$E$, temperature $\theta$ and pressure $p$ are defined by

$\rho=\langle
f\rangle,\quad\bm{U}=\frac{\langle\bm{\xi}f\rangle}{\rho}\in\mathbb{R}^{D},\quad
E=\frac{|\bm{U}|^{2}+(D+L)\theta}{2}=\frac{1}{\rho}\left\langle\frac{|\bm{\xi}|^{2}+|\bm{\zeta}|^{2}}{2}f\right\rangle,\quad
p=\rho\theta,$ (2.3)

where the bracket $\langle\cdot\rangle$ is defined as the integral $\langle
g(\bm{\xi},\bm{\zeta})\rangle=\int_{\mathbb{R}^{L}}\int_{\mathbb{R}^{D}}g(\bm{\xi},\bm{\zeta})d\bm{\xi}d\bm{\zeta}$
for any reasonable $g(\bm{\xi},\bm{\zeta})$. The equilibrium can be rewritten
in a concise form

$\mathcal{E}[f]=\exp\left(\bm{\alpha}_{eq}\cdot\bm{m}(\bm{\xi},\bm{\zeta})\right)$

with

$\bm{\alpha}_{eq}=\left(\ln\frac{\rho}{(2\pi\theta)^{(D+L)/2}}-\frac{|\bm{U}|^{2}}{2\theta},\
\frac{\bm{U}}{\theta},\
-\frac{1}{\theta}\right)^{T}\quad\text{and}\quad\bm{m}(\bm{\xi},\bm{\zeta})=\left(1,\
\bm{\xi},\ \frac{|\bm{\xi}|^{2}+|\bm{\zeta}|^{2}}{2}\right)^{T}$ (2.4)

both being $(D+2)$-dimensional real vectors. This form enlightens us on the
model development in later sections.

The equilibrium distribution $\mathcal{E}[f]$ satisfies two important
properties. First, $\mathcal{E}[f]$ reproduces the local macroscopic
quantities in the same manner as $f$:

$\left\langle\bm{m}(\bm{\xi},\bm{\zeta})\mathcal{E}[f]\right\rangle=\bm{\rho}:=\left(\rho,\
\rho\bm{U},\ \rho E\right)^{T}\in\mathbb{R}^{D+2},$ (2.5)

and thus the BGK equation respects the conservation laws of mass, momentum,
and energy. Then, given any $\bm{\rho}\in\mathbb{R}^{D+2}$ with positive
components $\rho$ and $E$, $\mathcal{E}[f]$ is the unique non-negative
solution that minimizes the following kinetic entropy

$H[f]=\langle f\ln f-f\rangle$ (2.6)

subject to the constraint
$\left\langle\bm{m}(\bm{\xi},\bm{\zeta})f\right\rangle=\bm{\rho}$.

This property has been adapted in both a discrete-velocity model [31] and our
previous discrete-velocity-direction model for the BGK equation without
internal degrees of freedom [23].

###  2.2. Discrete-velocity-direction models

A discrete-velocity-direction model (DVDM) based on the BGK equation has been
proposed in our previous work [23]. Our aim here is to enhance the model and
extend it to the case with internal molecular degrees of freedom.

The DVDM assumes that the molecule transport is limited to $N$ prescribed
directions denoted by $\left\\{\bm{l}_{m}\right\\}_{m=1}^{N}$ with each
$\bm{l}_{m}$ located on the unit sphere $\mathbb{S}^{D-1}$, but the velocity
magnitude $\xi\in\mathbb{R}$ in each direction remains continuous. The
directions are selected with the following two requirements.

  1. (A)

$(\bm{l}_{1},\dots,\bm{l}_{N})\in\mathbb{R}^{D\times N}$ is of rank $D$ and
therefore $N\geq D$.

  2. (B)

Each direction $\bm{l}_{m}$ and its opposite $-\bm{l}_{m}$ belong to
$S_{m}\subset\mathbb{S}^{D-1}$, where the $S_{m}$’s constitute a disjoint
partition of the unit sphere $\mathbb{S}^{D-1}=\bigcup_{m=1}^{N}S_{m}$ and
each $S_{m}$ has the same measure.

The equal measure means that the directions are ‘uniformly distributed’. For
$D=2$, such a partition on $\mathbb{S}^{1}$ can be realized by setting
$\bm{l}_{m}=(\cos\gamma_{m},\sin\gamma_{m})$ and
$\gamma_{m}=\frac{(m-1)\pi}{N}$ or $\frac{(2m-1)\pi}{2N}$, which will be
adopted for all numerical tests in this paper. For $D=3$, the algorithm in
[28] can help to yield such a partition on $\mathbb{S}^{2}$.

Once the directions are selected, the distribution
$f=f(t,\bm{x},\bm{\xi},\bm{\zeta})$ is replaced by $N$ distributions
$\left\\{f_{m}(t,\bm{x},\xi,\bm{\zeta})\right\\}_{m=1}^{N}$ with
$\xi\in\mathbb{R}$ and $\bm{\zeta}\in\mathbb{R}^{L}$. The transport velocity
for $f_{m}=f_{m}(t,\bm{x},\xi,\bm{\zeta})$ is $\xi\bm{l}_{m}$, and the
governing equation for $f_{m}$ becomes

$\partial_{t}f_{m}+\xi\bm{l}_{m}\cdot\nabla_{\bm{x}}f_{m}=\frac{1}{\tau}(\mathcal{E}_{m}-f_{m}),\quad
m=1,\dots,N,$ (2.7)

with the local equilibriums
$\mathcal{E}_{m}=\mathcal{E}_{m}(t,\bm{x},\xi,\bm{\zeta})$ yet to be
determined.

For Eq. (2.7), we use the weight function $|\xi|^{D-1}$ and define new fluid
quantities

$\displaystyle\rho$
$\displaystyle=s\sum_{m=1}^{N}\int_{\mathbb{R}}\int_{\mathbb{R}^{L}}f_{m}|\xi|^{D-1}d\bm{\zeta}d\xi,\quad\rho\bm{U}=s\sum_{m=1}^{N}\int_{\mathbb{R}}\int_{\mathbb{R}^{L}}\xi\bm{l}_{m}f_{m}|\xi|^{D-1}d\bm{\zeta}d\xi,$
(2.8) $\displaystyle\rho E$
$\displaystyle=s\sum_{m=1}^{N}\int_{\mathbb{R}}\int_{\mathbb{R}^{L}}\frac{\xi^{2}+|\bm{\zeta}|^{2}}{2}f_{m}|\xi|^{D-1}d\bm{\zeta}d\xi.$

Here $s$ is half of the measure of $S_{m}$.

######  Remark 2.1.

In contrast to our previous model in [23], the local equilibrium for the new
model Eq. (2.7) will be evaluated at the just defined fluid quantities
computed with the weight function $|\xi|^{D-1}$. This weight function is
inspired by changing variables from the Cartesian coordinate to polar or
spherical coordinates. Its introduction is independent of the internal degrees
of freedom. Our numerical tests show that this weight function is substantial
for correctly reconstructing macroscopic quantities.

As for the equilibrium states $\\{\mathcal{E}_{m}\\}_{m=1}^{N}$ on the RHS of
Eq. (2.7), we require that the following conservation property must be
satisfied:

$s\sum_{m=1}^{N}\int_{\mathbb{R}}\int_{\mathbb{R}^{L}}\left(1,\xi\bm{l}_{m},\frac{\xi^{2}+|\bm{\zeta}|^{2}}{2}\right)\mathcal{E}_{m}|\xi|^{D-1}d\bm{\zeta}d\xi=(\rho,\rho\bm{U},\rho
E).$ (2.9)

This can be viewed as a discrete-velocity-direction analogue of Eq. (2.5),
while $\rho$, $\bm{U}$, and $E$ are computed with the weighted integrals Eq.
(2.8) based on $f_{m}$. In this way, we can derive the classical Euler
equations by multiplying $1$, $\xi\bm{l}_{m}$ and
$\frac{\xi^{2}+|\bm{\zeta}|^{2}}{2}$ on both sides of Eq. (2.7) and taking the
weighted integrals; see details in [23].

Next we assume that the local equilibrium
$\mathcal{E}_{m}=\mathcal{E}_{m}(t,\bm{x},\xi,\bm{\zeta})$ has the variable-
separating form

$\mathcal{E}_{m}(t,\bm{x},\xi,\bm{\zeta})=\mathcal{E}_{tr,m}(t,\bm{x},\xi)\mathcal{E}_{in,m}(t,\bm{x},\bm{\zeta}),$

which is consistent with the Maxwellian Eq. (2.2) of the BGK equation. The
internal part $\mathcal{E}_{in,m}$ is taken to be the same as that in Eq.
(2.2):

$\mathcal{E}_{in,m}(t,\bm{x},\bm{\zeta})=\frac{1}{\sqrt{(2\pi\theta)^{L}}}\exp\left(-\frac{|\bm{\zeta}|^{2}}{2\theta}\right).$

Notice that the equilibrium temperature $\theta$ is

$\theta=\frac{2E-|\bm{U}|^{2}}{D+L}$ (2.10)

due to Eq. (2.3). Substituting such an $\mathcal{E}_{in,m}$ into Eq. (2.9), we
derive constraints for the transport part $\mathcal{E}_{tr,m}(t,\bm{x},\xi)$:

$s\sum_{m=1}^{N}\int_{\mathbb{R}}\bm{m}_{m}\mathcal{E}_{tr,m}|\xi|^{D-1}d\xi=\left(\rho,\rho\bm{U},\rho\left(E-\frac{L}{2}\theta\right)\right)^{T}=:\bm{\rho}_{tr},$
(2.11)

where
$\bm{m}_{m}(\xi)=\left(1,\xi\bm{l}_{m},\xi^{2}/2\right)^{T}\in\mathbb{R}^{D+2}$.
To determine the transport part, we refer to the minimum entropy property of
the Maxwellian in Eq. (2.6) and require that $\mathcal{E}_{tr,m}$ minimizes a
discrete analogue of the entropy

$H\left[\\{g_{m}\\}_{m=1}^{N}\right]:=s\sum_{m=1}^{N}\int_{\mathbb{R}}\int_{\mathbb{R}^{L}}(g_{m}\ln
g_{m}-g_{m})|\xi|^{D-1}d\bm{\zeta}d\xi$ (2.12)

among all possible 1-D distributions $\left\\{g_{m}(\xi)\geq
0\right\\}_{m=1}^{N}$ satisfying

$s\sum_{m=1}^{N}\int_{\mathbb{R}}\bm{m}_{m}g_{m}|\xi|^{D-1}d\xi=\bm{\rho}_{tr}.$

For the transport part, we have the following theorem which can be proved with
the same argument as that of Theorem 2.1 in our previous work [23].

######  Theorem 2.2.

Given $\bm{\rho}_{tr}\in\mathbb{R}^{D+2}$ satisfying
$0<|\bm{\rho}_{tr}|<\infty$, if there exists $\\{g_{m}(\xi)\geq
0\\}_{m=1}^{N}$ such that

$s\sum_{m=1}^{N}\int_{\mathbb{R}}\bm{m}_{m}g_{m}|\xi|^{D-1}d\xi=\bm{\rho}_{tr},$

then the discrete kinetic entropy Eq. (2.12) has a unique minimizer
$\\{\mathcal{E}_{tr,m}\\}_{m=1}^{N}$. Moreover, the minimizer has the
exponential form

$\mathcal{E}_{tr,m}=\exp(\bm{\alpha}\cdot\bm{m}_{m})$

and
$\bm{\alpha}=(\alpha_{0},\hat{\bm{\alpha}},\alpha_{D+1})\in\mathbb{R}^{D+1}\times\mathbb{R}^{-}$
is the unique minimizer of the following convex function

$J(\bm{\alpha}):=s\sum_{m=1}^{N}\int_{\mathbb{R}}\exp\left(\bm{\alpha}\cdot\bm{m}_{m}\right)|\xi|^{D-1}d\xi-\bm{\rho}_{tr}\cdot\bm{\alpha}.$
(2.13)

Thanks to this result, the computation of $\mathcal{E}_{tr,m}$ only requires
solving $\bm{\alpha}\in\mathbb{R}^{D+2}$ by minimizing $J(\bm{\alpha})$. This
is particularly beneficial when a large number of directions are used, i.e.,
$N\gg D+2$. We will present the algorithm for this optimization problem in
later section, which is shown to be highly efficient.

In some cases, it is convenient to rewrite $\mathcal{E}_{tr,m}$ in the form of
standard Gaussian distribution

$\mathcal{E}_{tr,m}=\exp(\bm{\alpha}\cdot\bm{m}_{m})=\frac{\rho_{m}}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{(\xi-
u_{m})^{2}}{2\sigma^{2}}\right),$ (2.14)

and the parameters $\rho_{m},\ u_{m},$ and $\sigma^{2}$ are related to
$\bm{\alpha}=(\alpha_{0},\hat{\bm{\alpha}},\alpha_{D+1})$ as follows:

$\sigma^{2}=-\frac{1}{\alpha_{D+1}},\quad
u_{m}=\left(\hat{\bm{\alpha}}\cdot\bm{l}_{m}\right)\sigma^{2},\quad\rho_{m}=\sqrt{2\pi\sigma^{2}}\exp\left(\alpha_{0}+\frac{u_{m}^{2}}{2\sigma^{2}}\right).$
(2.15)

######  Remark 2.3.

Due to the weighted integral in Eq. (2.8), $\bm{\rho}_{tr}$ cannot be
expressed, in general, by $\rho_{m}$, $u_{m}$ and $\sigma^{2}$ with the simple
algebraic relations

$\rho\neq s\sum_{m=1}^{N}\rho_{m},\quad\rho\bm{U}\neq
s\sum_{m=1}^{N}\rho_{m}u_{m}\bm{l}_{m},\quad\rho(|\bm{U}|^{2}+D\theta)\neq
s\sum_{m=1}^{N}\rho_{m}(u_{m}^{2}+\sigma^{2}).$

##  3\. Spatial-time models

Our BGK-DVDM model Eq. (2.7) contains continuous variables $\xi\in\mathbb{R}$
and $\bm{\zeta}\in\mathbb{R}^{L}$. In this section, we treat these variables
to derive spatial-time models with only $t$ and $\bm{x}$ as continuous
variables. As to $\bm{\zeta}$, we define

$g_{m}(t,\bm{x},\xi)=\int_{\mathbb{R}^{L}}f_{m}(t,\bm{x},\xi,\bm{\zeta})d\bm{\zeta},\quad
h_{m}(t,\bm{x},\xi)=\int_{\mathbb{R}^{L}}|\bm{\zeta}|^{2}f_{m}(t,\bm{x},\xi,\bm{\zeta})d\bm{\zeta}$

and derive from Eq. (2.7):

$\left\\{\begin{aligned}
&\partial_{t}g_{m}+\xi\bm{l}_{m}\cdot\nabla_{\bm{x}}g_{m}=\frac{1}{\tau}\left(\exp(\bm{\alpha}\cdot\bm{m}_{m})-g_{m}\right),\\\
&\partial_{t}h_{m}+\xi\bm{l}_{m}\cdot\nabla_{\bm{x}}h_{m}=\frac{1}{\tau}\left(L\theta\exp(\bm{\alpha}\cdot\bm{m}_{m})-h_{m}\right)\end{aligned}\right.$
(3.1)

for $m=1,...,N$.

Next we treat $\xi$ with the following three methods: the discrete-velocity
method (DVM) [31], the extended quadrature method of moment (EQMOM) [5], and
the Hermite spectral method (HSM) [20].

###  3.1. Discrete-velocity model

To derive this kind of model, we choose a positive integer $M$, a positive
real number $\Delta\xi$, and a real number $\xi_{0}$, which can vary for
different directions. Set $\xi_{mk}=k\Delta\xi+\xi_{0}$ for $k=1,...,M$. Based
on Eq. (3.1), the discrete-velocity model is

$\left\\{\begin{aligned}
&\partial_{t}g_{mk}+\xi_{mk}\bm{l}_{m}\cdot\nabla_{\bm{x}}g_{mk}=\frac{1}{\tau}(g^{eq}_{mk}-g_{mk}),\\\
&\partial_{t}h_{mk}+\xi_{mk}\bm{l}_{m}\cdot\nabla_{\bm{x}}h_{mk}=\frac{1}{\tau}(h^{eq}_{mk}-h_{mk}),\end{aligned}\right.$
(3.2)

for $k=1,...,M$ and $m=1,...,N$, where $g_{mk}^{eq}$ and $h_{mk}^{eq}$ need to
be determined.

To this end, we first compute

$\displaystyle\rho$
$\displaystyle=s\sum_{m=1}^{N}\sum_{k=1}^{M}g_{mk}|\xi_{mk}|^{D-1}\Delta\xi,\quad\rho\bm{U}=s\sum_{m=1}^{N}\sum_{k=1}^{M}\xi_{mk}\bm{l}_{m}g_{mk}|\xi_{mk}|^{D-1}\Delta\xi,\quad$
(3.3) $\displaystyle\rho E$
$\displaystyle=s\sum_{m=1}^{N}\sum_{k=1}^{M}\frac{\xi_{mk}^{2}g_{mk}+h_{mk}}{2}|\xi_{mk}|^{D-1}\Delta\xi,$

corresponding to the last model. With these fluid quantities,
$\mathcal{E}_{tr,m}$ in Eq. (2.14) can be derived by finding the minimizer of
the convex function $J(\bm{\alpha})$ in Eq. (2.13) (see detailed algorithms in
Section 5.1).

Having $\mathcal{E}_{tr,m}$, we determine the discretized equilibriums
$g_{mk}^{eq}$ as the minimizer of the discrete entropy

$H_{m}\left[\\{u_{mk}\\}_{k}\right]=\sum_{k=1}^{M}(u_{mk}\ln
u_{mk}-u_{mk})|\xi_{mk}|^{D-1}\Delta\xi$

among all $\\{u_{mk}\geq 0\\}_{k=1}^{M}$ satisfying the conservation
constraint in the $m$-th direction:

$\sum_{k=1}^{M}\bm{m}_{mk}u_{mk}|\xi_{mk}|^{D-1}\Delta\xi=\int_{\mathbb{R}}\bm{m}_{m}|\xi|^{D-1}\mathcal{E}_{tr,m}d\xi\in\mathbb{R}^{D+2}=:\bm{\rho}_{tr,m},$
(3.4)

where we have
$\bm{m}_{mk}=\left(1,\xi_{mk}\bm{l}_{m},\xi_{mk}^{2}/2\right)^{T}$. With the
argument in [31], we can prove that this discretized equilibrium has the form

$g_{mk}^{eq}=\exp(\bm{\alpha}_{m}\cdot\bm{m}_{mk}),$ (3.5)

where $\bm{\alpha}_{m}\in\mathbb{R}^{D+2}$ is the unique minimizer of the
convex function

$J_{m}(\bm{\alpha})=\sum_{k=1}^{M}\exp(\bm{\alpha}\cdot\bm{m}_{mk})|\xi_{mk}|^{D-1}\Delta\xi-\bm{\alpha}\cdot\bm{\rho}_{tr,m}.$

Then we set $h_{mk}^{eq}=L\theta g_{mk}^{eq}$ with
$\theta=\frac{2E-|\bm{U}|^{2}}{D+L}$ as defined in Eq. (2.10).

With $g^{eq}_{mk}$ and $h^{eq}_{mk}$ determined above, our discrete-velocity
spatial-time model reads as

$\left\\{\begin{aligned}
&\partial_{t}g_{mk}+\xi_{mk}\bm{l}_{m}\cdot\nabla_{\bm{x}}g_{mk}=\frac{1}{\tau}(\exp(\bm{\alpha}_{m}\cdot\bm{m}_{mk})-g_{mk}),\\\
&\partial_{t}h_{mk}+\xi_{mk}\bm{l}_{m}\cdot\nabla_{\bm{x}}h_{mk}=\frac{1}{\tau}(L\theta\exp(\bm{\alpha}_{m}\cdot\bm{m}_{mk})-h_{mk})\end{aligned}\right.$

for $m=1,...,N$, and $k=1,...,M$. We close this subsection with two remarks on
this model.

######  Remark 3.1 (Computation of $\bm{\rho}_{tr,m}$).

Since $\mathcal{E}_{tr,m}$ is of the Gaussian form Eq. (2.14), it is not
difficult to verify that
$\bm{\rho}_{tr,m}=\rho_{tr,m}(1,u_{tr,m}\bm{l}_{m},E_{tr,m})^{T}$ in Eq. (3.4)
can be computed with the following formulae

$\displaystyle\rho_{tr,m}$
$\displaystyle=\frac{\rho_{m}u_{m}}{2}P\left(\frac{u_{m}}{\sqrt{2\sigma^{2}}}\right)+\rho_{m}\sqrt{\frac{2\sigma^{2}}{\pi}}\exp\left(-\frac{u_{m}^{2}}{2\sigma^{2}}\right),$
$\displaystyle\rho_{tr,m}u_{tr,m}$
$\displaystyle=\frac{\rho_{m}(u_{m}^{2}+\sigma^{2})}{2}P\left(\frac{u_{m}}{\sqrt{2\sigma^{2}}}\right)+\rho_{m}u_{m}\sqrt{\frac{2\sigma^{2}}{\pi}}\exp\left(-\frac{u_{m}^{2}}{2\sigma^{2}}\right),$
$\displaystyle 2\rho_{tr,m}E_{tr,m}$
$\displaystyle=\frac{\rho_{m}u_{m}(u_{m}^{2}+3\sigma^{2})}{2}P\left(\frac{u_{m}}{\sqrt{2\sigma^{2}}}\right)+\rho_{m}(u_{m}^{2}+2\sigma^{2})\sqrt{\frac{2\sigma^{2}}{\pi}}\exp\left(-\frac{u_{m}^{2}}{2\sigma^{2}}\right),$

where $P(x):=\mathrm{erfc}(-x)-\mathrm{erfc}(x)$ and
$\mathrm{erfc}(x)=\frac{2}{\sqrt{\pi}}\int_{x}^{\infty}e^{-\eta^{2}}d\eta$ is
the complementary error function. $\rho_{m}$, $u_{m}$, and $\sigma^{2}$ are
defined in Eq. (2.15).

######  Remark 3.2.

Besides the above procedure in deriving a discrete-velocity model, there is
another way to close Eq. (3.2). Indeed, the equilibrium
$\left\\{g_{mk}^{eq}\right\\}_{m,k}$ in Eq. (3.2) can be taken as the
minimizer of the ‘total’ discrete entropy

$H\left[\\{u_{mk}\\}_{m,k}\right]=s\sum_{m=1}^{N}\sum_{k=1}^{M}(u_{mk}\ln
u_{mk}-u_{mk})|\xi_{mk}|^{D-1}\Delta\xi$

among all $\\{u_{mk}\geq 0\\}_{m,k}$ satisfying the conservation constraint

$s\sum_{m=1}^{N}\sum_{k=1}^{M}\left(1,\xi_{mk}\bm{l}_{m},\frac{|\xi_{mk}^{2}|}{2}\right)u_{mk}|\xi_{mk}|^{D-1}\Delta\xi=\left(\rho,\rho\bm{U},\rho\left(E-\frac{L}{2}\theta\right)\right).$

In this way, $g_{mk}^{eq}=\exp(\bm{\alpha}\cdot\bm{m}_{mk})$ and
$\bm{\alpha}\in\mathbb{R}^{D+1}\times\mathbb{R}^{-}$ minimizes the convex
function

$J(\bm{\alpha})=s\sum_{m=1}^{N}\sum_{k=1}^{M}\exp(\bm{\alpha}\cdot\bm{m}_{mk})|\xi_{mk}|^{D-1}\Delta\xi-\bm{\alpha}\cdot\bm{\rho}_{tr}.$

Here $\bm{\rho}_{tr}$ is defined as in Eq. (2.11). This treatment is a variant
of that for the DVM in [31], except that the discrete velocity nodes are
chosen radially with a weight function $|\xi_{mk}|^{D-1}$. By contrast, the
conventional DVM prefers discrete velocity nodes in a cubic lattice in
$\mathbb{R}^{D}$.

Compared with the DVM approach that solves one (larger-scale) optimization
problem, our DVD-DVM requires additional computation of $\mathcal{E}_{tr,m}$,
$\bm{\rho}_{tr,m}$ and $N$ optimization problems. But computing
$\mathcal{E}_{tr,m}$ and $\bm{\rho}_{tr,m}$ are numerically efficient, and
minimizing each $J_{m}(\bm{\alpha})$ has a smaller scale than
$J(\bm{\alpha})$. Therefore, the computational cost is acceptable.

###  3.2. Gaussian-EQMOM

In this subsection, we apply a method of moment to the BGK-DVDM Eq. (3.1). The
$k$-th velocity moment of $g_{m}(t,\bm{x},\xi)$ and $h_{m}(t,\bm{x},\xi)$ are
defined as

$M^{[g]}_{m,k}(t,\bm{x})=\int_{\mathbb{R}}\xi^{k}g_{m}(t,\bm{x},\xi)d\xi,\quad
M^{[h]}_{m,k}(t,\bm{x})=\int_{\mathbb{R}}\xi^{k}h_{m}(t,\bm{x},\xi)d\xi$

for $k=0,1,2,...$. To derive the evolution equations for $M^{[g]}_{m,k}$ and
$M^{[h]}_{m,k}$, we integrate the BGK-DVDM Eq. (3.1) to get

$\displaystyle\partial_{t}M^{[g]}_{m,k}+\bm{l}_{m}\cdot\nabla_{\bm{x}}M^{[g]}_{m,k+1}=\frac{1}{\tau}\left(\rho_{m}\Delta_{k}(u_{m},\sigma^{2})-M^{[g]}_{m,k}\right),$
(3.6)
$\displaystyle\partial_{t}M^{[h]}_{m,k}+\bm{l}_{m}\cdot\nabla_{\bm{x}}M^{[h]}_{m,k+1}=\frac{1}{\tau}\left(L\rho_{m}\theta\Delta_{k}(u_{m},\sigma^{2})-M^{[h]}_{m,k}\right)$

for $k=0,1,2,...$. Here $\Delta_{k}(u,\sigma^{2})$ denotes the $k$-th moment
of the normalized Gaussian function centered at $u$ with a variance
$\sigma^{2}$. Eq. (3.6) contains infinitely many equations.

To get a system with finite equations, we resort to the Gaussian-EQMOM method.
In this method, it is assumed that the 1-D distribution $g_{m}(\xi)$ (and
$h_{m}(\xi)$) is a sum of $M$ Gaussian functions [30]:

$\phi_{m}(\xi)=\sum_{\alpha=1}^{M}\frac{w^{[\phi]}_{m,\alpha}}{\sqrt{2\pi\vartheta^{[\phi]}_{m}}}\exp\left(-\frac{(\xi-v^{[\phi]}_{m,\alpha})^{2}}{2\vartheta^{[\phi]}_{m}}\right),\quad\text{for
}\phi=g\text{ or }h.$ (3.7)

The variance $\vartheta^{[\phi]}_{m}>0$ is independent on the index $\alpha$.
With this ansatz, the moments can be expressed as

$M^{[\phi]}_{m,k}=\sum_{\alpha=1}^{M}w^{[\phi]}_{m,\alpha}\Delta_{k}\left(v^{[\phi]}_{m,\alpha},\vartheta^{[\phi]}_{m}\right)\quad\text{for
}k=0,1,\dots.$ (3.8)

The ansatz above has $2M+1$ parameters
$\left(w^{[\phi]}_{m,\alpha},v^{[\phi]}_{m,\alpha},\vartheta^{[\phi]}_{m}\right)$
for $g_{m}$ or $h_{m}$.

To fix these parameters, we reserve the equations in Eq. (3.6) with
$k=0,...,2M$ and then solve the first $2M+1$ equations in Eq. (3.8) to express
the parameters in terms of the reserved lower moments $M^{[g]}_{m,k}$ and
$M^{[h]}_{m,k}$ with $k=0,...,2M$. An algorithm to solve this set of nonlinear
algebraic equations can be found in the literature [5, 30], which is uniquely
solvable in most practical situations [22]. In this way, the higher moments
$M^{[g]}_{m,2M+1}$ and $M^{[h]}_{m,2M+1}$ in the governing equation of
$M^{[g]}_{m,2M}$ and $M^{[h]}_{m,2M}$ can also be expressed in terms of the
lower moments

$M^{[\phi]}_{m,2M+1}=\sum_{\alpha=1}^{M}w^{[\phi]}_{m,\alpha}\Delta_{2M+1}\left(v^{[\phi]}_{m,\alpha},\vartheta^{[\phi]}_{m}\right).$

Consequently, the equations in Eq. (3.6) with $k=0,...,2M$ are closed.

With the ansatz Eq. (3.7), the macroscopic quantities are naturally computed
as

$\displaystyle\rho$
$\displaystyle=s\sum_{m,\alpha}\int_{\mathbb{R}}|\xi|^{D-1}\mathcal{N}\left(\xi;W^{[g]}_{m,\alpha}\right)d\xi,$
(3.9) $\displaystyle\rho\bm{U}$
$\displaystyle=s\sum_{m,\alpha}\bm{l}_{m}\int_{\mathbb{R}}\xi|\xi|^{D-1}\mathcal{N}\left(\xi;W^{[g]}_{m,\alpha}\right)d\xi,$
$\displaystyle\rho E$
$\displaystyle=\frac{s}{2}\sum_{m,\alpha}\int_{\mathbb{R}}|\xi|^{D-1}\left[\xi^{2}\mathcal{N}\left(\xi;W^{[g]}_{m,\alpha}\right)+\mathcal{N}\left(\xi;W^{[h]}_{m,\alpha}\right)\right]d\xi,$

where

$\mathcal{N}\left(\xi;W^{[\phi]}_{m,\alpha}\right)=\frac{w^{[\phi]}_{m,\alpha}}{\sqrt{2\pi\vartheta^{[\phi]}_{m}}}\exp\left(-\frac{(\xi-v^{[\phi]}_{m,\alpha})^{2}}{2\vartheta^{[\phi]}_{m}}\right)\quad\text{
and }\quad
W^{[\phi]}_{m,\alpha}=\left(w^{[\phi]}_{m,\alpha},v^{[\phi]}_{m,\alpha},\vartheta^{[\phi]}_{m}\right)\in\mathbb{R}^{3}.$

Notice that due to the weight function $|\xi|^{D-1}$ in Eq. (3.9), we
generally have

$\rho\neq s\sum_{m}M^{[g]}_{m,0},\quad\rho\bm{U}\neq
s\sum_{m}\bm{l}_{m}M^{[g]}_{m,1},\quad\rho E\neq
s\sum_{m}\frac{1}{2}\left(M^{[g]}_{m,2}+M^{[h]}_{m,0}\right).$

Eqs. (3.6,3.8-3.9) make up a spatial-time model by incorporating Gaussian-
EQMOM into the BGK-DVDM Eq. (3.1). This model, denoted as DVD-EQMOM, is a
convenient multidimensional version of quadrature-based method of moments,
which seems better understood than those in [5, 30]. Moreover, the moment
system is hyperbolic, indicating a well-posed extension of the EQMOM. The
proof is similar to that in our previous work [23] for the BGK equation
without internal degrees of freedom. It mainly relies on Ref. [22], where the
hyperbolicity of the 1-D EQMOM was thoroughly analyzed.

###  3.3. Hermite spectral method

In this subsection we treat the continuous variable $\xi$ with the Hermite
spectral method (HSM) proposed in [20]. In this method, it is assumed that the
distribution $\phi=\phi(t,\bm{x},\xi)$ is a truncation

$\phi(t,\bm{x},\xi)=\sum\limits_{k=0}^{M-1}\phi_{k}(t,\bm{x})\mathcal{H}_{k}^{[\bar{u},\bar{\theta}]}(\xi)$
(3.10)

of a series with the basis function

$\mathcal{H}_{n}^{[\bar{u},\bar{\theta}]}(\xi)=\bar{\theta}^{-n/2}H_{n}\left(\frac{\xi-\bar{u}}{\sqrt{\bar{\theta}}}\right)\frac{1}{\sqrt{2\pi\bar{\theta}}}e^{-\frac{(\xi-\bar{u})^{2}}{2\bar{\theta}}}.$

Here $M$ is a given integer,

$H_{n}(x)=(-1)^{n}e^{\frac{x^{2}}{2}}\left(\frac{d^{n}}{dx^{n}}e^{-\frac{x^{2}}{2}}\right)$

is the $n$th-order Hermite polynomial, and $\bar{u}$, $\bar{\theta}$ are two
constant parameters. In this paper, we always set $\bar{u}=0$ and determine
$\bar{\theta}$ by the initial flow condition.

Due to the orthogonality of the Hermite polynomials:

$\int_{\mathbb{R}}\mathcal{H}_{n}^{[\bar{u},\bar{\theta}]}(\xi)\frac{\bar{\theta}^{m/2}}{m!}H_{m}\left(\frac{\xi-\bar{u}}{\sqrt{\bar{\theta}}}\right)d\xi=\delta_{nm},$

the coefficient $\phi_{k}(t,\bm{x})$ in Eq. (3.10) can be uniquely determined
as

$\phi_{k}(t,\bm{x})=\int_{\mathbb{R}}\phi(t,\bm{x},\xi)\frac{\bar{\theta}^{k/2}}{k!}H_{k}\left(\frac{\xi-\bar{u}}{\sqrt{\bar{\theta}}}\right)d\xi.$

Thus the $M$-truncation of $\phi$ is fully determined.

To incorporate the HSM with the BGK-DVDM Eq. (3.1), we set

$\displaystyle\phi^{[g]}_{m,k}(t,\bm{x})=\int_{\mathbb{R}}g_{m}(t,\bm{x},\xi)|\xi|^{D-1}\frac{\bar{\theta}^{k/2}}{k!}H_{k}\left(\frac{\xi-\bar{u}}{\sqrt{\bar{\theta}}}\right)d\xi,$
$\displaystyle\phi^{[h]}_{m,k}(t,\bm{x})=\int_{\mathbb{R}}h_{m}(t,\bm{x},\xi)|\xi|^{D-1}\frac{\bar{\theta}^{k/2}}{k!}H_{k}\left(\frac{\xi-\bar{u}}{\sqrt{\bar{\theta}}}\right)d\xi.$

Then we multiply the both sides of Eq. (3.1) with
$|\xi|^{D-1}\frac{\bar{\theta}^{k/2}}{k!}H_{k}\left(\frac{\xi-\bar{u}}{\sqrt{\bar{\theta}}}\right)$
for $k=0,...,M-1$ and integrate over $\xi\in\mathbb{R}$ to obtain

$\partial_{t}\bm{\Phi}_{m}+\mathcal{A}\bm{l}_{m}\cdot\nabla_{\bm{x}}\bm{\Phi}_{m}=\frac{1}{\tau}(\bm{\Phi}_{m}^{eq}-\bm{\Phi}_{m}).$
(3.11)

Here $\bm{\Phi}_{m}=(\phi_{m,0},...,\phi_{m,M-1})^{T}\in\mathbb{R}^{M}$ with
$\phi_{m,k}=\phi_{m,k}^{[g]}$ or $\phi_{m,k}^{[h]}$, the constant matrix
$\mathcal{A}\in\mathbb{R}^{M\times M}$ is tridiagonal [20]:

$\mathcal{A}=\left(\begin{matrix}\bar{u}&1&&&\\\ \bar{\theta}&\bar{u}&2&&\\\
&\bar{\theta}&\ddots&\ddots&\\\ &&\ddots&\ddots&M-1\\\
&&&\bar{\theta}&\bar{u}\end{matrix}\right),$ (3.12)

and the equilibrium
$\bm{\Phi}_{m}^{eq}=(\phi^{eq}_{m,0},...,\phi^{eq}_{m,M-1})^{T}$ has
components

$\phi^{eq}_{m,k}=\int_{\mathbb{R}}\frac{\bar{\theta}^{k/2}}{k!}H_{k}\left(\frac{\xi-\bar{u}}{\sqrt{\bar{\theta}}}\right)\phi_{m}^{eq}d\xi\quad\text{for
}\ \phi_{m}^{eq}=\mathcal{E}_{tr,m}|\xi|^{D-1}\ \text{ or }\
L\theta\mathcal{E}_{tr,m}|\xi|^{D-1}.$ (3.13)

The corresponding macroscopic quantities are computed as

$\displaystyle\rho$
$\displaystyle=s\sum_{m=1}^{N}\phi^{[g]}_{m,0},\quad\rho\bm{U}=s\sum_{m=1}^{N}\bm{l}_{m}\left(\phi^{[g]}_{m,1}+\bar{u}\phi^{[g]}_{m,0}\right),$
$\displaystyle\rho E$
$\displaystyle=s\sum_{m=1}^{N}\frac{1}{2}\left(2\phi^{[g]}_{m,2}+2\bar{u}\phi^{[g]}_{m,1}+(\bar{\theta}+\bar{u}^{2})\phi^{[g]}_{m,0}+\phi^{[h]}_{m,0}\right).$

The equations in Eq. (3.11) constitute our third kind of models, denoted as
DVD-HSM.

We end this subsection with details on computing $\phi_{m,k}^{eq}$ in Eq.
(3.13) when $\bar{u}=0$. Clearly, we only need to consider
$\phi^{eq}_{m}=\mathcal{E}_{tr,m}|\xi|^{D-1}$ with $\mathcal{E}_{tr,m}$ given
in Eq. (2.14). Moreover, only the 2-D case is presented because for $D=3$ the
weight function has a simpler expression $|\xi|^{D-1}=\xi^{2}$ and therefore
the 3-D case is easier to handle. To simplify the notation, we set

$\displaystyle a_{k}$
$\displaystyle:=\phi_{m,k}^{eq}=\int_{\mathbb{R}}\frac{\bar{\theta}^{k/2}}{k!}H_{k}\left(\frac{\xi}{\sqrt{\bar{\theta}}}\right)\mathcal{E}_{tr,m}(\xi)|\xi|d\xi,$
$\displaystyle b_{k}$
$\displaystyle:=\int_{\mathbb{R}}\frac{\bar{\theta}^{k/2}}{k!}H_{k}\left(\frac{\xi}{\sqrt{\bar{\theta}}}\right)\mathcal{E}_{tr,m}(\xi)\mathrm{sgn}(\xi)d\xi,$

for $k=0,1,...$. It is not difficult to see from the recursive formula of
Hermite polynomials [20]

$H_{0}(x)=1,\quad H_{1}(x)=x,\quad H_{n+1}(x)=xH_{n}(x)-nH_{n-1}(x)$

and the relation $|\xi|=\xi\mathrm{sgn}(\xi)$ that

$a_{0}=b_{1}\quad\text{and}\quad
a_{k}=(k+1)b_{k+1}+\bar{\theta}b_{k-1}\quad\text{for}\quad k\geq 1.$

Thus, it suffices to compute $b_{k}$. Write

$b_{k}=\int_{0}^{+\infty}\frac{\bar{\theta}^{k/2}}{k!}H_{k}\left(\frac{\xi}{\sqrt{\bar{\theta}}}\right)\mathcal{E}_{tr,m}(\xi)d\xi-\int_{-\infty}^{0}\frac{\bar{\theta}^{k/2}}{k!}H_{k}\left(\frac{\xi}{\sqrt{\bar{\theta}}}\right)\mathcal{E}_{tr,m}(\xi)d\xi=:b_{k}^{+}-b_{k}^{-}.$

With the recursive formula, we can obtain

$\displaystyle b_{k+1}^{+}=$
$\displaystyle\int_{0}^{+\infty}\frac{\bar{\theta}^{(k+1)/2}}{(k+1)!}\left[\frac{\xi}{\sqrt{\bar{\theta}}}H_{k}\left(\frac{\xi}{\sqrt{\bar{\theta}}}\right)-kH_{k-1}\left(\frac{\xi}{\sqrt{\bar{\theta}}}\right)\right]\mathcal{E}_{tr,m}(\xi)d\xi$
$\displaystyle=$
$\displaystyle\int_{0}^{+\infty}\frac{\bar{\theta}^{k/2}}{(k+1)!}\xi
H_{k}\left(\frac{\xi}{\sqrt{\bar{\theta}}}\right)\mathcal{E}_{tr,m}(\xi)d\xi-\frac{\bar{\theta}}{k+1}b_{k-1}^{+}$
$\displaystyle=$
$\displaystyle\int_{0}^{+\infty}\frac{\bar{\theta}^{k/2}}{(k+1)!}(\xi-
u_{m})H_{k}\left(\frac{\xi}{\sqrt{\bar{\theta}}}\right)\mathcal{E}_{tr,m}(\xi)d\xi+\frac{u_{m}}{k+1}b_{k}^{+}-\frac{\bar{\theta}}{k+1}b_{k-1}^{+}.$

Note that $\frac{d}{d\xi}\mathcal{E}_{tr,m}(\xi)=-\frac{\xi-
u_{m}}{\sigma^{2}}\mathcal{E}_{tr,m}(\xi)$ and
$H_{n}^{\prime}(x)=nH_{n-1}3(x)$. Using the integration by parts gives

$\displaystyle b_{k+1}^{+}=$
$\displaystyle-\int_{0}^{+\infty}\frac{\sigma^{2}\bar{\theta}^{k/2}}{(k+1)!}H_{k}\left(\frac{\xi}{\sqrt{\bar{\theta}}}\right)d\mathcal{E}_{tr,m}(\xi)+\frac{u_{m}}{k+1}b_{k}^{+}-\frac{\bar{\theta}}{k+1}b_{k-1}^{+}$
$\displaystyle=$
$\displaystyle\frac{\sigma^{2}\bar{\theta}^{k/2}}{(k+1)!}H_{k}(0)\mathcal{E}_{tr,m}(0)+\int_{0}^{+\infty}\frac{\sigma^{2}\bar{\theta}^{(k-1)/2}}{(k+1)!}kH_{k-1}\left(\frac{\xi}{\sqrt{\bar{\theta}}}\right)\mathcal{E}_{tr,m}(\xi)d\xi+\frac{u_{m}b_{k}^{+}-\bar{\theta}b_{k-1}^{+}}{k+1}$
$\displaystyle=$
$\displaystyle\frac{\sigma^{2}\bar{\theta}^{k/2}}{(k+1)!}H_{k}(0)\mathcal{E}_{tr,m}(0)+\frac{u_{m}b_{k}^{+}+(\sigma^{2}-\bar{\theta})b_{k-1}^{+}}{k+1}.$

A similar computation can be done for $b_{k}^{-}$ and finally we get the
following recursive formula

$b_{k+1}=\frac{2\sigma^{2}\bar{\theta}^{k/2}}{(k+1)!}H_{k}(0)\mathcal{E}_{tr,m}(0)+\frac{u_{m}b_{k}+(\sigma^{2}-\bar{\theta})b_{k-1}}{k+1}.$

Additionally, a direct computation gives

$b_{0}=\frac{1}{2}\left(\mathrm{erfc}\left(\frac{u_{m}}{\sqrt{2\sigma^{2}}}\right)-\mathrm{erfc}\left(\frac{u_{m}}{\sqrt{2\sigma^{2}}}\right)\right),\quad
b_{1}=u_{m}b_{0}+2\sigma^{2}\mathcal{E}_{tr,m}(0).$

###  3.4. Brief summary of the models

Fig. 1 presents a brief summary and the hierarchy of the several models
proposed up to now. The original BGK equation with internal molecular degrees
of freedom is reviewed in Section 2.1. The DVDM assumes that the particles
move in $N$ fixed orientations, leading to the model Eq. (2.7) for $f_{m}$.
Then, Section 3 develops three spatial-time models by eliminating the
continuous variables $\xi$ and $\bm{\zeta}$, including DVD-DVM in Section 3.1,
DVD-EQMOM in Section 3.2, and DVD-HSM in Section 3.3. For these models, the
boundary conditions and numerical schemes need to be specified before
practical flow simulations.

Figure 1. Model hierarchy of the DVDM. The variable and the equilibrium state
are shown in each model (block).

##  4\. Boundary conditions

Let $\Omega\subset\mathbb{R}^{D}$ be the computational domain and denote by
$\bm{n}=\bm{n}(\bm{x})$ the outward unit normal vector of the boundary
$\partial\Omega$ at $\bm{x}$. Two types of boundary conditions are considered
in this paper. The first one is the Neumann condition
$\bm{n}\cdot\nabla_{\bm{x}}\phi=0$ with $\phi$ representing any unknown
variables in the DVDM submodels (see Fig. 1).

The second one is the solid wall conditions. For simplicity, let the boundary
velocity $\bm{U}_{w}$ at $\bm{x}_{w}\in\partial\Omega$ be perpendicular to
$\bm{n}=\bm{n}(\bm{x}_{w})$. For the original BGK equation, the boundary
distribution $f(t,\bm{x}_{w},\bm{\xi},\bm{\zeta})$ for reflecting particles,
i.e. $\bm{\xi}\cdot\bm{n}<0$, should be given by the distribution of outgoing
particles, i.e. $\bm{\xi}\cdot\bm{n}>0$. Two specific boundary conditions are
the diffuse-scattering law and the bounce-back rule (also termed specular-
reflection law) [34, 14]. The first one assumes that the distribution of
reflecting particles is a Maxwellian:

$f(t,\bm{x}_{w},\bm{\xi},\bm{\zeta})=\sqrt{\frac{2\pi}{\theta_{w}}}\frac{j(t,\bm{x}_{w})}{\sqrt{(2\pi\theta_{w})^{D+L}}}\exp\left(-\frac{|\bm{\xi}-\bm{U}_{w}|^{2}+|\bm{\zeta}|^{2}}{2\theta_{w}}\right),\quad\bm{\xi}\cdot\bm{n}<0,$
(4.1)

where $\theta_{w}$ is the boundary temperature at $\bm{x}_{w}$ and
$j(t,\bm{x}_{w})$ is the outward-flowing mass flux defined by

$j(t,\bm{x}_{w})=\int_{\mathbb{R}^{L}}\int_{\bm{\xi}\cdot\bm{n}>0}\bm{n}\cdot\bm{\xi}f(t,\bm{x}_{w},\bm{\xi},\bm{\zeta})d\bm{\xi}d\bm{\zeta}.$

This condition ensures no particle penetration through the boundary. The
bounce-back rule is widely used in the lattice Boltzmann method [26]. It reads
as

$f(t,\bm{x}_{w},\bm{\xi},\bm{\zeta})=f(t,\bm{x}_{w},-\bm{\xi},\bm{\zeta})+2\rho_{w}(t,\bm{x}_{w})\mathcal{E}(\bm{\xi},\bm{\zeta})\frac{\bm{\xi}\cdot\bm{U}_{w}}{\theta_{w}},\quad\bm{\xi}\cdot\bm{n}<0.$

Here

$\mathcal{E}(\bm{\xi},\bm{\zeta})=\frac{1}{\sqrt{(2\pi\theta_{w})^{D+L}}}\exp\left(-\frac{|\bm{\xi}|^{2}+|\bm{\zeta}|^{2}}{2\theta_{w}}\right)$

and

$\rho_{w}(t,\bm{x}_{w})=\frac{2\int_{\mathbb{R}^{L}}\int_{\bm{\xi}\cdot\bm{n}>0}f(t,\bm{x}_{w},\bm{\xi},\bm{\zeta})d\bm{\xi}d\bm{\zeta}}{1-\frac{2}{\theta_{w}}\int_{\mathbb{R}^{L}}\int_{\bm{\xi}\cdot\bm{n}<0}\bm{\xi}\cdot\bm{U}_{w}\mathcal{E}(\bm{\xi},\bm{\zeta})d\bm{\xi}d\bm{\zeta}}.$

This condition ensures that the macroscopic velocity $\bm{U}(t,\bm{x}_{w})$
equals $\bm{U}_{w}$. Since $\bm{U}_{w}$ is assumed to be perpendicular to
$\bm{n}$, $\rho_{w}$ is simplified as

$\rho_{w}(t,\bm{x}_{w})=2\int_{\mathbb{R}^{L}}\int_{\bm{\xi}\cdot\bm{n}>0}f(t,\bm{x}_{w},\bm{\xi},\bm{\zeta})d\bm{\xi}d\bm{\zeta}.$

We now illustrate how these kinetic boundary conditions are adapted to the new
DVDM submodels in Section 3. The main idea is to replace the integrals above
by proper discrete sums.

For the DVD-DVM in Subsection 3.1, the diffuse-scattering law is converted to

$\displaystyle g_{mk}(t,\bm{x}_{w})$
$\displaystyle=\sqrt{\frac{2\pi}{\theta_{w}}}j(t,\bm{x}_{w})\mathcal{E}_{tr,mk}[\bm{U}_{w},\theta_{w}],\quad$
$\displaystyle\xi_{k}\bm{l}_{m}\cdot\bm{n}<0,$ $\displaystyle
h_{mk}(t,\bm{x}_{w})$ $\displaystyle=L\theta_{w}g_{mk}(t,\bm{x}_{w}),\quad$
$\displaystyle\xi_{k}\bm{l}_{m}\cdot\bm{n}<0,$

with

$j(t,\bm{x}_{w})=s\sum_{m=1}^{N}\sum_{k=1}^{M}\xi_{k}\bm{l}_{m}\cdot\bm{n}g_{mk}(t,\bm{x}_{w})|\xi_{k}|^{D-1}\Delta\xi\bm{1}_{\xi_{k}\bm{l}_{m}\cdot\bm{n}>0},$

and $\mathcal{E}_{tr,mk}[\bm{U}_{w},\theta_{w}]$ the discrete equilibrium
defined in Eq. (3.5) with density $1$, velocity $\bm{U}_{w}$, and temperature
$\theta_{w}$. On the other hand, we assume that $\\{\xi_{k}\\}_{k=1}^{M}$
satisfies $\xi_{k}=-\xi_{M+1-k}$ for $k=1,...,M$ to apply the bounce-back rule
to the DVD-DVM. With this assumption, the discrete-velocity version of bounce-
back rule becomes

$\displaystyle g_{mk}(t,\bm{x}_{w})$
$\displaystyle=g_{m,M+1-k}(t,\bm{x}_{w})+2\rho_{w}(t,\bm{x}_{w})\mathcal{E}_{tr,mk}[\bm{0},\theta_{w}]\frac{\xi_{k}\bm{l}_{m}\cdot\bm{U}_{w}}{\theta_{w}},\quad$
$\displaystyle\xi_{k}\bm{l}_{m}\cdot\bm{n}<0,$ $\displaystyle
h_{mk}(t,\bm{x}_{w})$
$\displaystyle=h_{m,M+1-k}(t,\bm{x}_{w})+2\rho_{w}(t,\bm{x}_{w})L\theta_{w}\mathcal{E}_{tr,mk}[\bm{0},\theta_{w}]\frac{\xi_{k}\bm{l}_{m}\cdot\bm{U}_{w}}{\theta_{w}},\quad$
$\displaystyle\xi_{k}\bm{l}_{m}\cdot\bm{n}<0,$

where

$\rho_{w}(t,\bm{x}_{w})=\frac{2s\sum_{m=1}^{N}\sum_{k=1}^{M}g_{mk}(t,\bm{x}_{w})|\xi_{k}|^{D-1}\Delta\xi\bm{1}_{\xi_{k}\bm{l}_{m}\cdot\bm{n}>0}}{1-\frac{2}{\theta_{w}}s\sum_{m=1}^{N}\sum_{k=1}^{M}\xi_{k}\bm{l}\cdot\bm{U}_{w}\mathcal{E}_{tr,mk}[\bm{0},\theta_{w}]|\xi_{k}|^{D-1}\Delta\xi\bm{1}_{\xi_{k}\bm{l}_{m}\cdot\bm{n}<0}}.$

For the DVD-EQMOM in Subsection 3.2, only the diffuse-scattering law is used,
which reconstructs the velocity distributions of reflecting particles as

$\displaystyle g_{m}(t,\bm{x}_{w},\xi)$
$\displaystyle=\sqrt{\frac{2\pi}{\theta_{w}}}j(t,\bm{x}_{w})\mathcal{E}_{tr,m}[(1,\bm{U}_{w},\theta_{w})],\quad$
$\displaystyle\xi\bm{l}_{m}\cdot\bm{n}<0,$ $\displaystyle
h_{m}(t,\bm{x}_{w},\xi)$
$\displaystyle=L\theta_{w}g_{m}(t,\bm{x}_{w},\xi),\quad$
$\displaystyle\xi\bm{l}_{m}\cdot\bm{n}<0,$

with

$j(t,\bm{x}_{w})=s\sum_{m=1}^{N}\int_{\xi\bm{l}_{m}\cdot\bm{n}>0}\xi\bm{l}_{m}\cdot\bm{n}|\xi|^{D-1}\sum_{\alpha=1}^{M}\mathcal{N}\left(\xi;W^{[g]}_{m,\alpha}\right)d\xi,$

and $\mathcal{E}_{tr,m}[\bm{U}_{w},\theta_{w}]$ the discrete equilibrium
defined in Eq. (2.14) with density $1$, velocity $\bm{U}_{w}$, and temperature
$\theta_{w}$. Other notations follow the definitions in Section 3.2. The
moments on the boundary can then be evaluated as

$M_{m,k}^{[\phi]}(t,\bm{x}_{w})=\int_{\\{\xi\bm{l}_{m}\cdot\bm{n}<0\\}\bigcup\\{\xi\bm{l}_{m}\cdot\bm{n}>0\\}}\xi^{k}\phi_{m}(t,\bm{x}_{w},\xi)d\xi$

for $\phi_{m}=g_{m}$ or $h_{m}$. Note that the integrand takes different forms
in the two sets.

For the DVD-HSM in Subsection 3.3, further boundary conditions are left for
future work.

##  5\. Algorithms

###  5.1. Algorithm for the discrete equilibrium

Solving the discrete equilibrium $\mathcal{E}_{tr,m}$ defined in Eq. (2.11)
out of a known $\bm{\rho}_{tr}$ is necessary for all DVDM submodels in Section
3. Theorem 2.2 indicates that all we need is an
$\bm{\alpha}\in\mathbb{R}^{D+1}\times\mathbb{R}^{-}$ that minimizes the convex
function $J(\bm{\alpha})$ in Eq. (2.13). The gradient descent method was used
in our previous work [23], while we use the BFGS quasi-Newton method [4] in
this work.

Fig. 2 presents the performance of the BFGS quasi-Newton method and the
gradient descend (GD) method for $D=2$. Here we set $\rho=1$,
$\bm{U}=(0.4,0.8)^{T}$, and $\theta=0.8$. The discrete directions are chosen
as
$\left\\{\bm{l}_{m}=\left(\cos\frac{(m-1)\pi}{N},\sin\frac{(m-1)\pi}{N}\right)^{T}\right\\}_{m=1}^{N}$.
The BFGS method requires much less iteration steps to converge for $N\leq 5$.
Notably, when $N\geq 7$, the initial value $\bm{\alpha}_{eq}$ is so close to
the minimizer $\bm{\alpha}$ that only one step of iteration leads to
convergence. Therefore, the computation of discrete equilibrium in the DVDM is
numerically efficient.

Figure 2. Performance of the BFGS quasi-Newton method and the gradient descend
(GD) method for $D=2$, $\rho=1$, $\bm{U}=(0.4,0.8)^{T}$, and $\theta=0.8$.

###  5.2. Numerical schemes

In this subsection we present some numerical schemes to solve the DVDM
submodels proposed before. Recall that both the DVD-DVM and DVD-HSM can be
written in a unified form as

$\partial_{t}\bm{\Phi}_{m}+\mathcal{A}\bm{l}_{m}\cdot\nabla_{\bm{x}}\bm{\Phi}_{m}=\frac{1}{\tau}(\bm{\Phi}_{m}^{eq}-\bm{\Phi}_{m})=:\bm{\Omega}_{m},\quad
m=1,...,N.$ (5.1)

For the DVD-DVM, we have
$\bm{\Phi}_{m}=(\phi_{m1},...,\phi_{mM})^{T}\in\mathbb{R}^{M}$ with
$\phi_{mk}=g_{mk}$ or $h_{mk}$, and the matrix $\mathcal{A}={\rm
diag}\\{\xi_{1},...,\xi_{M}\\}$. For the DVD-HSM,
$\bm{\Phi}_{m}=(\phi_{m,0},...,\phi_{m,M-1})^{T}\in\mathbb{R}^{M}$ and
$\mathcal{A}$ are defined in Eqs. (3.11 & 3.12).

For a time discretization of Eq. (5.1), the implicit-explicit Runge-Kutta
(IMEX-RK) schemes [35] can be applied. Here we only use a second-order scheme
denoted by SSP2. It is characterised by a double tableau [35]

$\begin{array}[]{c|cc}0&0&0\\\ 1&1&0\\\
\hline\cr&1/2&1/2\end{array}\quad\begin{array}[]{c|cc}\gamma&\gamma&0\\\
1-\gamma&1-2\gamma&\gamma\\\
\hline\cr&1/2&1/2\end{array},\quad\gamma=1-\frac{1}{\sqrt{2}}.$

Although the source term $\bm{\Omega}_{m}$ is implicitly discretized, its
relaxation structure renders a well-known way to solve the equations
explicitly (see e.g. [14, 15]). The convection term is treated with the third-
order energy stable WENO (ES-WENO) scheme [43]. For the DVD-DVM, the Godunov
flux [35] is adopted, while the HLL flux [17, 20] is used for the DVD-HSM.

On the other hand, for the DVD-DVM, Eq. (5.1) can also be discretized with
upwind schemes of first-order accuracy, which renders a easier way to treat
the boundary conditions. An implicit discretization for the collision term can
be treated similarly as in the IMEX-RK scheme.

Finally for the DVD-EQMOM, the $M^{[g]}_{m,k}$-equation in Eq. (3.6) can be
approximated by the 2-D upwind scheme

$\displaystyle M^{[g],n+1}_{m,k,ij}=$ $\displaystyle
M^{[g],n}_{m,k,ij}-\frac{\Delta t}{\Delta
x}\bm{l}_{m}\cdot\bm{e}_{1}\left(\mathcal{G}^{n}_{m,k+1,i+\frac{1}{2},j}-\mathcal{G}^{n}_{m,k+1,i-\frac{1}{2},j}\right)$
(5.2) $\displaystyle-$ $\displaystyle\frac{\Delta t}{\Delta
y}\bm{l}_{m}\cdot\bm{e}_{2}\left(\mathcal{G}^{n}_{m,k+1,i,j+\frac{1}{2}}-\mathcal{G}^{n}_{m,k+1,i,j-\frac{1}{2}}\right)+\frac{\Delta
t}{\tau}\left(M^{[g],n}_{\mathcal{E}m,k,ij}-M^{[g],n+1}_{m,k,ij}\right)$

with a partially implicit collision term. The $M^{[h]}_{m,k}$-equation is
treated similarly. Here the fluxes

$\mathcal{G}^{n}_{m,k+1,i+\frac{1}{2}}=\left\\{\begin{aligned}
\int_{0}^{\infty}\xi^{k+1}g_{m,ij}^{n}d\xi+\int_{-\infty}^{0}\xi^{k+1}g_{m,i+1,j}^{n}d\xi,\quad&\text{if
}\bm{l}_{m}\cdot\bm{e}_{1}>0,\\\
\int_{0}^{\infty}\xi^{k+1}g_{m,i+1,j}^{n}d\xi+\int_{-\infty}^{0}\xi^{k+1}g_{m,ij}^{n}d\xi,\quad&\text{if
}\bm{l}_{m}\cdot\bm{e}_{1}<0\end{aligned}\right.$ (5.3)

are the same as those in [5, 30]. The moments

$M^{[g],n}_{\mathcal{E}_{m,k,ij}}=\rho^{n}_{m,ij}\Delta_{k}\left(u^{n}_{m,ij},(\sigma^{2})^{n}_{ij}\right)$

correspond to the equilibrium state, where $\Delta_{k}(u,\sigma^{2})$ is
defined in Section 3.2. The equilibrium state parameters $\rho^{n}_{m,ij},\
u^{n}_{m,ij}$ and $(\sigma^{2})^{n}_{ij}$ are obtained by solving the local
equilibrium Eq. (2.14).

##  6\. Numerical results

In this section, we present the results of some numerical tests based on the
discretizations of the previous DVDM submodels. The tests only involve planar
flows ($D=2$).

###  6.1. 1-D Riemann Problems

We start with 1-D Riemann problems. Assume no internal degrees of freedom
($L=0$). The Riemann initial data of the fluid quantities read as [10]:

$\rho(0,x)=\left\\{\begin{aligned} &3.093,&x<0,\\\
&1,&x>0,\end{aligned}\right.\quad\bm{U}(0,x)=\bm{0},\quad\theta(0,x)=1.$

Both the continuum (infinitely fast collision limit $\tau=0$) and free-
molecular (no collision limit $\tau=\infty$) regimes are considered. The
theoretical solutions for both cases can be found in [29] and [15]. The 1-D
physical domain $[-0.5,0.5]$ is divided into $200$ uniform cells. The Neumann
boundary condition $\frac{\partial f}{\partial\bm{n}}=0$ is applied by
extending the values on the boundary cells constantly along the outward-facing
unit normal vector $\bm{n}$. We test all three DVDM submodels with this
problem.

The continuum regime is characterised with $\tau=10^{-4}$. In all DVDM
submodels, we set $N=8$ and the directions
$\left\\{\bm{l}_{m}=\left(\cos\frac{(2m-1)\pi}{16},\sin\frac{(2m-1)\pi}{16}\right)^{T}\right\\}_{m=1}^{8}$.
In the DVD-DVM, the discrete velocity nodes in each direction are selected as
$\xi_{k}=0.4k-5$ for $k=1,...,24$. In the DVD-HSM, we choose the order $M=12$
for the truncated series in Eq. (3.10). The SSP2 scheme in Section 5.2 is
applied to both the DVD-DVM and DVD-HSM. For the DVD-EQMOM, we set $M=2$.

Fig. 3 shows the spatial distributions of the macroscopic quantities
$(\rho,u,E,p)$ at $t=0.2$. Both the simulated results and theoretical
solutions are plotted. The shock wave that goes right, the rarefraction wave
that goes left, and the discontinuity between them are all well captured. It
is seen that what produced by both the DVD-DVM and DVD-HSM agree well with the
analytical solutions except some oscillations near the discontinuities, while
the two-node EQMOM is less accurate. However, the DVD-DVM yields the worst
result for the heat flux
$\bm{q}=\frac{1}{2}\left\langle(\bm{\xi}-\bm{U})(|\bm{\xi}-\bm{U}|^{2}+|\bm{\zeta}|^{2})f\right\rangle$,
which should be zero since it is easy to verify that
$\left\langle(\bm{\xi}-\bm{U})(|\bm{\xi}-\bm{U}|^{2}+|\bm{\zeta}|^{2})\mathcal{E}[f]\right\rangle=0$
for $\mathcal{E}[f]$ in Eq. (2.2) (the bracket $\langle\cdot\rangle$ is
defined in Eq. (2.3)). More directions and discrete nodes may be needed to
reduce such a discrepancy.

Figure 3. 1-D Riemann problem with $\tau=10^{-4}$: profiles of density $\rho$,
velocity $u$, energy $E$, pressure $p$ and heat flux $q$ at $t=0.2$. In all
models, we set $N=8$ and the directions
$\left\\{\bm{l}_{m}=\left(\cos\frac{(2m-1)\pi}{16},\sin\frac{(2m-1)\pi}{16}\right)^{T}\right\\}_{m=1}^{8}$.
In the DVD-DVM, the discrete velocities in each direction are $\xi_{k}=0.4k-5$
for $k=1,...,24$. We set $M=12$ for the DVD-HSM and $M=2$ for the DVD-EQMOM.

We emphasize that the weighted integral in Eq. (2.8), with the weight function
$|\xi|^{D-1}$, is a key feature different from our previous model in [23].
This weight function has been carefully treated in all DVDM submodels in
Section 3. As a direct comparison, Fig. 4 shows that without this weight
function, the predicted heat flux $q$ deviates significantly from zero, which
contradicts the Euler limit solution. Other properties have larger errors as
well. Therefore, only with this weight function $|\xi|^{D-1}$, the resultant
DVDM can produce satisfactory results.

Figure 4. 1-D Riemann problem with $\tau=10^{-4}$: A comparison between the
DVD-DVM predictions with and without the weight function $|\xi|^{D-1}$. All
other setups are the same as in the previous case.

As for the free-molecular flow regime, we take $\tau=10^{4}$ to create a near-
zero collision term. In all DVDM submodels, we set $N=18$ and the directions
$\left\\{\bm{l}_{m}=\left(\cos\frac{(2m-1)\pi}{36},\sin\frac{(2m-1)\pi}{36}\right)^{T}\right\\}_{m=1}^{18}$.
In the DVD-DVM, the discrete velocity nodes in each direction are selected as
$\xi_{k}=0.4k-7.4$ for $k=1,...,36$. In the DVD-HSM, we still choose the order
$M=12$. The SSP2 scheme is applied to both the DVD-DVM and DVD-HSM. For the
DVD-EQMOM, we set $M=2$. Fig. 5 presents the resulting profiles of macroscopic
quantities at $t=0.2$. Obviously there is no shock in this case, and the DVD-
DVM shows the highest accuracy. The relatively large error of the DVD-EQMOM is
partly due to the small number of nodes ($M=2$) used in the simulation.

Figure 5. 1-D Riemann problem with $\tau=10^{4}$: profiles of density $\rho$,
velocity $u$, energy $E$, pressure $p$ and heat flux $q$ at $t=0.2$. In all
models, we set $N=18$ and the directions
$\left\\{\bm{l}_{m}=\left(\cos\frac{(2m-1)\pi}{36},\sin\frac{(2m-1)\pi}{36}\right)^{T}\right\\}_{m=1}^{18}$.
In the DVD-DVM, the discrete velocities in each direction are
$\xi_{k}=0.4k-7.4$ for $k=1,...,36$. We set $M=12$ for the DVD-HSM and $M=2$
for the DVD-EQMOM.

###  6.2. Couette Flow

The flow is confined between two infinite parallel walls located at $x=\pm
0.5H$. The left and right walls move with constant velocities $\pm
v_{w}\bm{e}_{y}$ to drive the fluid between them to a steady state. In this
way, the flow reduces to a spatially 1-D problem in $x$. Assume $D=2$ and
$L=0$ (no internal degrees of freedom). Let $H=1$, $v_{w}=0.1$ and the wall
temperature $\theta_{w}=2$. The initial values of the fluid are
$(\rho_{0},\bm{U}_{0},\theta_{0})=(1,\bm{0},2)$. These settings ensure a small
Mach number.

In the Couette flow, different flow regimes are characterized by the parameter
$\kappa:=(\sqrt{\pi}/2)\mathrm{Kn}$, where the Knudsen number $\mathrm{Kn}$ is
defined as [14]

$\mathrm{Kn}=\frac{\tau}{H}\sqrt{\frac{\pi\theta_{0}}{2}}.$

Thus, the flow regime can be tuned by varying the values of $\tau$.

Both the DVD-DVM and DVD-EQMOM are used with the first-order upwind scheme
(see Section 5.2). The 1-D physical domain $[-0.5,0.5]$ is divided into $200$
uniform cells. The diffuse-scattering law is applied as the wall boundary
condition. The computation stops when the $L^{2}$-norm of the difference of
$\bm{U}$ between two consecutive time steps is smaller than $10^{-6}$, which
indicates that the flow is in a steady state. We set $N=15$ and the directions
$\bm{l}_{\mathcal{L}}=\left\\{\left(\cos\frac{(m-1)\pi}{15},\sin\frac{(m-1)\pi}{15}\right)^{T}\right\\}_{m=1}^{15}$
in all computations. For the DVD-DVM, the velocity nodes in each direction are
chosen as $\xi_{k}=0.5k-5.75$ for $k=1,...,22$. For the DVD-EQMOM, we let
$M=2$.

Fig. 6 shows the steady-state vertical velocity profiles on the positive
domain $x>0$ for different values of $\kappa$. The velocity is normalized by
the wall velocity $v_{w}$. The DSMC results in [1] are included for a
comparison. Apparently, higher values of $\kappa$ correspond to more rarefied
gases and less momentum transfer from the moving wall to the fluids. Both the
DVD-DVM and DVD-EQMOM reproduce the velocity profiles quite close to the
reference data for all three values of $\kappa$. Fig. 6 further presents the
shear stress $\tau_{xy}$ defined by

$\tau_{xy}=\int_{\mathbb{R}^{2}}(\xi_{x}-u)(\xi_{y}-v)f(\bm{\xi})d\bm{\xi}=\int_{\mathbb{R}^{2}}\xi_{x}\xi_{y}f(\bm{\xi})d\bm{\xi}-\rho
uv$

for Kn ranging from 0.01 to 100. Here we denote
$\bm{\xi}=(\xi_{x},\xi_{y})^{T}$ and $\bm{U}=(u,v)^{T}\in\mathbb{R}^{2}$. The
shear stress is normalized by the free-molecular stress $\tau_{\infty}=-\rho
u_{w}\sqrt{2\theta/\pi}$. Our DVDM results are generally in good agreement
with the DSMC results [1]. It is seen that the two-node DVD-DVDM has more
significant errors at larger Kn (rarefied flow) conditions, as compared with
the DVD-DVM made up by more velocity nodes.

Figure 6. Couette flow: (a) Steady-state vertical velocity profiles, and (b)
Shear stress for different values of Kn. In all models, $N=15$ and the
directions are
$\\{\bm{l}_{m}=(\cos\frac{(m-1)\pi}{15},\sin\frac{(m-1)\pi}{15})\\}_{m=1}^{15}$.
In the DVD-DVM, the discrete velocities in each direction are
$\xi_{k}=0.5k-5.75$ for $k=1,...,22$. The DSMC data is from [1].

###  6.3. 2-D Riemann Problems

Two-dimensional Riemann problems have been studied in [27]. Here we consider
the following initial data

$(\rho,u,v,p)=\left\\{\begin{array}[]{lll}(\rho_{1},u_{1},v_{1},p_{1})=(0.5313,\
0,\ 0,\ 0.4),&x>0,&y>0,\\\ (\rho_{2},u_{2},v_{2},p_{2})=(1,\ 0.7276,\ 0,\
1),&x\leq 0,&y>0,\\\ (\rho_{3},u_{3},v_{3},p_{3})=(0.8,\ 0,\ 0,\ 1),&x\leq
0,&y\leq 0,\\\ (\rho_{4},u_{4},v_{4},p_{4})=(1,\ 0,\ 0.7276,\ 1),&x>0,&y\leq
0,\end{array}\right.$

which was also studied in [15]. In contrast to the previous subsections, the
internal degrees of freedom is involved here. Thus we set $L=3$ and the
specific heat ratio $\gamma=(2+D+L)/(D+L)=1.4$. The computational domain is
$[-0.5,0.5]^{2}$. The Neumann condition is applied on the boundary. Like in
Subsection 6.1, only the continuum and collisionless limits are considered.

The continuum regime is characterised again by $\tau=10^{-4}$. Both the DVD-
DVM and DVD-HSM are tested in this case. We set $N=8$ and the directions
$\bm{l}_{\mathcal{L}}=\left\\{\left(\cos\frac{(2m-1)\pi}{16},\sin\frac{(2m-1)\pi}{16}\right)^{T}\right\\}_{m=1}^{8}$.
For the DVD-DVM, the discrete velocity nodes in each direction are taken as
$\xi_{k}=k-8.5$ for $k=1,...,16$. For the DVD-HSM, we set $M=12$. The SSP2
scheme is applied for both models. The physical domain $[-0.5,0.5]^{2}$ is
divided into a $400\times 400$ uniform mesh. Fig. 7 shows the density contours
at $t=0.25$ simulated by the both models. The shock waves and contact
discontinuities are clearly manifested, which agree reasonably well with the
solutions of kinetic equation in [15] and Euler equation in [27]. We again
remark that the weight $|\bm{\xi}|^{D-1}$ in the DVDM Eq. (2.8) is necessary.
As is revealed in Fig. 8, if such a weight is absent, neither the DVD-DVM nor
DVD-HSM correctly predicts the density contour for the Riemann problem in the
continuum limit.

Figure 7. 2-D Riemann problem with $\tau=10^{-4}$: density contours at
$t=0.25$ simulated by (a) DVD-DVM and (b) DVD-HSM. Let $N=8$ and
$\bm{l}_{m}=\left(\cos\frac{(2m-1)\pi}{16},\sin\frac{(2m-1)\pi}{16}\right)^{T}$
for $m=1,...,8$ in both models. For the DVD-DVM, the discrete velocity nodes
in each direction are taken as $\xi_{k}=k-8.5$ for $k=1,...,16$. For the DVD-
HSM, we have $M=12$.

Figure 8. 2-D Riemann problem with $\tau=10^{-4}$: density contours at
$t=0.25$ simulated by (a) DVD-DVM and (b) DVD-HSM without the weight function
$|\bm{\xi}|^{D-1}$. All other setups are the same as in the previous case.

The collisionless free-molecular regime is characterised with $\tau=10^{4}$.
The analytical results can be found in [15]. All three DVDM submodels are used
in the simulation. The non-equilibrium flow generally requires more elaborated
discretization of the velocity space than the continuum case, while the
absence of shock or discontinuity allows greater sizes of the spatial cells.
Thus, the physical domain $[-0.5,0.5]^{2}$ is discretized into a $80\times 80$
uniform mesh. For the DVD-DVM, we set $N=24$ and the directions
$\bm{l}_{m}=\left(\cos\frac{(2m-1)\pi}{48},\sin\frac{(2m-1)\pi}{48}\right)^{T}$.
The discrete velocity nodes in each direction are taken as $\xi_{k}=0.4k-9.8$
for $k=1,...,48$. Remark that the total number of velocity nodes 1152 is much
smaller than that used in [15] (over 40 000). For the DVD-HSM, we set $N=30$
and $M=14$. The SSP2 scheme is used for both the DVD-DVM and DVD-HSM. For the
DVD-EQMOM, we set $N=30$ and $M=2$. The upwind scheme is employed.

The contours of density, temperature and velocity magnitude at $t=0.15$ are
presented in Figs. 9-11 by using different models. Also plotted are the
analytical solutions (black dashed line). It is clearly seen that the DVD-DVM
yields accurate predictions. In contrast, the DVD-EQMOM and DVD-HSM exhibit
greater errors, especially in the temperature profiles. This may be partly
attributed to the lower-order approximation in the velocity space (i.e., small
values of $M$) or lower-order discretization scheme (i.e., the first-order
upwind scheme for the DVD-EQMOM). Future work is needed to address these
issues.

Figure 9. 2-D Riemann problem with $\tau=10^{4}$ by the DVD-DVM: contours of
(a) density, (b) velocity magnitude and (c) temperature at $t=0.15$. The
analytical solutions are shown as black dashed lines. We set $N=24$ and the
directions are of the similar form as before. Discrete velocities are
$\xi_{k}=0.4k-9.8$ for $k=1,...,48$.

Figure 10. 2-D Riemann problem with $\tau=10^{4}$ by the DVD-EQMOM: contours
of (a) density, (b) velocity magnitude and (c) temperature at $t=0.15$. The
analytical solutions are shown as black dashed lines. $N$ is set to be $30$.

Figure 11. 2-D Riemann problem with $\tau=10^{4}$ by the DVD-HSM: contours of
(a) density, (b) velocity magnitude and (c) temperature at $t=0.15$. The
analytical solutions are shown as black dashed lines. We set $N=30$ and
$M=14$.

###  6.4. Lid-Driven Cavity Flow

Our last case is the two-dimensional lid-driven flow in a square cavity
$[0,H]^{2}$. The upper wall moves horizontally with a constant speed $u_{w}$
to drive the fluid while the other three walls are fixed. There are two types
of lid-driven cavity flows. The first type is also termed as the microcavity
flow, where the Reynolds number $\rm Re$ is so small that the flow is mainly
characterized by the Knudsen number [34, 14, 24]. The other type with ${\rm
Re}\gg 1$ has been widely studied by either solving the Navier-Stokes equation
[12] or employing the lattice Boltzmann method [19].

In this case, the internal degrees of freedom are neglected, that is, $D=2$
and $L=1$. Our aim is to derive the steady-state flow field and the
simulations start from a static flow ($\bm{U}=\bm{0}$) in equilibrium with a
constant density $\rho=1$ at $t=0$. Let the initial temperatures of both the
fluid and the walls be $\theta_{0}$ and assume that the walls keep this
temperature. Then the upper wall starts to move and drive the fluid in the
cavity. The computation lasts until the flow becomes steady when the
$L^{2}$-norm of the difference of $\bm{U}$ between two consecutive time steps
is smaller than $10^{-6}$.

We first consider the microcavity flows where the Knudsen number ${\rm
Kn}=\frac{\tau}{H}\sqrt{\frac{\pi\theta_{0}}{2}}$. In this case, we set $H=1$,
$\theta_{0}=2.4$ and $u_{w}=0.32$, resulting in a Mach number of 0.16. We thus
tune $\rm Kn$ by taking different values of $\tau$. Only the DVD-DVM is used
to simulate the flow, for which we set $N=30$ and the directions
$\bm{l}_{m}=\left(\cos\frac{(2m-1)\pi}{60},\sin\frac{(2m-1)\pi}{60}\right)^{T}$
($m=1,\dots,30$). The discrete velocity nodes are taken as $\xi_{k}=0.3k-8.85$
for $k=1,...,60$. The cavity $[0,H]^{2}$ is divided into $100\times 100$
uniform cells. The upwind scheme incorporated with diffuse-scattering boundary
laws is applied here.

Fig. 12 depicts the streamlines and the flow vector field for microcavity
flows with various ${\rm Kn}$. A bulk vortex is clearly observed and the
streamlines are almost axisymmetric about the horizontal center $x=0.5H$. As
Kn increases, the height ($y$-value) of the vortex center reduces. These
features were also presented in previous works [34, 14]. Fig. 13 gives a
comparison of the velocity profiles $\bm{U}=(u,v)^{T}$ across the cavity
center with the reference data [24]. Both $u(y)|_{x=0.5H}$ and
$v(x)|_{y=0.5H}$ are plotted together for each $\rm Kn$. It is seen that the
DVD-DVM results are in good agreement with the reference data.

Figure 12. Microcavity flow: velocity streamlines of (a) Kn = 0.1, (b) Kn = 1
and (c) Kn = 8. For discrete directions and velocities, we set $N=30$ and
$\\{\bm{l}_{m}=(\cos\frac{(2m-1)\pi}{60},\sin\frac{(2m-1)\pi}{60})\\}_{m=1}^{30}$
while $\xi_{k}=0.3k-8.85$ for $k=1,...,60$.

Figure 13. Microcavity flow: profiles of $u(y)|_{x=0.5H}$ and $v(x)|_{y=0.5H}$
for various Knudsen numbers. The red solid lines are our DVD-DVM simulations,
and the black dashed lines are the reference data in [24]. For discrete
directions and velocities, we set $N=30$ and
$\\{\bm{l}_{m}=(\cos\frac{(2m-1)\pi}{60},\sin\frac{(2m-1)\pi}{60})\\}_{m=1}^{30}$
while $\xi_{k}=0.3k-8.85$ for $k=1,...,60$.

We next consider the flow with high Reynolds numbers $\rm Re$, where

${\rm Re}=\frac{u_{w}H}{\theta_{0}\tau}.$

In this case, we set ${\rm Re}=1000$ by taking $u_{w}=0.2$, $H=1$,
$\theta_{0}=1$, and $\tau=2\times 10^{-4}$. This set of parameters
characterises a nearly incompressible flow. For the DVD-DVM, we set $N=6$ and
the directions
$\bm{l}_{m}=\left(\cos\frac{(2m-1)\pi}{12},\sin\frac{(2m-1)\pi}{12}\right)^{T}$.
The SSP2 scheme is adopted with the bounce-back boundary condition for no-slip
walls. The discrete velocity nodes in each direction are taken as
$\xi_{k}=k-8.5$ for $k=1,...,16$.

Fig. 14 shows the steady-state velocity profiles across the cavity center. The
benchmark data are from [12]. The physical domain $[0,H]^{2}$ is discretized
to uniform cells. It is seen that when the uniform grids get finer (from
$80\times 80$ to $160\times 160$), the simulation results become more accurate
and well captures the highly nonlinear boundary profiles.

Figure 14. Lid-driven cavity flow: (a) profiles of $u(y)|_{x=0.5H}$ and
$v(x)|_{y=0.5H}$ and (b) the streamlines for ${\rm Re}=1000$. Red circles are
benchmark data [12]. The lines are the DVD-DVM results. The green dashed lines
are from a $80\times 80$ uniform grid and the blue solid lines are from a
$160\times 160$ uniform grid. We set $N=6$ and choose
$\\{\bm{l}_{m}\\}_{m=1}^{6}$ as above. The discrete velocities in each
direction are $\xi_{k}=k-8.5$ for $k=1,...,16$. The streamlines are based on
data from the $160\times 160$ uniform grid.

##  7\. Conclusions

In this article, we have proposed a discrete-velocity-direction model (DVDM)
based on the BGK equation with the internal molecular degrees of freedom.
Assuming that the molecule velocity is restricted to a few prescribed
directions but the velocity magnitude is still continuous, a semi-continuous
DVDM is obtained, where the local discrete equilibrium in each direction is
derived by the minimum entropy principle subject to the conservation laws. A
key feature of the new model is the introduction of the weight function
$|\xi|^{D-1}$ in the evaluation of the macroscopic fluid quantities.

This DVDM can be combined with various treatments of 1-D velocity distribution
functions to develop multidimensional spatial-time approximations of the
original BGK equation. Specifically, three spatial-time DVD-submodels are
derived by incorporating the discrete-velocity model (DVM), the 1-D Gaussian-
EQMOM and a Hermite spectral method (HSM). We remark that the DVD-DVM allows
radially-positioned discrete velocity nodes, whereas the DVD-EQMOM and DVD-HSM
can be regarded as alternative multidimensional versions of EQMOM and HSM,
respectively.

The feasibility of three spatial-time models have been verified numerically.
For the numerical tests, the DVD-DVM and DVD-HSM are discretized with the
second-order implicit-explicit Runge-Kutta scheme, while only the first-order
upwind scheme is used for the DVD-EQMOM. Two widely-used limiting gas-solid
boundary conditions, including the diffuse-scattering law and the bounce-back
rule, are properly specified for the DVD-DVM and DVD-EQMOM. Only the Neumann
condition is applied for the DVD-HSM. The numerical results for 1-D and 2-D
Riemann problems, especially in both the hydrodynamic and rarefied limits,
illustrate the ability of the DVDM submodels to capture flow discontinuities.
Furthermore, the simulations of the planar Couette flow and lid-driven cavity
flow agree reasonably well with the benchmark data in a wide range of flow
regimes.

The numerical tests suggest that the DVD-DVM should be used for the rarefied
flows. On the other hand, our numerical results are just preliminary. Better
results are expected by using higher-order numerical schemes for spatial-time
models or by enlarging the order $M$ of the DVD-EQMOM. These and the
simulation of 3-D flows are our ongoing projects.

##  8\. Acknowledgments

This work is supported by the National Key Research and Development Program of
China (Grant no. 2021YFA0719200) and the National Natural Science Foundation
of China (Grant no. 51906122 and 12071246).

## References

  * [1] P. Bahukudumbi, J.H. Park, A. Beskok, A Unified Engineering Model for Steady and Quasi-Steady Shear-Driven Gas Microflows, Microscale Therm. Eng., 7 (2003) 291–315.
  * [2] N. Bellomo, D. Burini, G. Dosi, L. Gibelli, D. Knopoff, N. Outada, P. Terna, M.-E. Virgillito, What is life? A perspective of the mathematical kinetic theory of active particles, Math. Mod. Meth. Appl. Sci. 31 (2021) 1821–1866.
  * [3] P. Bhatnagar, E. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Phys. Rev., 94 (1954) 511–525.
  * [4] C. Broyden, The Convergence of a Class of Double-Rank Minimization Algorithms, IMA J. Appl. Math., 6 (1970) 76–90.
  * [5] C. Chalons, R. Fox, F. Laurent, M. Massot, A. Vi$\acute{e}$, Multivariate Gaussian extended quadrature method of moments for turbulent disperse multiphase flow, Multiscale Model. Simul. 15 (2017) 1553–1583.
  * [6] C.-K. Chu, Kinetic-theoretic description of the formation of a shock wave, Phys. Fluids 8 (1965) 12–22.
  * [7] G. Dimarco, L. Pareschi, Numerical methods for kinetic equations, Acta Numerica (2014) 369–520.
  * [8] F. Fei, P. Jenny, A hybrid particle approach based on the unified stochastic particle Bhatnagar-Gross-Krook and DSMC methods, J. Comput. Phys. 424 (2021) 109858.
  * [9] W.A. Fiveland, Discrete ordinate methods for radiative heat transfer in isotropically and anisotropically scattering media, J. Heat Transfer, 109 (1987) 809–812.
  * [10] R. Fox, A quadrature-based third-order moment method for dilute gas-particle flows, J. Comput. Phys. 227 (2008) 6313–6350.
  * [11] S. Friedlander, Smoke, Dust, and Haze: Fundamentals of Aerosol Dynamics, 2nd ed., Oxford University Press, 2000.
  * [12] U. Ghia, K. N. Ghia, C. T. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method, J. Comput. Phys., 48 (1982) 387–411.
  * [13] H. Grad, On the kinetic theory of rarefied gases, Comm. Pure Appl. Math., 2 (1949), 331–407.
  * [14] 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.
  * [15] Z. Guo, R. Wang, K. Xu, Discrete unified gas kinetic scheme for all Knudsen number flows. II. Thermal compressible case, Phys. Rev. E 91 (2015) 033313.
  * [16] S. Harris, An Introduction to the Theory of the Boltzmann Equation, Dover Publications, New York, 2004.
  * [17] A. Harten, P. Lax, B. Vanleer, On upstream differencing and Godunov-type schemes for hyperbolic conservation laws, SIAM Review 25 (1982) 35–61.
  * [18] Lowell H. Holway, New Statistical Models for Kinetic Theory: Methods of Construction, Phys. Fluids, 9 (1966) 1658–1673.
  * [19] S. Hou, Q. Zou, S. Chen, G. Doolen, A.C. Cogley, Simulation of Cavity Flow by the Lattice Boltzmann Method, J. Comput. Phys., 118 (1995) 329–347.
  * [20] Z. Hu, Z. Cai, Y. Wang, Numerical Simulation of Microflows Using Hermite Spectral Methods, SIAM J. Sci. Comput. 42 (2020) B105–B134.
  * [21] Q. Huang, S. Li, G. Li, Q. Yao. Mechanisms on the size partitioning of sodium in particulate matter from pulverized coal combustion, Combust. Flame, 182 (2017) 313–323.
  * [22] Q. Huang, S.Q. Li, W.-A. Yong, Stability analysis of quadrature-based moment methods for kinetic equations, SIAM J. Appl. Math. 80 (2020) 206–231.
  * [23] Q. Huang, Y. Chen, W.-A. Yong, Discrete-velocity-direction models of BGK-type with minimum entropy: I. basic idea, arxiv: 2206.00572.
  * [24] B. John, X.-J. Gu, D. Emerson, Investigation of Heat and Mass Transfer in a Lid-Driven Cavity Under Nonequilibrium Flow Conditions, Numerical Heat Transfer, Part B: Fundamentals, 58 (2010) 287–303.
  * [25] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E.M. Viggen, The Lattice Boltzmann Method: Principles and Practice, Springer, 2017.
  * [26] A.J.C. Ladd, Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation, J. Fluid Mech., 271 (1994) 285–309.
  * [27] P.D. Lax, X.-D. Liu, Solution of two-dimensional Riemann problems of gas dynamics by positive schemes, SIAM J. Sci. Comput. 19 (1998) 319–340.
  * [28] P. Leopardi, A partition of the unit sphere into regions of equal area and small diameter, Electronic Transactions on Numerical Analysis 25 (2006) 309–327.
  * [29] F. D. Lora-Clavijo, J.P. Cruz-Pérez, F.S. Guzmán, J.A. González, Exact solution of the 1D Riemann problem in Newtonian and relativistic hydrodynamics, Revista Mexicana de Física E, 59 (2013) 28–50.
  * [30] D. Marchisio, R. Fox, Computational Models for Polydisperse Particulate and Multiphase Systems, Cambridge University Press, Cambridge, 2013.
  * [31] L. Mieussens, Discrete velocity model and implicit scheme for the BGK equation of rarefied gas dynamics, Math. Mod. Meth. Appl. Sci. 10 (2000) 1121–1149.
  * [32] L. Mieussens, Discrete-velocity models and numerical schemes for the Boltzmann-BGK equation in plane and axisymmetric geometries, J. Comput. Phys. bf 162 (2000) 429–466.
  * [33] A.B. Morris, D.B. Goldstein, P.L. Varghese, L.M. Trafton, Plume impingement on a dusty lunar surface, AIP Conf. Proc. 1333 (2011) 1187–1192.
  * [34] S. Naris, D. Valougeorgis, The driven cavity flow over the whole range of the Knudsen number, Phys. Fluids 17 (2005) 097106.
  * [35] L. Pareschi, G. Russo, Implicit–Explicit Runge–Kutta Schemes and Applications to Hyperbolic Systems with Relaxation, J. Sci. Comput., 25 (2005) 129–155.
  * [36] B. Perthame, Global existence to the BGK model of Boltzmann equation, J. Differ. Equ. 82 (1989) 191–205.
  * [37] W.F. Phillips, Drag on a small sphere moving through a gas, Phys. Fluids 18 (1975) 1089–1093.
  * [38] E.M. Shakhov, Generalization of the Krook kinetic relaxation equation, Fluid Dyn., 3 (1968) 95–96.
  * [39] F. Sharipov, Rarefied gas dynamics and its applications to vacuum technology, CAS - CERN Accelerator School: Vacuum in Accelerators (2007) 1–14. https://doi.org/10.5170/CERN-2007-003.1.
  * [40] F. Sharipov, Rarefied Gas Dynamics: Fundamentals for Research and Practice, Wiley-VCH, Weinheim, 2016.
  * [41] C.-W. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes, Acta Numerica, 29 (2020) 701–762.
  * [42] P. Welander, On the temperature jump in a rarefied gas, Arkiv. Fysik. 7 (1954) 507–553.
  * [43] N.K. Yamaleev, M.H. Carpenter, Third-order Energy Stable WENO scheme, J. Comput. Phys., 228 (2009) 3025–3047.
  * [44] Z. Zhang, J. Xu, Z. Qi, G. Xi, A discrete velocity direction model for the Boltzmann equation and applications to micro gas flows, J. Comput. Phys. 227 (2008) 5256–5271.