See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/342387574

Comparison of the Shakhov and ellipsoidal models for the Boltzmann equation
and DSMC for ab initio-based particle interactions

Article  in  Computers & Fluids · June 2020

DOI: 10.1016/j.compfluid.2020.104637

CITATIONS READS

8 153

3 authors:

Victor E. Ambruş Felix Sharipov

Goethe-Universität Frankfurt am Main Universidade Federal do Paraná
74 PUBLICATIONS   522 CITATIONS    200 PUBLICATIONS   5,606 CITATIONS   

SEE PROFILE SEE PROFILE

Victor Sofonea
Center for Fundamental and Advanced Technical Research, Romanian Academy, Tim…
76 PUBLICATIONS   1,407 CITATIONS   

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

primary pressure standard View project

Large vacuum system modelling View project

All content following this page was uploaded by Victor Sofonea on 30 September 2020.

The user has requested enhancement of the downloaded file.

 Comparison of the Shakhov and ellipsoidal models for
the Boltzmann equation and DSMC for ab initio-based
arXiv:1909.11552v2 [physics.comp-ph] 18 Jul 2020

particle interactions

Victor E. Ambrus, a,b , Felix Sharipovc , Victor Sofonead

a Department of Physics, West University of Timis, oara,
Bd. Vasile Pârvan 4, Timis, oara 300223, Romania†
b Department of Mathematics and Statistics, Old Dominion University,

Norfolk, VA 23529, USA

c Departamento de Fı́sica, Universidade Federal do Paraná, Curitiba, 81531-980 Brazil
d Center for Fundamental and Advanced Technical Research, Romanian Academy,

Bd. Mihai Viteazu 24, Timis, oara 300223, Romania

Abstract

In this paper, we consider the capabilities of the Boltzmann equation with the
Shakhov and ellipsoidal models for the collision term to capture the characteris-
tics of rarefied gas flows. The benchmark is performed by comparing the results
obtained using these kinetic model equations with direct simulation Monte Carlo
(DSMC) results for particles interacting via ab initio potentials. The analysis
is restricted to channel flows between parallel plates and we consider three flow
problems, namely: the heat transfer between stationary plates, the Couette flow
and the heat transfer under shear. The simulations are performed in the non-
linear regime for the 3 He, 4 He, and Ne gases. The reference temperature ranges
between 1 K and 3000 K for 3 He and 4 He and between 20 K and 5000 K for
Ne. While good agreement is seen up to the transition regime for the direct
phenomena (shear stress, heat flux driven by temperature gradient), the rela-
tive errors in the cross phenomena (heat flux perpendicular to the temperature
gradient) exceed 10% even in the slip-flow regime. The kinetic model equations
are solved using the finite difference lattice Boltzmann algorithm based on half-

Email addresses: [email protected] (Victor E. Ambrus, ),

[email protected] (Felix Sharipov), [email protected],
[email protected] (Victor Sofonea)
† Permanent address.

Preprint submitted to Computers & Fluids July 21, 2020

 range Gauss-Hermite quadratures with the third order upwind method used for
the implementation of the advection.
Keywords: Ab initio, DSMC, Ellipsoidal model, Shakhov model, Half-range
Gauss-Hermite quadrature

1. Introduction

Finding accurate solutions of the kinetic equations governing rarefied gas

flows is a challenging task because of their complexity [1, 2]. In the case of
channel flows, it has been shown under quite general assumptions that the ve-
5 locity field in the vicinity of solid boundaries is non-analytic, its normal deriva-
tive presenting a logarithmic singularity with respect to the distance to the wall
[3]. Understanding the main properties of such flows is crucial when devising
micro/nano-electro-mechanical systems (MEMS/NEMS) [4].
Since the kinetic equation is difficult to solve analytically, numerical meth-
10 ods remain the primary tools available for its investigation. It has been es-
tablished in the research community that the direct simulation Monte Carlo
(DSMC) method [5] can provide solutions to realistic systems in a wide range
of flow regimes. The main ingredient controlling the relevance of the DSMC
formulation lies in specifying the interparticle interactions. Recently, ab initio
15 potentials have been implemented into the DSMC method [6, 7, 8, 9, 10]. A
quantum consideration of interatomic collisions [11, 12] allowed to extend an
application of ab initio potentials to low temperatures. To reduce the com-
putational effort, lookup tables for the deflection angle of binary collisions of
helium-3 (3 He), helium-4 (4 He), and neon (Ne) atoms have been calculated and
20 reported in the Supplementary material to Ref. [12]. The lookup tables can be
used for any flow of these gases over a wide range of temperature. Due to the
stochastic nature of the DSMC method, its results often exhibit steady-state
fluctuations, which are especially significant in the slip-flow regime and at small
Mach numbers. Filtering out these fluctuations is a computationally demanding
25 part of the algorithm, making this method computationally convenient only in

2
 the transition and free molecular flow regimes.
Another approach for the description of rarefied gas flows starts from the
Boltzmann equation, where the collision integral takes into account the details of
the interparticle interactions. While recent years have seen significant progress
30 in the development of numerical methods for evaluating the Boltzmann collision
integral [13, 14, 15, 16, 17], this operation still remains the most expensive part
of the solver, making the application of such methods for complex systems
computationally prohibitive.
As argued in the early ’50s, the features of the collision integral can be
35 preserved, at least for small Knudsen numbers and mildly non-linear systems,
by replacing the collision term through a relaxation time approach. The BGK
model, introduced by Bhatnagar, Gross and Krook [18], employed a single re-
laxation time τ to control the departure of the Boltzmann distribution function
f from local thermal equilibrium. This parameter could be used to match re-
40 alistic flows by ensuring the correct recovery of the dynamic viscosity µ in the
hydrodynamic regime, however it could not allow the heat conductivity κ to be
controlled independently. This difficulty was later alleviated through two exten-
sions, known as the ellipsoidal-BGK (ES) and Shakhov (S) models, proposed in
the late ’60s by Holway [19] and Shakhov [20, 21], respectively. The accuracy
45 of these models has been tested by considering the comparison to experimental
[22, 23, 24] or DSMC [25, 26, 27] results. In the following, we refer to these two
models (the ES and S models) as the model equations.
Various methods have been developed over the years to solve the model
equations and their variations. Amongst these, we mention the discrete velocity
50 method (DVM) [2, 28, 29, 30], the discrete unified gas kinetic scheme (DUGKS)
[31, 32, 33], the discrete Boltzmann method (DBM) [34, 35, 36] (generally re-
stricted to the Navier-Stokes regime due to the small velocity set size) and the
lattice Boltzmann (LB) method [37, 38, 39, 40] with its finite difference (FDLB)
version [25, 41, 42, 43, 44].
55 In the LB approach, the kinetic equation is employed to obtain an accurate
account of the evolution of the macroscopic moments of f [45, 46, 47, 48]. Less

3
 attention is directed to the distribution f itself. This allows the momentum
space to be sampled in a manner optimized for the recovery of the moments of
f [49]. Since the moments are defined as integrals of f , the momentum space
60 discretization can be viewed as a quadrature method [42]. Our implementation
is based on the idea of Gauss quadratures [50, 51], which provide a prescription
of choosing optimal quadrature points for the recovery of polynomial integrals,
given a certain domain and integration weight.
In this paper, we consider the systematic comparison between the numerical
65 solutions of the Boltzmann equation with the S and ES models for the collision
term, obtained using the FDLB algorithm, and the numerical results obtained
using DSMC. The comparison is made in the frame of channel flows between
parallel plates, where the fluid is assumed to be homogeneous with respect to
the directions parallel to the plates. Specifically, we address three flow prob-
70 lems. The first one is the heat transfer between stationary plates at differing
temperatures. The second is the Couette flow between parallel plates at equal
temperatures. The third problem refers to the heat transfer between plates at
differing temperatures undergoing parallel motion. In future studies, it may be
interesting to perform the comparison in more complex configurations, such as
75 the thermal transpiration through a long channel attached to two vessels with
different temperatures considered in Refs. [52] and [53], or the pressure-driven
flow through a long rectangular channel setup considered in Ref. [54]; as well as
in the highly nonlinear context of shock wave structure [8]. Since the focus in
this paper is on introducing kinetic models for the simulation of various kinds of
80 gas particles interacting via ab initio interparticle potentials, the present study
is restricted only to the channel flows mentioned above.
In channel flows, it is known that the particle-wall interaction induces a dis-
continuity in the distribution function [3, 55]. This discontinuity is responsible
for microfluidics effects, such as the development of a slip velocity and temper-
85 ature jump near the walls. Another important consequence of the discontinuity
of f is that the velocity profile becomes non-analytic in the vicinity of the wall,
where its derivative diverges logarithmically with respect to the distance to the

4
 wall [3, 56, 57].
As highlighted already in the late ’50s by Gross and his collaborators [55,
90 58, 59, 60], taking into account the discontinuity of the distribution function by
considering separately its moments with respect to the vectors pointing towards
and away from the wall (px > 0 and px < 0, respectively) can give a dramatic in-
crease in the accuracy of the Knudsen layer representation, compared to the full
momentum space projection approach. Recent works have focused on employ-
95 ing half-range quadratures [61, 62, 63] for the (semi-)analytical analysis of the
solutions of the (linearised or non-linear) Boltzmann equation in the relaxation
time approximation,
An important step in employing the idea of treating separately the dis-
tribution function for incoming and outgoing particles with respect to solid
100 walls in the numerical simulation of rarefied gas flows was taken in the ’60s
by Huang and Giddens [64], who computed the quadrature points and weights
for the one-dimensional half-range Gauss-Hermite quadrature with the weight
2
function ω(x) = e−x , up to 8th order. The extension of the procedure to
higher orders through a recurrence relation was discussed by Ball in Ref. [65]
105 and the algorithm was adapted in Ref. [44] to the case of the weight function
2 √
ω(x) = e−x /2 / 2π. A half-range (or modified) Gauss-Hermite quadrature was
used in the early 2000’s by Li and his collaborators [66, 67] for kinetic theory
simulations in the context of unbounded flows. Recently, the half-range Gauss-
Hermite quadrature was shown to offer significantly more accurate solutions of
110 the kinetic model equations than the full-range Gauss-Hermite quadrature with
the same number of quadrature points for the moderate and highly rarefied
regimes [43, 44, 68]. As a side note, similarly accurate results can be obtained
when the Gauss-Laguerre quadrature is used on the semi-axis, instead of the
Gauss-Hermite quadrature [69, 70].
115 In order to take advantage of the geometry of the channel flows considered in
this paper, we solve the kinetic model equations by employing the mixed quadra-
tures concept, according to which the quadrature is controlled separately on each
axis [44, 62]. This approach allows the half-range Gauss-Hermite quadrature to

5
 be employed on the x axis, which is perpendicular to the channel walls. On the
120 axes parallel to the walls, the full-range Gauss-Hermite quadrature can be em-
ployed. Details regarding Gauss quadratures can be found in various textbooks,
of which we remind Refs. [50, 51].
Furthermore, in the channel flows considered in this paper, the dynamics is
non-trivial only along d < D degrees of freedom (DOFs), where D = 3 is the
125 number of DOFs for an ideal monatomic gas. In particular, we consider d = 1
when the walls are stationary and d = 2 when the plates are in motion. We then
introduce two reduced distributions, φ and χ, which are obtained by integrating
the distribution function f multiplied by 1 and [p2d+1 + . . . p2D ]/m with respect
to dpd+1 · · · dpD [24]. Thus, φ can be seen to describe the mass and momen-
130 tum evolution and χ contributes to the energy evolution [27, 71]. When d = 2,
we employ the mixed quadrature paradigm [44, 61] and discretize the momen-
tum along the direction parallel to the wall using the full range Gauss-Hermite
quadrature. Furthemore, the homogeneity of the fluid along these directions
allows the system to be exactly described (i.e., without introducing any errors)
135 using a relatively low order quadrature [27, 44]. In this paper, we introduce
a novel expansion of the Shakhov and ellipsoidal collision terms with respect
to the full-range (standard) Hermite polynomials which allows the quadrature
orders along the y axis (which is parallel to the walls) to be set to Qφy = 4
and Qχy = 2 for the φ and χ distributions, respectively. The resulting expan-
140 sion coefficients remain dependent on the momentum component px which is
perpendicular to the boundary. While replacing the distributions with their
truncated polynomial expansions is inherited from the standard LB algorithm
[49], the expansion coefficients are evaluated directly (without resorting to poly-
nomial expansions), which is closer to the standard DVM practice [28]. In this
145 sense, our approach is a hybrid FDLB-DVM method (we refer to it as the hybrid
method), combining the advantages of both LB and DVM. For brevity, we use
the notation FDLB to refer to the numerical scheme that we employ to solve
the kinetic model equations. For completeness, the subsequent projection with
respect to px onto the half-range Hermite polynomials is also discussed (we refer

6
150 to this latter approach as the projection method). The latter approach is more
efficient than the hybrid approach in the hydrodynamic (small Kn) regime.
For the analysis presented in this paper, only the stationary state is of in-
terest. Since the transient solution is not important, iterative schemes can be
employed to solve the kinetic model equation, as described, e.g., in Refs. [72, 73,
155 74, 75]. However, since the computations in the one-dimensional settings that we
consider in this paper are not very demanding, we compute the stationary solu-
tion using explicit time marching, implemented using the third order total vari-
ation diminishing Runge-Kutta (RK-3) method introduced in Refs. [76, 77, 78].
For the advection operator, we introduce a third order upwind scheme which
160 preserves the order of accuracy in the presence of diffuse reflecting boundaries
which extends the one considered in Ref. [79] for the linearised Boltzmann-BGK
equation. We further increase the resolution inside the Knudsen layer by em-
ploying a grid stretching procedure [80, 81, 82]. The upwind method is known
to introduce numerical dissipation [83]. The numerical errors due to this spuri-
165 ous dissipation can be controlled by refining the grid. This becomes especially
important in the inviscid regime, where the numerical viscosity can dominate
over the physical one [83, 84]. In the slip-flow and transition regimes, we find
that 2S = 64 points per channel width are sufficient to obtain results which
have errors less than 0.1% (for more details, see Sec. 4.1), which is acceptable
170 from a computational cost point of view, e.g., compared to 100 [24], 101 [85],
500 [74] or 5000 [86] grid points employed in previous studies.
For simplicity, in this paper we only consider the Maxwell diffuse reflection
model with complete accommodation at the bounding walls. The methodology
can easily be extended to the case of more complex boundary conditions, such
175 as the diffuse-specular [29] and the Cercignani-Lampis [87] boundary models.
This paper is organised as follows. The kinetic models and the connection
to the DSMC simulations via the transport coefficients is presented in Sec. 2.
The FDLB algorithm is summarized in Sec. 3 and the simulation methodology
employed in the frame of the FDLB and DSMC approaches is summarized in
180 Sec. 4. Appendix A discusses the application of the FDLB method to the

7
 hydrodynamic regime. Sections 5, 6 and 7 present the numerical results for the
heat transfer between stationary plates, the Couette flow and the heat transfer
between moving plates problems, respectively. Section 8 concludes this paper.

2. Kinetic models and connection to DSMC

185 Subsection 2.1 introduces briefly the Shakhov and ellipsoidal-BGK models.
Subsection 2.2 introduces the reduced distribution functions employed in the
context of the channel flows discussed in this paper. The implementation of the
transport coefficients using the numerical data obtained from ab initio potentials
at the level of the model equations is discussed in Subsec. 2.3. Finally, our non-
190 dimensionalization conventions are summarized in Subsec. 2.4.

2.1. Model equations in the relaxation time approximation

In this paper, we focus on the study of channel flows between parallel plates.
The coordinate system is chosen such that the x
e axis is perpendicular to the
walls. The discussion in this section is presented at the level of dimensional
quantities, which are denoted explicitly via an overhead tilde. The origin of
the coordinate system is taken to be on the channel centerline, such that the
left and right walls are located at x e and x
e = −L/2 e
e = L/2, respectively. The
flow is studied in the Galilean frame where the left and right plates move with
velocities −e
uw and u
ew , respectively (e
uw = 0 for the heat transfer problem
between stationary plates discussed in Sec. 5). The temperatures of the left and
right plates are set to Teleft = Teref − ∆T
g /2 and Teright = Teref + ∆T
g /2, respectively
g = 0 for the Couette flow problem discussed in Sec. 6). In this case, the
(∆T
Boltzmann equation in the relaxation time approximation for the collision term
can be written as follows:

∂ fe pex ∂ fe 1
+ = − (fe − fe∗ ), (1)
∂t e m
e ∂ex τ
e∗

where fe is the particle distribution function, pex is the particle momentum along
the direction perpendicular to the walls, m
e is the particle mass and τe∗ is the

8
relaxation time. The collision term governs the relaxation of fe towards the local
equilibrium distribution function fe∗ . The star subscript in Eq. (1) distinguishes
between the two models that we consider in this paper, namely the Shakhov
model (∗ = S) and the ellipsoidal-BGK (∗ = ES) model. We consider in this
paper only monatomic ideal gases, for which fe∗ reduces at global thermodynamic
equilibrium to the Maxwell-Boltzmann distribution function feMB :

feMB (e e , Te) =e
n, u neg(e ex , Te)e
px , u g (e ey , Te)e
py , u g (e ez , Te),
pz , u
 
1 (ep − mee u)2
g(e
e e, Te) = q
p, u exp − . (2)
e B Te 2meK e B Te
2π m
eK

e is the particle number density, Te is the temperature and u

Here n eα (α ∈
{x, y, z}) are the components of the macroscopic velocity. These quantities
are obtained as moments of fe and fe∗ via the following relations:
     
n
e 1 1
  Z   Z  
  3   3  
 ρeu e  = d pe  pe  fe = d pe  pe  fe∗ , (3)
     
3 e e
n
e K B T e2 /2m
ξ e e2 /2m
ξ e
2

where ξe = pe − me
e u is the peculiar momentum. The last equality above is
a statement that the model equations preserve the collision invariants, ψ ∈
{1, p, p2 /2m}, of the original Boltzmann collision operator.
In the case of the Shakhov (S) model, the local equilibrium can be written
as [20, 21, 25, 27, 88]:
!
1 − Pr ξe2
feS = feMB (1 + S), S= e
− 1 qe · ξ, (4)
n
eKe 2 Te2 5m
eKe B Te
B

where the heat flux qe is obtained via

Z
ξe2 ξe
qe = d3 pe fe . (5)
2m e m
e

In the S model, the dynamic viscosity and the heat conductivity are controlled
by the relaxation time τeS and the Prandtl number Pr through

e
cp µ
eS e B τeS Pe
5K
eS = τeS Pe ,
µ κ
eS = = , (6)
Pr 2mPr
e

9
195 where e e B /2m
cp = 5 K e is the specific heat at constant pressure for an ideal
monatomic gas.
In the ellipsoidal-BGK (ES) model, the equilibrium distribution feES can be
written as [19, 26, 29, 36]:
!
B−1 e e
n
e αβ ξα ξβ
feES = √ exp − , (7)
(2π m
eKe B Te)3/2 detB 2m
eKe B Te

where B is an invertible 3 × 3 matrix (1 ≤ α, β ≤ D = 3) having the following

components: " #
1 Teαβ
Bαβ = δαβ − (1 − Pr) . (8)
Pr Pe

In the above, Pe = n
eKe B Te is the ideal gas pressure, while the Cartesian com-

ponents Teαβ of the pressure tensor are obtained as second order moments of fe:

Z
ξeα ξeβ
Teαβ = d3 pe fe . (9)
me
In the ES model, the transport coefficients are retrieved through:

e
cp µ
eES e B τeES Pe
5K
eES = τeES Pr Pe ,
µ κES =
e = . (10)
Pr 2m e

Eq. (1) is supplemented by boundary conditions. In this paper, we restrict

the analysis to the case of diffuse reflection with complete accommodation, such
that the distribution of the particles emerging from the wall back into the fluid
is described by the Maxwell-Boltzmann distribution [29]:

fe(−L/2,
e t) =feMB (e
pex > 0, e e w , Teleft),
nleft , −u

fe(L/2,
e t) =feMB (e
pex < 0, e e w , Teright ),
nright , u (11)

where n
eleft and n
eright are determined by imposing zero mass flux through the
walls: Z
d3 pe fe(±L/2,
e p, e t)e
px = 0. (12)

10
 Substituting Eq. (11) into Eq. (12) gives [24]:
s Z
2π
n
eleft = − d3 pe fe(−L/2,
e p, ee t)epx ,
meKe B Teleft pex <0
s Z
2π
n
eright = d3 pe fe(L/2,
e ee
p, t)e
px . (13)
m e e
e KB Tright pex >0

2.2. Reduced distributions

In the context of the channel flows considered in this paper, the dynamics
along the z direction is trivial. Moreover, in the heat transfer problem without
200 shear, the dynamics along the y axis also becomes trivial. In this context, it is
convenient to integrate out the trivial momentum space degrees of freedom at
the level of the model equation.
For notational convenience, let D = 3 represent the total number of degrees
of freedom of the momentum space. Denoting by d the number of non-trivial
momentum space degrees of freedom, the D − d degrees of freedom can be inte-
grated out and two reduced distribution functions, φe and χ
e, can be introduced
as follows [24, 27, 67, 71, 82, 89]:
Z Z
pe2d+1 + · · · pe2D e
φe = D−d
d pefe, χ
e= dD−d pe f. (14)
me

The evolution equations for φe and χ

e can be obtained by multiplying Eq. (1)
with the appropriate factors and integrating with respect to the D − d trivial
momentum space degrees of freedom:
     
e e e e
∂ φ pex ∂ φ 1 φ − φ∗ 
+ =− . (15)
∂et χe me ∂e
x χ e τe∗ χe−χe∗

Denoting using latin indices i and j the components corresponding to the

non-trivial directions (1 ≤ i, j ≤ d), the macroscopic moments given in Eqs. (3),

11
(9) and (5) can be obtained through:
   
n
e 1
  Z  
   
ρeu  d  e
 ei  = d pe  pei  φ,
   
Teij ξei ξej /m
e
    !
3 e e Z
e KB T
2n 1 ee
  = dd pe   ξj ξj φe + 1 χ
e , (16)
qei ξei /m e 2m
e 2

where the summation over the repeated index j is implied.

For the Shakhov model, φeS and χ
eS are given by [27]:
!
1 − Pr ξej ξej
φeS =φeMB (1 + Sφ ), Sφ = − d − 2 qei ξei ,
(D + 2)e
nKe 2 Te2 m
eK e B Te
B
!
e B TeφeMB (1 + Sχ ), 1 − Pr ξej ξej
χ
eS =(D − d)K Sχ = − d qei ξei ,
(D + 2)e
nKe 2 Te2 m
eK e B Te
B

(17)
where again the summation over the repeated index j is implied. The reduced
Maxwell-Boltzmann distribution φeMB is:

φeMB = n
egex (e ex , Te) · · · ged (e
px , u ed , Te).
pd , u (18)

Before discussing the ES model, we first mention that the representation as

a D × D matrix of the pressure tensor Teαβ admits the following block decom-
position:  
Teij 0ib
Teαβ =  , (19)
0aj Pered δab
where the latin indices at the beginning of the alphabet run over the trivial
degrees of freedom, i.e. d < a, b ≤ D. With this convention, the top left and
bottom right blocks are d × d and (D − d) × (D − d) matrices with components
Teij and Pered δab , respectively, while the top right and bottom left blocks are
d × (D − d) and (D − d) × d null matrices, respectively. The Kronecker delta
δab takes the value 1 when a = b and 0 otherwise. The scalar quantity Pered is
obtained from Eq. (16):
Pd
DPe − j=1 Tejj
Pered = . (20)
D−d

12
 Using the same decomposition as in Eq. (19), the matrix Bαβ can be written as:
 
Bij 0ib
Bαβ =  , (21)
0aj Bred δab ,

where
1 1 − Pr Teij
Bij = δij − . (22)
Pr Pr Pe
The scalar quantity Bred is given by:

1 1 − Pr Pered
Bred = − . (23)
Pr Pr Pe
It can be seen that the determinant of B can be written as:

D−d
det B = Bred det B. (24)

This allows the integral of feES over the D − d trivial degrees of freedom to be
performed analytically, giving:
−1 e e
!
n
e Bij ξi ξj
φeES = √ exp − , (25)
(2π m
eKe B Te)d/2 det B 2m
eK e B Te

while χ e B TeredφeES and K

eES = (D − d)K e B Tered = Pered /e
n.

205 2.3. Ab initio transport coefficients

In this paper, we consider a series of comparisons between the results ob-

tained in the frame of the model equations introduced in the previous subsec-
tions and the results obtained using the DSMC method with ab initio particle
interactions. The connection between these two formulations can be made at
the level of the transport coefficients. The basis for the approach that we take
in this paper is to note that in the variable hard spheres model, the viscosity
has a temperature dependence of the form [5]

µ eref (Te/Teref )ω ,
e=µ (26)

where the tilde denotes dimensionful quantities, as discussed in the previous

subsection. The viscosity index ω introduced above takes the values 1/2 and

13
 1 for hard sphere and Maxwell molecules, respectively. For real gases, ω is in
general temperature-dependent. This temperature dependence is not known
210 analytically, however the values of µ
e and κ
e corresponding to a gas comprised
of molecules interacting via ab initio potentials can be computed numerically.
The supplementary material in Ref. [90] contains the data corresponding to 3 He
and 4 He consistuents in the temperature ranges 1 K ≤ Te ≤ 10000 K, while
the data for Ne covering the range 20 K ≤ Te ≤ 10000 K can be found in the
215 tables reported in Ref. [91]. In order to perform simulations of the heat transfer
problem (discussed in Sec. 7) at Teref = 1 K (for He constituents) and 20 K
(for Ne constituents), we require data for the transport coefficients also in the
temperature range 0.25 K ≤ Te ≤ 1 K and 5 K ≤ Te ≤ 20 K, respectively. These
data were obtained by the method described in Ref. [91].
The temperature dependence of the viscosity index ω is accounted for by
employing Eq. (26) in a piecewise fashion. Let n (1 ≤ n ≤ N ) represent the
index of the tabulated values Te1 < Te2 < . . . TeN of the temperature, where N
is the total number of available entries. Considering a temperature interval
Ten ≤ Te ≤ Ten+1 , we define
!ωn
Te ln(e
µn+1 /eµn )
µ
e (n)
(Te) = µ
en , ωn = , (27)
Ten ln(Ten+1 /Ten )
220 where µ
en and µ
en+1 are the tabulated values of the viscosity corresponding to
the temperatures Ten and Ten+1 , respectively. The above formula ensures that
the function µ e(n) (Ten ) = µ
e(n) satisfies µ e(n) (Ten+1 ) = µ
en and µ en+1 .
The Prandtl number Pr is also defined in a piecewise fashion. For the tem-
perature range Ten ≤ Te < Ten+1 , we define Prn as
cep µ
en
Prn = , (28)
κ
en
κn is the heat conductivity corresponding to the temperature Te = Ten ,
where e
retrieved from the tabulated data mentioned above.
225 In general, the temperatures encountered in our simulations are within the
bounds of the temperature range for which data is available for interpolation.
For completeness, we present a possible extension of the above procedure for

14
 values of the temperature which are outside the range spanned by the tabulated
data. In the case when Te < Te2 , we propose to use µ
e(Te) = µ
e(1) (Te) and Pr(T ) =
230 Pr1 . For Te > TeN , where TeN is the highest available temperature in the tabulated
e(Te) = µ
data, we propose to use µ e(N −1) (Te) and Pr(Te) = PrN .
The algorithm described in this section can be summarized through [88]:


 e(1) (Te),
µ Te < Te2 ,



e(Te) = µ
µ e(n) (Te), Ten < Te < Ten+1 , ,




µe(N −1) (Te), TeN < Te,


 Pr1 , Te < Te2 ,



Pr(Te) = Prn , Ten < Te < Ten+1 , , (29)




Pr , Te < Te,
N N

where n = 2, 3, . . . N − 1 refers to the index of the tabulated data.

2.4. Non-dimensionalization convention

All simulation results reported in this paper are based on the nondimen-
sionalization conventions employed in Ref. [27], which are summarized here for
completeness. In general, the dimensionless form A of a dimensional quantity
e is obtained by dividing the latter with respect to its reference value, A
A eref :

e
A
A= . (30)
eref
A
We employ the convention that dimensionless quantities are denoted without the
overhead tilde encountered for their dimensionful counterparts. The reference
temperature is taken as the average of the wall temperatures:

Teleft + Teright
Teref = . (31)
2

The reference speed is defined through:

s
e B Teref
K
e
cref = , (32)
me

15
 e takes the values 5.0082373 × 10−27 kg, 6.6464764 ×
where the particle mass m
235 10−27 kg, and 3.3509177 × 10−26 kg for 3 He, 4 He, and Ne, respectively.
The reference particle number density is taken as the average particle number
density over the channel:
Z e
L/2
1
n
eref = de
xne. (33)
e
L e
−L/2

The reference length is taken to be the channel width:

e ref = L.
L e (34)

Finally, the reference time is

s
e ref
L m
e
e
tref = e
=L . (35)
cref
e KB Teref
e

The dimensionless relaxation time τ∗ = τe∗ /e

tref in the S and ES models
becomes:
µ(T ) µ(T )
τS = √ , τES = √ , (36)
Pδ 2 Pr P δ 2
where the rarefaction parameter δ is defined through [11]:

Le Peref
δ= √ , (37)
µ
eref e
cref 2

In the above, µ e(Teref ) and Peref = n

eref = µ e B Teref .
eref K
The distribution function fe is nondimensionalized via

fepeD
ref
f= , (38)
n
eref
q
where peref = m e B Teref . The reduced distributions can be nondimensionalized
eK
in a similar fashion:
epd
φe epedref
χ
ref
φ= , χ= . (39)
n
eref Peref
This allows Eq. (15) to be written as:
     
∂  φ  px ∂  φ  1 φ − φ∗ 
+ =− . (40)
∂t χ m ∂x χ τ∗ χ − χ
∗

16
 3. Finite difference lattice Boltzmann models with mixed quadratures

In this section, the FDLB algorithm employed to solve Eq. (40) is briefly
described. There are three pieces to the algorithm, which will be outlined in
240 the following subsections. The first piece concerns the implementation of both
the time stepping and the advection, which will be addressed in Subsec. 3.1.
The second concerns the discretization of the momentum space using the full-
range and half-range Gauss-Hermite quadratures. This will be discussed in
Subsec. 3.2. The third and final piece is the projection of the collision term in the
245 model equation on the space generated by the full-range Hermite polynomials
for the direction parallel to the wall. Details will be given in Subsec. 3.3.

3.1. Time stepping and advection

In order to describe the time stepping algorithm, Eq. (40) is written as:

∂t F = L[F ], (41)

where F ∈ {φ, χ} represents the reduced distributions. Considering the equidis-

tant discretization of the time variable using intervals δt and using tn = nδt to
denote the time coordinate after n iterations, we employ the third order total
variation diminishing Runge-Kutta scheme to obtain Fn+1 at time tn+1 through
two intermediate steps [76, 77, 78]:

Fn(1) =Fn + δtL[Fn ],

3 1 (1) 1
Fn(2) = Fn + F + δtL[Fn(1) ],
4 4 n 4
1 2 (2) 2
Fn+1 = Fn + F + δtL[Fn(2) ]. (42)
3 3 n 3

As pointed out by various authors [80, 81, 82], an accurate account for the
Knudsen layer phenomena requires a sufficiently fine grid close to the wall.
This can be achieved by performing a coordinate change from x = x e to the
e/L
coordinate η, defined through [27, 71, 82, 88]:

tanh η
x= , (43)
2A

17
 where the stretching parameter A controls the grid refinement. When A →
0, the grid becomes equidistant, while as A → 1, the grid points accumulate
250 towards the boundaries at x = ±1/2. The channel walls are located at η =
±arctanh A. In this paper, we always use the stretching corresponding to A =
0.98.
The η coordinate is discretized symmetrically with respect to the channel
centerline (where x = 0 and η = 0). On the right half of the channel, S
equidistant intervals of size δη = arctanh A/S are employed. In the case of the
Couette flow, which is symmetric with respect to the channel centerline, the
simulation setup contains only the domain 0 ≤ η ≤ A and the total number of
grid points is equal to S [27, 71]. The center of cell s (1 ≤ s ≤ S for the right
half of the channel and −S < s ≤ 0 for its left half) is located at ηs = (s − 12 )δη.
At each node s, the advection term is computed using the third order upwind
method, implemented using a flux-based approach:
     
px ∂F px ∂η ∂F Fs+1/2 − Fs−1/2
= = 2A cosh2 ηs + O(δη 3 ).
m ∂x s m ∂x s ∂η s δη
(44)
The stencil employed for the flux Fs+1/2 is chosen depending on the sign of the
advection velocity px /m:

1
 5 1
px  3 Fs+1 + 6 Fs − 6
Fs−1 , px > 0,
Fs+1/2 = (45)
m  1 5 1
 Fs + Fs+1 − Fs+2 , px < 0.
3 6 6
The diffuse reflection boundary conditions in Eq. (11) specify the distribu-
tions φ and χ on the channel walls. For definiteness, we will refer to the right
boundary, which is located at ηS+1/2 = arctanhA. In order to perform the ad-
vection at node S for the particles traveling towards the wall (having px > 0),
the value of the distribution function in the node s = S + 1 is required. This
value can be obtained using a third order extrapolation from the fluid nodes:
px >0
FS+1 = 4FS − 6FS−1 + 4FS−2 − FS−3 . (46)

It can be shown that the third order accuracy in the sense of Eq. (44) is preserved
when the fluxes FS+1/2 and FS−1/2 are computed using Eq. (45). For the

18
 particles traveling towards the fluid (px < 0), the nodes at S + 1 and S + 2
must be populated. According to the diffuse reflection concept, summarized in
Eq. (11), the reduced distributions at s = S + 1/2 are set to:

φS+1/2 =φMB (nright , uw , Tright ),

χS+1/2 =(D − d)Tright φS+1/2 , (47)

where Tright = 1 + ∆T /2 is the temperature on the right wall (∆T = 0 in the

case of Couette flow). The distributions in the ghost nodes at S + 1 and S + 2
can be set to [79]:

px <0 16 1
FS+1 = F − 3FS + FS−1 − FS−2 ,
5 S+1/2 5
px <0
FS+2 =4FS+1 − 6FS + 4FS−1 − FS−2 . (48)

The expression for FS+2 can be seen to represent a third order extrapolation
from the nodes with S − 2 ≤ s ≤ S + 1. In the expression for FS+1 , the distri-
255 bution at the wall, FS+1/2 is employed. It can be checked by direct substitution
in Eq. (44) that the third order accuracy is preserved when the ghost nodes are
populated as indicated above.
The density nright in Eq. (47) can be computed using the discrete equivalent
of Eq. (12): Z
dd p ΦS+1/2 = 0, (49)

where ΦS+1/2 is the flux corresponding to the reduced distribution φ, computed

using Eq. (45). Using Eq. (48), the flux for outgoing particles is:
 
px <0 px 8 5 1 2
ΦS+1/2 = φS+1/2 + φS − φS−1 + φS−2 ,
m 15 6 2 15
 
px <0 p x 8 4 1 1
ΦS−1/2 = − φS+1/2 + φS + φS−1 + φS−2 , (50)
m 15 3 6 30

where the flux ΦS−1/2 is given above for completeness. Due to the above expres-
px <0
sion for ΦS+1/2 , it can be seen that the unknown density, nright , enters Eq. (49)
through the distribution φS+1/2 , which is fixed by boundary conditions for mo-
menta pointing towards the fluid (px < 0), according to Eq. (47). Splitting the

19
integration domain in Eq. (49) in two domains, corresponding to positive and
negative values of px , the integral for px < 0 of φS+1/2 can be computed as
follows: Z r
d px Tright
d p φS+1/2 = −nright . (51)
px <0 m 2πm
Taking into account Eqs. (50) and (51), the following expression is obtained for
nright :
s Z
15 2πm
nright = dd p ΦS+1/2
8 Tright px >0
Z  
px 5 1 2
− dd p φS − φS−1 + φS−2 . (52)
px <0 m 6 2 15

For completeness, we also give below the expressions for ΦS+1/2 when px > 0:
 
px >0 px 13 13 4 1
ΦS+1/2 = φS − φS−1 + φS−2 − φS−3 . (53)
m 6 6 3 3

In the case of the Couette flow, only the nodes with 1 ≤ s ≤ S comprise
the fluid domain, while bounce-back boundary conditions are imposed on the
channel centerline (s = 1/2). The nodes with s < 1 become ghost nodes, which
are populated according to:


F0 (p) = −F1 (−p),
px <0 : F0 (p) = −F1 (−p), px >0 : (54)

F (p) = −F (−p).
−1 2

3.2. Momentum space discretization

Through the discretization of the momentum space, the integrals defining

the macroscopic moments in Eq. (16) are replaced by quadrature sums, i.e.:
   
n 1
   
  X 
ρu  ≃  p  φκ ,
 i  κ;i 
  κ  
Tij ξκ;i ξκ;j /m
     
3
nT X 1 X 1
2  ≃   ξκ;j ξκ;j φκ + 1   χσ , (55)
q ξ /m 2m 2 σ ξ /m
i
κ κ;i σ;i

20
 where κ and σ collectively denote the indices labeling the momenta correspond-
260 ing to the discrete populations φκ and χσ .
The discretization on the axis perpendicular to the walls (the x axis), is
performed using the half-range Gauss-Hermite quadrature prescription [44]. In
principle, px can be discretized separately for the φ and χ distributions. We
take advantage of this freedom when deriving the FDLB models for the hydro-
dynamic regime, which is discussed in Appendix A. For the flows considered
in Sections 5, 6 and 7, we consider the same quadrature order Qφx = Qχx ≡ Qx
on each x semiaxis. For definiteness, we assume Qφx = Qχx henceforth, unless
otherwise stated. Focusing on the distribution φ, the discrete momentum com-
ponents px,kx (1 ≤ kx ≤ 2Qx) are linked to the roots of the half-range Hermite
polynomial hQx (z) of order Qx via:


p0,x zkx , 1 ≤ kx ≤ Qx ,
px,kx = (56)

−p z
0,x kx −Qx , Qx < kx ≤ 2Qx ,

where hQx (zkx ) = 0 for 1 ≤ kx ≤ Qx , while p0,x represents a constant momen-

tum scale (we set p0,x = 1 for the rest of this paper). The same considerations
apply for the distribution χ, after replacing kx with the index sx (1 ≤ sx ≤ 2Qx ).
When d = 1, the populations φκ ≡ φkx and χσ ≡ χsx are linked to the
continuum distributions φ and χ through:

p0,x wkhx (Qx ) p0,x wshx (Qx )

φkx = φ(px,kx ), χsx = χ(px,sx ), (57)
ω(px,kx ) ω(px,sx )

where px ≡ px /p0,x , while the weight function ω(z) is defined through:

1 2
ω(z) = √ e−z /2 . (58)
2π

The quadrature weights wkh (Q) can be computed using [44, 92]

px,k a2Q
wkh (Q) = √
h2Q+1 (px,k )[px,k + h2Q (0)/ 2π]
px,k a2Q−1
= √ . (59)
h2Q−1 (px,k )[px,k + h2Q (0)/ 2π]

21
In the above, aQ = hQ+1,Q+1 /hQ,Q represents the ratio of the coefficients of
the leading power of px in hQ+1 (px ) and hQ (px ). Specifically, the notation hℓ,s
refers to the coefficient of xs appearing in hℓ (px ), namely:
ℓ
X
hℓ (px ) = hℓ,s psx . (60)
s=0

In the case when the boundaries are moving, d = 2 and the momentum
component py is discretized using the full-range Gauss-Hermite quadrature pre-
scription. As remarked in Refs. [27, 44], a small quadrature order is sufficient
to ensure the exact recovery of the dynamics along this axis. To assess the
quadrature orders for the φ and χ distributions, we consider the expansions of
φ and χ with respect to the full-range Hermite polynomials for the py degree of
freedom:
       
X∞ Z ∞
φ ω(py ) 1 Φℓ Φℓ φ
 = Hℓ (py )   ,  = dpy Hℓ (py )   , (61)
χ p0,y ℓ! Xℓ Xℓ −∞ χ
ℓ=0

where Φℓ and Xℓ are expansion coefficients [not to be confused with the fluxes
in Eq. (49)]. Substituting the above expansions in Eq. (40) gives:
   
 
∂ p x ∂  Φℓ  1  Φℓ − Φ∗ℓ 
+ =− . (62)
∂t m ∂x Xℓ τ∗ Xℓ − X ∗
ℓ

It can be seen that the moments Φℓ and Xℓ of order ℓ are coupled with those
of order ℓ′ 6= ℓ only through the collision term. However, the equilibrium pop-
ulations φ∗ and χ∗ are determined exclusively by the macroscopic quantities
corresponding to the collision invariants, n, u and T , as well as Tij (for the ES
model) and qi (for the S model). These quantities can be written in terms of

22
 the coefficients Φℓ′ with 0 ≤ ℓ′ ≤ 3 and Xℓ′ with 0 ≤ ℓ′ ≤ 1, as follows:
   
n Z Φ0
  ∞  
   
ρux  = dpx  Φ0 px  ,
  −∞  
ρuy Φ1 p0,y
   
2
Txx ξ Φ
  Z ∞ dp  x 0

  x  
Txy  =  ξx (Φ1 p0,y − muy Φ0 ), ,
  −∞ m  
Tyy Φ2 p20,y − 2mp0,y uy Φ1 + (p20,y + m2 u2y )Φ0
    !
3 Z ∞
2 nT 1 ξx2 + p20,y + m2 u2y p20,y 1
 = dpx   Φ0 − p0,y uy Φ1 + Φ2 + X 0 ,
qx −∞ ξx /m 2m 2m 2
Z ∞ "
p30,y 3p20,y uy p0,y 2
qy = dpx 2
Φ 3 − Φ2 + (ξ + 3p20,y + 3m2 u2y )Φ1
2 x
−∞ 2m 2m 2m
uy 2 p0,y uy i
− (ξx + 3p20,y + m2 u2y )Φ0 + X1 − X0 . (63)
2m 2m 2

It can be seen that for a given value of ℓ, Eq. (62) involves only terms with ℓ′
265 such that 0 ≤ ℓ′ ≤ max(ℓ, 3) for Φℓ and 0 ≤ ℓ′ ≤ max(ℓ, 1) for Xℓ . Thus, it can
be concluded that the moment system with respect to the py degree of freedom
is closed when the terms up to ℓ = 3 and 1 in the series expansions of φ and
χ, respectively, are included. Moreover, the dynamics (and therefore stationary
state properties) of the moments in Eq. (63) is recovered exactly when the
270 series for φ and χ in Eq. (61) are truncated at ℓ = 3 and 1, respectively. This
truncation is equivalent to considering the quadrature orders Qφy = 4 and Qχy =
2, in the sense that employing higher order quadratures yields results which are
exactly equivalent (up to numerical errors due to finite machine precision) to
those obtained using Qφy = 4 and Qχy = 2. We discuss below the discretization
275 corresponding to these quadrature orders.
The roots of the Hermite polynomial H4 (z) = z 4 − 6z 2 + 3 of order 4 are
known analytically [49, 93, 94, 95, 96, 97]:
q q
√ √
pφy,1 = − 3 + 6, pφy,2 = − 3 − 6,
q q
√ √
pφy,3 = 3 − 6, pφy,4 = 3 + 6, (64)

23
 where pφy,ky ≡ pφy,ky /pφ0,y is normalized with respect to an arbitrary scaling
factor pφ0,y , which we set to 1 in this paper. For the χ populations, the discrete
momentum components along the y axis can be found via the roots of H2 (z) =
z 2 − 1:
pχy,1 = −1, pχy,2 = 1, (65)

where pχy,sy ≡ pχy,sy /pχ0,y and pχ0,y = 1.

The connection between the discrete populations φκ and χσ and their con-
tinuous counterparts is given by the 2D extension of Eq. (57):
φ
p0,x wkhx (Qx ) p0,y wkHy (Qφy )
φκ = φ(px,kx , pφy,ky ),
ω(px,kx ) ω(pφy,ky )
χ
p0,x wshx (Qx ) p0,y wsHy (Qχy )
χσ = χ(px,sx , pχy,sy ), (66)
ω(px,sx ) ω(pχy,sy )

where ω(z) is defined in Eq. (58). The quadrature weights for the full-range
Gauss-Hermite quadrature can be computed via [44, 50, 51]:

Q∗y !
wkH (Q∗y ) = , (67)
[HQ∗y +1 (zk )]2

where zk (1 ≤ k ≤ Q∗y ) is the k’th root of HQ∗y (z). In particular, the weights
for Qφy = 4 and Qyχ = 2 are given by:
√
H H 3− 6
wy,1 (4) =wy,4 (4)
= ,
12√
H H 3+ 6
wy,2 (4) =wy,3 (4) = ,
12
H H 1
wy,1 (2) =wy,2 (2) = . (68)
2

We now summarize the procedure described above. In the case of the heat
transfer problem, the one-dimensional momentum space is discretized following
the half-range Gauss-Hermite quadrature prescription using Qφx = Qχx = Qx
280 quadrature points on each semiaxis for both φ and χ.
For the shear flow problems, the y axis of the momentum space is discretized
separately for φ and χ. The total number of quadrature points used to dis-
cretize the momentum space for φ is 2Qx Qφy = 8Qx , while for χ, 2Qx Qχy = 4Qx

24
 quadrature points are required, resulting in a total number of 12Qx discrete
285 populations.

3.3. Projection of the collision term

Part of the lattice Boltzmann paradigm is to replace the local equilibrium

distribution by a polynomial expansion, such that the collision invariants ψ ∈
{1, pi , p2 /2m} are exactly preserved. This requires that, after the discretization
of the momentum space, the folllowing quadrature sums are exact:
   
X 1 n X ξκ;i ξκ;i 1X 3
  φ∗;κ =   , φ∗;κ + χ∗;σ = nT. (69)
κ pκ;i ρui κ
2m 2 σ 2

The above relations can be exactly ensured by first expanding φ∗ and χ∗ with
respect to the Hermite polynomials (half-range on the x and full-range on the
φ/χ
y axes, if required), followed by a truncation of the sums at orders 0 ≤ Ni <
φ/χ
290 Qi . This approach is followed in Appendix A in order to tackle flows in the
hydrodynamic regime.
For the flows with 0.1 ≤ δ ≤ 10 considered in Sections 5, 6 and 7, we
follow a hybrid approach. Namely, the equilibrium distributions φ∗ and χ∗
are projected onto the set of full-range Hermite polynomials with respect to
295 the axis parallel to the walls (no projection is required in the case of the heat
transfer between stationary plates problem). Then, the expansion coefficients
are evaluated directly, following the standard DVM approach. This hybrid
approach is motivated as follows.
On the x axis, the quadrature order Qx is considered to be equal for both
300 φ and χ. Since we are interested in performing simulations in the slip flow and
transition regime, we need in general high values of Qx (i.e., Qx ≥ 7 will be
required [27]). Let us now assume the equilibrium distributions are expanded
with respect to the half-range Hermite polynomials up to order Nx = Qx − 1. It
is expected that the coefficients of the expansion grow with Nx as ∼ Nx !MaNx .
305 Since the simulations that we are considering are performed in the non-linear
regime, where Ma > 1, high expansion orders may be required (we use Qx = 50

25
 at δ = 0.1), such that the individual terms in the series expansion can be large.
The addition and subtraction of these terms typically leads to a significantly
smaller remainder, which can easily be poluted by numerical errors due to finite
310 numerical precision. It is a well-known limitation of the (FD)LB algorithm that
the polynomial expansion of the equilibrium distribution is not well suited for
high-Mach number flows. On the other hand, directly evaluating the equilibrium
distributions discussed in Sec. 2.1 when computing the equilibrium moments
in Eq. (69) at Qx ≥ 7 is already quite accurate for δ . 10 when the half-
315 range Gauss-Hermite quadrature is employed (in this case, 2Qx ≥ 14 quadrature
points are employed on the px axis). Thus, we find the loss in precision due to the
integration using Gauss quadratures of non-polynomial functions via Eq. (69)
to be irrelevant.
At larger values of δ, the physical time to reach the steady state increases.
320 Over a longer time interval, the errors in the recovery of the conservation laws
due to the inaccurate integration of the equilibrium distribution accumulate,
affecting the accuracy of the properties of the stationary state. This prob-
lem can be alleviated by projecting the equilibrium distribution onto the space
of half-range Hermite polynomials also on the px direction, as discussed in
325 Appendix A.
We further discuss in detail the implementation of the S and ES collision
terms for the d = 1 case encountered in the heat transfer between stationary
plates problem (Subsec. 3.3.1). In the d = 2 case, encountered for the Couette
flow and heat transfer between moving plates problem, the implementation of
330 the ES and S models is discussed separately in Subsecs. 3.3.2 and 3.3.3, respec-
tively.

3.3.1. d = 1 case
In the case of the ES model, the equilibrium distribution functions can be
found from Eq. (25). When d = 1, the equilibrium distribution function is
 
n (px − mux )2
φES = √ exp − , (70)
2πmT Bxx 2mT Bxx

26
while χES = 2TredφES , where Tred = Pred /n. In the above, Bxx and Pred are
given by:
1 1 − Pr Txx 3 1
Bxx = − , Pred = P − Txx . (71)
Pr Pr P 2 2
The transition to the discrete system is made via Eq. (57):

p0,x wkhx (Qx ) p0,x wshx (Qx )

φES;kx = φES (px;kx ), χES;sx = χES (px;sx ), (72)
ω(px;kx ) ω(px;sx )
where 1 ≤ kx , sx ≤ 2Qx and ω(z) is defined in Eq. (58).
For the S model, the equilibrium distributions φS and χS can be obtained
from Eq. (17):
 
1 − Pr ξx2
φS =ngx (1 + Sφ ), Sφ = − 3 qx ξx ,
5nT 2 mT
 2 
1 − Pr ξx
χS =2P gx (1 + Sχ ), Sχ = − 1 qx ξx , (73)
5nT 2 mT
where P = nT and gx ≡ g(px , ux , T ) is given through Eq. (2), while ξx =
px − mux . As in Eq. (72), the equilibrum distributions after discretization are
computed using:
p0,x wkhx (Qx ) p0,x wshx (Qx )
φS;kx = φS (px;kx ), χS;sx = χS (px;sx ). (74)
ω(px;kx ) ω(px;sx )

3.3.2. d = 2 case: ES model

−1
In the d = 2 case, the exponent Bij ξi ξj in Eq. (25) can be written as:
!2
−1
−1 −1 Bxy ξx2 −1 −1 −1 2


Bij ξi ξj = Byy ξy + −1 ξx + −1 Bxx Byy − (Bxy ) . (75)
Byy Byy

Noting that the inverse of Bij is given by:

 
1 B −Bxy
−1
Bij =  yy , (76)
detB −Bxy Bxx

φES can be factorized as follows:

 
ξx Bxy T detB
φES = ng(px , ux , T Bxx)g py , uy + , . (77)
mBxx Bxx
A similar factorization holds for χES = Tred φES , where Tred = Pred /n and

Pred = 3P − Txx − Tyy . (78)

27
 (N φ ) (N χ )
We now seek to replace φES and χES with the expansions φESy and χESy
with respect to the Hermite polynomials Hℓ (py ) containing only terms up to
orders Nyφ = Qφy − 1 = 3 and Nyχ = Qχy − 1 = 1, respectively. Defining:
ξx Bxy detB
ζy = uy + , Ty = T, (79)
mBxx Bxx
Eq. (77) reduces to φES = ng(px , ux , T Bxx)g(py , ζy , Ty ). The trailing function
g(py , ζy , Ty ) is expanded with respect to Hℓ (py ) up to order Ny∗ ∈ {Nyφ , Nyχ }, as
follows:
N∗
(Ny∗ )
ω(py ) Xy
1
g (py , ζy , Ty ) = Hℓ (py )GℓH (ζy , Ty ). (80)
p0,y ℓ!
ℓ=0

The expansion coefficients GℓH (ζy , Ty ) were obtained analytically in Eq. (C.13)
in Ref. [44]. Below, we reproduce the coefficients for 0 ≤ ℓ ≤ 3:

G0H = 1, G1H = U, G2H = U2 + I, G3H = U3 + 3UI. (81)

Identifying U and I from Eq. (C.16) of Ref. [44] with the following expressions,
mζy mTy
U(ζy ) = , I(Ty ) = − 1, (82)
p0,y p20,y

g (1) (py , ζy , Ty ) necessary for the construction of χES can be written as:
ω(py )  
g (1) (py , ζy , Ty ) = H0 (py ) + H1 (py )U(ζy ) . (83)
p0,y
The function g (3) (py , ζy , Ty ) required for φES , is given by:

ω(py ) 1  
g (3) (py , ζy , Ty ) = H0 (py ) + H1 (py )U(ζy ) + H2 (py ) U2 (ζy ) + I(Ty )
p0,y 2!

1  3 
+ H3 (py ) U (ζy ) + 3U(ζy )I(Ty ) . (84)
3!
With the above ingredients, after discretization, φES
κ can be evaluated using:
φ
p0,x wkhx (Qx ) p0,y wkHy (Qφy )
φES
κ =n φ
g(px,kx , ux , T Bxx)g (3) (pφy,ky , ζy;kx , Ty ), (85)
ω(px,kx ) ω(py,ky )
Bxy
where ζy;kx = uy + mBxx ξx,kx . Similarly, χES
σ is:
χ
p0,x wshx (Qx ) p0,y wsHy (Qχy )
χES = Pred g(px,sx , ux , T Bxx)g (1) (pχy,sy , ζy;sx , Ty ).
σ
ω(px,sx ) ω(pχy,sy )
(86)

28
 In Eqs. (85) and (86), the function g(px , ux , T Bxx) is evaluated directly. Its
expression is reproduced below for convenience:
 
1 (px,kx − mux )2
g(px,kx , ux , T Bxx) = p exp − . (87)
2πmT Bxx 2mT Bxx

335 3.3.3. d = 2 case: S model

In the case of the Shakhov model, φS and χS can be written as:
!
1 − Pr ξx2 + ξy2
φS =ngx gy (1 + Sφ ), Sφ = − 4 (qx ξx + qy ξy ),
5nT 2 mT
!
1 − Pr ξx2 + ξy2
χS =nT gxgy (1 + Sχ ), Sχ = − 2 (qx ξx + qy ξy ), (88)
5nT 2 mT

where gx ≡ g(px , ux , T ) and gy ≡ g(py , uy , T ) are one-dimensional Maxwell-

Boltzmann distributions introduced in Eq. (2). The functions φS and χS can be
expanded with respect to the full-range Hermite polynomials Hℓ (py ), as follows:
   
∞ φ;H
ω(py ) X
φS 1 GS;ℓ
 = Hℓ (py )  . (89)
χS p0,y ℓ! G χ;H
ℓ=0 S;ℓ

φ/χ;H
The expansion coefficients GS;ℓ can be written as:
     
φ;H φ;H
GS;ℓ n 1 − Pr G
  = gx   GℓH (uy , T ) +  S;ℓ  . (90)
χ;H
GS;ℓ P 5nT 2 Gχ;H
S;ℓ

The coefficients GℓH have the same form as in Eq. (81), where the factors U ≡
U(uy ) and I ≡ I(T ) are given by Eq. (C.16) in Ref. [44]:

muy mT
U(uy ) = , I(T ) = − 1. (91)
p0,y p20,y

Denoting:
    
Z ∞
Isφ ξx2 + ξy2 4
 = dpy g(py , uy , T )(qx ξx + qy ξy )  −   ξys , (92)
Isχ −∞ mT 2

29
 the coefficients G∗;H
S;ℓ (∗ ∈ {φ, χ}) in Eq. (90) can be obtained as:

1
G∗;H ∗
S;0 =I0 , G∗;H
S;1 = (I1∗ + muy I0∗ ),
p0,y
1
G∗;H
S;2 = [I ∗ + 2muy I1∗ + (m2 u2y − p20,y )I0∗ ],
p20,y 2
1
G∗;H
S;3 = 3 [I3∗ + 3muy I2∗ + 3(m2 u2y − p20,y )I1∗
p0,y

+ muy (m2 u2y − 3p20,y )I0∗ ]. (93)

Finally, the terms Is∗ can be obtained by direct integration in Eq. (92):
         
I0φ ξ 2 3 I φ
ξ 2 −1
  = qx ξx  x −   ,  1  = qy mT  x +   ,
I0χ mT 1 I1χ mT 1
         
I2φ ξ 2 −1 I φ
ξ 2 1
  = qx ξx mT  x +   ,  3  = 3qy (mT )2  x +   .
I χ mT 1 I χ mT 3
2 3

(94)

Putting the pieces together, the discrete populations φS;κ and χS;σ can be
computed using:
3  
p0,x wkhx (Qx ) H
X 1 φ H 1 − Pr φ;H
φS;κ =n g(px,kx , ux , T )wky (4) Hℓ (py,ky ) Gℓ + G ,
ω(px,kx ) ℓ! 5nT 2 S;ℓ
ℓ=0
X 1  
p0,x wshx (Qx ) H 1 χ H 1 − Pr χ;H
χS;σ =nT g(px,sx , ux , T )wsy (2) Hℓ (py,sy ) Gℓ + G .
ω(px,sx ) ℓ! 5nT 2 S;ℓ
ℓ=0

(95)

4. Simulation methodology

This section briefly summarizes the methodology employed for obtaining

the numerical results discussed in the next sections. Three applications are
considered in this paper, namely the heat transfer between stationary plates
340 (Sec. 5), the Couette flow between plates at the same temperature (Sec. 6), and
the heat transfer between moving plates (Sec. 7).
In all cases, the simulation results are presented for three values of the rar-
efaction parameter, namely δ = 10, 1 and 0.1. For all applications, we take the

30
 working gas to be comprised of 3 He or 4 He molecules. Additionally, in the case
345 of the heat transfer between moving plates problem, we also report results for
Ne. The reference temperature Teref , defined in Eq. (31), varies between 1 K and
3000 K for the 3 He and 4 He constituents and between 20 K and 5000 K for the
Ne constituents.
Quantitative comparisons are performed by considering a set of dimension-
less numbers. In the context of the flows between moving walls (discussed in
Sections 6 and 7), the shear stress is used to define the quantity [11]

Texy e
cref
Π=− √ . (96)
Peref u
ew 2
It can be shown that, in the stationary state, Π is constant throughout the
channel. In order to access the non-linear regime, we set the wall velocities to
√ q
u cref 2 = 2K
ew = e e B Teref /m,
e such that the Mach number is

2e
uw
Ma = ≃ 2.19, (97)
e
cs
q
where e e B Teref /m
cs = γ K e is the speed of sound and γ = 5/3 is the adiabatic
index for a monatomic ideal gas. After non-dimensionalization, Π is computed
through
1
Π = − Txy . (98)
2
In the heat transfer problems, discussed in Sections 5 and 7, the longitudinal
heat flux qex (perpendicular to the x axis) is used to introduce

qx + Texy u
(e ey )Teref
Q=− √ , (99)
Peref e g 2
cref ∆T
which is again constant throughout the channel. The second term in the nu-
merator vanishes when the walls are stationary (i.e., in Sec. 5). We consider
the nonlinear regime, in which the ratio between the temperature difference
g = Teright − Teleft and Teref , defined in Eq. (31), is
∆T

g
∆T Teright − Teleft
=2 = 1.5. (100)
Teref Teright + Teleft

31
 After non-dimensionalization, the wall temperatures are Tleft = 0.25 and Tright =
1.75, while Q is obtained via:
√ √
2 2 2
Q= Πuy − qx . (101)
3 3
In the context of the Couette flow, we further consider two more quantities.
The first is the dimensionless half-channel heat flow rate, defined through
Z e
L/2
2 qey
Qy = de
x . (102)
e
L 0 e
Pref u
ew
The second is related to the heat transfer through the domain wall, and is
defined through:
e cref
qex (L/2)e u
ey (1/2)
Qw = √ = Π, (103)
e
Pref u 2
ew 2 u
ew

where the second equality follows after noting that qex + Texy u
ey = 0 in the sta-
350 tionary state of the Couette flow.
In practice, the quantities Π and Q exhibit a mild coordinate dependence
in the stationary state due to the errors of the numerical scheme. The values
reported in the applications sections are obtained by averaging Π and Q over
the simulation domain, as follows:
   
Z L/2
e
Π 1 Π(e
x)
 = x
de . (104)
Q e −L/2
L e Q(e
x)

In the case of the Couette flow, Π(−e

x) = Π(e
x) is used to reduce the integration
domain to 0 ≤ x e
e ≤ L/2.
The FDLB methodology is discussed in Subsec. 4.1 and the DSMC method-
ology is summarized in Subsec. 4.2.

355 4.1. FDLB methodology

For the heat transfer problems, the FDLB simulations are performed on grids
comprised of 2S cells. For the Couette flow simulations, we take advantage of
the symmetry and use only S cells on the half-channel. Each cell has the width
δη = arctanhA/S with respect to the η coordinate and the stretching parameter
360 entering Eq. (43) is set to A = 0.98.

32
 δ 1000 100 10 1 0.1
Nodes per half-channel S 32 32 32 32 16
−4 −4 −4 −4
Time step δt 10 2.5 × 10 5 × 10 5 × 10 5 × 10−4
Qx 11 8 7/8 11 50
Discrete populations (d = 1) 44 32 28/32 44 200
Discrete populations (d = 2) 132 96 84/96 132 600

Table 1: Discretisation details for the FDLB method for various values of the rarefaction
parameter δ. For δ = 10, both Qx = 7 (for the heat transfer between stationary plates and
the Couette flow) and Qx = 8 (for the heat transfer between moving plates) are employed.

At δ = 10, the quadrature order on the x axis is set to Qx = 7 for the

Couette flow and heat transfer between stationary plates problems, while for
the heat transfer between moving plates, Qx = 8 is used. For δ = 1 and 0.1,
the quadrature order is increased to Qx = 11 and 50, respectively, in order
365 to capture the rarefaction effects. The quadrature orders and total number
of discrete populations are shown in Table 1. For completeness, Table 1 also
includes information for the δ = 1000 and 100 cases. For these larger values
of δ, lower quadrature orders can be employed if the equilibrium distributions
are projected onto the space of half-range Hermite polynomials, as discussed in
370 Appendix A.
The simulation is performed until the stationary state is achieved. The time
steps δt = 10−4 , 2.5 × 10−4 and 5 × 10−4 were employed for δ = 1000, δ = 100
and δ ≤ 10, respectively. The number of points on the half-channel is set to
S = 32 and 16 for δ ≥ 1 and δ = 0.1, respectively. The number of iterations
375 performed to reach the stationary state is 5 × 106 , 2 × 105 , 6 × 104 , 4 × 104 and
2 × 105 for δ = 1000, 100, 10, 1 and 0.1, respectively.
In order to assess the accuracy of the simulation results, another set of
simulations is performed using Qx = 16 for δ = 1000 and 100, Qx = 40 for
δ = 10 and 1, while for δ = 0.1, Qx = 200 is employed. The spatial grid is
380 refined by a factor of 4, such that S = 128 is used for δ ≥ 1, while for δ = 0.1,

33
 S = 64 points are used on the half-channel. The time step in this case is set
to δt = 5 × 10−5 for δ ≥ 1, and δt = 4 × 10−5 for δ = 0.1. We compared the
results obtained for Q, Π, Qw and Qy and found that the relative differences
between the results obtained within the two sets of simulations were below 0.1%
385 for δ ≤ 10. At δ = 100 and 1000, the relative error has a larger magnitude, since
it is computed at the level of quantities that tend to zero in the large δ limit. It
can be seen in Figs. A.15 and A.16 that the absolute error is confortably under
0.1% at S = 32.
In order to compute the integrals over the discretized domain, a fourth order
rectangle method is used, summarized below:
Z e
L/2 Z arctanh A
1 1 dη
de
x M (e
x) = M (η)
e
L e
−L/2 A −arctanh A cosh2 η
S
X
arctanh A Ms
= fs , (105)
AS
s=−S+1
cosh2 ηs

where Ms ≡ M (ηs ) and



 13
 , s = 4i or 4i + 1,
fs = 12 (106)

 11
 , s = 4i + 2 or 4i + 3.
12
A comparison with various approaches described in the literature confirms
390 that our method is efficient. For example, in Ref. [85], the heat transfer between
stationary plates problem is simulated in the linear regime using a full-range
Gauss-Hermite discretisation of order 128 (only one distribution is required in
this case), which is above the number of velocities that we employ at δ = 0.1,
when 2Qx = 100 (we note that this discretisation allows us to access the non-
395 linear regime). Another example can be seen in Ref. [24], where the 2D velocity
space comprised of px and py is discretised using polar coordinates (p, θ). In
the transition regime, a number of 16 × 101 = 1616 velocities are employed [24],
compared to only 100 with our approach (when d = 2, we employ 400 and 200
velocities for the φ and χ distributions, respectively).

34
 HT SH HT-SH

δ TS (s) TES (s) TS (s) TES (s) TS (s) TES (s)

1000 2094 1556 4621 3523 8900 6621
100 62 47 133 99 272 194
10 16 12 35 26 78 58
1 16 12 36 26 72 52
0.1 170 120 400 300 790 570

Table 2: Computational times (in seconds) required to reach the steady state using the FDLB
method for the heat transfer between stationary plates (HT), Couette flow (SH) and heat
transfer under shear (HT-SH) problems, considered in Sections 5, 6 and 7. The data for
δ = 100 and 1000 is added for completeness.

400 4.2. DSMC methodology

e ≤x
The DSMC calculations were carried out dividing the space −L/2 e
e ≤ L/2
into 800 cells, considering 200 particles per cell in average, and using the time
√ q
step δe e 2e
t equal to 0.002L/ cref = K
cref , where e e B Teref /m
e is defined in Eq. (32).
The shear stress Π and heat flux Q, defined in Eqs. (96) and (99), were calculated
405 by counting the momentum and energy brought and taken away by all particles
on both surfaces. To reduce the statistical scattering, the macroscopic quantities
were calculated by averaging over 5 × 105 samples. These parameters of the
numerical scheme provide the total numerical error of Q and Π less than 0.1%,
estimated by carrying out test calculations with the double number of cells,
410 the double number of particles and reducing the time step by a factor of 2.
The relative divergence of Π and Q, calculated on the difference surface using
an additional accuracy criterion, does not exceed 0.01%. The details of the
numerical scheme and the method used to calculate the look-up tables can be
found in Ref. [11].

415 4.3. Computational time analysis

It is known that the DSMC method suffers from stochastic noise, which
persists after the steady state is reached. This noise can be eliminated through

35
 100
(a) HT HT 100
-1 (b) HT HT
10 SH SH SH SH
HTSH HTSH HTSH HTSH
10-2
-2
(S) (ES) 10
(S) (ES)
-3
10
10-4
10-4
L2 [ T ]

L2 [ T ]
10-5
-6
10
10-6

10-7
10-8
-8
10
(δ = 10) (δ = 1)
10-9 10-10
0 5 10 15 20 25 30 2 4 6 8 10 12 14 16 18 20
t t

0
10
(c) HT HT
SH SH
-2 HTSH HTSH
10
(S) (ES)

-4
10
L2 [ T ]

-6
10

-8
10

(δ = 0.1)
-10
10
10 20 30 40 50 60 70 80 90 100
t

Figure 1: Approach to steady state, assessed at the level of the temperature profile, at (a)
δ = 10, (b) δ = 1 and (c) δ = 0.1. The results for heat transfer under stationary walls (HT),
Couette flow (SH) and heat transfer under shear (HTSH) are represented using squares, circles
and rhombi, respectively. The results for the S and ES models are represented using black
lines with filled symbols and red lines with empty symbols, respectively.

36
 averaging over a large number of time steps, which can be time consuming espe-
cially at large values of δ. The time required to complete the DSMC simulations
420 in this paper is about 20 hours using an MPI parallel code which runs on 32
processor cores.
In the case of the kinetic solver, we estimate the computational efficiency by
considering simulations on a single core of an i7-4790K processor, running at a
frequency of 4.0 GHz. The simulation time is very short at δ = 10 and δ = 1
425 – of the order of one minute. This is because the quadrature order employed
can be very small. At δ = 0.1, the quadrature must be increased, leading to
computational times of the order of 5–10 minutes on a single processor core.
The exact figures are summarised in Table 2. For completeness, in this table
we included also the simulation times required to reach the stationary states
430 when using the hybrid approach at δ = 100 and 1000, corresponding to the
hydrodynamics regime. It can be seen that the hybrid approach described in this
section becomes inefficient when δ increases. This happens because, for δ & 10,
the number of iterations required to reach the steady state increases dramatically
with δ, being around two orders of magnitude larger at δ = 1000 than at δ = 10.
435 It is noteworthy that the projection method discussed in Appendix A performs
better at larger values of δ, as can be seen in Table A.3. In particular, the
computing times required by the projection method at δ = 1000 are shorter
than those required by the hybrid method by factors of about 4 and 5 for the S
and ES models, respectively. Further decreases in computational time at large
440 values of δ can be expected when the explicit time-stepping method employed in
this paper is replaced by, e.g., the implicit-explicit (IMEX) method that treats
the collision term implicitly [98, 99, 100], or the iterative methods discussed in
Refs. [72, 73, 74, 75].
Before ending this section, we briefly mention the procedure employed to
judge the approach to steady state. We consider 10 batches of NT iterations
each, with a time step of δt. The time interval corresponding to one batch is
∆t = NT δt. Denoting via T n (x) the temperature profile after the nth batch, we
compute the L2 norm of the relative difference between two successive batches,

37
 periodic

Tright = 1 + ∆T /2
Tleft = 1 − ∆T /2

x = −1/2 periodic x = 1/2

Figure 2: The simulation setup for the heat transfer problem. The vertical dashed lines show
a sample grid employing S = 4 points on each half of the channel, stretched according to
Eq. (43) with A = 0.95.

as follows: 
Z e
!2 1/2
1 L/2
Ten+1 (e
x)
L2 [T n+1 ] =  de
x −1  , (107)
e
L e
−L/2 e n
T (ex)

where the integration is performed as indicated in Eqs. (105) and (106). Figure 1
445 shows that L2 [T ] steadily decreases with the time t (tn = n∆t).

5. Heat transfer

The first application considered in this paper concerns the heat transfer
between stationary parallel plates problem. The simulation setup is represented
schematically in Fig. 2. In our simulations, the reference temperature, Teref =
450 (Teleft + Teright )/2, is varied between 1 K and 3000 K.
Representative profiles of the density n and temperature T are shown for
3
He constituents at Teref = 100 K in Fig. 3. The DSMC results are shown
using solid lines. The FDLB results obtained with the S model are shown using
red dashed lines with empty symbols. The FDLB results obtained with the
455 ES model are shown using black dotted lines with filled symbols. The FDLB
data corresponding to δ = 10, 1 and 0.1 are shown with squares, circles and

38
 2.7 1.8
S:δ=10 S:δ=10
2.4 δ=1 δ=1
δ=0.1 1.5 δ=0.1
2.1 ES:δ=10 ES:δ=10
δ=1 δ=1
1.8 δ=0.1 1.2 δ=0.1
DSMC DSMC
1.5

T
n

1.2 0.9

0.9
0.6
0.6
~ ~
(3He,Tref=100 K) (a) (3He,Tref=100 K) (b)
0.3 0.3
-0.5 -0.25 0 0.25 0.5 -0.5 -0.25 0 0.25 0.5
x x

Figure 3: Comparison between the S (red dashed lines and empty symbols) and ES (black
dotted lines and filled symbols) results and the DSMC (continuous lines) results in the context
of the heat transfer between stationary plates problem for the profiles of (a) n and (b) T
through the channel (−1/2 ≤ x ≤ 1/2), for 3 He gas constituents. The reference temperature
is set to Teref = 100 K, while the temperature difference between the two walls is ∆T
g = 1.5Teref .

triangles, respectively. Very good agreement can be seen between the results
obtained using the ES model and the DSMC data. There is a visible discrepancy
in the temperature profile obtained with the Shakohv model at δ = 1.
460 A more quantitative analysis is performed at the level of the quantity Q,
introduced in Eq. (99), with u
ey set to 0. Figure 4 compares the FDLB and
DSMC results for Q with respect to Teref for 1 K ≤ Teref ≤ 3000 K, at δ = 10
(top line), 1 (middle line) and 0.1 (bottom line). The value of Q is represented
in the left column of Fig. 4, while the relative error QFDLB /QDSMC − 1 is shown
465 in the right column of Fig. 4. These results were obtained using the S (empty
symbols) and the ES (filled symbols) models, for both the 3 He (red lines with
squares) and the 4 He (black lines with circles) constituents. At δ = 10, the S
model overestimates the DSMC results. Contrary to the S model, these DSMC
results are underestimated by the ES model. The relative errors are roughly the
470 same in absolute values. At smaller values of δ, the ES model provides results
which are more accurate than those obtained using the S model. The highest

39
 0.13 1.6
(δ = 10)
1.2
0.128
0.8

[QFDLB / QDSMC - 1] (%)

0.126 0.4
S:3He S:4He
0 ES:3He ES:4He
Q

0.124
-0.4

0.122 -0.8
3 4
S: He S: He
3 4 -1.2
0.12 ES: He ES: He
(a) DSMC DSMC (d) (δ = 10)
-1.6
1 10 100 1000 1 10 100 1000 10000
~ ~
Tref(K) Tref(K)

5
0.325
(δ = 1) S:3He S:4He
3 4
4 ES: He ES: He
0.32
[QFDLB / QDSMC - 1] (%)
3
0.315
2
Q

0.31
1

0.305 3 4
0
S: He S: He
3 4
ES: He ES: He
DSMC DSMC (e) (δ = 1)
(b)
0.3 -1
1 10 100 1000 1 10 100 1000 10000
~ ~
Tref(K) Tref(K)

1.8
0.4 (δ = 0.1)
3 4
S: He S: He
1.6 ES:3He ES:4He
0.398 1.4
[QFDLB / QDSMC - 1] (%)

1.2
0.396
Q

0.8
0.394
0.6
3 4
S: He S: He
3 4 0.4
0.392 ES: He ES: He
(c) DSMC DSMC (f) (δ = 0.1)
0.2
1 10 100 1000 1 10 100 1000 10000
~ ~
Tref(K) Tref(K)

Figure 4: (Left) Dependence of the constant Q, computed for the heat transfer between
ey = 0, on the average wall temperature Teref .
stationary plates problem using Eq. (99) with u
(Right) Relative error QFDLB /QDSMC − 1 of the FDLB results with respect to the DSMC
results. Both 3 He (red dashed lines with squares) and 4 He (black dotted lines with circles)
are considered within the S (empty symbols) and ES (filled symbols) models and the results
are represented at δ = 10 (top), 1 (middle) and 0.1 (bottom).
40
 periodic

Bounce-back
uw

x=0 periodic x = 1/2

Figure 5: The simulation setup for the Couette flow problem. The vertical dashed lines show
a sample grid employing S = 8 points, stretched according to Eq. (43) with A = 0.95.

relative discrepancy with respect to the DSMC data can be observed at δ = 1,

when the relative error of the S model reaches almost 5%, while for the ES
model, it stays below 3%.

475 6. Couette flow

The second application concerns the Couette flow between parallel plates.
Due to the symmetry of the flow, only the right half of the channel (0 ≤ x ≤ 1/2)
is considered in the simulation setup, as shown in Fig. 5. The walls are kept at
constant temperatures Teleft = Teright = Teref and Teref is varied bewteen 1 K and
q √
480 3000 K. The wall velocity u ew = 2K e B Teref /m
e takes the value uw = 2 after
non-dimensionalization.
Aside from the transversal component qx of the heat flux, which can be
related at large δ to the temperature variations with respect to the coordinate
x via Fourier’s law, qx = −κ∂x T , the Couette flow exhibits a non-vanishing
485 longitudinal heat flux, qy , which is a purely microfluidics effect. Figure 6 shows
a comparison between the FDLB results for the S (dashed red lines and empty
symbols) and ES (dotted black lines and filled symbols) models and the DSMC
results (solid purple lines). The wall temperature is set to Teref = 300 K and

41
 1.16 0.9
S: δ=10 S: δ=10 (b)
(a)
1.14 ES: δ=10 0.8 ES: δ=10
S: δ=1 S: δ=1
1.12
ES: δ=1 0.7 ES: δ=1
1.1 S:δ=0.1 S:δ=0.1
0.6 ES:δ=0.1
1.08 ES:δ=0.1
DSMC DSMC
0.5

uy / uw
1.06 4 4
( He,T=300K) ( He,T=300K)
n

1.04 0.4
1.02 0.3
1
0.2
0.98
0.96 0.1

0.94 0
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5
x x

0
1.7 (4He,T=300K) (c) (d)

1.6 -0.1

1.5 -0.2
4
( He,T=300K)
qy / uw

-0.3
T

1.4

1.3 -0.4

(S model) (ES model) (S model) (ES model)

1.2 δ=10 δ=10 -0.5 δ=10 δ=10
δ=1 δ=1 δ=1 δ=1
δ=0.1 δ=0.1 δ=0.1 δ=0.1
1.1 -0.6
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5
x x

Figure 6: Comparison between the FDLB results for the S model (dashed red lines and
empty symbols) and ES model (dotted black lines and filled symbols) and the DSMC results
(continuous lines) for the profiles of (a) n, (b) uy , (c) T and (d) qy through the half-channel
(0 ≤ x ≤ 1/2), for 4 He gas constituents, in the context of the Couette flow. The wall
√
temperature is set to Teref = 300 K, while the wall velocity is uw = 2.

42
 0.41
3
S: He
4
S: He (b) S: 3He S:4He
(a)
3
ES: He
4
ES: He ES: 3He ES:4He
0.105 0.4 DSMC DSMC
DSMC DSMC
0.39
0.1
0.38

Π
Π

0.095 0.37

0.36
0.09
0.35
(δ = 10) (δ = 1)
0.085 0.34
1 10 100 1000 1 10 100 1000
~ ~
Tref(K) Tref(K)

0.55 3 4
(c) S: He S: He
3 4
ES: He ES: He
DSMC DSMC
0.545

0.54
Π

0.535

0.53

(δ = 0.1)
0.525
1 10 100 1000
~
Tref(K)

Figure 7: Dependence of Π, computed using Eq. (98) in the context of the Couette flow, on
the wall temperature Teref for both 3 He (squares) and 4 He (circles), at (a) δ = 10, (b) δ = 1
and (c) δ = 0.1.

43
 0.096 3 4
0.21 3 4
(a) S: He S: He (b) S: He S: He
3 4
ES: He ES: He 0.208 ES: 3He ES:4He
DSMC DSMC DSMC DSMC
0.092 0.206
0.204
0.088 0.202
Qw

Qw
0.2
0.084 0.198
0.196
0.08 0.194
0.192
(δ = 10) (δ = 1)
0.076 0.19
1 10 100 1000 1 10 100 1000
~ ~
Tref(K) Tref(K)

0.12 3 4
(c) S: He S: He
ES: 3He ES:4He
0.11 DSMC DSMC

0.1
Qw

0.09

0.08

0.07

(δ = 0.1)
0.06
1 10 100 1000
~
Tref(K)

Figure 8: Dependence of Qw , computed in the context of the Couette flow using Eq. (103), on
the wall temperature Teref for both 3 He (squares) and 4 He (circles), at (a) δ = 10, (b) δ = 1
and (c) δ = 0.1.

44
 -0.014 3 4
-0.115 3 4
(a) S: He S: He (b) S: He S: He
3 4 3 4
ES: He ES: He -0.12 ES: He ES: He
-0.016 DSMC DSMC DSMC DSMC
-0.125
-0.018
-0.13
Qy

Qy
-0.02 -0.135

-0.14
-0.022
-0.145
-0.024
-0.15
(δ = 10) (δ = 1)
-0.026 -0.155
1 10 100 1000 1 10 100 1000
~ ~
Tref(K) Tref(K)

-0.055 3 4
(c) S: He S: He
3 4
-0.06 ES: He ES: He
DSMC DSMC
-0.065

-0.07
Qy

-0.075

-0.08

-0.085

-0.09
(δ = 0.1)
-0.095
1 10 100 1000
~
Tref(K)

Figure 9: Dependence of Qy , computed in the context of the Couette flow using Eq. (102), on
the wall temperature Teref for both 3 He (squares) and 4 He (circles), at (a) δ = 10, (b) δ = 1
and (c) δ = 0.1.

45
 2.5 5 (ES model)
4 3
He:δ=10 He:δ=10
δ=1 δ=1
2 4
δ=0.1 δ=0.1

[ΠFDLB / ΠDSMC - 1] (%)

[ΠFDLB / ΠDSMC - 1] (%)

1.5 4 3
3 He:δ=10 He:δ=10
δ=1 δ=1
1 δ=0.1 δ=0.1
2
0.5
1
0

0
-0.5
(S model) (a) (d)
-1 -1
1 10 100 1000 1 10 100 1000
~ ~
Tref(K) Tref(K)

5
0
0
[Qw;FDLB / Qw;DSMC - 1] (%)

[Qw;FDLB / Qw;DSMC - 1] (%)

-5
-5 4 3
He:δ=10 He:δ=10
δ=1 δ=1 -10
δ=0.1 δ=0.1 4 3
-10 -15 He:δ=10 He:δ=10
δ=1 δ=1
-20 δ=0.1 δ=0.1

-15 -25

-30

-20 -35
(S model) (b) (ES model) (e)
-40
1 10 100 1000 1 10 100 1000
~ ~
Tref(K) Tref(K)

20
20 4 3 4 3
He:δ=10 He:δ=10 He:δ=10 He:δ=10
δ=1 δ=1 δ=1 δ=1
δ=0.1 δ=0.1 10 δ=0.1 δ=0.1
[Qy;FDLB / Qy;DSMC - 1] (%)

[Qy;FDLB / Qy;DSMC - 1] (%)

10
0

0 -10

-20
-10

-30

-20 (S model) (c) (ES model) (f)

-40
1 10 100 1000 1 10 100 1000
~ ~
Tref(K) Tref(K)

Figure 10: Relative errors ΠFDLB /ΠDSMC − 1 (top), Qw;FDLB /Qw;DSMC − 1 (middle) and
Qy;FDLB/Qy;DSMC − 1 (bottom) between the DSMC and FDLB results for the S model (left)
and ES model (right), at δ = 10 (squares), 1 (circles) and 0.1 (triangles) for 1 K ≤ Teref ≤
3000 K, computed in the context of the Couette flow.
46
 4
He gas constituents are considered for δ = 10, 1 and 0.1. Both the S and
490 ES models are in good agreement with the DSMC data at δ = 10. When δ
decreases, the agreement deteriorates, being slightly worse in the case of the ES
model. Remarkably, the density profiles are well recovered with both models at
all tested values of δ.
We now consider a more quantitative analysis at the level of Π, Qw and
495 Qy , computed via Eqs. (96), (103) and (102), respectively. The variations with
the plate temperature Teref of Π, Qw and Qy for 3 He and 4 He are shown in
Figs. 7, 8 and 9 for (a) δ = 10, (b) δ = 1 and (c) δ = 0.1. Each plot shows
curves corresponding to the S model (dashed lines with empty symbols), ES
model (dotted lines with filled symbols) and DSMC (solid lines). The data
500 corresponding to 3 He is shown using red squares, while the data for 4 He is shown
with black circles. It can be seen that in general, the agreement between the
results obtained with the model equations and the DSMC results deteriorates as
δ is decreased. Contrary to the results obtained in the case of the heat transfer
problem, the S model gives more accurate results compared to the ES model,
505 confirming the results reported in Ref. [89]. Figure 10 shows the relative errors
computed with respect to the DSMC results, obtained with the S (left column)
and ES (right column) models. The results for 4 He are shown with solid lines
and filled symbols, while those for 3 He are shown with dashed lines and empty
symbols. The data corresponding to δ = 10, 1 and 0.1 are shown with red
510 squares, green circles and amber triangles, respectively. In the case of Π, the
relative error of the ES model is roughly twice that of the S model.
It is remarkable that the relative errors for both Qw and Qy (shown in
Figs. 8 and 9) reach values around 20% for δ = 0.1. This can be explained since
the heat fluxes decrease to 0 as δ is decreased, while Π, for which the relative
515 error is below 5%, attains a finite value as the ballistic regime is approached
(limδ→0 Π = π −1/2 ). Thus, the relative errors for Qw and Qy are computed by
dividing the FDLB values by small numbers. However, in the case of Qy , the
errors are around 20% even when δ = 10, whereas for both Qw and Π, the error
at δ = 10 is less than 1%. This disagreement between the model equations and

47
 periodic

Tright = 1 + ∆T /2
Tleft = 1 − ∆T /2
−uw uw

x = −1/2 periodic x = 1/2

Figure 11: The simulation setup for the heat transfer under shear problem. The vertical
dashed lines show a sample grid employing S = 4 points on each half of the channel, stretched
according to Eq. (43) with A = 0.95.

520 the DSMC data can be attributed to the nature of Qy . Since the longitudinal
heat flux, qy , is not generated by a temperature gradient (through the so-called
direct phenomenon), its characteristics must depend on higher order transport
coefficients, which are visible only at the Burnett level [101]. Since the model
equations are constructed to ensure consistency only at the Navier-Stokes level
525 (corresponding to the first order in the Chapman-Enskog expansion), it is not
surprising that such cross phenomena are not accurately recovered.

7. Heat transfer under shear

The final example considered in this paper is the heat transfer between
parallel plates in motion. The simulation setup is represented in Fig. 11. This
530 example combines the features of the heat transfer between stationary plates dis-
cussed in Sec. 5 and those of the Couette flow discussed in Sec. 6. The reference
temperature Teref = (Teleft + Teright )/2, is varied between 1 K and 3000 K for 3 He
and 4 He constituents, while for Ne, the range for Teref is 20 K ≤ Teref ≤ 5000 K.
g = Teright − Teleft obeys Eq. (100).
As in Sec. 5, the temperature difference ∆T
535 e left = −e
Furthermore, the plates have velocities u e right = u
uw j and u ew j, where
q
ew = 2K
u e B Teref /m,
e such that the Mach number is given by Eq. (97).

48
 3 ~ 1
(Ne,Tref=300 K)
S:δ=10 0.8 S:δ=10
δ=1 δ=1
2.5 δ=0.1 0.6 δ=0.1
ES:δ=10 ES:δ=10
δ=1 0.4 δ=1
2 δ=0.1 0.2 δ=0.1
DSMC DSMC

uy / uw
0
n

1.5
-0.2
-0.4
1
-0.6
-0.8 ~
0.5 (a) (b) (Ne,Tref=300 K)
-1
-0.5 -0.25 0 0.25 0.5 -0.5 -0.25 0 0.25 0.5
x x

1.8 ~
(Ne,Tref=300 K)
1.6

1.4

1.2
T

1 S:δ=10
δ=1
0.8 δ=0.1
ES:δ=10
δ=1
0.6
δ=0.1
DSMC
0.4 (c)
-0.5 -0.25 0 0.25 0.5
x

Figure 12: Comparison between the FDLB results (dotted lines and points) obtained using
the S (red empty symbols) and ES (black filled symbols) models and the DSMC (continuous
lines) results for the profiles of n (a), uy (b) and T (c) through the channel (−1/2 ≤ x ≤ 1/2),
for Ne gas constituents, in the context of the heat transfer between moving plates problem.
The reference temperature is set to Teref = 300 K, the temperature difference between the two
q
g = 1.5Teref and the wall velocity is u
walls is ∆T ew = 2K e B Teref /m.
e

49
 0.1 0.165
S:3He
4
S: He
3 4 S:3He
ES: He ES: He
ES:3He
DSMC DSMC 0.16 DSMC
0.095 4
S: Ne S: He
ES: Ne ES:4He
0.155 DSMC
DSMC

Q
0.09
Π

0.15

0.145 S: Ne
0.085 ES: Ne
DSMC
0.14 (f) (ES-BGK) (δ = 10)
(δ = 10)
0.08 1 10 100 1000 10000
~ ~
Tref(K) Tref(K)

(δ = 1) 3
S: He
0.38 3
0.31 ES: He
DSMC

0.37
0.3
3 4 S: Ne
S: He S: He
3 4
Q
Π

ES: He 0.36 ES: He ES: Ne

0.29 DSMC DSMC DSMC

0.35
0.28 S:4He S: Ne
ES:4He ES: Ne
DSMC DSMC (δ = 1)
0.34
1 10 100 1000 10000 1 10 100 1000 10000
~ ~
Tref(K) Tref(K)

0.404 0.416

0.402
0.414
0.4 3 4 S:3He S:4He
S: He S: He
3
ES: He
4
ES: He 0.412 ES:3He ES:4He
0.398 DSMC DSMC DSMC DSMC
Π

S: Ne 0.41 S: Ne
0.396 ES: Ne
ES: Ne
DSMC DSMC
0.394 0.408

0.392
0.406
(δ = 0.1) (δ = 0.1)
0.39
1 10 100 1000 10000 1 10 100 1000 10000
~ ~
Tref(K) Tref(K)

Figure 13: Dependence of Π (left column) and Q (right column), defined in Eqs. (96) and
(99) for the heat transfer between moving plates problem, on the average wall temperature
Teref = (Teleft + Teright )/2 for 3 He (red squares), 4 He (black circles) and Ne (blue triangles), at
δ = 10 (top line), 1 (middle line) and 0.1 (bottom line).

50
 Figure 12 shows the profiles of the density (a), velocity (b) and temperature
(c) for the case of Ne constituents at T = 300 K. In general, good agreement
can be seen between the results corresponding to the model equations and the
540 DSMC results. A larger discrepancy can be seen between the ES model and the
DSMC results, especially in the temperature profile at δ = 1 and 0.1.
A quantitative analysis can be made at the level of the nondimensional quan-
tities Π and Q, computed using Eqs. (96) and (99). Figure 13 shows a compari-
son between the FDLB results for the S (dashed lines with empty symbols) and
545 the ES (dotted lines with filled symbols) models and the DSMC results (solid
lines), obtained for 3 He (squares), 4 He (circles), and Ne (triangles) constituents.
Figure 14 shows the relative errors in Q (dashed lines and empty symbols) and
Π (dotted lines and filled symbols) computed for the S model (left column) and
ES model (right column) with respect to the DSMC results for 3 He (squares),
4
550 He (circles) and Ne (triangles). At δ = 10 (top line), the results obtained using
the ES model seem to be in better agreement with the DSMC results than those
obtained using the S model. At δ = 1 (middle line) and 0.1 (bottom line), the
two models give results with similar accuracy. As noticed in the case of the heat
transfer between stationary plates and in the case of the direct phenomena in
555 the Couette flow, the relative erros are highest at δ = 1, where they take values
between 6 − 8% (about 1% higher for Q than for Π).

8. Conclusion

In this paper, we presented a systematic comparison between the results

obtained using the Boltzmann equation with the Shakhov (S) and Ellipsoidal-
560 BGK (ES) models for the collision term and those obtained using the direct
simulation Monte Carlo (DSMC) method for three benchmark channel flows
between parallel plates, namely: heat transfer between static walls, Couette
flow and heat transfer under shear. The results were obtained numerically
in the nonlinear regime [Ma ≃ 2.19 for the case when the parallel plates are
565 moving and 2(Teright − Teleft)/(Teright + Teleft) = 1.5 for the heat transfer problems],

51
 1
3
He
2 0.8 4
Q: He
Ne
0.6
1.6
Relative error (%)

Relative error (%)

0.4
3 3
1.2 He He
3
Q:4He Π:4He 0.2 He
Ne Ne Π:4He
0.8 0 Ne

-0.2
0.4
(a) (S-model) (δ = 10) -0.4 (d) (ES-BGK) (δ = 10)
1 10 100 1000 10000 1 10 100 1000 10000
~ ~
Tref(K) Tref(K)

7 7

6 6
Relative error (%)

Relative error (%)

5 5

Q: 3He Π: 3He Q: 3He Π: 3He

4 4 4 4
4 He He 4 He He
Ne Ne Ne Ne

3 (b) (S-model) (δ = 1) 3 (e) (ES-BGK) (δ = 1)

1 10 100 1000 10000 1 10 100 1000 10000
~ ~
Tref(K) Tref(K)

2.6
2.6
2.4
2.4
2.2
2.2
Relative error (%)

Relative error (%)

2
2
1.8
1.8
1.6
1.6
1.4
3 3 3 3
1.2 Q: He Π: He 1.4 Q: He Π: He
4 4 4 4
He He He He
1 Ne Ne 1.2 Ne Ne
(c) (S-model) (δ = 0.1) (f) (ES-BGK) (δ = 0.1)
0.8 1
1 10 100 1000 10000 1 10 100 1000 10000
~ ~
Tref(K) Tref(K)

Figure 14: Dependence of the relative errors QFDLB /QDSMC − 1 (dashed lines and empty
symbols) and ΠFDLB /ΠDSMC − 1 (dotted lines and filled symbols), expressed in percentages,
where the FDLB results are obtained using the S (left column) and ES (right column) models,
for the heat transfer between moving plates problem, on the average wall temperature Tref =
(Tleft + Tright )/2 for 3 He (squares), 4 He (circles) and Ne (triangles), at δ = 10 (top line), 1
(middle line) and 0.1 (bottom line).
52
 by considering 3 He and 4 He constituents interacting via ab initio potentials. We
also considered Ne constituents for the heat transfer under shear problem.
In the kinetic theory setup, the connection with the DSMC simulations was
established at the level of the transport coefficients (dynamic viscosity µ
e and
570 κ). For 3 He and 4 He, the range of values for the reference
heat conductivity e
temperature Teref = (Teright + Teleft )/2 was 1 K ≤ Teref ≤ 3000 K, while for the
Ne constituents, it was 20 K ≤ Teref ≤ 5000 K. We considered three values for
the rarefaction parameter, namely δ = 10 (slip flow regime), δ = 1 (transition
regime) and δ = 0.1 (early free molecular flow regime).
575 We first conducted a qualitative comparison at the level of the profiles of
the density, temperature, velocity and heat flux. In all cases considered, the
density profile was well recovered with both kinetic models, for all values of the
rarefation parameter. In the context of the heat transfer problem, the results
obtained using the ES model were in better agreement with the DSMC results for
580 the temperature profile. In the Couette and heat transfer with shear problems,
the S model seemed to give results which were closer to the DSMC predictions
for all quantities (temperature, velocity and heat flux).
We next considered a quantitative comparison of the performance of the
kinetic models with respect to the DSMC data by comparing the numerical
585 values for non-dimensional quantities derived from the longitudinal heat flux
(in the case of heat transfer between stationary and moving plates, denoted Q),
shear stress (in the case of Couette flow and heat transfer between moving plates,
denoted Π), as well as the half-channel heat flow rate, Qy , and heat transfer
rate through the boundary, Qw (in the case of the Couette flow). Among these
590 quantities, we can distinguish two categories. The first category (containing Q,
Π and Qw ) refers to quantities related to “direct phenomena,” which are driven
by, e.g., shear rate ∂x uy for Π and temperature gradient ∂x T for Q, as predicted
by the Navier-Stokes-Fourier theory. The second category (containing Qy ) refers
to quantities related to “cross phenomena,” visible at the level of the Burnett
595 equations, in which the usual thermodynamic forces driving the non-equilibrium
quantity are absent (i.e., non-vanishing qy when ∂y T = 0).

53
 For the quantities in the first category (corresponding to direct phenomena),
the agreement between the kinetic models and the DSMC results was within a
few percent at δ = 10, which confirms the validity of these models in the slip
600 flow regime. At δ = 1, the errors seem to be bounded within 8% for both
models, with the ES model giving better results in the heat transfer between
stationary plates problem, while the S model performs better for the Couette
flow and heat transfer under shear problems. When δ = 0.1, the free molecular
flow regime is approached. For the quantities that attain a finite value in this
605 regime (Q in the heat transfer problems and Π in the Couette flow problem),
the relative errors drop compared to δ = 1, to within 2% − 3%. On the contrary,
the relative errors for the heat flux Qw measured at the wall in the Couette flow
grow to around 20% for the S model and 30% for the ES model. This can be
attributed to the fact that Qw decreases towards 0 as the free molecular flow
610 regime is approached, such that the relative errors are computed by dividing
the results obtained within the model equations by a small quantity.
When considering the quantity Qy from the second category, which is gener-
ated through the cross-phenomena, the results of the kinetic models had relative
errors of the order of 20% even at δ = 10, highlighting that the model equations
615 do not accurately take into account such phenomena. At δ = 1, the relative
errors decrease to around 10% for the S model and 15% for the ES model. At
δ = 0.1, they increase again to around 20% and 35% for the S and ES models,
respectively. As was the case for Qw , the large values of the relative errors of Qy
encountered at δ = 10 and δ = 0.1 may be caused by the fact that Qy vanishes
620 in the inviscid (δ → ∞) and free molecular flow (δ → 0) regimes.
In conclusion, our results demonstrate that even in the strongly non-linear
regime, the model equations can give reasonably accurate results, with errors of
up to 10% for quantities related to direct phenomena throughout the rarefaction
spectrum (provided they remain finite in the free molecular flow regime), while
625 the errors for the cross phenomena-related quantities seem to be within 35%.
Due to the computational efficiency of the finite difference lattice Boltzmann
(FDLB) algorithm employed in this paper, solving the kinetic model equations

54
 can provide a cheap and reasonably accurate solution for the flow properties in
the case of realistic monatomic gases under rarefied conditions.
630 Acknowledgments. VEA gratefully acknowledges the generous support of
the Romanian-U.S. Fulbright Commission through The Fulbright Senior Post-
doctoral Program for Visiting Scholars 2017-2018, Grant number 678/2018. FS
acknowledges the Brazilian Agency CNPq, Brazil, for the support of his research,
grant 304831/2018-2. VEA is grateful to Professor P. Dellar (Oxford University,
635 UK) for preliminary discussions regarding the projection of the Shakhov collision
term onto orthogonal polynomials. The numerical simulations were performed
on the Turing High Performance Computing cluster and the Computing cluster
at the Computer Science Department of the Old Dominion University (Norfolk,
VA, USA). The computer simulations reported in this paper were done using the
640 Portable Extensible Toolkit for Scientific Computation (PETSc 3.6) developed
at Argonne National Laboratory, Argonne, Illinois [102, 103].

Appendix A. FDLB models for the hydrodynamic regime

When δ & 100, the flow enters the hydrodynamic regime, where it can be
approximately described by the Navier-Stokes equations [104],

Dρ
+ ρ∇ · u =0,
Dt
Dui
ρ + ∂i P = − ∂j σij ,
Dt
De
ρ + P ∇ · u = − σij ∂i uj − ∇ · q, (A.1)
Dt

where D/Dt = ∂t + u ·∇ is the convective (material) derivative, σij = Tij − P δij

3
is the shear stress tensor, q is the heat flux, e = 2m T is the specific energy for
an ideal monatomic gas. The Newtonian fluid model and Fourier’s law give the
following constitutive equations for τij and q:
 
2
τij = −µ ∂i uj + ∂j ui − δij ∇ · u , q = −κ∇T, (A.2)
3

where µ and κ are the dynamic viscosity and heat conductivity, respectively.

55
 According to the Chapman-Enskog expansion, briefly discussed in Appendix A.1,
645 low order moments of the reduced distribution functions φ and χ are required to
ensure the relations in Eq. (A.2). The evolution and stationary state properties
of these moments can be obtained by employing similarly low order quadratures
(i.e., Qφx = 5 and Qφy = 4 for the φ distribution; and Qχx = 3 and Qχy = 2 for the
χ distribution).
650 As the quadrature order is lowered and δ is increased, the recovery of the
conservation equations becomes increasingly challenging when the distributions
are evaluated directly (i.e., using the hybrid method described in Sec. 3.3). In
the traditional lattice Boltzmann framework, the key to employing the low order
quadratures is to project the local equilibrium distribution on a set of orthogonal
655 polynomials, which is subsequently truncated at an order Nx∗ (∗ ∈ {φ, χ}).

Appendix A.1. Chapman-Enskog analysis

To derive the hydrodynamic regime from the kinetic model equation,

∂f p 1
+ · ∇f = − (f − f∗ ), (A.3)
∂t m τ∗

the fluid can be assumed to be very close to isotropic thermal equilibrium,

described by f = f∗ = f (eq) , where f (eq) is the Maxwell-Boltzmann distribution.
The deviations of f and f∗ from f (eq) , denoted by δf = f − f (eq) and δf∗ =
f∗ − f (eq) , can be assumed to be of the same order as the relaxation time τ∗ ,
which is considered to be small. To first order with respect to τ∗ , the deviation
δf can be written as
 
∂f (eq) p (eq)
δf = δf∗ − τ∗ + · ∇f , (A.4)
∂t m

where f (eq) is determined by the density ρ, velocity u and temperature T .

The constitutive relations for σij = Tij − P δij and qi given in Eq. (A.2) can
be obtained by taking the second and third order moments of Eq. (A.4) with

56
respect to the momentum space:
 Z Z 
∗ ∂ 3 pi pj (eq) 3 pi pj pk (eq)
σij =σij − τ∗ d p f + ∂k d p f ,
∂t m m2
 Z Z 
∂ p2 pi (eq) 2
3 p pi pk (eq)
qi + uj σij =qi∗ + uj σij
∗
− τ∗ d3 p f + ∂k d p f ,
∂t 2m2 2m3
(A.5)
∗
where σij and qi∗ are obtained by taking moments of δf∗ ,
Z Z
3 pi pj ∗ p2 pi
d p δf∗ = σij , d3 p δf∗ = qi∗ + σij
∗
uj . (A.6)
m 2m2

For completeness, the details of the Chapman-Enskog procedure for the S and
ES models employed in this paper are briefly presented. The integrals of f (eq)
entering Eq. (A.5) are
Z
pi pj (eq)
d3 p f =P δij + ρui uj ,
m
Z
pi pj pk (eq)
d3 p f =P (δij uk + δjk ui + δki uj ) + ρui uj uk ,
m2
Z
p2 pi (eq) 5P ρu2
d3 p f = ui + ui ,
2m2 2 2
Z    
p2 pi pk (eq) 5T u2 7P ρu2
d3 p f =P + δik + + ui uk . (A.7)
2m3 2m 2 2 2

Using the above relations and replacing the time derivatives using the Euler
(inviscid) form of the Navier-Stokes equations (A.1),

Dρ Du De
+ ρ∇ · u = O(τ∗ ), ρ + ∇P = O(τ∗ ), ρ + P ∇ · u = O(τ∗ ),
Dt Dt Dt
(A.8)
it can be seen that
 
∗ 2 5τ∗ P
σij =σij − τ∗ P ∂i uj + ∂j ui − ∇ · uδij , q =q∗ − ∇T. (A.9)
3 2m

In the Shakhov (S) and ellipsoidal (ES) models, when f∗ are given by fS and
∗
fES introduced in Eqs. (4) and (7), σij and q∗ are given by:

S
σij =0, qS =(1 − Pr)q,
ES 1 − Pr
σij =− σij , qES =0. (A.10)
Pr

57
 Substituting Eq. (A.10) into Eq. (A.9) and comparing the result to Eq. (A.2),
the relations given in Eqs. (6) and (10) between the relaxation times, τS and
τES , and the transport coefficients µ and κ can be obtained.
660 The recovery of the constitutive equations (A.2) is conditioned by the correct
recovery of the integrals (A.6) and (A.7) of f∗ and f (eq) .
In the case when the flow is trivial along the z direction (d = 2), the pz
degree of freedom can be integrated automatically and Eqs. (A.7) and (A.6) can
be written in terms of the reduced distributions φ and χ introduced in Eq. (14).
The highest moments required for the reduced equilibrium distributions φ∗ and
χ∗ are derived from the fourth order moment on the last line of Eq. (A.7).
Substituting f (eq) → f∗ , the following equation is obtained:
Z Z " #
2 2 2
p p i p k (p x + p y )p i p k p i p k
d3 p f∗ = dpx dpy φ∗ + χ∗ . (A.11)
2m3 2m3 2m2

It can be seen that the above integrals require the correct recovery of the mo-
ments with respect to px of order 4 for φ∗ and of order 2 for χ∗ . This can be
achieved using the half-range Gauss-Hermite quadrature employing Qφx = 5 and
665 Qχx = 3 points on each of the px > 0 and px < 0 semiaxes. Furthermore, φ∗ and
χ∗ must be replaced by a truncated expansion with respect to the half-range
Hermite polynomials of orders Q∗x − 1. More details for each of the collision
models are given below.

Appendix A.2. Maxwell-Boltzmann distribution

670 The projection of the Maxwell-Boltzmann distribution f (eq) with respect

to the half-range Gauss-Hermite polynomials was derived in Ref. [44]. Here,
we only summarise the details. Considering the factorisation f (eq) = ngx gy gz
introduced in Eq. (2), the expansion of f (eq) can be performed at the level
of each gi factor individually. Specifically, the gz factor is integrated out when
675 introducing the reduced distributions. The gy factor is expanded up to order Ny∗
(∗ ∈ {φ, χ}) with respect to the full-range Hermite polynomials, as summarised
in Eqs. (80)–(81) and (90)–(91) in the contexts of the ES and S models. For the

58
 Maxwell-Boltzmann distribution, Eqs. (80) and (81) can be used by substituting
ζy and Ty by uy and T , respectively.
The gx factor is expanded with respect to the half-range Hermite polynomials
up to orders Nxφ = Qφx − 1 = 4 and Nxχ = Qχx − 1 = 2:
N∗
(Nx∗ ) ω(px ) Xx

g (px , ux , T ) = hr (|px |)[θ(px )Gr+ (ux , T ) + θ(−px )Gr− (ux , T )],

p0,x r=0
(A.12)
where θ(x) is the Heaviside step function. The coefficients Gr± are given by

r
!s/2  
1X mT 2 −ζ 2 ∗
Gr± (ux , T ) = (±1)s hr,s (1 ± erf ζ)Ps+ (ζ) ± √ e Ps (ζ) ,
2 s=0 2p20,x π
(A.13)
p
where ζ ≡ ζ(ux , T ) = ux m/2T and hr,s represent the coefficients of xs in the
expression for hr (x), as shown in Eq. (60). The polynomial Ps∗ (ζ),
s−1  
X s
Ps∗ (ζ) = Pj+ (ζ)Ps−j−1
−
(ζ), (A.14)
j=0
j

is defined with the help of the polynomials Ps± (ζ), which satisfy:

2 ds ±ζ 2
Ps± (ζ) = e∓ζ e . (A.15)
dζ s

680 Appendix A.3. d = 1: ES model

The mass equilibrium distribution φES , given by Eq. (70), can be rewritten
in the language of the previous subsection in the following form:

φES = ng(px , ux , T Bxx), (A.16)

such that its truncated version is

Nφ
(N φ ) ω(px ) Xx

φESx = hr (|px |)[θ(px )Φ+ −

ES;r (ux , T ) + θ(−px )ΦES;r (ux , T )],
p0,x r=0

Φ± ±
ES;r (ux , T ) =nGr (ux , T Bxx ), (A.17)

where Bxx is introduced in Eq. (71). The expansion coefficients Φ±

ES;r (ux , T )

can be written entirely in terms of the expansion coefficients Gr± corresponding

59
to the Maxwell-Boltzmann distribution. The energy equilibrium distribution
satisfying χES = 2TredφES admits a similar decomposition:
Nχ
(N χ ) ω(px ) Xx
+ −
χESx = hr (|px |)[θ(px )XES;r (ux , T ) + θ(−px )XES;r (ux , T )],
p0,x r=0

XES;r (ux , T ) =2Pred Gr± (ux , T Bxx). (A.18)

where Pred = nTred = 32 P − 21 Txx .

Appendix A.4. d = 1: S model

In the S model, the distributions φS and χS , given in Eq. (73), can be
expanded as:
      
Nx∗ + −
φS ω(p ) X Φ S;r Φ S;r
 = x
hr (|px |) θ(px )   + θ(−px )   ,
χS p0,x r=0 +
XS;r −
XS;r
 
1 − Pr φ;±
Φ±
S;r (u x , T ) =n Gr
±
(u x , T ) + G (u x , T ) ,
5nT 2 S;r
 
± 1 − Pr χ;±
XS;r (ux , T ) =2nT Gr± (ux , T ) + G (u x , T ) , (A.19)
5nT 2 S;r
∗;±
where Gr± (ux , T ) is given in Eq. (A.13), while GS;r (ux , T ) can be found as fol-
lows:
    
φ;± Z ∞ 2
GS;r (±p x − mu x ) 3
  =qx dpx gx (±px , ux , T )  −  
χ;±
GS;r 0 mT 1

× (±px − mux )h(px ). (A.20)

∗;±
Employing the expansion for gx (px , ux , T ) in Eq. (A.12), the coefficients GS;r (ux , T )
can be obtained as follows:
∞
X Z ∞ 3
X
∗;±
GS;r (ux , T ) = qx Gr±′ (ux , T ) dpx ω(px )hr′ (px ) (±1)k A∗k [pkx h(px )],
r ′ =0 0 k=0
(A.21)
where the coefficients A∗k are given by
        
Aφ0 3 mu 2 Aφ1mu 2
p
x 0,x
3
  =   − x
mux ,   =3 −   p0,x ,
χ
A0 1 T Aχ T 1
1

3p20,x ux p30,x
Aφ2 =Aχ2 = − , Aφ3 =Aχ3 = . (A.22)
T mT

60
 In order to evaluate Eq. (A.21) using the orthogonality relation for the half-
range Hermite polynomials [44],
Z ∞
dpx ω(px )hr (px )hr′ (px ) = δr,r′ , (A.23)
0

the following recurrence relation can be employed to eliminate the factors of px :

1 br cr
px hr (px ) = hr+1 (px ) − hr (px ) − hr−1 (px ), (A.24)
ar ar ar

where the recurrence coefficients ar , br and cr can be obtained by the procedure

described in Ref. [44]. There is an easy relation allowing cr to be eliminated in
favour of ar :
ar
cr = − . (A.25)
ar−1
The recurrence in Eq. (A.24) can be used to obtain the following relation:
n
X (n)
pnx hr (px ) = Br,r′ hr′ (px ), (A.26)
r ′ =−n

(n) (n)
where it is understood that Br,r′ = 0 when r + r′ < 0. The coefficients Br,r′
depend only on the properties of the half-range Hermite polynomials and can
685 thus be computed automatically at runtime. Their explicit values are given for
0 ≤ n ≤ 3 at the end of this subsection.
After applying the recurrence relations to eliminate all factors of px , Eq. (A.21)
becomes
3
X i
X
∗;± (i) ±
GS;r (ux , T ) = qx (±1)i A∗i Br;r′ Gr+r ′ (ux , T ). (A.27)
i=0 r ′ =−i

(i)
The coefficients Br;r′ entering the above expression can be computed as follows.
For i = 0, we have
(0)
Br;0 = 1. (A.28)

At i = 1, there are three non-vanishing coefficients:

(1) 1 (1) br (1) 1

Br;−1 = , Br;0 = − , Br;1 = . (A.29)
ar−1 ar ar

61
 (i) (1)
We remind the reader that Br;r′ = 0 whenever r + r′ < 0, e.g. B0;−1 = 0. At
i = 2, we find
 
(2) 1 (2) 1 br br−1 (2) 1 + b2r 1
Br;−2 = , Br;−1 = − + , Br;0 = + 2 ,
ar−2 ar−1 ar−1 ar ar−1 a2r ar−1
 
(2) 1 br+1 br (2) 1
Br;1 =− + , Br;2 = . (A.30)
ar ar+1 ar ar ar+1
(3)
Finally, Br;r′ is given by
 
(3) 1 (3) 1 br−2 br−1 br
Br;−3 = , =− Br;−2 + + ,
ar−3 ar−2 ar−1 ar−2 ar−1 ar−2 ar−1 ar
 
(3) 1 1 1 + b2r−1 br−1 br 1 + b2r
Br;−1 = + + + ,
ar−1 a2r−2 a2r−1 ar−1 ar a2r
 
(3) br−1 2br br (2 + b2r ) br+1
Br;0 = − 3 + 2 + 3
+ 2 ,
ar−1 ar−1 ar ar ar ar+1
 
(3) 1 1 1 + b2r br br+1 1 + b2r+1
Br;1 = + + + ,
ar a2r−1 a2r ar ar+1 a2r+1
 
(3) 1 br br+1 br+2 (3) 1
Br;2 =− + + , Br;3 = . (A.31)
ar ar+1 ar ar+1 ar+2 ar ar+1 ar+2

Appendix A.5. d = 2: ES model

When d = 2 and φES is given by Eq. (77), we seek the expansion coefficients
±,H ±,H
ΦES;r,ℓ and XES;r,ℓ defined through
φ
Nφ N
(N φ ,N φ ) ω(px )ω(py ) Xx Xy
1
φESx y = hr (|px |)Hℓ (py )[θ(px )Φ+,H −,H
ES;r,ℓ + θ(−px )ΦES;r,ℓ ],
p0,x p0,y r=0 ℓ!
ℓ=0
χ
Nxχ Ny
(N χ ,Nyχ ) ω(px )ω(py ) X X 1 +,H −,H
χESx = hr (|px |)Hℓ (py )[θ(px )XES;r,ℓ + θ(−px )XES;r,ℓ ].
p0,x p0,y r=0 ℓ!
ℓ=0

(A.32)

Due to the relation χES = TredφES , where Tred = Pred /n and Pred is defined in
±,H
Eq. (78), the coefficients XES;r,ℓ and Φ±,H
ES;r,ℓ can be related via

±,H
XES;r,ℓ = TredΦ±,H
ES;r,ℓ . (A.33)

Inverting Eq. (A.32) and using Eq. (77) to replace φES , we find
∞
X Z ∞
±,H ± ω(px ) H ±
ΦES;r,ℓ = n Gr′ (ux , T Bxx) dpx G (ζ , Ty )hr (px )hr′ (px ), (A.34)
0 p0,x ℓ y
′ r =0

62
where the expansions in Eqs. (A.12) and (80) were used to replace the functions
g(px , ux , T Bxx) and g(py , ζy , Ty ), respectively. The summation range for r′ was
extended to ∞ to allow coefficients Gr±′ of orders r′ > Nxφ to be taken into
account. The superscript ± in ζy± indicates the sign of px in Eq. (79), i.e.
Bxy p0,x Bxy
ζy± = uy − ux ± p . (A.35)
Bxx mBxx x
The coefficients GℓH (ζy± , Ty ) are given explicitly for 0 ≤ ℓ ≤ 3 in Eq. (81). For
larger values of ℓ, the following formula can be employed [44]:
⌊ℓ/2⌋
X ℓ!
GℓH (ζy± , Ty ) = Uℓ−2s (ζy± )Is (Ty ). (A.36)
s=0
2s s!(ℓ − 2s)!

It can be seen that GℓH (ζy± , Ty ) is a polynomial of order ℓ with respect to ζy± ,
and therefore with respect to px . Thus, it can be expanded as follows:
ℓ
X
GℓH (ζy± , Ty ) = H k
(±1)k Gℓ;k px . (A.37)
k=0

Using now the recurrsion relation for the half-range Hermite polynomials given
in Eq. (A.26), it can be shown that Φ±,H
ES;r,ℓ reduces to

ℓ
X k
X (k)
Φ±,H
ES;r,ℓ = n (±1)k Gℓ;k
H ±
Br,r′ Gr+r ′ (ux , T Bxx ). (A.38)
k=0 r=−k
H
The first few functions Gℓ;k can be obtained by inspection:
H H H
G0;0 = 1, G1;0 = Uy , G1;1 = B,
H
G2;0 = Uy2 + I(Ty ), H
G2;1 = 2BUy , H
G2;2 = B2,
H
G3;0 = Uy3 + 3Uy I(Ty ), H
G3;1 = 3B[Uy2 + I(Ty )],
H
G3;2 = 3B 2 Uy , H
G3;3 = B3, (A.39)
mTy detB
where I(Ty ) = p20,y
− 1 and Ty = Bxx T are defined in Eqs. (82) and (79),
respectively, while Uy and B are introduced below:
 
m Bxy p0,x Bxy
Uy = uy − ux , B= . (A.40)
p0,y Bxx p0,y Bxx
For larger values of ℓ, the following formula can be used:
ℓ−k
⌊ 2 ⌋
H Bk X ℓ!
Gℓ;k = Is (Ty )Uyℓ−2s−k . (A.41)
k! s=0 2s s!(ℓ − 2s − k)!

63
Appendix A.6. d = 2: S model

For the d = 2 case, the same strategy employed in Subsec. Appendix A.4
can be employed. Expanding φS and χS via
 
Nx ∗ N∗
ω(px )ω(py ) X
φS X y
1
 = hr (|px |) Hℓ (py )
χS p0,x p0,y r=0 ℓ!
ℓ=0
    
Φ+;H Φ −;H
× θ(px )  S;r,ℓ  + θ(−px )  S;r,ℓ  , (A.42)
+ −
XS;r,ℓ XS;r,ℓ

the coefficients Φ± ±
S;r,ℓ and XS;r,ℓ can be obtained by taking into account Eq. (90):
 
   
Φ±,H
S;r,ℓ n 1 − Pr Gφ;±,H
  =   GℓH (uy , T )Gr± (ux , T ) +  S;r,ℓ  . (A.43)
±,H
XS;r,ℓ P 5nT 2 Gχ;±,H
S;r,ℓ

The coefficients G∗;±,H

S;r,ℓ are obtained by computing the following integrals:

∞
X Z ∞
ω(px ) ∗;H
G∗;±,H
S;r,ℓ = Gr±′ (ux , T ) dpx G (±px )hr (px )hr′ (px ), (A.44)
0 p0,x S;ℓ
r ′ =0

where G∗;H
S;ℓ (±px ) was introduced in Eq. (90). As in Eq. (A.34), the summation

with respect to r′ was extended to ∞. Using now an expansion of G∗;H

S;ℓ (±px )

similar to that introduced in Eq. (A.37),

3
X
G∗;H
S;ℓ (±px ) = (±1)k G∗;H k
S;ℓ;k px , (A.45)
k=0

the integral in Eq. (A.44) can be performed using the recurrence relation in
Eq. (A.26):
3
X k
X (k)
G∗;±,H
S;r,ℓ = (±1)k G∗;H
S;ℓ;k
±
Br,r′ Gr+r ′ (ux , T ). (A.46)
k=0 r ′ =−k

We close this subsection by giving explicitly the expressions for G∗;H

S;ℓ;k for

the values of ℓ and k that are relevant to this paper. From Eq. (93), it can be
seen that px enters G∗;H ∗
S;ℓ only through the terms Is , defined in Eq. (92) and

given explicitly for 0 ≤ s ≤ 3 in Eq. (94). Thus, the expression for G∗;H
S;ℓ;k is the

64
 HT SH HT-SH

δ TS (s) TES (s) TS (s) TES (s) TS (s) TES (s)

1000 674 320 992 699 2106 1368
100 67 32 108 68 197 135

Table A.3: Computational times (in seconds) required to reach the steady state using the
projection method described in Appendix A at δ = 100 and δ = 1000, in the context of the
heat transfer between stationary plates (HT), Couette flow (SH) and heat transfer under shear
(HT-SH) problems, considered in Sections 5, 6 and 7.

same as that for G∗;H ∗ ∗ ∗

S;ℓ , but with Is replaced by Is;k , where Is;k represents the

coefficient of pkx in Is∗ . Specifically, we find

         
φ 2 φ 2
I0;0 3 mu I0;1 3mu 3
  =mux qx   − x
,   =p0,x qx  x
−   ,
χ
I0;0 1 T χ
I0;1 T 1

∗
3p20,x ∗
p30,x
I0;2 =− ux qx , I0;3 = qx ,
  T   mT
φ 2
I1;0 −1
  =mT qy  mux +   , ∗
I1;1 = − 2mux p0,x qy ,
χ
I1;0 T 1
∗
I1;2 =p20,x qy , ∗
I1;3 =0,
         
φ 2 φ 2
I2;0 −1 mu I2;1 3mu −1
  = − m2 ux qx T   + x
,   =mp0,x qx T  x
+   ,
χ
I2;0 1 T χ
I2;1 T 1
∗
I2;2 = − 3mp20,x ux qx , ∗
I2;3 =p30,x qx ,
    
φ 2
I3;0 mu x 1
  =3(mT )2 qy  +   , ∗
I3;1 = − 6m2 T p0,xux qy ,
χ
I3;0 T 3
∗
I3;2 =3mT p20,xqy , ∗
I3;3 =0. (A.47)

Appendix A.7. Numerical results

690 In order to demonstrate the capabilities of the model introduced in the
previous subsections, we performed simulations in the context of the problems
introduced in Sections 5, 6 and 7 for δ = 100 and 1000. For definiteness, we
considered the 4 He gas at Teref = 300 K. The simulations were performed using

65
 10-1 10-3
S: Q S: Q
ES: Q ES: Q
Slope -2.1 Slope -2.1

10-2 10-4

Absolute error
Relative error

10-3 10-5

(a) (δ = 100) (b) (δ = 100)

10-4 10-6
10 100 10 100
S S

0 -2
10 10
S: Π (δ = 100) S: Π
Qy Qy
-1 Qw -3 Qw
10 10
ES: Π ES: Π
Qy Qy
-2 -4
10 Qw 10 Qw
Absolute error
Relative error

-3 -5
10 10

-4 -6
10 10

-5 Slope -2.9 -7
Slope -2.9
10 10

(c) (δ = 100) (d)

-6 -8
10 10
10 100 10 100
S S

-1 -2
10 10
S: Π S: Π
Q Q
ES: Π ES: Π
Q Q
10-3
10-2 Slope -2.1 Slope -2.1
Absolute error
Relative error

10-4

10-3

10-5

-4
10
(e) (δ = 100) (f) (δ = 100)
10-6
10 100 10 100
S S

Figure A.15: Convergence test for the results obtained using the projection method with
respect to the number S of nodes in the half channel. The heat transfer between stationary
plates (top line), Couette flow (middle line) and heat transfer between moving plates (bottom
line) are considered at δ = 100. The results66obtained using the S model are shown using
filled symbols, while the ES model results are represented using empty symbols. The relative
(left) and absolute (right) errors are computed at the level of the various quantities introduced
in Sec. 4 by taking the results obtained using the hybrid method with Qφ χ
x = Qx = 16 and

S = 128 as a reference. The slope of the dotted line indicates the convergence order.
 -3
10
-2 S: Q S: Q
10
ES: Q ES: Q
Slope -3.6 -4 Slope -3.6
10

10-3
10-5

Absolute error
Relative error

10-6
10-4

10-7
-5
10
-8
10

(a) (δ = 1000) (b) (δ = 1000)

10-6 10-9
10 100 10 100
S S

1 -2
10 10
S: Π (δ = 1000) S: Π
Qy Qy
0 Qw -3 Qw
10 10
ES: Π ES: Π
Qy Qy
-1 -4
10 Qw 10 Qw
Absolute error
Relative error

-2 -5
10 10

-3 -6 Slope -2.2
10 10
Slope -2.2
-4 -7
10 10

(c) (δ = 1000) (d)

-5 -8
10 10
10 100 10 100
S S

10-1 10-2
S: Π S: Π
Q Q
ES: Π 10
-3 ES: Π
-2 Q Q
10
Slope -2.6 Slope -2.6
10-4
Absolute error
Relative error

10-3 10-5

10-6
10-4
10-7

(e) (δ = 1000) (f) (δ = 1000)

-5 -8
10 10
10 100 10 100
S S

Figure A.16: Same as Fig. A.15, for δ = 1000.

67
 the quadrature orders Qφx = 5 and Qχx = 3 for the half-range Gauss-Hermite
695 quadrature employed on the px axis. In the case of the heat transfer between
stationary plates (discussed in Sec. 5), there are no other non-trivial degrees of
freedom and the total number of distinct populations is 2(Qφx + Qχx ) = 16. In
the Couette flow and heat transfer between moving plates problems, discussed
in Sections 6 and 7, the py axis was discretised using the full-range Gauss-
700 Hermite quadrature of orders Qφy = 4 and Qχy = 2, as discussed in Sec. 3.2. The
total number of distinct population in this case is 2(Qφx Qφy + Qχx Qχy ) = 52. In
the context of this subsection, we refer to the method employing the models
described above as the “projection method.”
In order to validate the results obtained with the projection method, we also
705 performed simulations using the hybrid models introduced in Sec. 3, which differ
from the former since the equilibrium distributions are replaced by truncated
expansions only with respect to the py axis. The quadrature orders in the hybrid
approach are set to Qφx = Qχx = 16, while the quadrature orders Qφy and Qχy
remain unchanged with respect to those employed within the projection method.
710 Figures A.15 and A.16 show convergence tests performed at δ = 100 and
1000, respectively, comparing the results obtained using the projection method
for various values of S with those obtained using the hybrid method using
Qx = 16 and S = 128, in the context of the heat transfer between stationary
plates (top line), Couette flow (middle line) and heat transfer between moving
715 plates (bottom line) problems. The left columns of Figs. A.15 and A.16 show
the relative errors, while the right columns show the corresponding absolute
errors, computed at the level of the Q, Π, Qy and Qw quantities introduced
in Sec. 4. The resulting convergence orders take values between 2.1 and 3.6,
depending on δ and on the quantity being analysed.
720 While the relative errors in Q and Π quickly approach 0.1% as S is increased,
it can be seen from the left columns of Figs. A.15 and A.16 that achieving the
same level of relative error for Qw and Qy in the context of the Couette flow
becomes more challenging at large δ. Taking into account that all quantities
decrease in absolute value as δ −1 (except Qy , which decreases as δ −2 ), the

68
725 relative error is correspondingly amplified and thus becomes less relevant. By
comparison, the right columns of the same figures show that the absolute errors
are several orders of magnitude below the relative errors. It can be seen that
reasonable results are obtained using the projection method with S = 32, which
gives absolute errors that are below 10−4 for al quantities under consideration.
730 The corresponding runtimes for S = 32 are summarised in Table A.3.

References

[1] C. Cercignani. Rarefied Gas Dynamics: From Basic Concepts to Actual

Calculations. Cambridge University Press, Cambridge, 2000.

[2] Y. Sone. Molecular Gas Dynamics: Theory, Techniques and Applications.

735 Birkhäuser, Boston, 2007.

[3] S. Takata and H. Funagane. Singular behaviour of a rarefied gas on a

planar boundary. J. Fluid Mech., 717:30–47, 2013.

[4] M. Gad-el-Haq. MEMS Handbook. CRC Press, Boca Raton, 2006.

[5] G. A. Bird. Molecular Gas Dynamics and the Direct Simulation of Gas
740 Flows. Oxford University Press, Oxford, 1994.

[6] F. Sharipov and J. L. Strapasson. Benchmark problems for mixtures of

rarefied gases. I. Couette flow. Phys. Fluids, 25:027101, 2013.

[7] F. Sharipov and J. L. Strapasson. Ab initio simulation of rarefied gas flow

through a thin orifice. Vacuum, 109:246–252, 2014.

745 [8] F. Sharipov and C. F. Dias. Ab initio simulation of planar shock waves.
Computers and Fluids, 150:115–122, 2017.

[9] A. N. Volkov and F. Sharipov. Flow of a monatomic rarefied gas over a

circular cylinder: Calculations based on the ab initio potential method.
Int. J. Heat Mass Transfer., 114:47–61, 2017.

69
750 [10] L. Zhu, L. Wu, Y. Zhang, and F. Sharipov. Ab initio calculation of rarefied
flows of helium-neon mixture: Classical vs quantum scatterings. Int. J.
Heat Mass Tran., 145:118765, 2019.

[11] F. Sharipov. Modeling of transport phenomena in gases based on quantum

scattering. Physica A, 508:797–805, 2018.

755 [12] F. Sharipov and C. F. Dias. Temperature dependence of shock wave

structure in helium and neon. Phys. Fluids, 31:037109, 2019.

[13] C. Mouhot and L. Pareschi. Fast algorithms for computing the Boltzmann
collision operator. Math. Comput., 75:1833–1852, 2006.

[14] F. Filbet. On deterministic approximation of the Boltzmann equation in

760 a bounded domain. Multiscale Model. Simul., 10:792–817, 2012.

[15] L. Wu, C. White, T. J. Scanlon, J. M. Reese, and Y. Zhang. Deterministic

numerical solutions of the Boltzmann equation using the fast spectral
method. J. Comput. Phys., 250:27–52, 2013.

[16] L. Wu, H. Liu, Y. Zhang, and J. M. Reese Influence of intermolecular

765 potentials on rarefied gas flows: Fast spectral solutions of the Boltzmann
equation. Phys. Fluids, 27:082002, 2015.

[17] I. M. Gamba, J. R. Haack, C. D. Hauck, and J. Hu. A fast spectral method

for the Boltzmann collision operator with general collision kernels. SIAM
J. Sci. Comput., 39:B658–B674, 2017.

770 [18] P. L. Bhatnagar, E. P. Gross, and M. Krook. A model for collision pro-
cesses in gases. I. Small amplitude processes in charged and neutral one-
component systems. Phys. Rev., 94:511–525, 1954.

[19] L. H. Holway, Jr. New statistical models for kinetic theory: methods of
construction. Phys. Fluids, 9:1658–1673, 1966.

775 [20] E. M. Shakhov. Generalization of the Krook kinetic relaxation equation.

Fluid Dyn., 3:95–96, 1968.

70
 [21] E. M. Shakhov. Approximate kinetic equations in rarefied gas theory.
Fluid Dyn., 3:112–115, 1968.

[22] F. Sharipov. Application of the Cercignani-Lampis scattering kernel to

780 calculations of rarefied gas flows. I. Plane flow between two parallel plates.
Eur. J. Mech. B-Fluid, 21:113–123, 2002.

[23] F. Sharipov. Application of the Cercignani-Lampis scattering kernel to

calculations of rarefied gas flows. II. Slip and jump coefficients. Eur. J.
Mech. B-Fluid, 22:133–143, 2003.

785 [24] I. A. Graur and A. P. Polikarpov. Comparison of different kinetic models

for the heat transfer problem. Heat Mass Transfer, 46:237–244, 2009.

[25] V. E. Ambrus, and V. Sofonea. High-order thermal lattice Boltzmann

models derived by means of Gauss quadrature in the spherical coordinate
system. Phys. Rev. E, 86:016708, 2012.

790 [26] J. P. Meng, Y. H. Zhang, N. G. Hadjiconstantinou, G. A. Radtke, and

X. W. Shan. Lattice ellipsoidal statistical BGK model for thermal non-
equilibrium flows. J. Fluid. Mech., 718:347–370, 2013.

[27] V. E. Ambrus, and V. Sofonea. Half-range lattice Boltzmann models for

the simulation of Couette flow using the Shakhov collision term. Phys.
795 Rev. E, 98:063311, 2018.

[28] J. E. Broadwell. Study of rarefied shear flow by the discrete velocity

method. J. Fluid Mech., 19:401–414, 1964.

[29] F. Sharipov. Rarefied gas dynamics: Fundamentals for research and prac-
tice. Wiley-VCH, Weinheim, 2016.

800 [30] M. T. Ho and I. Graur. Heat transfer through rarefied gas confined be-
tween two concentric spheres. Int. J. Heat Mass Transfer, 90:58–71, 2015.

71
 [31] Z. Guo, K. Xu, and R. Wang. Discrete unified gas kinetic scheme for
all Knudsen number flows: Low-speed isothermal case. Phys. Rev. E,
88:033305, 9 2013.

805 [32] Z. Guo, R. Wang, and K. Xu. Discrete unified gas kinetic scheme for
all Knudsen number flows. II. thermal compressible case. Phys. Rev. E,
91:033313, 3 2015.

[33] L. Zhu, P. Wang, and Z. Guo. Performance evaluation of the general

characteristics based off-lattice Boltzmann scheme and DUGKS for low
810 speed continuum flows. J. Comput. Phys., 333:227–246, 2017.

[34] X. He, X. Shan, and G. D. Doolen. Discrete Boltzmann equation model

for nonideal gases. Phys. Rev. E, 57:R13–R16, 1998.

[35] C. Lin, K. H. Luo, L. Fei, and S. Succi. A multi-component discrete

Boltzmann model for nonequilibrium reactive flows. Sci. Rep., 7:14580,
815 2017.

[36] Y.-D. Zhang, A.-G. Xu, G.-C. Zhang, Z.-H. Chen, and P. Wang. Discrete
ellipsoidal statistical BGK model and Burnett equations. Front. Phys.,
13:135101, 2018.

[37] W. P. Yudistiawan, S. Ansumali, and I. V. Karlin. Hydrodynamics beyond

820 Navier-Stokes: The slip flow model. Phys. Rev. E, 78:016705, 2008.

[38] W. P. Yudistiawan, S. K. Kwak, D. V. Patil, and S. Ansumali. Higher-

order Galilean-invariant lattice Boltzmann model for microflows: Single-
component gas. Phys. Rev. E, 82:046701, 2010.

[39] C. Feuchter and W. Schleifenbaum. High-order lattice Boltzmann mod-

825 els for wall-bounded flows at finite Knudsen numbers. Phys. Rev. E,
94:013304, 2016.

[40] M. Atif, M. Namburi, and S. Ansumali. Higher-order lattice Boltzmann

model for thermohydrodynamics. Phys. Rev. E, 98:053311, 2018.

72
 [41] C. K. Aidun and J. R. Clausen. Lattice-Boltzmann method for complex
830 flows. Annu. Rev. Fluid Mech., 42:439–472, 2010.

[42] J. P. Meng and Y. H. Zhang. Gauss-Hermite quadratures and accuracy

of lattice Boltzmann models for nonequilibrium gas flows. Phys. Rev. E,
83:036704, 2011.

[43] Y. Shi, Y. W. Yap, and J. E. Sader. Linearized lattice Boltzmann method

835 for micro- and nanoscale flow and heat transfer. Phys. Rev. E, 92:013307,
2015.

[44] V. E. Ambrus, and V. Sofonea. Lattice Boltzmann models based on half-

range Gauss-Hermite quadratures. J. Comput. Phys., 316:760–788, 2016.

[45] S. Succi. The Lattice Boltzmann Equation for Fluid Dynamics and Be-
840 yond. Clarendon Press, Oxford, 2001.

[46] P. Fede, V. Sofonea, R. Fournier, S. Blanco, O. Simonin, G. Lepoutère,

and V. E. Ambrus, . Lattice Boltzmann model for predicting the deposition
of inertial particles transported by a turbulent flow. Int. J. Multiph. Flow,
76:187–197, 2015.

845 [47] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. M.

Viggen. Lattice Boltzmann Method: Principles and Practice. Springer,
2017.

[48] S. Succi. The Lattice Boltzmann Equation: For Complex States of Flowing
Matter. Oxford University Press, 2018.

850 [49] X. W. Shan, X. F. Yuan, and H. D. Chen. Kinetic theory representation

of hydrodynamics: a way beyond the Navier-Stokes equation. J. Fluid.
Mech., 550:413–441, 2006.

[50] F. B. Hildebrand. Introduction to Numerical Analysis. Dover Publications,

second edition edition, 1987.

73
855 [51] B. Shizgal. Spectral Methods in Chemistry and Physics: Applications
to Kinetic Theory and Quantum Mechanics (Scientific Computation).
Springer, 2015.

[52] F. Sharipov. Rarefied gas flow through a long tube at any temperature
ratio. J. Vac. Sci. Technol. A, 14:2627–2635, 1996.

860 [53] L. Wu, J. M. Reese, and Y. Zhang. Solving the Boltzmann equation de-
terministically by the fast spectral method: application to gas microflows.
J. Fluid Mech., 746:53–84, 2014.

[54] F. Sharipov. Rarefied gas flow through a long rectangular channel. J.

Vac. Sci. Technol. A, 17:3062–3066, 1999.

865 [55] E. P. Gross, E. A. Jackson, and S. Ziering. Boundary value problems in

kinetic theory of gases. Ann. Phys., 1:141–167, 1957.

[56] Y. Sone. Kinetic theory analysis of linearized Rayleigh problem. J. Phys.

Soc. Jpn., 19:1463–1473, 1964.

[57] S. Jiang and L.-S. Luo. Analysis and accurate numerical solutions of
870 the integral equation derived from the linearized BGKW equation for the
steady Couette flow. J. Comput. Phys, 316:416–434, 2016.

[58] E. P. Gross and S. Ziering. Kinetic theory of linear shear flow. Phys.
Fluids, 1:215–224, 1958.

[59] S. Ziering. Shear and heat flow for Maxwellian molecules. Phys. Fluids,
875 3:503–509, 1960.

[60] P. L. Bhatnagar and M. P. Srivastava. Heat transfer in plane Couette

flow of a rarefied gas using Bhatnagar-Gross-Krook model. Phys. Fluids,
12:938–940, 1969.

[61] A. Frezzotti, L. Gibelli, and B. Franzelli. A moment method for low speed
880 microflows. Continuum Mech. Thermodyn., 21:495–509, 2009.

74
 [62] L. Gibelli. Velocity slip coefficients based on the hard-sphere Boltzmann
equation. Phys. Fluids, 24:022001, 2012.

[63] G. P. Ghiroldi and L. Gibelli. A direct method for the Boltzmann equation
based on a pseudo-spectral velocity space discretization. J. Comput. Phys.,
885 258:568–584, 2014.

[64] A. B. Huang and D. P. Giddens. A new table for a modified (half

range) GaussHermite quadrature with an evaluation of the integral
R ∞ −u2 −(z/u)
0 e du. J. Math. Phys., 47:213–218, 1968.

[65] J. S. Ball. Half-range generalized Hermite polynomials and the related

890 Gaussian quadratures. SIAM J. Numer. Anal., 40:2311–2317, 2003.

[66] Z.-H. Li and H.-X. Zhang. Numerical investigation from rarefied flow to
continuum by solving the Boltzmann model equation. Int. J. Numer.
Meth. Fluids, 42:361–382, 2003.

[67] Z.-H. Li and H.-X. Zhang. Study on gas kinetic unified algorithm for
895 flows from rarefied transition to continuum. J. Comput. Phys, 193:708–
738, 2004.

[68] G. P. Ghiroldi and L. Gibelli. A finite-difference lattice Boltzmann ap-

proach for gas microflows. Commun. Comput. Phys., 17:1007–1018, 2015.

[69] V. E. Ambrus, and V. Sofonea. Lattice Boltzmann models based on Gauss

900 quadratures. Int. J. Mod. Phys. C, 25:1441011, 2014.

[70] V. E. Ambrus, and V. Sofonea. Implementation of diffuse-reflection bound-

ary conditions using lattice Boltzmann models based on half-space Gauss-
Laguerre quadratures. Phys. Rev. E, 89:041301(R), 4 2014.

[71] V. E. Ambrus, and V. Sofonea. Quadrature-based lattice Boltzmann mod-

905 els for rarefied gas flow. In: F. Toschi and M. Sega, editors, Flowing
Matter. Soft and Biological Matter, chapter 9, 271–299. Springer, Cham,
2019.

75
 [72] D. Valougeorgis and S. Naris. Acceleration schemes of the discrete veloc-
ity method: Gaseous flows in rectangular microchannels. SIAM J. Sci.
910 Comput., 25:534–552, 2003.

[73] L. Wu, J. Zhang, H. Liu, Y. Zhang, and J. M. Reese. A fast iterative

scheme for the linearized Boltzmann equation. J. Comput. Phys., 338:431–
451, 2017.

[74] W. Su, P. Wang, H. Liu, and L. Wu. Accurate and efficient computation of
915 the Boltzmann equation for Couette flow: Influence of intermolecular po-
tentials on Knudsen layer function and viscous slip coefficient. J. Comput.
Phys., 378:573–590, 2019.

[75] L. Zhu, X. Pi, W. Su, Z.-H. Li, Y. Zhang, and L. Wu. General Synthetic
Iteration Scheme for Non-linear Gas Kinetic Simulation of Multi-scale
920 Rarefied Gas flows. arXiv:2004.10530 [physics.comp-ph].

[76] Chi-Wang Shu and Stanley Osher. Efficient implementation of essentially

non-oscillatory shock-capturing schemes. J. Comput. Phys., 77:439–471,
1988.

[77] S. Gottlieb and C.-W. Shu. Total variation diminishing Runge-Kutta

925 schemes. Math. Comput., 67:73–85, 1998.

[78] J. A. Trangenstein. Numerical solution of hyperbolic partial differential

equations. Cambridge University Press, New York, 2007.

[79] V. E. Ambrus, and L.-S. Luo. Analysis of Knudsen layer phenomena using
half-range quadratures. 2019. In preparation.

930 [80] R. Mei and W. Shyy. On the finite difference-based lattice Boltzmann
method in curvilinear coordinates. J. Comput. Phys., 143:426–448, 1998.

[81] Z. Guo and T. S. Zhao. Explicit finite-difference lattice Boltzmann method

for curvilinear coordinates. Phys. Rev. E, 67:066709, 2003.

76
 [82] S. Busuioc and V. E. Ambrus, . Lattice Boltzmann models based on the
935 vielbein formalism for the simulation of flows in curvilinear geometries.
Phys. Rev. E, 99:033304, 2019.

[83] V. Sofonea and R. F. Sekerka. Viscosity of finite difference lattice Boltz-

mann models. J. Comput. Phys., 183:422–434, 2003.

[84] V. Sofonea, A. Lamura, G. Gonnella, and A. Cristea. Finite-difference

940 lattice Boltzmann model with flux limiters for liquid-vapor systems. Phys.
Rev. E, 70:046702, 2004.

[85] S. Naris, D. Valougeorgis, D. Kalempa, and F. Sharipov. Gaseous mixture

flow between two parallel plates in the whole range of the gas rarefaction.
Physica A, 336:294–318, 2004.

945 [86] C. Tantos and D. Valougeorgis. Conductive heat transfer in rarefied binary
gas mixtures confined between parallel plates based on kinetic modeling.
Int. J. Heat Mass Tran., 117:846–860, 2018.

[87] C. Cercignani and M. Lampis. Kinetic model for gas-surface interaction.

Transp. Theory Stat. Phys., 1:101–114, 1971.

950 [88] V. E. Ambrus, , F. Sharipov, and V. Sofonea. Lattice Boltzmann approach

to rarefied gas flows using half-range Gauss-Hermite quadratures: Com-
parison to DSMC results based on ab initio potentials. AIP Conf. Proc.,
2132:060012, 2019.

[89] J. Meng, L. Wu, J. M. Reese, and Y. Zhang. Assessment of the ellipsoidal-

955 statistical Bhatnagar-Gross-Krook model for force-driven Poiseuille flow.
J. Comput. Phys., 251:383–395, 2013.

[90] W. Cencek, M. Przybytek, J. Komasa, J. B. Mehl, B. Jeziorski, and

K. Szalewicz. Effects of adiabatic, relativistic and quantum electrody-
namics interactions on the pair potential and thermophysical properties
960 of helium. J. Chem. Phys., 136:224303, 2012.

77
 [91] F. Sharipov and V. J. Benites. Transport coefficients of helium-neon mix-
tures at low density computed from ab initio potentials. J. Chem. Phys.,
147:224302, 2017.

[92] V. E. Ambrus, and V. Sofonea. Application of mixed quadrature lattice

965 Boltzmann models for the simulation of Poiseuille flow at non-negligible
values of the Knudsen number. J. Comput. Sci., 17:403–417, 2016.

[93] S. Ansumali, I. V. Karlin, and H. C. Öttinger. Minimal entropic kinetic

models for hydrodynamics. Europhys. Lett., 63:798, 2003.

[94] A. Bardow, I. V. Karlin, and A. A. Gusev. General characteristic-based

970 algorithm for off-lattice Boltzmann simulations. Europhys. Lett., 75:434,
2006.

[95] A. Bardow, I. V. Karlin, and A. A. Gusev. Multispeed models in off-lattice

Boltzmann simulations. Phys. Rev. E, 77:025701(R), 2008.

[96] T. Bicius, că, A. Horga, and V. Sofonea. Simulation of liquid-vapour phase

975 separation on GPUs using Lattice Boltzmann models with off-lattice ve-
locity sets. C. R. Mecanique, 343:580–588, 2015.

[97] V. Sofonea, T. Bicius, că, S. Busuioc, V. E. Ambrus, , G. Gonnella, and

A. Lamura. Corner-transport-upwind lattice Boltzmann model for bubble
cavitation. Phys. Rev. E, 97:023309, 2018.

980 [98] L. Pareschi and G. Russo. Implicit-Explicit Runge-Kutta schemes and ap-
plications to hyperbolic systems with relaxation. J. Sci. Comput., 25:129–
155, 2005.

[99] Y. Wand, Y. L. He, T. S. Zhao, G. H. Tang, and W. Q. Tao. Implicit-

Explicit finite-difference lattice Boltzmann method for compressible flows.
985 Int. J. Mod. Phys. C, 18:1961–1983, 2007.

[100] S. T. Kis and V. E. Ambrus, . Implicit-explicit finite-difference lattice Boltz-

mann model with varying adiabatic index. AIP Conf. Proc., 2218:050008,
2020.

78
 [101] Jr. W. Marques, G. M. Kremer, and F. M. Sharipov. Couette flow with slip
990 and jump boundary conditions. Continuum Mech. Thermodyn., 12:379–
386, 2000.

[102] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Jed Brown, Peter
Brune, Kris Buschelman, Lisandro Dalcin, Victor Eijkhout, William D.
Gropp, Dinesh Kaushik, Matthew G. Knepley, Lois Curfman McInnes,
995 Karl Rupp, Barry F. Smith, Stefano Zampini, and Hong Zhang. PETSc
users manual. Technical Report ANL-95/11 - Revision 3.6, Argonne Na-
tional Laboratory, 2015.

[103] Satish Balay, William D. Gropp, Lois Curfman McInnes, and Barry F.
Smith. Efficient management of parallelism in object oriented numerical
1000 software libraries. In E. Arge, A. M. Bruaset, and H. P. Langtangen,
editors, Modern Software Tools in Scientific Computing, pages 163–202.
Birkhäuser Press, 1997.

[104] P. K. Kundu, I. M. Cohen, D. R. Dowling, Fluid Mechanics, 6th Ed.

(Academic Press, 2016).

79

View publication stats