# Bounded flows of dense gases

Sergiu Busuioc<EMAIL_ADDRESS>Department of Physics, West University
of Timişoara
Bd. Vasile Pârvan 4, 300223 Timişoara, Romania  Institute for Advanced
Environmental Research, West University of Timişoara
Bd. Vasile Pârvan 4, 300223 Timişoara, Romania  Victor Sofonea
<EMAIL_ADDRESS>Department of Physics, West University of Timişoara
Bd. Vasile Pârvan 4, 300223 Timişoara, Romania  Center for Fundamental and
Advanced Technical Research, Romanian Academy
Bd. Mihai Viteazul 24, 300223, Timişoara, Romania

###### Abstract

Numerical solutions of the Enskog equation obtained employing a Finite-
Difference Lattice Boltzmann (FDLB) and a Direct Simulation Monte Carlo
(DSMC)-like particle method (PM) are systematically compared to determine the
range of applicability of the simplified Enskog collision operator implemented
in the Lattice Boltzmann framework. Three types of bounded flows of dense
gases – namely the Fourier, the Couette, and the Poiseuille flows – are
investigated for a wide range of input parameters. For low to moderate reduced
density, the proposed FDLB model exhibits commendable accuracy for all bounded
flows tested in this study, with substantially lower computational cost than
the PM method.

Lattice Boltzmann; Simplified Enskog collision operator; Dense gases; Bounded
flows; Discrete Simulation Monte Carlo.

##  I Introduction

In recent decades, the study of rarefied gas flows characterized by non-
negligible values of the Knudsen number ($Kn$), which represents the ratio
between the mean free path of fluid particles in a gas and the characteristic
length of the flow domain, has yielded significant progress. These flows have
been numerically studied using the Boltzmann equation, considering the fluid
constituents as point particles. However, when the mean free path of the fluid
particles becomes comparable to their size, the influence of the finite
molecular size becomes crucial [1]. This scenario arises in various practical
applications, such as gas extraction in unconventional reservoirs [2, 3],
high-pressure shock tubes [4], flows through microfabricated nanomembranes [5]
and single-bubble sonoluminescence [6].

The Enskog equation offers a means to extend the kinetic theory description of
fluids beyond the dilute-gas Boltzmann limit [7, 8, 1, 9, 10]. Unlike the
Boltzmann approach, the Enskog equation considers the finite size of gas
molecules and incorporates the space correlations between colliding molecules,
the molecular mutual shielding, as well as the reduction of the available
volume. This equation can be solved numerically using probabilistic or
deterministic methods, similar to the Boltzmann equation. Deterministic
approaches, such as the Monte Carlo quadrature method [11], the fast spectral
method [12, 2] and the Fokker-Planck approximation [13, 14], have been
employed to solve the Enskog equation in recent years. Furthermore,
probabilistic methods have emerged following the success of the Direct
Simulation Monte Carlo method (DSMC) [15] developed by Alexander et al. [16],
Montanero et al. [17] and Frezzotti [18].

Over time, the Enskog equation has been utilized to investigate the properties
of dense gases composed of hard spheres near solid walls in micro- and nano-
channels [19, 20, 18, 21, 22, 23, 24]. Its extension to weakly attracting
hard-sphere systems has proven successful in describing liquid-vapor flows of
monoatomic [25, 26, 27, 28] and polyatomic fluids [29, 30] or mixtures [31],
the formation and the rupture of liquid menisci in nanochannels [32], as well
as the growth/ collapse of spherical nano-droplets/bubbles [33].

Although the aforementioned methods are known for their reliability and
accuracy, their high computational costs often make them impractical for
various applications. To address this limitation, a convenient approach is to
simplify the non-local Enskog collision integral by expanding it into a Taylor
series around a specific point $\bm{x}$ in the coordinate space. The first
term of this expansion corresponds to the conventional Boltzmann collision
operator. Additionally, the second term can be further simplified by replacing
the distribution function with the local equilibrium distribution function, a
procedure which holds when the fluid is close to equilibrium [8, 9]. This
simplification has been employed in Lattice Boltzmann (LB) models to
investigate non-ideal gases [34, 35, 36] and multiphase flows by incorporating
long-range attractive forces [37].

More recently, the simplified Enskog collision operator has been successfully
implemented in various solvers, including the discrete velocity method [38],
the discrete unified gas kinetic scheme (DUGKS) [39], the double-distribution
LB model [40] and the discrete Boltzmann method [41, 42, 43]. By employing the
simplified Enskog collision operator, these solvers offer computationally
efficient alternatives for investigating micro-scale flow phenomena while
maintaining reasonable accuracy.

In this paper, the recently introduced Lattice Boltzmann model for dense gas
flows [43], based on the full-range Gauss-Hermite quadrature, is further
developed to account for the bounded flow of the Enskog gas between parallel
walls. In order to tackle the wall-induced discontinuity, which becomes
effective at non-negligible values of the Knudsen number, the half-range
Gauss-Hermite quadrature method [44, 45, 46, 47] is used to reduce the
numerical errors, as well as the computational costs.

Throughout the paper, we use the non-dimensionalization procedure based on the
reference quantities described in Ref. [46]: $L_{\text{ref}}$ (length),
$n_{\text{ref}}$ ( particle number density) and $T_{\text{ref}}$
(temperature). Accordingly, the reference momentum is defined as
$p_{\text{ref}}=\sqrt{m_{\text{ref}}k_{B}T_{\text{ref}}}$ and the reference
time is given by $t_{\text{ref}}=m_{\text{ref}}L_{\text{ref}}/p_{\text{ref}}$,
where $m_{\text{ref}}$ represents the mass of a fluid particle.

This paper is organized as follows. In sec. II, the simplified Enskog equation
is presented along with the Shakhov collision term. The Finite Difference
Lattice Boltzmann (FDLB) model used to numerically solve the simplified Enskog
equation when the flow domain is bounded by parallel walls is introduced in
Sec. III. This model relies on half-range Gauss-Hermite quadratures in order
to account for the boundary-induced discontinuities. The computer simulation
results are reported in Sec. IV, which has three subsections dedicated to the
Fourier, the Couette, and the Poiseuille flow, respectively. We conclude the
paper in Sec.V.

For the convenience of the reader, further information about the full-range
Gauss-Hermite quadrature and a comparison to the half-range quadrature results
is presented in Appendix A, the heat flux components evaluation in the PM and
the local maxima in the total heat flux are discussed in Appendix B, while
Appendix C briefly presents the particle method (PM) of solving the Enskog
equation [18], which is systematically used to validate the FDLB results for
the bounded flows listed above.

##  II The Enskog equation

The Enskog equation, proposed in 1922 [7], describes the evolution of a system
consisting of rigid spherical molecules. Unlike Boltzmann’s equation, which
assumes molecules as point-like particles subjected to local collisions,
Enskog’s equation considers the volume of fluid particles (i.e., molecules).
This volume restricts the free movement space available to each particle,
leading to an increased number of collisions. Additionally, the collisions
between particles are non-local, occurring when the centers of the two
colliding molecules are separated by one molecular diameter. The Enskog
equation can be written as follows [8, 9, 10]:

$\frac{\partial f}{\partial
t}+\frac{\bm{p}}{m}\cdot\bm{\nabla}_{\bm{x}}f+\bm{F}\cdot\bm{\nabla}_{\bm{p}}f=J_{E}$
(1)

where $m$ denotes the particle mass, $\bm{F}=m\bm{a}$ represents the external
body force and $f(\bm{x},\bm{p},t)$ is the single-particle distribution
function. At time $t$, the distribution function $f$ provides the number of
particles located within the phase space volume $d\bm{x}d\bm{p}$ centered at
the point $(\bm{x},\bm{p})$. The right-hand side of the equation is given by
the Enskog collision operator $J_{E}$, expressed as:

$J_{E}=\int\left\\{\chi\left({\bm{x}}+\frac{\sigma}{2}{\bm{k}}\right)f({\bm{x}},\bm{p^{*}})f({\bm{x}}+\sigma{\bm{k}},\bm{p_{1}^{*}})\right.\\\
-\left.\chi\left({\bm{x}}-\frac{\sigma}{2}{\bm{k}}\right)f({\bm{x}},\bm{p})f({\bm{x}}-\sigma{\bm{k}},\bm{p_{1}})\right\\}\sigma^{2}({\bm{p_{r}}}\cdot{\bm{k}})d{\bm{k}}d{\bm{p_{1}}}$
(2)

In the above equation, $\sigma$ represents the molecular diameter,
$\bm{p_{r}}=\bm{p_{1}}-\bm{p}$ is the relative momentum, and $\bm{k}$ is the
unit vector specifying the relative position of the two colliding particles.
The time dependence of the distribution function is omitted for brevity.

The contact value of the pair correlation function $\chi$ incorporates the
effect of the molecular diameter $\sigma$ on the collision frequency. In the
standard Enskog theory (SET), $\chi\equiv\chi_{\text{\tiny SET}}$ is evaluated
at the contact point of two colliding particles in a fluid assumed to be in
uniform equilibrium [9]. An approximate but accurate expression for
$\chi_{\text{\tiny SET}}$, namely:

$\chi_{\text{\tiny
SET}}[n]=\frac{1}{nb}\left(\frac{P^{hs}}{nk_{B}T}-1\right)=\frac{1}{2}\frac{2-\eta}{(1-\eta)^{3}},$
(3)

is derived from the equation of state (EOS) for the hard-sphere fluid

$P^{hs}=nk_{B}T\frac{1+\eta+\eta^{2}-\eta^{3}}{(1-\eta)^{3}}$ (4)

proposed by Carnahan and Starling [48]. Here $n$ represents the particle
number density, $\eta=b\rho/4$ is the reduced particle density
($b=2\pi\sigma^{3}/3m$), $P^{hs}$ is the pressure of the hard-sphere fluid,
$k_{B}$ is the Boltzmann constant and $T$ is the temperature. The square
brackets in Eq. (3) indicate a functional dependence.

In the revised Enskog theory (RET), the fluid is considered to be in a non-
uniform equilibrium state [9, 10, 49], hence the particle number density is
position dependent. In this case, an effective approximation for the radial
distribution function is obtained using the Fischer-Methfessel (FM)
prescription [50]. This approach involves the replacement in Eq.(3) of the
actual value of the particle number density $n$ with the average particle
density $\overline{n}$ computed over a spherical volume of radius $\sigma$,
centered at $\bm{x}-\frac{\sigma}{2}{\bm{k}}$:

$\chi_{\mbox{\tiny RET-
FM}}\left[n\Big{(}\bm{x}-\frac{\sigma}{2}{\bm{k}}\Big{)}\right]=\chi_{\mbox{\tiny
SET}}\left[\overline{n}\Big{(}\bm{x}-\frac{\sigma}{2}{\bm{k}}\Big{)}\right].$
(5)

The average particle density $\overline{n}$ is given by:

$\overline{n}(\bm{x})=\frac{3}{4\pi\sigma^{3}}\int_{\mathbb{R}^{3}}n(\bm{x}^{\prime})w(\bm{x},\bm{x}^{\prime})\,d\bm{x}^{\prime},\quad
w(\bm{x},\bm{x}^{\prime})=\left\\{\begin{array}[]{cc}1,&\qquad\|\bm{x}^{\prime}-\bm{x}\|<\sigma,\\\
0,&\qquad\|\bm{x}^{\prime}-\bm{x}\|\geq\sigma.\end{array}\right.$ (6)

In the rest of the paper the subscript RET-FM and the functional dependence of
$\chi_{\mbox{\tiny RET-FM}}[n(\bm{x}-\frac{\sigma}{2}{\bm{k}})]$ will be
omitted for brevity.

The Enskog collision operator $J_{E}$ in Eq. (2) can be seen as a
generalization of the Boltzmann collision operator to account for particles
with spatial extent. When the molecular diameter $\sigma$ approaches zero, the
contact value of the pair correlation function approaches unity
($\chi\rightarrow 1$), which recovers the Boltzmann collision operator [8, 9].

###  II.1 Simplified Enskog collision operator

Assuming that the contact value of the pair correlation function $\chi$ and
the distribution functions $\\{f^{*}\equiv
f({\bm{x}},\bm{p^{*}}),f_{1}^{*}\equiv
f({\bm{x}}+\sigma{\bm{k}},\bm{p_{1}^{*}}),f\equiv
f({\bm{x}},\bm{p}),f_{1}\equiv f({\bm{x}}-\sigma{\bm{k}},\bm{p_{1}})\\}$,
which appear in the Enskog collision integral $J_{E}$ given in Eq.(2), are
smooth around the contact point $\bm{x}$, we can approximate these functions
using a Taylor series expansion around $\bm{x}$. The simplified Enskog
collision operator is obtained after retaining the resulting expansion of
$J_{E}\approx\,J_{0}+J_{1}$ up to first-order gradients, namely [8, 9]:

$\displaystyle J_{0}\equiv J_{0}[f]$ $\displaystyle=$
$\displaystyle\chi\int(f^{*}f_{1}^{*}-ff_{1})\sigma^{2}({\bm{p_{r}}}\cdot{\bm{k}})d{\bm{k}}d{\bm{p_{1}}}$
(7) $\displaystyle J_{1}\equiv J_{1}[f]$ $\displaystyle=$
$\displaystyle\chi\sigma\int\bm{k}(f^{*}\bm{\nabla}f_{1}^{*}+f\bm{\nabla}f_{1})\sigma^{2}({\bm{p_{r}}}\cdot{\bm{k}})d{\bm{k}}d{\bm{p_{1}}}$
(8) $\displaystyle+$
$\displaystyle\frac{\sigma}{2}\int\bm{k}\bm{\nabla}\chi(f^{*}f_{1}^{*}+ff_{1})\sigma^{2}({\bm{p_{r}}}\cdot{\bm{k}})d{\bm{k}}d{\bm{p_{1}}}$

The functions $f^{*},f_{1}^{*},f,f_{1}$ and $\chi$ in the two equations above
are evaluated at the point ${\bm{x}}$.

The term $J_{0}[f]$ corresponds to the conventional collision term of the
Boltzmann equation multiplied by $\chi$ and is treated as such by applying the
relaxation time approximation. In this study, we employ the Shakhov collision
term [51, 52]:

$J_{0}[f]=-\frac{1}{\tau}(f-f^{S}),$ (9)

where $\tau$ represents the relaxation time and $f^{S}$ is the equilibrium
Maxwell-Boltzmann distribution multiplied by a correction factor [51, 52, 53,
54]:

$f^{S}=f^{\text{\tiny
MB}}\left[1+\frac{1-\text{Pr}}{P_{i}k_{B}T}\left(\frac{\bm{\xi}^{2}}{5mk_{B}T}-1\right)\bm{\xi}\cdot\bm{q}\right]$
(10)

The Maxwell-Boltzmann distribution $f^{\text{\tiny MB}}$ is defined as

$f^{\text{\tiny MB}}=\frac{n}{(2m\pi
k_{B}T)^{3/2}}\exp{\left(-\frac{\bm{\xi}^{2}}{2mk_{B}T}\right)}$ (11)

and the heat flux $\bm{q}$ is obtained using:

$\bm{q}=\int d^{3}pf\frac{\bm{\xi}^{2}}{2m}\frac{\bm{\xi}}{m},$ (12)

where $\bm{\xi}=\bm{p}-m\bm{u}$ represents the peculiar momentum,
$\text{Pr}=c_{P}\mu/\lambda$ denotes the Prandtl number, $c_{P}=5k_{B}/2m$ is
the specific heat at constant pressure, $\mu$ is the shear viscosity,
$\lambda$ is the thermal conductivity and $P_{i}=\rho RT=nk_{B}T$ is the ideal
gas equation of state, with $R$ being the specific gas constant. It is
important to note that although the Shakhov model does not guarantee non-
negativity of the correction factor and the H-theorem has not been proven, the
model has been successfully implemented and its accuracy has been tested
through comparisons with experimental [55, 56, 53] or DSMC [46, 57, 54, 58]
results.

The second term of $J_{E}$, denoted as $J_{1}[f]$, can be approximated by
replacing the distribution functions ($f^{*},f_{1}^{*},f,f_{1}$) with their
corresponding equilibrium distribution functions. By using $f^{*,{\text{\tiny
MB}}}f_{1}^{*,\text{\tiny MB}}=f^{\text{\tiny MB}}f^{\text{\tiny MB}}_{1}$ and
integrating over $\bm{k}$ and $\bm{p_{1}}$, we obtain [8, 9]:

$J_{1}[f]\approx J_{1}[f^{\text{\tiny MB}}]=-b\rho\chi f^{\text{\tiny
MB}}\left\\{\bm{\xi}\cdot\left[\bm{\nabla}\ln(\rho^{2}\chi
T)+\frac{3}{5}\left(\zeta^{2}-\frac{5}{2}\right)\bm{\nabla}\ln
T\right]\right.\\\ \left.+\frac{2}{5}\left[2\bm{\zeta}\bm{\zeta}\bm{:\nabla
u}+\left(\zeta^{2}-\frac{5}{2}\right)\bm{\nabla\cdot u}\right]\right\\}$ (13)

where $\bm{\zeta}=\bm{\xi}/\sqrt{2RT}$. By incorporating these approximations,
the Enskog equation Eq. (1) can be expressed as:

$\frac{\partial f}{\partial
t}+\frac{\bm{p}}{m}\bm{\nabla}_{\bm{x}}f+\bm{F}\cdot\bm{\nabla}_{\bm{p}}f=-\frac{1}{\tau}(f-f^{S})+J_{1}[f^{\text{\tiny
MB}}]$ (14)

The macroscopic quantities can be determined by evaluating the corresponding
moments of the distribution function:

$\begin{pmatrix}n\\\ \rho\bm{u}\\\ \frac{3}{2}nk_{B}T\\\ {\bm{\Pi}}^{kin}\\\
{\bm{q}}^{kin}\end{pmatrix}=\int d^{3}p\begin{pmatrix}1\\\ \bm{p}\\\
\frac{\bm{\xi}^{2}}{2m}\\\ \frac{\bm{\xi}\bm{\xi}}{m}\\\
\frac{\bm{\xi}\bm{\xi}^{2}}{2m^{2}}\end{pmatrix}f$ (15)

where $\rho=mn$, while ${\bm{\Pi}}^{kin}$ and ${\bm{q}}^{kin}$ denote the
kinetic stress tensor and the kinetic heat flux, respectively. Multiplying the
Enskog equation Eq.(1) by the collision invariants $1,{\bm{p}}$, and
${\bm{p}}^{2}/2m$, and integrating over the momentum space, we obtain the
following conservation equations for mass, momentum, and energy [9]:

$\displaystyle\frac{D\rho}{Dt}+\rho\nabla\cdot\bm{u}$ $\displaystyle=0$ (16a)
$\displaystyle\rho\frac{D{\bm{u}}}{Dt}+\nabla P$
$\displaystyle=-\nabla\cdot{\bm{\Pi}}$ (16b)
$\displaystyle\rho\frac{De}{Dt}+P\nabla\cdot\bm{u}$
$\displaystyle=-\nabla\cdot\bm{q}+{\bm{\Pi}}\bm{:}\nabla\bm{u}$ (16c)

Here, $D/Dt=\partial_{t}+\bm{u}\cdot\nabla$ represents the material
derivative, and $P=P_{i}(1+b\rho\chi)$ denotes the equation of state for a
non-ideal gas. The heat flux $\bm{q}$ and the viscous part of the stress
tensor $\bm{\Pi}$ are given by:

$\displaystyle\bm{q}=-\lambda\nabla T,$ (17a)
$\displaystyle\bm{\Pi}=-\mu_{v}\mathcal{I}\bm{\nabla}\cdot\bm{u}-\mu\left(\bm{\nabla
u}+(\bm{\nabla u})^{T}-\frac{2}{3}\mathcal{I}\bm{\nabla}\cdot\bm{u}\right)$
(17b)

where $\mathcal{I}$ represents the identity matrix. The bulk viscosity
$\mu_{v}$, the shear viscosity $\mu$ and the thermal conductivity $\lambda$,
which appear in Eqs. (17), are given by [9]:

$\displaystyle\mu_{v}$ $\displaystyle=$
$\displaystyle\frac{16}{5\pi}\mu_{0}b^{2}\rho^{2}\chi\,,$ (18)
$\displaystyle\mu=\tau P_{i}$ $\displaystyle=$
$\displaystyle\mu_{0}\left[\frac{1}{\chi}+\frac{4}{5}(b\rho)+\frac{4}{25}\left(1+\frac{12}{\pi}\right)(b\rho)^{2}\chi\right]\,,$
(19) $\displaystyle\lambda=\frac{5k_{B}}{2m}\frac{\tau P_{i}}{\text{Pr}}$
$\displaystyle=$
$\displaystyle\lambda_{0}\left[\frac{1}{\chi}+\frac{6}{5}(b\rho)+\frac{9}{25}\left(1+\frac{32}{9\pi}\right)(b\rho)^{2}\chi\right]\,.$
(20)

In these equations, $\mu_{0}$ and $\lambda_{0}$ represent the viscosity
coefficient and thermal conductivity for hard-sphere molecules at temperature
$T$, namely [9]:

$\mu_{0}=\frac{5}{16\sigma^{2}}\sqrt{\frac{mk_{B}T}{\pi}},\quad\lambda_{0}=\frac{75k_{B}}{64m\sigma^{2}}\sqrt{\frac{mk_{B}T}{\pi}}$
(21)

For a dense gas, the Prandtl number Pr is expressed as [9]:

$\text{Pr}=\frac{2}{3}\,\frac{1+\frac{4}{5}b\rho\chi+\frac{4}{25}\left(1+\frac{12}{\pi}\right)(b\rho\chi)^{2}}{1+\frac{6}{5}b\rho\chi+\frac{9}{25}\left(1+\frac{32}{9\pi}\right)(b\rho\chi)^{2}}$
(22)

The dilute limit corresponds to $\text{Pr}=2/3$. The Chapman-Enskog expansion
of Eq. (14) provides relationships between the relaxation time $\tau$ and the
transport coefficients. In this context, the relaxation time $\tau$ is
expressed as:

$\tau=\frac{\mu}{P_{i}}$ (23)

The quantity $\mu$ encompasses both kinetic and potential contributions, which
account for the flow of molecules and the collisional effects on the transfer
of momentum and energy in the gas [8, 9]. The relaxation time approximation
effectively captures the collisional transfer resulting from non-local
molecular collisions. It is worth noting that the viscosity of a dense gas
with a fixed reduced density $\eta$ can be adjusted by varying the molecular
diameter $\sigma$ and the number density $n$.

The Knudsen number is defined as the ratio of the mean free path and a
characteristic length (in our case the channel width):

$Kn=\frac{\lambda}{L}=\frac{1}{\sqrt{2}\pi\sigma^{2}n_{0}\chi(n_{0})L}=\frac{1}{6\sqrt{2}\eta_{0}\chi(\eta_{0})R}$
(24)

where $R=L/\sigma$ is the confinement ratio [12, 23, 24].

This paper primarily concentrates on benchmarking the Fourier, Couette, and
Poiseuille flow cases. In these scenarios, the steady flow either lacks bulk
motion or exhibits motion perpendicular to the direction in which gas density
varies. Consequently, the bulk viscosity does not have an impact in these
cases.

The model equations involving the simplified Enskog collision operator $J_{1}$
Eq. (13) are formulated using only a limited number of low-order derivatives,
resulting in the omission of information contained in higher-order terms. As a
consequence, the high-order information excluded in $J_{1}$, which is not
present in the collisional momentum and energy transfer, is reintroduced in
the kinetic transfer of momentum and energy through the relaxation time (23)
and the Prandtl number (22) in the collision term $J_{0}[f]$ (9). This ensures
that the total stress tensor and heat flux derived from the current kinetic
model align with those obtained from the Enskog equation, at least up to the
first-order approximation[9, 24]. As such, when we will compare the heat
fluxes in the Fourier, Couette, and Poiseuille setups, the simulation results
obtained using the PM will contain the total heat flux $q_{x}$, the kinetic as
well as the potential contributions defined in Appendix B.

###  II.2 Reduced distributions

In the channel flows examined in this paper, the dynamics along the $z$
direction is straightforward. Furthermore, in the heat transfer (i.e.,
Fourier) problem, also the dynamics along the y-axis is straightforward. In
this context, it is advantageous to integrate out the trivial degrees of
freedom in the momentum space at the level of the model equation.

####  II.2.1 1D flows

In the Fourier flow, the dynamics along the $y$ and $z$ directions is trivial.
After integrating along these Cartesian axes, two reduced distribution
functions, namely $\phi$ and $\theta$, are introduced [59, 53, 60, 46, 47,
61]:

$\displaystyle\phi_{\text{1D}}({\bm{x}},p_{x},t)$ $\displaystyle=\int
dp_{y}dp_{z}f({\bm{x}},{\bm{p}},t),$ (25)
$\displaystyle\theta_{\text{1D}}({\bm{x}},p_{x},t)$ $\displaystyle=\int
dp_{y}dp_{z}\frac{p_{y}^{2}+p_{z}^{2}}{m}f({\bm{x}},{\bm{p}},t)$ (26)

In this case, the macroscopic quantities are given by:

$\displaystyle\begin{pmatrix}n\\\ \rho u_{x}\\\ \Pi_{xx}\end{pmatrix}$
$\displaystyle=\int dp_{x}\begin{pmatrix}1\\\ p_{x}\\\
\frac{\xi_{x}^{2}}{m}\end{pmatrix}\phi_{\text{1D}},$ (27a)
$\displaystyle\begin{pmatrix}\frac{3}{2}nk_{B}T\\\ q_{x}\end{pmatrix}$
$\displaystyle=\int dp_{x}\begin{pmatrix}1\\\
\frac{\xi_{x}}{m}\end{pmatrix}\left(\frac{\xi_{x}^{2}}{2m}\phi_{\text{1D}}+\frac{1}{2}\theta_{\text{1D}}\right)$
(27b)

and the evolution equations for the reduced distribution functions become:

$\frac{\partial}{\partial t}\begin{pmatrix}\phi_{\text{1D}}\\\
\theta_{\text{1D}}\end{pmatrix}+\frac{p_{x}}{m}\frac{\partial}{\partial
x}\begin{pmatrix}\phi_{\text{1D}}\\\
\theta_{\text{1D}}\end{pmatrix}=-\frac{1}{\tau}\begin{pmatrix}\phi_{\text{1D}}-\phi_{\text{1D}}^{S}\\\
\theta_{\text{1D}}-\theta_{\text{1D}}^{S}\end{pmatrix}+\begin{pmatrix}J_{1}^{\phi_{\text{1D}}}\\\
J_{1}^{\theta_{\text{1D}}}\end{pmatrix}$ (28)

In the above equations, $\phi_{\text{1D}}^{S}$ and $\theta_{\text{1D}}^{S}$
are given by:

$\displaystyle\phi^{S}=f_{x}^{\text{\tiny
MB}}\left[1+\frac{1-\text{Pr}}{5P_{i}mk_{B}T}\left(\frac{\xi_{x}^{2}}{mk_{B}T}-3\right)\xi_{x}q_{x}\right],$
(29) $\displaystyle\theta^{S}=2k_{B}Tf_{x}^{\text{\tiny
MB}}\left[1+\frac{1-\text{Pr}}{5P_{i}mk_{B}T}\left(\frac{\xi_{x}^{2}}{mk_{B}T}-1\right)\xi_{x}q_{x}\right]$
(30)

where

$f_{x}^{\text{\tiny MB}}=\frac{n}{(2m\pi
k_{B}T)^{1/2}}\exp{\left(-\frac{\xi_{x}^{2}}{2mk_{B}T}\right)}$ (31)

while the first order corrections $J_{1}^{\phi}$ and $J_{1}^{\theta}$ are:

$J_{1}^{\phi_{\text{1D}}}=-\left[\xi_{x}\partial_{x}\ln\chi+2\xi_{x}\partial_{x}\ln\rho+{5}\left(\frac{\xi_{x}^{2}}{mk_{B}T}-1\right)\partial_{x}u_{x}\right.\\\
\left.+\frac{3}{10}\left(\frac{\xi_{x}^{3}}{m^{2}k_{B}T}+\frac{\xi_{x}}{3m}\right)\partial_{x}\ln
T\right]f_{x}^{\text{\tiny MB}}b\rho\chi$ (32a)
$J_{1}^{\theta_{\text{1D}}}=-\left[\xi_{x}\partial_{x}\ln\chi+2\xi_{x}\partial_{x}\ln\rho+\frac{3}{5}\left(\frac{\xi_{x}^{2}}{mK_{B}T}-\frac{1}{3}\right)\partial_{x}u_{x}\right.\\\
\left.+\frac{3}{10}\left(\frac{\xi_{x}^{3}}{m^{2}k_{B}T}+\frac{7\xi_{x}}{3m}\right)\partial_{x}\ln
T\right]2mk_{B}Tf_{x}^{\text{\tiny MB}}b\rho\chi$ (32b)

####  II.2.2 2D flows

To reduce the computational costs when simulating the Couette and the
Poiseuille flows, where only the dynamics along the $z$ direction is trivial,
we integrate along this direction in the momentum space and get the reduced
distribution functions [59, 53, 60, 46, 47, 61]:

$\displaystyle\phi_{\text{2D}}({\bm{x}},p_{x},p_{y},t)$ $\displaystyle=\int
dp_{z}f({\bm{x}},{\bm{p}},t),$ (33)
$\displaystyle\theta_{\text{2D}}({\bm{x}},p_{x},p_{y},t)$ $\displaystyle=\int
dp_{z}\frac{p_{z}^{2}}{m}f({\bm{x}},{\bm{p}},t)$ (34)

This way, the macroscopic quantities are given by:

$\displaystyle\begin{pmatrix}n\\\ \rho u_{i}\\\ \Pi_{ij}\end{pmatrix}$
$\displaystyle=\int dp_{x}dp_{y}\begin{pmatrix}1\\\ p_{i}\\\
\xi_{i}\xi_{j}/m\end{pmatrix}\phi_{\text{2D}},$ (35a)
$\displaystyle\begin{pmatrix}\frac{3}{2}nk_{B}T\\\ q_{i}\end{pmatrix}$
$\displaystyle=\int dp_{x}dp_{y}\begin{pmatrix}1\\\
\xi_{i}/m\end{pmatrix}\left(\frac{\xi_{j}\xi_{j}}{2m}\phi_{\text{2D}}+\frac{1}{2}\theta_{\text{2D}}\right)$
(35b)

where the sum over repeated Cartesian indices $j\in\\{x,\,y\\}$ is implicitely
understood. In this case, the evolution equations for the reduced distribution
functions are:

$\frac{\partial}{\partial t}\begin{pmatrix}\phi_{\text{2D}}\\\
\theta_{\text{2D}}\end{pmatrix}+\frac{p_{x}}{m}\frac{\partial}{\partial
x}\begin{pmatrix}\phi_{\text{2D}}\\\
\theta_{\text{2D}}\end{pmatrix}+F_{y}\frac{\partial}{\partial
p_{y}}\begin{pmatrix}\phi_{\text{2D}}\\\
\theta_{\text{2D}}\end{pmatrix}=-\frac{1}{\tau}\begin{pmatrix}\phi_{\text{2D}}-\phi_{\text{2D}}^{S}\\\
\theta_{\text{2D}}-\theta_{\text{2D}}^{S}\end{pmatrix}+\begin{pmatrix}J_{1}^{\phi_{\text{2D}}}\\\
J_{1}^{\theta_{\text{2D}}}\end{pmatrix}$ (36)

where, for brevity, we included only the body force term necessary for the
Poiseuille flow.

In the evolution equations above, $\phi_{\text{2D}}^{S}$ and
$\theta_{\text{2D}}^{S}$ are given by:

$\displaystyle\phi_{\text{2D}}^{S}=f_{2D}^{\text{\tiny
MB}}\left[1+\frac{1-\text{Pr}}{5P_{i}mk_{B}T}\left(\frac{\xi_{x}^{2}+\xi_{y}^{2}}{mk_{B}T}-4\right)(\xi_{x}q_{x}+\xi_{y}q_{y})\right],$
(37a) $\displaystyle\theta_{\text{2D}}^{S}=k_{B}Tf_{2D}^{\text{\tiny
MB}}\left[1+\frac{1-\text{Pr}}{5P_{i}mk_{B}T}\left(\frac{\xi_{x}^{2}+\xi_{y}^{2}}{mk_{B}T}-2\right)(\xi_{x}q_{x}+\xi_{y}q_{y})\right]$
(37b)

where

$f_{2D}^{\text{\tiny MB}}=\frac{n}{(2m\pi
k_{B}T)}\exp{\left(-\frac{\xi_{x}^{2}+\xi_{y}^{2}}{2mk_{B}T}\right).}$ (38)

The first order corrections $J_{1}^{\phi_{2D}}$ and $J_{1}^{\theta_{2D}}$ in
Eqs. (36) are:

$J_{1}^{\phi_{\text{2D}}}=-\left[\xi_{x}\partial_{x}\ln\chi+2\xi_{x}\partial_{x}\ln\rho+\frac{2}{5}\left(\frac{\xi_{x}^{2}}{mk_{B}T}+\frac{\xi_{x}^{2}+\xi_{y}^{2}}{2mk_{B}T}-2\right)\partial_{x}u_{x}\right.\\\
\left.+\frac{3\xi_{x}}{10m}\left(\frac{\xi_{x}^{2}+\xi_{y}^{2}}{mk_{B}T}-\frac{2}{3}\right)\partial_{x}\ln
T\right]f_{2D}^{\text{\tiny MB}}b\rho\chi$ (39a)
$J_{1}^{\theta_{\text{2D}}}=-\left[\xi_{x}\partial_{x}\ln\chi+2\xi_{x}\partial_{x}\ln\rho+\frac{2}{5}\left(\frac{\xi_{x}^{2}}{mk_{B}T}+\frac{\xi_{x}^{2}+\xi_{y}^{2}}{2mk_{B}T}-1\right)\partial_{x}u_{x}\right.\\\
\left.+\frac{3\xi_{x}}{10m}\left(\frac{\xi_{x}^{2}+\xi_{y}^{2}}{mk_{B}T}+\frac{4}{3}\right)\partial_{x}\ln
T\right]2mk_{B}Tf_{2D}^{\text{\tiny MB}}b\rho\chi$ (39b)

##  III Finite-difference Enskog Lattice Boltzmann model

In this section, we outline the FDLB algorithm used to solve Eqs. (28) and
(36) numerically. This algorithm involves two main stages. The first one is
the discretization of the momentum space and the second one is the choice of
the numerical scheme used to handle the advection term in the evolution
equations.

###  III.1 Discretization of the momentum space

In our work, we use Gauss-Hermite quadrature methods of various orders for the
discretization of the momentum space [62, 63, 44, 45, 46, 47]. Unlike Ref.
[43], where the full-range Gauss-Hermite quadrature was adequate for the use
in the Enskog FDLB model for the investigation of problems involving only
periodic boundary conditions, namely the propagation of the sound and the
shock waves in a one-dimensional domain, in this paper we investigate fluid
flow problems bounded by two parallel walls perpendicular to the $x$ axis of
the Cartesian system. As known in the kinetic theory of gases, the presence of
the walls induces a discontinuity of the distribution function, which becomes
more effective at higher values of the Knudsen number. In such cases, the
half-range Gauss-Hermite quadrature method in confined fluid flow was already
proved to be more accurate and more efficient to capture the effects of the
wall-induced discontinuity, when compared to the corresponding full-range
quadrature method [44, 45, 46, 47].

After discretization of the momentum space, the integrals in Eqs.(27) and (35)
are replaced with sums over the $\mathcal{K}$ elements of the discrete
momentum set, i.e. the momentum vectors ${\mathbf{p}}_{\kappa}$ or the
corresponding peculiar momenta ${\bm{\xi}}_{\kappa}$,
$1\leq\kappa\leq\mathcal{K}$:

$\displaystyle\begin{pmatrix}n\\\ \rho u_{i}\\\ \Pi_{i,j}\end{pmatrix}$
$\displaystyle=\sum_{\kappa=1}^{\mathcal{K}}\begin{pmatrix}1\\\
p_{\kappa;i}\\\
\frac{\xi_{\kappa;j}\xi_{\kappa;j}}{m}\end{pmatrix}\phi_{\gamma}({\mathbf{p}}_{\kappa})\,,$
(40a) $\displaystyle\begin{pmatrix}\frac{3}{2}nk_{B}T\\\ q_{i}\end{pmatrix}$
$\displaystyle=\sum_{\kappa=1}^{\mathcal{K}}dp_{\kappa;i}\begin{pmatrix}1\\\
\frac{\xi_{\kappa;i}}{m}\end{pmatrix}\left(\frac{\xi_{\kappa;j}\xi_{\kappa;j}}{2m}\phi_{\gamma}({\mathbf{p}}_{\kappa})+\frac{1}{2}\theta_{\gamma}({\mathbf{p}}_{\kappa})\right)$
(40b)

The indices $i,\,j\,\in{\\{x,\,y\\}}$ denote the Cartesian components of the
vectors $\mathbf{p_{\kappa}}$ and $\bm{\xi_{\kappa}}$, while
$\gamma\in\\{\text{1D,\,2D}\\}$.

In order to take advantage of the geometry of the channel flows considered in
this paper, when investigating the Couette and the Poiseuille flow we solve
the evolution equations of the reduced distribution functions by employing the
mixed quadratures concept, according to which the quadrature is controlled
separately on each axis [44, 61, 54]. The details regarding the 2D LB model
with mixed quadratures are found in Subsec. III.1.1 below.

####  III.1.1 Mixed quadrature LB models with reduced distribution functions

In mixed quadrature lattice Boltzmann (LB) models, the momentum space is
constructed using a direct product rule. This allows for independent
construction of quadratures on each axis, by taking into account the flow
characteristics.

In the Couette and Poiseuille flows considered in this paper, the presence of
diffuse-reflective walls introduces a significant discontinuity in the
distribution functions $\phi_{2D}$ and $\theta_{2D}$ when the Knudsen number
(${Kn}$) becomes sufficiently large. For such values of $Kn$, the efficiency
of the full-range Gauss-Hermite quadrature on the axis perpendicular to the
wall diminishes compared to the half-range Gauss-Hermite quadrature [45, 44,
61, 54]. On the other hand, for axes not bounded by walls, a relatively low-
order quadrature (usually a full-range Gauss-Hermite quadrature) is adequate.
The detailed information about the full-range Hermite polynomials is provided
in Appendix A, which also includes a comparison of simulation results obtained
using both full and half range approaches.

When using reduced distribution functions to address 3D flow problems that are
homogeneous along the $z$ axis and bounded by walls perpendicular to the $x$
axis, as encountered in the subsequent sections, appropriate LB models employ
a half-range Gauss-Hermite quadrature of order $Q_{x}^{h}$ along the $x$ axis
and a full-range Gauss-Hermite quadrature of order $Q_{y}$ along the $y$ axis.
The number of momentum vectors used in the FDLB model is
${\mathcal{K}}=\mathcal{Q}_{x}\times\mathcal{Q}_{y}$, where
$\mathcal{Q}_{x}=2Q_{x}^{h}$ and $\mathcal{Q}_{y}=Q_{y}$. This technique
allows one to reduce the computational costs [45, 46, 54]. Note that, when
simulating the Fourier flow there is no variation of the reduced distribution
functions along the $y$ axis. In this case, the sole use of the half-range
Gauss-Hermite quadrature on the $x$ axis is sufficient and the mixed
quadratures are no longer needed.

####  III.1.2 Half-range Gauss-Hermite quadrature

To ensure the recovery of the half-range moments $M_{s}^{\pm}$ and
$M_{s}^{(\mathrm{eq}),\pm}$ of finite order $s$ for the distribution functions
$f$, half-range quadratures are employed on the $x$ axis. The moments are
defined as follows:

$\begin{pmatrix}M_{s}^{+}\\\
M_{s}^{{(\mathrm{eq})},+}\rule{0.0pt}{19.91692pt}\end{pmatrix}=\int_{0}^{\infty}dp\,\begin{pmatrix}f(p)\\\
f^{\mathrm{eq}}(p)\rule{0.0pt}{19.91692pt}\end{pmatrix}p^{s},\qquad\begin{pmatrix}M_{s}^{-}\\\
M_{s}^{{(\mathrm{eq})},-}\rule{0.0pt}{19.91692pt}\end{pmatrix}=\int_{-\infty}^{0}dp\,\begin{pmatrix}f(p)\\\
f^{\mathrm{eq}}(p)\rule{0.0pt}{19.91692pt}\end{pmatrix}p^{s}.$ (41)

The half-range Gauss-Hermite quadrature on the $x$ axis is defined with
respect to the weight function $\omega(p)$:

$\displaystyle\omega(p_{x})$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\pi}}e^{\bar{p}_{x}^{2}/2},\quad\bar{p}_{x}\equiv
p_{x}/p_{0}$ (42a)
$\displaystyle\int_{0}^{\infty}dp\,\omega(\bar{p}_{x})P_{s}(\bar{p}_{x})$
$\displaystyle\simeq$
$\displaystyle\sum_{k_{x}=1}^{Q_{x}^{h}}w(p_{k_{x}})P_{s}(p_{k_{x}}),$ (42b)

where $P_{s}$ represents a polynomial of order $s$, $p_{0}$ is a
characteristic momentum scale and the equality (42b) is exact if the number of
quadrature points $Q_{x}^{h}$ satisfies $2Q_{x}^{h}>s$ [62, 64]. In this paper
we set $p_{0}=\sqrt{m_{ref}k_{B}T_{ref}}$. The quadrature points $p_{k_{x}}$
($k_{x}=1,2,\dots,Q_{x}^{h}$) are the positive roots of the half-range Hermite
polynomial $\mathfrak{h}_{Q}(p)$ of order Q. The quadrature weights
$w_{k_{x}}^{h}$ are given by [62, 64, 45, 44]:

$w_{k_{x}}^{h}=\frac{p_{k_{x}}a_{Q}^{2}}{\mathfrak{h}_{Q+1}^{2}(p_{k_{x}})\left[p_{k_{x}}+\mathfrak{h}_{Q}^{2}(0)/\sqrt{2\pi}\right]},$
(43)

where $a_{Q}=\mathfrak{h}_{Q+1,Q+1}/\mathfrak{h}_{Q,Q}$ and
$\mathfrak{h}_{\ell,s}$ represents the coefficient of $p^{s}$ in
$\mathfrak{h}_{\ell}(p)$:

$\mathfrak{h}_{\ell}(p)=\sum_{s=0}^{\ell}\mathfrak{h}_{\ell,s}\bar{p}^{s}.$
(44)

The half-range Hermite polynomials are normalized according to:

$\int_{0}^{\infty}dp\,\omega(p)\mathfrak{h}_{\ell}(p)\mathfrak{h}_{\ell^{\prime}}(p)=\delta_{\ell,\ell^{\prime}}.$
(45)

To apply the half-range Gauss-Hermite quadrature, the distribution function
$f$ and the equilibrium distribution function $f^{\mathrm{eq}}$ must be
expanded using the half-range Hermite polynomials. Since these polynomials are
defined only on half of the momentum axis, $f$ is split using the Heaviside
step function $\theta(p)$:

$f(p)=\theta(p)f^{+}(p)+\theta(-p)f^{-}(p),\qquad\theta(p)=\begin{cases}1,&p>0,\\\
0,&p<0.\end{cases}$ (46)

The functions $f^{+}(p)$ and $f^{-}(p)$ are defined on the positive and
negative momentum semi-axis, respectively, and can be expanded as:

$f^{+}=\frac{\omega(p)}{p_{0}}\sum_{\ell=0}^{\infty}\mathcal{F}^{+}\mathfrak{h}_{\ell}(p),\qquad
f^{-}=\frac{\omega(-p)}{p_{0}}\sum_{\ell=0}^{\infty}\mathcal{F}^{-}\mathfrak{h}_{\ell}(-p).$
(47)

The coefficients $\mathcal{F}^{\pm}$ can be obtained using the orthogonality
of the half-range Hermite polynomials:

$\mathcal{F}^{+}_{\ell}=\int_{0}^{\infty}dp\,f(p)\mathfrak{h}_{\ell}(p),\qquad\mathcal{F}^{-}_{\ell}=\int_{-\infty}^{0}dp\,f(p)\mathfrak{h}_{\ell}(-p).$
(48)

Expanding $f^{\mathrm{eq}}$ in a similar manner, with expansion coefficients
$\mathcal{G}^{\pm}_{\ell}$, the half-range moments can be recovered:

$\begin{pmatrix}M_{s}^{+}\\\
M_{s}^{-}\rule{0.0pt}{19.91692pt}\end{pmatrix}=\sum_{\ell=0}^{\infty}\frac{1}{\ell!}\begin{pmatrix}\mathcal{F}_{\ell}^{+}\\\
\mathcal{F}_{\ell}^{-}\rule{0.0pt}{19.91692pt}\end{pmatrix}\int_{0}^{\infty}dp\,\omega(p)\,\mathfrak{h}_{\ell}(p)(\pm
p)^{s}.$ (49)

Truncating the expansion at $\ell=Q-1$ ensures that a quadrature of order $Q$
can recover the moments for $0\leq s\leq Q$. The discrete momentum set of the
one-dimensional half-range Gauss-Hermite lattice Boltzmann model consists of
$\mathcal{Q}=2Q$ elements, defined as:

$\bar{p}_{k_{x}}=p_{0}p_{k_{x}},\qquad\bar{p}_{k_{x}+Q}=-p_{k_{x}}\qquad(1\leq
k\leq Q_{x}^{h}).$ (50)

Thus, the half-range moments are recovered as:

$M_{s}^{+}=\sum_{k=1}^{Q}f_{k}p_{k_{x}}^{s},\qquad
M_{s}^{-}=\sum_{k=Q+1}^{2Q}f_{k}p_{k_{x}}^{s},$ (51)

where

$f_{k}=\frac{w_{k}^{h}p_{0}}{\omega(p_{k})}f(p_{k}),\qquad
f_{k+Q}=\frac{w_{k}^{h}p_{0}}{\omega(p_{k})}f(-p_{k})\qquad(1\leq k\leq
Q_{x}^{h}).$ (52)

Now let us consider the expansion of $f^{\mathrm{eq}}$ in terms of the half-
range Hermite polynomials. We can write
$f_{x}^{\text{MB}}=\theta(p)g_{+}(p)+\theta(-p)g_{-}(p)$, where

$g_{\pm}=\frac{\omega(\left|p\right|)}{p_{0}}\sum_{\ell=0}^{\infty}\mathcal{G}^{\pm}_{\ell}\mathfrak{h}_{\ell}(\left|p\right|).$
(53)

The expansion coefficients $\mathcal{G}^{\pm}_{\ell}$ can be obtained in a
similar way to Eq. (48).

For more information on one-dimensional lattice Boltzmann models based on
half-range Gauss-Hermite quadratures, we refer the reader to Refs. [44, 61,
54].

###  III.2 Numerical schemes

For reader’s convenience, in this subsection we detail the numerical schemes
employed to solve the evolution equations (28) and (36), namely the third-
order total variation diminishing (TVD) Runge-Kutta method for time-stepping,
the fifth-order WENO-5 advection scheme, and the $6$th order central
difference scheme used for gradient evaluation.

Diffuse boundary conditions are applied on both walls/boundaries. This implies
that the molecules striking the walls are re-emitted according to the
Maxwellian distribution with parameters $T_{w}$ and $U_{w}$, where $T_{w}$ and
$U_{w}$ represent the predetermined wall temperature and velocity,
respectively. The value of $n_{w}$ is determined by satisfying the
impenetrable condition [44, 65].

####  III.2.1 Third-order TVD Runge-Kutta method

After the momentum space discretization, it is convenient to cast the Enskog
equation (14) in the following form:

$\partial_{t}f_{\kappa}=L[f_{\kappa}],\qquad
L[f_{\kappa}]=-\frac{{p}_{\kappa;\,x}}{\,m\,}\cdot\nabla_{x}f_{\kappa}+(\partial_{p_{y}}f)_{\kappa}-\frac{1}{\tau}[f_{{\kappa}}-f^{S}_{{\kappa}}]+J_{1,\kappa}.$
(54)

in order to implement the time-stepping algorithm. In the above the subscript
$\kappa$ refers to the discretized functions corresponding to $\xi_{\kappa}$.

The third-order Runge-Kutta integrator gives the following three-step
algorithm for computing the values of $f_{\kappa}$ at time $t+\delta t$ [66,
67, 68]:

$\displaystyle f_{\kappa}^{(1)}(t)=$ $\displaystyle f_{\kappa}(t)+\delta
t\,L[f_{\kappa}(t)],$ $\displaystyle f_{\kappa}^{(2)}(t)=$
$\displaystyle\frac{3}{4}f_{\kappa}(t)+\frac{1}{4}f_{\kappa}^{(1)}(t)+\frac{1}{4}\delta
t\,L[f_{\kappa}^{(1)}(t)],$ $\displaystyle f_{\kappa}(t+\delta t)=$
$\displaystyle\frac{1}{3}f_{\kappa}(t)+\frac{2}{3}f_{\kappa}^{(2)}(t)+\frac{2}{3}\delta
t\,L[f_{\kappa}^{(2)}(t)].$ (55)

####  III.2.2 WENO-5 advection scheme

The advection term which appears in Eq. (54) above, namely
$p_{\kappa;x}\cdot\nabla_{x}f_{\kappa}/m$ is computed using the Weighted
Essentially Non-Oscillatory scheme of order $5$ (WENO-5) along each coordinate
[69, 70]. We will describe in the following the one-dimensional case. Assuming
that the flow domain is discretized using $1\leq i\leq N$ nodes on the $x$
axis, the advection term becomes:

$\left(\frac{p_{\kappa}}{m}\partial_{x}f_{\kappa}\right)_{\kappa;i}=\frac{\mathcal{F}_{\kappa;i+1/2}-\mathcal{F}_{k;i-1/2}}{\delta
s}$ (56)

where $\mathcal{F}_{k;i+1/2}$ represents the flux of $f$ advected with
velocity $p_{\kappa}/m$ through the interface between the cells centered on
$\bm{x}_{i}$ and $\bm{x}_{i+1}$. The construction of these fluxes is
summarized below, under the assumption of a positive advection velocity
$p_{\kappa}/m>0$. In this case, the flux $\mathcal{F}_{k;i+1/2}$ can be
computed using the following expression [69]:

$\mathcal{F}_{\kappa,i+1/2}=\overline{\omega}_{1}\mathcal{F}^{1}_{\kappa,i+1/2}+\overline{\omega}_{2}\mathcal{F}^{2}_{\kappa,i+1/2}+\overline{\omega}_{3}\mathcal{F}^{3}_{\kappa,i+1/2}.$
(57)

The interpolating functions $\mathcal{F}^{q}_{\kappa,i+1/2}$ ($q=1,2,3$) are
given by:

$\displaystyle\mathcal{F}^{1}_{\kappa,i+1/2}=$
$\displaystyle\frac{p_{\kappa}}{m}\left(\frac{1}{3}f_{\kappa,i-2}-\frac{7}{6}f_{\kappa,i-1}+\frac{11}{6}f_{\kappa,i}\right),$
$\displaystyle\mathcal{F}^{2}_{\kappa,i+1/2}=$
$\displaystyle\frac{p_{\kappa}}{m}\left(-\frac{1}{6}f_{\kappa,i-1}+\frac{5}{6}f_{\kappa,i}+\frac{1}{3}f_{\kappa,i+1}\right),$
$\displaystyle\mathcal{F}^{3}_{\kappa,i+1/2}=$
$\displaystyle\frac{p_{\kappa}}{m}\left(\frac{1}{3}f_{\kappa,i}+\frac{5}{6}f_{\kappa,i+1}-\frac{1}{6}f_{\kappa,i+2}\right).$
(58)

The weighting factors $\overline{\omega}_{q}$ appearing in Eq. (57) are given
by:

$\overline{\omega}_{q}=\frac{\widetilde{\omega}_{q}}{\widetilde{\omega}_{1}+\widetilde{\omega}_{2}+\widetilde{\omega}_{3}},\qquad\widetilde{\omega}_{q}=\frac{\delta_{q}}{\varphi^{2}_{q}}.$
(59)

The ideal weights $\delta_{q}$ are:
$\delta_{1}=\frac{1}{10},\,\delta_{2}=\frac{6}{10},\,\delta_{3}=\frac{3}{10},$
while the smoothness indicators $\varphi_{q}$ can be computed as follows:

$\displaystyle\varphi_{1}=$
$\displaystyle\frac{13}{12}\left(f_{i-2}-2f_{i-1}+f_{i}\right)^{2}+\frac{1}{4}\left(f_{i-2}-4f_{i-1}+3f_{i}\right)^{2},$
$\displaystyle\varphi_{2}=$
$\displaystyle\frac{13}{12}\left(f_{i-1}-2f_{i}+f_{i+1}\right)^{2}+\frac{1}{4}\left(f_{i-1}-f_{i+1}\right)^{2},$
$\displaystyle\varphi_{3}=$
$\displaystyle\frac{13}{12}\left(f_{i}-2f_{i+1}+f_{i+2}\right)^{2}+\frac{1}{4}\left(3f_{i}-4f_{i+1}+f_{i+2}\right)^{2}.$
(60)

where the index $\kappa$ was omitted for brevity. The computation of the
weighting factors $\overline{\omega}_{q}$ (59) implies the division between
the ideal weights $\delta_{q}$ and the smoothness indicators $\varphi_{q}$
(60). To avoid division by $0$ when either one, two, or all three of the
indicators of smoothness vanish, we follow Refs. [71, 61] and compute the
weighting factors $\overline{\omega}_{q}$ directly in the limiting cases when
any of the smoothness indicators vanishes.

####  III.2.3  The 6th order central difference scheme for the gradient
operator

For evaluating the gradients we employ the $6$th order central difference
scheme [72]:

$\partial_{x}Q(x)=\frac{1}{\Delta x}\left[-\frac{1}{60}Q(x-3\Delta
x)+\frac{3}{20}Q(x-2\Delta x)-\frac{3}{4}Q(x-\Delta x)\right.\\\
\left.+\frac{3}{4}Q(x+\Delta x)-\frac{3}{20}Q(x+2\Delta
x)+\frac{1}{60}Q(x+3\Delta x)\right]$ (61)

where $Q\in\\{\ln\rho,u,\ln T\\}$.

##  IV  Numerical results

Our attention is directed towards the study of dense gas flows confined
between two infinite parallel plates which are perpendicular to the $x$
direction. The plates are located at $x=-L/2$ and $L/2$, respectively, while
the diffuse reflection boundary conditions are applied at $x=\pm(L-\sigma)/2$.
The confinement ratio is defined as $R=L/\sigma$, where $L$ is the physical
domain width and $L_{c}=L-\sigma$ is the width of the computational domain.
The FDLB results are compared to the simulation results obtained with the
particle method (DSMC) which is briefly presented in Appendix C. If not stated
otherwise, the time step was set to $\Delta t=10^{-3}$ and the lattice spacing
(cell length for particle method) at $\Delta x=\sigma/100$. A number of $1000$
particles per cell was used in the particle method in order to obtain smooth
profiles of macroscopic quantities.

For our simulations, we have chosen the following value sets of the
confinement ratio $R$ and the mean (initial) reduced density $\eta_{0}$,
namely $R\in\\{4,\,10\\}$ and $\eta_{0}\in\\{0.01,\,0.10,\,0.20\\}$. We have
chosen these values to encompass a wide range of Knudsen numbers, spanning
from the slip to the early transition regime. The objective is to specifically
emphasize the unique characteristics of fluid flow when dense gas effects and
confinement are involved.

In the following subsections, we examine various scenarios to assess the
effectiveness of the proposed model in recovering the flow of a dense gas
bounded by parallel plates. When not stated otherwise, the orders $Q_{x}^{h}$
of the half-range quadratures used in this study for the three flow problems
(Fourier, Couette, and Poiseuille) are listed in Table 1. These values were
obtained by performing a convergence test, for which the following error
evaluation was used [44]:

$\epsilon_{M}^{\delta
x}=\frac{\max\left[M(x)-M_{ref}(x)\right]}{\max[M_{ref}(x)]-\min[M_{ref}(x)]}$
(62)

where $M\in\\{T,u_{y},u_{y}\\}$, for the Fourier, Couette, and Poiseuille
flow, respectively, and $\delta x$ is the lattice spacing which is kept fixed
for the convergence test. The reference profile $M_{ref}$ is obtained using
$Q_{x}^{h}=200$ quadrature points. The convergence test is satisfied when
$\epsilon_{m}^{\delta x=\sigma/100}<0.01$, which represents a $1\%$ error with
respect to the macroscopic quantity maximum variation in the channel. In the
case of the 2D flows a quadrature of order $Q_{y}=5$ was imposed in the
longitudinal direction for both reduced distributions [54].

The typical runtime for a PM simulation is about $3.5\times 10^{3}$ s for
$\eta_{0}=0.01$ and $1.85\times 10^{4}$ s for $\eta_{0}=0.2$, irrespective of
the flow studied, while the FDLB, using a quadrature order of $Q_{x}^{h}=8$,
takes around 17s for the 1D case of Fourier flow, and 100s for the 2D case of
Couette and Poisseuille flow. The running time for FDLB is independent of the
reduced density employed but it is directly proportional to the quadrature
order, resulting in a minimum runtime ratio of $\approx 150$ in the Fourier
flow case ($Q_{x}^{h}=11$ for $\eta_{0}=0.01$), and a minimum runtime ratio
$\approx 10$ in the Poisseuille case with $Q_{x}^{h}=29$ at $eta_{0}=0.01$.
These times were recorded using $R=4$ and a single core of an Intel(R) Xeon(R)
Gold 6330 CPU running at 2.0GHz.

| Fourier |  | Couette |  | Poiseuille |
---|---|---|---|---|---|---
$\eta_{0}\backslash R$ | 4 | 10 | 4 | 10 | 4 | 10
0.01 | 11 | 8 | 8 | 8 | 29 | 15
0.1 | 8 | 8 | 8 | 8 | 8 | 8
0.2 | 8 | 8 | 8 | 8 | 8 | 8
Table 1:  Quadrature order $Q^{h}_{x}$ for the parameters used in this study.

###  IV.1 Dense gas at rest near a reflective wall

Figure 1:  Dense gas at rest near a reflective wall: normnalized density
profiles $\eta/\eta_{0}$ for three values of the mean reduced density
$\eta_{0}\in\\{0.01,\,0.10,\,0.20\\}$ and two values of the confinement ratios
$R=10$ (a) and $R=4$ (b).

To begin, we consider a hard-sphere gas confined between two stationary
parallel walls, kept at the temperature $T_{w}=1$. We first assess the
accuracy of the proposed model in reproducing the reduced density profile near
a wall. Due to the symmetry of the simulated problem, only the half-channel
plot is presented ($x\in[0:L_{c}]$), as it will be also the case in the
Couette and the Poiseuille flows to be discussed further. In Fig. 1 one may
observe that the stationary profile of the reduced density $\eta$ is non-
monotonic near the wall, unlike in the case of a dilute gas. Indeed, it is
important to consider that when a fluid particle is located at a distance less
than $\sigma$ from the wall, a portion of its surface remains protected from
collisions. This occurs because there is insufficient space available for a
second molecule to occupy that region. As a result, the particle is pushed
toward the wall. It is important to note that an oscillating density profile
near the wall is a characteristic feature of dense gases. These density
variations emerge within a region approximately equivalent to the molecular
diameter, and their intensity diminishes as the $\eta_{0}\rightarrow 0$.
Consequently, in the Boltzmann limit, the density becomes spatially uniform.
To assess the accuracy of the FDLB model, we compare the density profiles to
the profiles obtained using the particle method (PM). As expected, the
difference between the FDLB and the PM profiles becomes significantly larger
when the mean reduced density $\eta_{0}$ increases. As seen in Fig.1, a fairly
good agreement between the PM and the FDLB results is obtained up to the mean
value $\eta_{0}=0.1$. For larger values of $\eta_{0}$, substantial differences
can be observed regardless of the value of the confinement ratio $R$. This
value $\eta_{0}$ can be assumed to be the limit of the simplified Enskog
collision model, for a reasonable accuracy.

###  IV.2  Fourier flow

In this subsection, we will analyze the thermal transfer in a dense gas
confined between two infinite parallel plates. The left and the right plate
temperatures are fixed at $T_{L}=T_{0}-\Delta T$ and $T_{R}=T_{0}+\Delta T$
respectively. We will explore two values of the temperature difference $\Delta
T\in\\{0.1,0.5\\}$, corresponding to $T_{R}/T_{L}\in\\{11/9,\,3\\}$. For each
value of the temperature difference $\Delta T$, the simulations were conducted
using three values of the mean reduced density $\eta_{0}\in\\{0.01,0.1,0.2\\}$
and two values of the confinement ratio $R\in\\{4,10\\}$($L_{c}\in\\{3,9\\}$).

Fig.2 presents the results for the smaller temperature difference $\Delta
T=0.1$. The profiles of the normalized reduced density $\eta/\eta_{0}$ and the
temperature $T$ are shown in the first row and the second row, respectively.
Excellent agreement is obtained for both confinement ratios at small values of
density $\eta_{0}$, while at higher values of $\eta_{0}$ the density and
temperature present discrepancies between the two simulation methods,
especially near the walls. Nevertheless, these results are in good relative
agreement.

Figure 3 depicts the outcomes for a significantly larger temperature
difference $\Delta T=0.5$, corresponding to the wall temperature ratio
$T_{R}/T_{L}=3$. The density profiles resemble those observed for smaller
temperature differences. However, in the case of temperature profiles for high
values of the mean reduced density $\eta_{0}$, noticeable discrepancies are
observed also in the bulk of the fluid (the relative error with respect to the
PM results is approximatively $5\%$).

Figure 2:  Fourier flow: Heat transfer at wall temperature difference $\Delta
T=0.1$ and three values of the mean reduced density $\eta_{0}$: Transversal
profiles of the ratio $\eta/\eta_{0}$ (first row) and of the temperature $T$
(second row) for $R=10$ (left column) and $R=4$ (right column).

Figure 3: Fourier flow: Heat transfer at wall temperature difference $\Delta
T=0.5$ and three values of the mean reduced density $\eta_{0}$ : Transversal
profiles of the ratio $\eta/\eta_{0}$ (first row) and of the temperature $T$
(second row) for $R=10$ (left column) and $R=4$ (right column).

According to the conservation laws, the transversal heat flux $q_{x}$ is
constant through the channel. The profiles of the transversal heat flux
$q_{x}$, obtained with both the FDLB and the PM models are presented in Fig.
4. The transversal heat flux is plotted in the left panel for $R\in\\{4,10\\}$
and $\Delta T\in\\{0.1,0.5\\}$. As the dilute gas limit is approached, the two
methods agree perfectly, and the departure of the FDLB results from the PM
results for increased values of $\eta_{0}$ is observed. On the other hand, in
the right panel, we plot the ratio $q_{x}/\eta_{0}$ with respect to various
values of the Knudsen number $Kn$ obtained by varying the width of the
computational domain at a constant value of the reduced density $\eta_{0}$.
The monotonic increase of the heat flux with respect to the Knudsen number
encountered for the dilute gas [73, 74] is no longer present for dense gases
at relatively high values of the reduced density $\eta_{0}$. Indeed, one may
observe a local maximum of the heat flux which emerges in both the FDLB and
the PM models. This maximum is further investigated using the PM in Appendix
B, where the two components of the heat flux (kinetic and potential) are
evaluated.

Figure 4:  Fourier flow: (a) Heat flux $q_{x}$ profiles at temperature
difference of $\Delta T=\\{0.1,0.5\\}$, confinement ratios of $R=\\{4,10\\}$.
(b) Dependence of the ratio $q_{x}/n_{0}$ with respect to the Knudsen number
$Kn$ for 3 values of the reduced density $\eta\in\\{0.01,0.1,0.2\\}$.

###  IV.3  Couette flow

In this subsection, we will analyze the Couette flow of a gas confined between
two infinite parallel plates. The plates move in opposite directions along the
$y$ axis with fixed velocity $U_{0}=\sqrt{k_{B}T/m}$, i.e. $U_{L}=-U_{0}$ and
$U_{R}=U_{0}$. The simulations were conducted for three values of the mean
reduced density, $\eta_{0}\in\\{0.01,0.1,0.2\\}$ and two values of the
confinement ratio $R\in\\{4,10\\}$, corresponding to $L_{c}\in\\{3,9\\}$.

In Fig.5 we plot the transversal density profile on the right half
$[0,\,L_{c}]$ of the computational domain. One can easily recognize the same
feature as in the case of stationary dense gas, namely the formation of a
layer in the proximity of each wall. The layer width is of the order of
$\sigma$ (the molecular diameter). Excellent agreement between the FDLB and
the PM profiles is observed up to $\eta_{0}=0.1$. The corresponding velocity
and temperature profiles are presented in Fig. 6, for the two values of the
confinement ratio considered in this paper. Very good agreement between the
FDLB and the PM profiles is obtained throughout the considered range of the
reduced density $\eta_{0}$, with small discrepancies observed near the wall
when the value of the confinement ratio is large enough ($R=10$). For the
smaller value of the confinement ratio ($R=4$), the discrepancies are larger,
but this is expected behavior due to the severe approximations involved in the
derivation of the simplified Enskog collision operator.

In addition to the transverse component $q_{x}$ of the heat flux, the Couette
flow exhibits a non-zero longitudinal heat flux $q_{y}$, which is a distinct
microfluidics phenomenon. The corresponding results are presented in Figure 7,
where the left and the right panels represent the simulation results conducted
with the values $R=10$ and $R=4$ of the confinement ratio, respectively.
Excellent agreement between the FDLB and the PM results is observed when the
reduced density is small enough ($\eta_{0}=0.01$). For larger values of the
reduced density, the discrepancies become noticeable. The plot specifically
showcases the right half of the channel, where the transverse heat flux values
are positive, and the longitudinal component consistently exhibits negative
values across the parameter range. Due to the fact that the longitudinal heat
flux, $q_{y}$, does not arise from a temperature gradient (known as the direct
phenomenon), its behavior is expected to depend on higher-order transport
coefficients. Hence, it is unsurprising that these cross phenomena are not
accurately captured by the model when $\eta_{0}$ is large enough.

Figure 5:  Couette flow: Normalised reduced density profile $\eta/\eta_{0}$ at
a wall velocity of $U_{w}=1$, three values of the mean reduced density
$\eta_{0}=\\{0.01,0.1,0.2\\}$ and (a) $R=10$ and (b) $R=4$.

Figure 6:  Couette flow: velocity and temperature profiles at a wall velocity
of $U_{w}=1$, three values of the mean reduced density
$\eta_{0}=\\{0.01,0.1,0.2\\}$ and $R=10$ (left column) and $R=4$ (right
column).

Figure 7:  Couette flow: Profiles of both the transversal ($q_{x}>0$) and the
longitudinal ($q_{y}<0$) heat fluxes for the wall velocity $U_{w}=1$, three
values of the mean reduced density $\eta_{0}=\\{0.01,0.1,0.2\\}$ and two
values of the confinement ratio (a) $R=10$ and (b) $R=4$.

###  IV.4  Poiseuille flow

In this section, we examine the Poiseuille flow, which is generated by an
external acceleration $a_{y}$ acting parallel to the plates. The results are
grouped into linear and non-linear flows, based on the driving force
magnitude. Nonlinearity arises when the heat generated by viscous dissipation
cannot be adequately dissipated through the diffusive boundary conditions. To
ensure that the response to the driving force remains within the linear
regime, a specific driving force of $a_{y}=0.001$ has been chosen, as detailed
below.

####  IV.4.1  Linear regime

Figure 8:  Poiseuille flow (linear regime): Density and velocity at an
external acceleration of $a_{y}=0.001$, three values of the mean reduced
density $\eta_{0}=\\{0.01,0.1,0.2\\}$ and $R=10$ (left column) and $R=4$
(right column).

At first, we will look at a flow with a small acceleration $a_{y}=0.001$. In
this case, the flow deviates slightly from the equilibrium, hence only the
velocity field has a statistically significant variation. In Fig. 8 we take
advantage of the symmetry of the Poiseuille flow between parallel walls and
restrict the plot of both the normalized reduced density $\eta/\eta_{0}$ and
the velocity $u_{y}$ on the right half $[0,\,L_{c}]$ of the computational
domain. In this figure, we can observe again the very good agreement between
the FDLB and the PM results for $\eta_{0}=0.01$. For the larger value of
$\eta_{0}$, the density profile has a similar behavior as in the case of dense
gas at rest between two parallel walls in Sec. IV.1, while the FDLB velocity
profile is within a few percent from the PM one. The results are similar for
both values of the confinement ratio used $R\in\\{4,10\\}$.

####  IV.4.2  Non-linear regime

Figure 9:  Poiseuille flow (non-linear regime): Velocity and temperature
profiles at an external acceleration of $a_{y}=0.1$, three values of the mean
reduced density $\eta_{0}=\\{0.01,0.1,0.2\\}$ and $R=10$ (left column) and
$R=4$ (right column).

Figure 10:  Poiseuille flow (non-linear regime): Transversal and longitudinal
heat flux, $q_{x}$ and $q_{y}$, at an external acceleration of $a_{y}=0.1$,
three values of the mean reduced density $\eta_{0}=\\{0.01,0.1,0.2\\}$ and
$R=10$ (left column) and $R=4$ (right column).

Moving to high acceleration of $a_{y}=0.1$, we also encounter transversal
variations in both the temperature and the heat fluxes. As in Couette flow, in
addition to the transverse component $q_{x}$ of the heat flux, the Poiseuille
flow displays a non-zero longitudinal heat flux $q_{y}$, which arises as a
distinct microfluidics phenomenon.

Fig. 9 plots the velocity and the temperature profiles for the acceleration
$a_{y}=0.1$. The density profile has insignificant variations with respect to
the small acceleration case (see Fig. 8). Concerning the velocity profiles,
the FDLB results are in good agreement with the PM curves, with larger
deviations observed in the center of the channel as the reduced density
$\eta_{0}$ is increased. Instead, a surprisingly good match between the FDLB
and the PM profiles is obtained near the walls. Also, the temperature profiles
are again in good agreement: the FDLB profiles exhibit a deviation less than
$5\%$ with respect to the PM profiles throughout the channel.

In Fig. 10 we present the profiles of both the transversal and the
longitudinal heat fluxes. An excellent agreement between FDLB and PM profiles
is seen for the transversal heat flux at low values of the reduced density
$\eta_{0}$. As seen also in the previous plots, at larger values of $\eta_{0}$
the deviations are huge, but this is an expected outcome due to the crude
approximations used in the simplified Enskog collision operator. In the case
of the longitudinal heat flux, the discrepancies are even more pronounced.
Nevertheless, the model still provides a reasonably accurate description of
the flow at low fluid densities.

####  IV.4.3  Mass flow rate

The mass flow rate (MFR) is defined as:

$\dot{m}=\int_{-\frac{L_{c}}{2}}^{\frac{L_{c}}{2}}\frac{n(x)u_{y}(x)}{a}dx$
(63)

and it was normalized by $\dot{m}_{0}=n_{0}maL_{c}^{2}/v_{m}$, where
$v_{m}=\sqrt{k_{B}T_{w}/m}$. The results are plotted in Fig. 11(a), where we
plot the normalized mass flow rate with respect to the Knudsen number for four
values of the confinement ratio $R\in\\{2,5,10,20\\}$. The tightness of wall
confinement significantly affects the mass flow rate in the Enskog equation.
When the channel width is large ($R>20$), the MFR at Knudsen numbers ($Kn$)
larger than 1 mostly matches that of the Boltzmann equation. However, for
$Kn<1$, the MFR in the Enskog equation is smaller than predicted by the
Boltzmann equation, and this difference becomes more pronounced as $Kn$
decreases. As the channel width decreases, the MFR becomes smaller, and the
deviation from the Boltzmann equation extends to larger Knudsen numbers. For
$R>5$, the MFR does not monotonically decrease with $Kn$ in the slip flow
regime ($Kn<0.1$). Instead, there exists a specific value of $Kn$ at which the
MFR locally reaches a maximum. When $R$ is less than or equal to 5, the
Knudsen minimum disappears, and the MFR only increases with the Knudsen
number. Following the study in Ref. [75], the Knudsen minimum disappearance
can be explained as follows. Under tight geometries, the combined contribution
of viscosity and density weakens, while the slip term remains constant as the
confinement ratio R decreases. As a result, the relative importance of the
slip term increases in this context. Consequently, the disappearance of the
Knudsen minimum under tight confinement can be attributed to the more
pronounced significance of fluid slippage at the wall.

In Fig. 11(b), the normalized mass flow rate $\dot{m}/\dot{m}_{0}$ is
evaluated for different quadrature orders $Q_{x}^{h}$ in the FDLB model with
respect to the Knudsen number. Specifically, four values of the quadrature
order $Q_{x}^{h}\in\\{8,16,32,64\\}$ were considered. This plot illustrates
that, as the Knudsen number $Kn$ increases, one has to employ a larger
quadrature order in order to approach the PM results. The results presented in
Fig. 11(a) are obtained using a varying quadrature order over the interval of
the Knudsen number. More specifically, a quadrature order of order
$Q_{x}^{h}=8$ is sufficient for $Kn<1$ and then we switched to a quadrature of
$Q_{x}^{h}=100$ to ensure the accuracy.

Figure 11:  Poiseuille flow (mass flow rate): (a) Normalised mass flow rate
$\dot{m}/\dot{m}_{0}$ at an external acceleration of $a_{y}=0.001$ with
respect to Knudsen number $Kn$, four values of the confinement ratio
$R=\\{2,5,10,20\\}$. (b) Normalized mass flow rate $\dot{m}/\dot{m}_{0}$ with
respect to the Knudsen number for four values of the quadrature order
$Q_{x}\in\\{8,16,32,64\\}$. One can observe that as the Knudsen number
increases one has to employ a higher order quadrature in order to recover the
mean flow rate, as is the case for dilute gases.

##  V Conclusions

In this work, a series of dense gas flows bounded by parallel plates were
simulated in order to validate the proposed finite-difference Lattice
Boltzmann model employing the simplified Enskog collision integral. In this
model, the Enskog collision integral is approximated using a Taylor expansion
and retaining the first-order gradients. We benchmarked the model in the
following setups: Fourier flow, Couette flow, and Poiseuille flow. The
simulation parameters range from a low reduced density value ($\eta_{0}=0.01$)
to a relatively high value ($\eta_{0}=0.2$) and include the confinement ratios
of $R\in\\{4,10\\}$.

The FDLB results obtained for the Couette flow, Fourier flow, and Poiseuille
flow were validated against the corresponding PM results. Reasonable agreement
was observed throughout the parameter range. More specifically, our kinetic
model adequately captures the effects of denseness, density inhomogeneity, as
well as nonequilibrium phenomena within the range of flow parameters
investigated.

In the case of the Fourier flow, we further examined the transversal heat flux
with respect to the Knudsen number $Kn$ where we observed a significant
deviation from the dilute gases results. Specifically, the heat flux is no
longer monotonic and moreover exhibits a local maximum, whose value for
$\eta_{0}=0.2$ is even larger than the ballistic limit. This is attributed to
the potential contributions to the heat flux, as detailed in Appendix B.

The Couette flow results are well captured by the FDLB model, even if the
deviations with respect to the PM results become larger as the fluid is more
confined (i.e. smaller $R$). Besides the density, velocity, and temperature
profiles, we have also presented the transversal and longitudinal heat flux
profiles. Since the simplified Enskog collision operator approximation
discards the higher-order contributions to the collisional momentum and energy
transfer, the FDLB results for the transversal and the longitudinal heat
fluxes agree with the PM results only at very low reduced density
($\eta_{0}=0.01$).

As in the case of the Fourier and Couette flow, also the Poiseuille flow
results are in good agreement with the PM results up to a mean reduced density
of $\eta_{0}=0.1$. As shown before [2, 23, 75], the Knudsen minimum in the
Poiseuille setup disappears for ultra-tight confinement ($R<5$). This effect
is well captured by the FDLB model.

In conclusion, the presented model demonstrates its capability in handling
moderately dense gases. Additionally, we examined the model’s performance in
dealing with flows characterized by sharp gradients in macroscopic quantities
arising in the gas-surface interaction. This was achieved with a much lower
computational cost compared to the particle method simulations. Moving
forward, our future plans include incorporating attractive forces between
molecules to address multiphase flows.

###### Acknowledgements.

This work was supported through a grant from the Ministry of Research,
Innovation and Digitization, CNCS - UEFISCDI, project number PN-
III-P1-1.1-PD-2021-0216, within PNCDI III. The authors thank dr. V.E. Ambrus
for useful discussions regarding the present study.

##  Appendix A Full-range Gauss Hermite quadratures: construction and
comparison to half-range GH quadrature results

Let’s examine integrals of $f$ and $f^{\mathrm{eq}}$ along the entire axis of
the 1D momentum space:

$\begin{pmatrix}M_{s}\\\
M_{s}^{(\rm{eq})}\rule{0.0pt}{19.91692pt}\end{pmatrix}=\int_{-\infty}^{\infty}dp\,\begin{pmatrix}f\\\
f^{\rm{eq}}\rule{0.0pt}{19.91692pt}\end{pmatrix}p^{s}.\qquad$ (64)

where $f^{eq}=ng$ and

$g\equiv g(x,p,t)=\frac{1}{\sqrt{2\pi
mT}}\exp\left[-\frac{(p-mu)^{2}}{2T}\right].$ (65)

The function $g$ can be expanded with respect to the full-range Hermite
polynomials ${H_{\ell}(p),\ell=0,1,\dots}$ as follows [45, 44]:

$g=\frac{\omega(\bar{p})}{p_{0}}\sum_{\ell=0}^{\infty}\frac{1}{\ell!}\mathcal{G}_{\ell}H_{\ell}(\bar{p}),\qquad\mathcal{G}_{\ell}=\sum_{s=0}^{\lfloor\ell/2\rfloor}\frac{\ell!}{2^{s}s!(\ell-2s)!}\left(\frac{mT}{p_{0}^{2}}-1\right)^{s}\left(\frac{mu}{p_{0}}\right)^{\ell-2s},$
(66)

where $\bar{p}\equiv p/p_{0}$ represents the particle momentum with respect to
an arbitrary momentum scale $p_{0}$, and $\lfloor\cdot\rfloor$ denotes the
floor function.

The full-range Hermite polynomials [44, 62, 63] satisfy the following
orthogonality relation with respect to the weight function $\omega(p)$:

$\int_{-\infty}^{\infty}dp\,\omega(p)H_{\ell}(\bar{p})H_{\ell^{\prime}}(\bar{p})=\ell!\,\delta_{\ell,\ell^{\prime}},\qquad\omega(p)=\frac{1}{\sqrt{2\pi}}e^{-\bar{p}^{2}/2}.$
(67)

The expansion coefficients $\mathcal{G}_{\ell}$ given in Eq. (66) are obtained
as follows:

$\mathcal{G}_{\ell}=\int_{-\infty}^{\infty}dp\,g\,H_{\ell}(p).$ (68)

By substituting Eq. (66) into Eq. (64), we obtain:

$M^{(\mathrm{eq})}_{s}=p_{0}^{s}\sum_{\ell=0}^{\infty}\frac{1}{\ell!}\,\mathcal{G}_{\ell}\int_{-\infty}^{\infty}dp\,\omega(p)\,H_{\ell}(p)\,p^{s}.$
(69)

For finite values of $s$ and $\ell$, the Gauss-Hermite quadrature can be
applied to compute the integral over $p$ across the entire momentum axis using
the following approach:

$\int_{-\infty}^{\infty}dp\,\omega(p)P_{s}(p)\simeq\sum_{k=1}^{Q}w_{k}^{H}P_{s}(p_{k}),$
(70)

where $P_{s}(p)$ is a polynomial of order $s$ in $p$, and the $Q$ quadrature
points $p_{k}$ ($k=1,2,\dots Q$) are the roots of the Hermite polynomial of
order $Q$, i.e., $H_{Q}(p_{k})=0$. The quadrature weights $w_{k}^{H}$ are
given by:

$w_{k}^{H}=\frac{Q!}{[H_{Q+1}(p_{k})]^{2}}.$ (71)

The equality in Eq. (70) is exact if $2Q>s$. In an LB simulation, $Q$ is fixed
at runtime. Thus, in order to ensure the exact recovery of
$M^{(\mathrm{eq})}_{s}$ (69), the sum over $\ell$ in Eq. (66) must be
truncated at a finite value $\ell=N$. Setting $Q>N$ ensures the exact recovery
of the first $N+1$ moments (i.e. $s=0,1,\dots N$) of $f^{eq}$, since the terms
of higher order in the expansion of $g$ are orthogonal to all polynomials
$P_{s}(p)$ of orders $0\leq s\leq N$, due to the orthogonality relation (67).
This allows us to obtain $M^{\mathrm{eq}}_{s}$ as follows:

$M^{(\mathrm{eq})}_{s}=\sum_{k=1}^{Q}f^{eq}_{k}\bar{p}_{k}^{s},\qquad
f^{eq}_{k}=ng_{k}^{H},\qquad
g_{k}^{H}=\frac{w_{k}^{H}p_{0}}{\omega(\bar{p}_{k})}g^{H,(N)}(\bar{p}_{k}),$
(72)

where $\bar{p}_{k}=p_{0}p_{k}$ represents the discrete momenta, and the
notation $g^{H,(N)}(p)$ indicates that the polynomial expansion (66) of $g(p)$
is truncated at order $\ell=N$ using the full-range Hermite polynomials. For
clarity, we provide the expression for $g^{H}_{k}$ [45, 44] below:

$g^{H,(N)}_{k}=w^{H}_{k}\sum_{\ell=0}^{N}H_{\ell}(\bar{p}_{k})\sum_{s=0}^{\lfloor\ell/2\rfloor}\frac{1}{2^{s}s!(\ell-2s)!}\left(\frac{mT}{p_{0}^{2}}-1\right)^{s}\left(\frac{mu}{p_{0}}\right)^{\ell-2s}.$
(73)

In the case of the Poiseuille flow, the momentum derivative
$\partial_{p_{y}}f$ can be written as:

$(\partial_{p_{y}}f)_{k}=\sum_{k,k^{\prime}}\mathcal{K}_{k,k^{\prime}}f_{k^{\prime}},$
(74)

where the kernel $\mathcal{K}_{k,k^{\prime}}$ has the following components
[76, 61]:

$\mathcal{K}_{k,k^{\prime}}=-\frac{w_{k}^{H}}{p_{0}}\sum_{\ell=0}^{Q-2}\frac{1}{\ell!}H_{\ell+1}(\bar{p}_{k})H_{\ell}(\bar{p}_{k^{\prime}}).$
(75)

In Fig.12 we plot the results obtained using the Full-range Hermite
polynomials against the half-range results for $R=10$ and $\Delta T=0.5$,
$u_{w}=1$ and $a_{y}=0.001$, for Fourier, Couette and Poiseuille flows,
respectively. The reduced density was set to $\eta_{0}=0.01$. We can easily
see that a significantly larger velocity set needs to be employed in the case
of Full-range polynomials in order to obtain similar results to those obtained
using a half-range quadrature of $Q^{h}_{x}=8$, with 16 velocities. In most
cases, a quadrature of $Q^{H}_{x}=200$ is needed in order to have relative
errors of less than $1\%$.

Figure 12:  Comparison of half and fll-range Hermite polynomials: (a) Fourier
flow (b) Couette flow and (c) Poiseuille flow at $R=10$ and $\Delta T=0.5$,
$u_{w}=1$ and $a_{y}=0.001$, while the reduced density was set to
$\eta_{0}=0.01$. These plots indicate a convergence of the results obtained
using full-range quadrature towards those obtained using half-range
quadrature.

##  Appendix B  Local maximum in heat flux for the Fourier flow

Further looking at the heat flux in the context of the Fourier flow, one can
evaluate the kinetic and potential contributions to the total heat flux. While
the kinetic component is evaluated as usual, the potential part is given by
[21]:

$q_{x}^{pot}=-m\frac{\sigma^{2}}{4}\int
d\bm{\xi}_{1}d\bm{\xi}_{2}d^{2}\bm{k}\int_{0}^{\sigma}d\alpha({\xi_{1}^{*}}^{2}-\xi_{1}^{2})\bm{k}\chi\left({\bm{x}}+\alpha\bm{k}-\frac{\sigma}{2}{\bm{k}}\right)\\\
f({\bm{x}}+\alpha\bm{k}-\sigma{\bm{k}},\bm{\xi}_{1})f({\bm{x}+\alpha\bm{k}},\bm{\xi}_{2})({\bm{p_{r}}}\cdot{\bm{k}})$
(76)

where $\alpha$ is a dummy variable.

Using the particle method described below one may assess the individual
contributions to the total heat flux. In Fig. 13 we plot the components of the
kinetic
$\widetilde{q_{x}^{kin}}=\int_{-\frac{L_{c}}{2}}^{\frac{L_{c}}{2}}q_{x}^{kin}dx$
and potential
$\widetilde{q_{x}^{pot}}=\int_{-\frac{L_{c}}{2}}^{\frac{L_{c}}{2}}q_{x}^{pot}dx$
heat flux contributions integrated over the channel length with respect to the
Knudsen number. The total heat flux $q_{x}=q_{x}^{kin}+q_{x}^{pot}$ is also
plotted in this figure, being the same curve corresponding to $\eta_{0}=0.2$
as in Fig. 4(b). One can easily observe the monotonic increase in the kinetic
component while the potential contribution has two local maxima in the
interval $Kn\in\\{1:10\\}$. The total heat flux has a value larger than the
free molecular limit due to the potential contribution to the heat flux.

Figure 13: Heat flux: The components of the heat flux, kinetic $q_{x}^{kin}$
and potential $q_{x}^{pot}$, with respect to the Knudsen number. The kinetic
component is monotonic with respect to the Knudsen number, while the potential
part has two local maxima between $Kn=1$ and 10, leading to a total heat flux
that exceeds the free molecular limit.

##  Appendix C  Particle method for the Enskog equation

The Enskog equation (Eq.(1)) is numerically solved using a particle method
that extends the original Direct Simulation Monte-Carlo (DSMC) method to
handle the nonlocal nature of the Enskog collision integral [18]. For a
comprehensive explanation of the numerical scheme and an analysis of its
computational complexity, please refer to Ref. [27]. Here, we provide a brief
overview of the scheme.

The DSMC framework used to solve the Boltzmann equation is maintained, with
modifications made to the collision algorithm to accommodate the nonlocal
structure of the Enskog collision operator. The distribution function is
represented by $N$ computational particles:

$f(\bm{x},\bm{p},t)=\frac{1}{m}\sum_{i=1}^{N}\delta{\left(\bm{x}-\bm{x}_{i}(t)\right)}\delta(\bm{p}-\bm{p}_{i}(t)),$
(77)

where $\bm{x}_{i}$ and $\bm{p}_{i}$ are the positions and momenta of the $i$th
particle at time $t$, respectively. The distribution function
$f(\bm{x},\bm{p},t)$ is updated using a fractional-step method that splits the
evolution operator into two sub-steps: free streaming and collision. In the
first stage, the distribution function is advanced from $t$ to $t+\Delta t$ by
neglecting particle collisions, i.e., by solving the equation:

$\frac{\partial f}{\partial
t}+\frac{\bm{p}}{m}\cdot\nabla_{\bm{x}}f+\bm{F}\cdot\bm{\nabla}_{\bm{p}}f=0,$
(78)

which leads to updating the positions of the computational particles as
follows:

$\displaystyle\bm{x}_{i}(t+\Delta t)$
$\displaystyle=\bm{x}_{i}(t)+\frac{\bm{p}_{i}}{m}\Delta
t+\frac{F}{m}\frac{(\Delta t)^{2}}{2},$ (79a)
$\displaystyle\bm{v}_{i}(t+\Delta t)$
$\displaystyle=\bm{v}_{i}(t)+\frac{F}{m}\Delta t.$ (79b)

resulting in the updated distribution function denoted
$\tilde{f}(\bm{x},\bm{p},t+\Delta t)$.

In the second stage, the short-range hard-sphere interactions are evaluated,
and the distribution function is updated according to:

$f(\bm{x},\bm{p},t+\Delta t)=\tilde{f}(\bm{x},\bm{p},t+\Delta
t)+J_{E}[\tilde{f}]\Delta t.$ (80)

During this stage, the $N$ particle positions $\bm{x}_{i}$ remain unchanged,
while their momenta $\bm{p}_{i}/m$ are modified based on stochastic rules that
essentially involve the Monte Carlo evaluation of the collision integral given
by Eq. (2) by selecting collision pairs accordingly.

The macroscopic quantities are obtained by averaging the microscopic states of
the particles over time and performing phase averaging by running
statistically independent simulations with identical macroscopic initial
conditions but different random seeds.

## References

  * Ferziger and Kaper [1972] J. Ferziger and H. Kaper,  _Mathematical Theory of Transport Processes in Gases._ (North-Holland Publishing Company, Amsterdam, London, 1972).
  * Wu  _et al._ [2016] L. Wu, H. Liu, J. M. Reese, and Y. Zhang, Non-equilibrium dynamics of dense gas under tight confinement, Journal of Fluid Mechanics 794, 252–266 (2016).
  * Sander  _et al._ [2017] R. Sander, Z. Pan, and L. D. Connell, Laboratory measurement of low permeability unconventional gas reservoir rocks: A review of experimental methods, Journal of Natural Gas Science and Engineering 37, 248 (2017).
  * Petersen and Hanson [2001] E. L. Petersen and R. K. Hanson, Nonideal effects behind reflected shock waves in a high-pressure shock tube, Shock Waves 10, 405 (2001).
  * Holt  _et al._ [2006] J. K. Holt, H. G. Park, Y. Wang, M. Stadermann, A. B. Artyukhin, C. P. Grigoropoulos, A. Noy, and O. Bakajin, Fast mass transport through sub-2-nanometer carbon nanotubes, Science 312, 1034 (2006).
  * Brenner  _et al._ [2002] M. P. Brenner, S. Hilgenfeldt, and D. Lohse, Single-bubble sonoluminescence, Rev. Mod. Phys. 74, 425 (2002).
  * Enskog [1922] D. Enskog,  _Kinetische Theorie der Wärmeleitung: Reibung und Selbst-diffusion in Gewissen verdichteten gasen und flüssigkeiten._ (Almqvist & Wiksells, 1922).
  * Chapman and Cowling [1970] S. Chapman and T. G. Cowling,  _The Mathematical Theory of Non-uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases._ (Cambridge University Press, 1970).
  * Kremer [2010] G. M. Kremer,  _An introduction to the Boltzmann equation and transport processes in gases_ (Springer-Verlag, Berlin Heidelberg, 2010).
  * Dorfman  _et al._ [2021] J. R. Dorfman, H. van Beijeren, and T. Kirkpatrick,  _Contemporary Kinetic Theory of Matter_ (Cambridge University Press, 2021).
  * Frezzotti and Sgarra [1993] A. Frezzotti and C. Sgarra, Numerical analysis of a shock-wave solution of the Enskog equation obtained via a Monte Carlo method, J. Stat. Phys. 73, 193 (1993).
  * Wu  _et al._ [2015] L. Wu, Y. Zhang, and J. M. Reese, Fast spectral solution of the generalized Enskog equation for dense gases, Journal of Computational Physics 303, 66 (2015).
  * Sadr and Gorji [2017] M. Sadr and M. H. Gorji, A continuous stochastic model for non-equilibrium dense gases, Physics of Fluids 29, 122007 (2017).
  * Sadr and Gorji [2019] M. Sadr and M. Gorji, Treatment of long-range interactions arising in the Enskog–Vlasov description of dense fluids, J. Comput. Phys. 378, 129 (2019).
  * Bird [1976] G. A. Bird,  _Molecular Gas Dynamics_ (Oxford Univ. Press, Oxford, England, UK, 1976).
  * Alexander  _et al._ [1995] F. J. Alexander, A. L. Garcia, and B. J. Alder, A consistent Boltzmann algorithm, Phys. Rev. Lett. 74, 5212 (1995).
  * Montanero and Santos [1996] J. M. Montanero and A. Santos, Monte Carlo simulation method for the Enskog equation, Phys. Rev. E 54, 438 (1996).
  * Frezzotti [1997] A. Frezzotti, A particle scheme for the numerical solution of the Enskog equation, Phys. Fluids 9, 1329 (1997).
  * Davis [1987] H. T. Davis, Kinetic theory of inhomogeneous fluid: Tracer diffusion, J. Chem. Phys. 86, 1474 (1987).
  * Din and Michaelides [1997] X.-D. Din and E. E. Michaelides, Kinetic theory and molecular dynamics simulations of microscopic flows, Physics of Fluids 9, 3915 (1997).
  * Frezzotti [1999] A. Frezzotti, Monte Carlo simulation of the heat flow in a dense hard sphere gas, European Journal of Mechanics - B/Fluids 18, 103 (1999).
  * Nedea  _et al._ [2006] S. Nedea, A. Frijns, A. van Steenhoven, A. Jansen, A. Markvoort, and P. Hilbers, Density distribution for a dense hard-sphere gas in micro/nano-channels: Analytical and simulation results, Journal of Computational Physics 219, 532 (2006).
  * Sheng  _et al._ [2020] Q. Sheng, L. Gibelli, J. Li, M. K. Borg, and Y. Zhang, Dense gas flow simulations in ultra-tight confinement, Phys. Fluids 32, 092003 (2020).
  * Su  _et al._ [2023] W. Su, L. Gibelli, J. Li, M. K. Borg, and Y. Zhang, Kinetic modeling of nonequilibrium flow of hard-sphere dense gases, Phys. Rev. Fluids 8, 013401 (2023).
  * Frezzotti  _et al._ [2005] A. Frezzotti, L. Gibelli, and S. Lorenzani, Mean field kinetic theory description of evaporation of a fluid into vacuum, Phys. Fluids 17, 012102 (2005).
  * Kon  _et al._ [2014] M. Kon, K. Kobayashi, and M. Watanabe, Method of determining kinetic boundary conditions in net evaporation/condensation, Phys. Fluids 26, 072003 (2014).
  * Frezzotti  _et al._ [2019] A. Frezzotti, P. Barbante, and L. Gibelli, Direct simulation Monte Carlo applications to liquid-vapor flows, Phys. Fluids 31, 062103 (2019).
  * Busuioc  _et al._ [2020] S. Busuioc, L. Gibelli, D. A. Lockerby, and J. E. Sprittles, Velocity distribution function of spontaneously evaporating atoms, Phys. Rev. Fluids 5, 103401 (2020).
  * Bruno and Frezzotti [2019] D. Bruno and A. Frezzotti, Dense gas effects in the Rayleigh-Brillouin scattering spectra of SF6, Chem. Phys. Lett. 731, 136595 (2019).
  * Busuioc and Gibelli [2020] S. Busuioc and L. Gibelli, Mean-field kinetic theory approach to Langmuir evaporation of polyatomic liquids, Physics of Fluids 32, 093314 (2020).
  * Kobayashi  _et al._ [2017] K. Kobayashi, K. Sasaki, M. Kon, H. Fujii, and M. Watanabe, Kinetic boundary conditions for vapor–gas binary mixture, Microfluid. Nanofluid. 21, 53 (2017).
  * Barbante  _et al._ [2015] P. Barbante, A. Frezzotti, and L. Gibelli, A kinetic theory description of liquid menisci at the microscale, Kinet. Relat. Mod. 8, 235 (2015).
  * Busuioc  _et al._ [2023] S. Busuioc, A. Frezzotti, and L. Gibelli, A weighted particle scheme for Enskog-Vlasov equation to simulate spherical nano-droplets/bubbles, Journal of Computational Physics 475, 111873 (2023).
  * Luo [1998] L.-S. Luo, Unified theory of lattice Boltzmann models for nonideal gases, Phys. Rev. Lett. 81, 1618 (1998).
  * Luo [2000] L.-S. Luo, Theory of the lattice Boltzmann method: Lattice Boltzmann models for nonideal gases, Phys. Rev. E 62, 4982 (2000).
  * Melchionna and Marconi [2007] S. Melchionna and U. M. B. Marconi, Lattice Boltzmann method for inhomogeneous fluids, Europhysics Letters 81, 34001 (2007).
  * He and Doolen [2002] X. He and G. Doolen, Thermodynamic foundations of kinetic theory and lattice Boltzmann models for multiphase flows., J. Stat. Phys. 107, 309 (2002).
  * Wang  _et al._ [2020] P. Wang, L. Wu, M. T. Ho, J. Li, Z.-H. Li, and Y. Zhang, The kinetic Shakhov–Enskog model for non-equilibrium flow of dense gases, Journal of Fluid Mechanics 883, A48 (2020).
  * Chen  _et al._ [2023] T. Chen, L. Wu, L.-P. Wang, and S. Chen, Rarefaction effects in head-on collision of two near-critical droplets, International Journal of Multiphase Flow 164, 104451 (2023).
  * Huang  _et al._ [2021] R. Huang, H. Wu, and N. A. Adams, Mesoscopic lattice Boltzmann modeling of the liquid-vapor phase transition, Phys. Rev. Lett. 126, 244501 (2021).
  * Zhang  _et al._ [2020] Y.-D. Zhang, A.-G. Xu, J.-J. Qiu, H.-T. Wei, and Z.-H. Wei, Kinetic modeling of multiphase flow based on simplified Enskog equation., Front. Phys. 15, 62503 (2020).
  * Gan  _et al._ [2022] Y. Gan, A. Xu, H. Lai, W. Li, G. Sun, and S. Succi, Discrete Boltzmann multi-scale modelling of non-equilibrium multiphase flows, Journal of Fluid Mechanics 951, A8 (2022).
  * Busuioc [2023] S. Busuioc, Quadrature-based lattice Boltzmann model for non-equilibrium dense gas flows, Physics of Fluids 35, 016112 (2023).
  * Ambruş and Sofonea [2016a] V. Ambruş and V. Sofonea, Lattice Boltzmann models based on half-range Gauss-Hermite quadratures, J. Comput. Phys. 316, 760 (2016a).
  * Ambruş and Sofonea [2016b] V. Ambruş and V. Sofonea, Application of mixed quadrature lattice Boltzmann models for the simulation of Poiseuille flow at non-negligible values of the Knudsen number, J. Comput. Science 17, 403 (2016b).
  * Ambruş and Sofonea [2018] V. E. Ambruş and V. Sofonea, Half-range lattice Boltzmann models for the simulation of Couette flow using the Shakhov collision term, Phys. Rev. E 98, 063311 (2018).
  * Ambruș and Sofonea [2019] V. E. Ambruș and V. Sofonea, Quadrature-based lattice Boltzmann models for rarefied gas flow, in  _Flowing Matter_ , edited by F. Toschi and M. Sega (Springer International Publishing, Cham, 2019) pp. 271–299.
  * Carnahan and Starling [1969] N. F. Carnahan and K. E. Starling, Equation of state for nonattracting rigid spheres, J. Chem. Phys. 51, 635 (1969).
  * Van Beijeren and Ernst [1973] H. Van Beijeren and M. Ernst, The modified Enskog equation, Physica 68, 437 (1973).
  * Fischer and Methfessel [1980] J. Fischer and M. Methfessel, Born-Green-Yvon approach to the local densities of a fluid at interfaces, Phys. Rev. A 22, 2836 (1980).
  * Shakhov [1968a] E. Shakhov, Approximate kinetic equations in rarefied gas theory, Fluid Dynamics 3, 95 (1968a).
  * Shakhov [1968b] E. Shakhov, Approximate kinetic equations in rarefied gas theory, Fluid Dynamics 3, 112 – 115 (1968b).
  * Graur and Polikarpov [2009] I. Graur and A. Polikarpov, Comparison of different kinetic models for the heat transfer problem, Heat Mass Transfer 46, 237 (2009).
  * Ambruş  _et al._ [2020] V. E. Ambruş, F. Sharipov, and V. Sofonea, Comparison of the Shakhov and ellipsoidal models for the Boltzmann equation and DSMC for ab initio-based particle interactions, Computers & Fluids 211, 104637 (2020).
  * Sharipov [2002] F. Sharipov, Application of the Cercignani–Lampis scattering kernel to calculations of rarefied gas flows. i. plane flow between two parallel plates, European Journal of Mechanics - B/Fluids 21, 113 (2002).
  * Sharipov [2003] F. Sharipov, Application of the Cercignani–Lampis scattering kernel to calculations of rarefied gas flows. ii. slip and jump coefficients, European Journal of Mechanics - B/Fluids 22, 133 (2003).
  * Zhang  _et al._ [2019] Y. Zhang, A. Xu, G. Zhang, Z. Chen, and P. Wang, Discrete Boltzmann method for non-equilibrium flows: Based on Shakhov model, Computer Physics Communications 238, 50 (2019).
  * Todorova  _et al._ [2020] B. N. Todorova, C. White, and R. Steijl, Modeling of nitrogen and oxygen gas mixture with a novel diatomic kinetic model, AIP Advances 10, 095218 (2020).
  * Li and Zhang [2004] Z.-H. Li and H.-X. Zhang, Study on gas kinetic unified algorithm for flows from rarefied transition to continuum, J. of Comput. Phys. 193, 708 (2004).
  * Meng  _et al._ [2013] J. Meng, L. Wu, J. M. Reese, and Y. Zhang, Assessment of the ellipsoidal-statistical Bhatnagar–Gross–Krook model for force-driven Poiseuille flows, Journal of Computational Physics 251, 383 (2013).
  * Busuioc and Ambruş [2019] S. Busuioc and V. E. Ambruş, Lattice Boltzmann models based on the vielbein formalism for the simulation of flows in curvilinear geometries, Phys. Rev. E 99, 033304 (2019).
  * Hildebrand [1987] F. B. Hildebrand,  _Introduction to Numerical Analysis, 2nd edition_ (Dover Publications, 1987).
  * Shan  _et al._ [2006] X. Shan, X.-F. Yuan, and H. Chen, Kinetic theory representation of hydrodynamics: a way beyond the Navier–Stokes equation, Journal of Fluid Mechanics 550, 413 (2006).
  * [64] B. Shizgal,  _Spectral methods in chemistry and physics_.
  * Ansumali and V. Karlin [2002] S. Ansumali and I. V. Karlin, Kinetic boundary conditions in the lattice Boltzmann method, Phys. Rev. E 66, 026311 (2002).
  * Shu and Osher [1988] C.-W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, Journal of Computational Physics 77, 439 (1988).
  * Gottlieb and Shu [1998] S. Gottlieb and C.-W. Shu, Total variation diminishing Runge-Kutta schemes, Mathematics of computation 67, 73 (1998).
  * Rezzolla and Zanotti [2013] L. Rezzolla and O. Zanotti,  _Relativistic Hydrodynamics_ (Oxford University Press, 2013).
  * Gan  _et al._ [2011] Y. Gan, A. Xu, G. Zhang, and Y. Li, Lattice Boltzmann study on Kelvin-Helmholtz instability: Roles of velocity and density gradients, Phys. Rev. E 83, 056704 (2011).
  * Jiang and Shu [1996] G.-S. Jiang and C.-W. Shu, Efficient implementation of Weighted ENO schemes, Journal of Computational Physics 126, 202 (1996).
  * Ambruş and Blaga [2018] V. E. Ambruş and R. Blaga, High-order quadrature-based lattice Boltzmann models for the flow of ultrarelativistic rarefied gases, Phys. Rev. C 98, 035201 (2018).
  * Fornberg [1988] B. Fornberg, Generation of finite difference formulas on arbitrarily spaced grids, Mathematics of computation 51, 699 (1988).
  * Bassanini  _et al._ [1968] P. Bassanini, C. Cercignani, and C. Pagani, Influence of the accommodation coefficient on the heat transfer in a rarefied gas, International Journal of Heat and Mass Transfer 11, 1359 (1968).
  * Sharipov and Seleznev [1998] F. Sharipov and V. Seleznev, Data on Internal Rarefied Gas Flows, Journal of Physical and Chemical Reference Data 27, 657 (1998).
  * Corral-Casas  _et al._ [2022] C. Corral-Casas, J. Li, M. K. Borg, and L. Gibelli, Knudsen minimum disappearance in molecular-confined flows, Journal of Fluid Mechanics 945, A28 (2022).
  * Ambrus  _et al._ [2017] V. E. Ambrus, V. Sofonea, R. Fournier, and S. Blanco, Implementation of the force term in half-range lattice Boltzmann models, (2017), arXiv:1708.03249 [physics.flu-dyn] .