^†^†thanks: Corresponding author

Probing double distribution function models in the lattice Boltzmann method for highly compressible flows

S. A. Hosseini Department of Mechanical and Process Engineering, ETH Zurich, 8092 Zurich, Switzerland A. Bhadauria Department of Mechanical and Process Engineering, ETH Zurich, 8092 Zurich, Switzerland I. V. Karlin [email protected] Department of Mechanical and Process Engineering, ETH Zurich, 8092 Zurich, Switzerland

Abstract

The double distribution function approach is an efficient route towards extension of kinetic solvers to compressible flows. With a number of realizations available, an overview and comparative study in the context of high speed compressible flows is presented. We discuss the different variants of the energy partition, analyses of hydrodynamic limits and a numerical study of accuracy and performance with the particles on demand realization. Out of three considered energy partition strategies, it is shown that the non-translational energy split requires a higher-order quadrature for proper recovery of the Navier–Stokes–Fourier equations. The internal energy split on the other hand, while recovering the correct hydrodynamic limit with fourth-order quadrature, comes with a non-local –both in space and time– source term which contributes to higher computational cost and memory overhead. Based on our analysis, the total energy split demonstrates the optimal overall performance.

I Introduction

Extension of the lattice Boltzmann method (LBM), and other discrete velocity methods for the Boltzmann equation, to compressible flow simulations has been a topic of interest over the past decade [1]. A number of different strategies have been devised to extend the original isothermal weakly compressible LBM to compressible flows with energy balance. One straightforward approach is the use of larger lattices allowing the model to properly recover more moments of the equilibrium distribution function -necessary for the energy balance up to Navier–Stokes–Fourier level [2, 3, 4, 5]. This strategy, as it has been observed over the past few years, has a number of limitations such as an increase in non-locality of the streaming operator, larger computational load and memory footprint. It also has a limited range of operation in terms of the temperature.

An alternative, allowing the use of smaller lattices, is the so-called double distribution function approach [6, 7, 8, 9, 10] inspired by Rykov’s original work [11]. In this approach, the distribution function is split into two reduced distribution functions where the second one carries a form of energy. This idea has become a popular approach in discrete kinetic theory of gases, in combination with the Bhatnagar–Gross–Krook (BGK) and similar collision operators in order to extend isothermal solvers to compressible flows. While in Rykov’s original model, the second distribution function was intended to solely transport internal energy due to non-translational degrees of freedom for polyatomic molecules and did not make any assumptions as to equilibrium between internal and translational degrees of freedom, in the context of discrete solvers the choice of the variable transported by the second distribution function is not unique and usually assumes equipartition, though recent proposals with non-equilibrium temperatures can be found [12]. Different choices of variables have been used and presented over the years. The aim of the present contribution is to discuss the non-uniqueness in the choice of the variable carried by the second population by looking at and analyzing the effects of each one of the popular choices in the literature, with a focus on applications involving larger Mach numbers. The focus of the discussion being on the double distribution function approach and effect of the choice of variable carried by the second distribution we will only consider BGK collision operators. Nonetheless, the results can readily be extended to variable Prandtl numbers and/or independent bulk viscosity via introduction of quasi-equilibrium states [13, 14], for instance. The article will begin with a brief overview of the double distribution function formalism and then discuss the consequences of three major choices, namely non-translational internal energy, internal energy and total energy, all in the context of a single relaxation time BGK collision operator. Finally numerical simulation results, conducted using a particles-on-demand realization [15], probing all hydrodynamic level properties will be presented and discussed.

II Background and theory

In the context of the present contribution we are interested in the Boltzmann equation with a BGK approximation for the collision operator [16], i.e.

\partial_{t}f+\bm{v}\cdot\bm{\nabla}f=\frac{1}{\tau}\left(f-f^{\rm eq}\right).

(1)

Here, $f$ represents the first reduced probability distribution function with $\bm{v}$ the particle velocity, $\bm{x}$ the position in physical space and $t$ time. The definition of the first reduced distribution function, for molecules endowed with internal degrees of freedom follows that of Rykov in [11]. The reduced distribution function $f$ is subject to the balance equation (1), with the equilibrium defined as,

f^{\rm eq}=\frac{\rho}{{(2\pi rT)}^{D/2}}\exp{\left(-\frac{{(\bm{v}-\bm{u})}^{2}}{2rT}\right)},

(2)

where $T$ is the temperature, $\rho$ density, $\bm{u}$ velocity and $r$ the gas constant. Moreover, $\tau$ is the relaxation time, tied to the fluid dynamic viscosity, $\mu$ , as,

\tau=\frac{\mu}{p},

(3)

where $p=\rho rT$ is the thermodynamic pressure. While the zeroth- and first-order conserved moments are computed with the first reduced distribution function,


$\displaystyle\int_{\mathbb{R}^{D}}fd\bm{v}$	$\displaystyle=\rho,$	(4a)
$\displaystyle\int_{\mathbb{R}^{D}}\bm{v}fd\bm{v}$	$\displaystyle=\rho\bm{u},$	(4b)

the second reduced distribution function carries a form of energy. Let us introduce the kinetic energy,

K=\frac{1}{2}\rho\bm{u}^{2},

(5)

and the internal energy,

U=\rho\int_{0}^{T}c_{v}dT,

(6)

where $c_{v}$ is specific heat at constant volume. The total energy $E$ is

E=U+K.

(7)

Furthermore, let $U^{\prime}$ be the part of the internal energy associated with the non-translational degrees of freedom,

U^{\prime}=U-\rho\frac{DrT}{2}.

(8)

The second distribution function can now be defined with respect to the total energy in several ways,


$\displaystyle E$	$\displaystyle=\int_{\mathbb{R}^{D}}gd\bm{v},$	(9a)
$\displaystyle E$	$\displaystyle=\int_{\mathbb{R}^{D}}gd\bm{v}+K,$	(9b)
$\displaystyle E$	$\displaystyle=\int_{\mathbb{R}^{D}}\left(g+\frac{\bm{v}^{2}}{2}f\right)d\bm{v}.$	(9c)

Note that the approach of Eq. (9a) was first proposed in [7, 8] in the context of the LBM while Eqs. (9b) and (9c) were discussed in [6] and [5], respectively. The model discussed here is slightly different from [6] in that it also accounts for internal degrees of freedom. For the sake of readability, in the remainder of the paper we will refer to Eqs.(9a), (9b) and (9c), respectively as the total, internal and non-translational energy splits. It is interesting to note that, while replacing the temperature $T$ in the Maxwell–Boltzmann equilibrium (2) with internal energy through Eq. (6), the reduced $f$ -equilibrium becomes dependent on the second $g$ -distribution function which carries a part or all of the internal energy. The following kinetic equations for the different realizations (9) of the second reduced distribution functions can be derived,


$\displaystyle\partial_{t}g+\bm{v}\cdot\bm{\nabla}g$	$\displaystyle=-\frac{1}{\tau}\left(g-g^{\rm eq}\right),$	(10a)
$\displaystyle\partial_{t}g+\bm{v}\cdot\bm{\nabla}g$	$\displaystyle=-\frac{1}{\tau}\left(g-g^{\rm eq}\right)-f(\bm{v}-\bm{u})\cdot\left(\partial_{t}\bm{u}+\bm{v}\cdot\bm{\nabla}\bm{u}\right),$	(10b)
$\displaystyle\partial_{t}g+\bm{v}\cdot\bm{\nabla}g$	$\displaystyle=-\frac{1}{\tau}\left(g-g^{\rm eq}\right),$	(10c)

where the reduced equilibrium distributions $g^{\rm eq}$ are respectively defined as,


$\displaystyle g^{\rm eq}$	$\displaystyle=\left(\frac{1}{\rho}U^{\prime}+\frac{{\bm{v}}^{2}}{2}\right)f^{\rm eq},$	(11a)
$\displaystyle g^{\rm eq}$	$\displaystyle=\left(\frac{1}{\rho}U^{\prime}+\frac{{(\bm{v}-\bm{u})}^{2}}{2}\right)f^{\rm eq},$	(11b)
$\displaystyle g^{\rm eq}$	$\displaystyle=\frac{1}{\rho}U^{\prime}f^{\rm eq}.$	(11c)

Details of the multi-scale Chapman-Enskog analysis are presented in the Appendix A. All realizations (9) lead to the Navier–Stokes equations for the momentum balance,

\partial_{t}\rho\bm{u}+\bm{\nabla}\cdot\rho\bm{u}\otimes\bm{u}+\bm{\nabla}p-\bm{\nabla}\cdot\bm{T}_{\rm NS}=0,

(12)

where the viscous stress tensor $\bm{T}_{\rm NS}$ is defined as,

\bm{T}_{\rm NS}=\mu\left[\bm{\nabla}\bm{u}+\bm{\nabla}\bm{u}^{\dagger}-\frac{2}{D}\bm{\nabla}\cdot\bm{u}\bm{I}\right]+\eta\bm{\nabla}\cdot\bm{u}\bm{I}.

(13)

Here we have introduced the bulk viscosity, $\eta$ . A closer look at the results show that all double distribution realizations maintain consistent bulk viscosity, i.e.

\eta=\left(\frac{2}{D}-\frac{r}{c_{v}}\right)\mu.

(14)

Finally, at the energy level all realizations recover the following,

\partial_{t}E+\bm{\nabla}\cdot\left(E+p\right)\bm{u}+\bm{\nabla}\cdot\bm{q}_{\rm NS}=0,

(15)

with,

\bm{q}_{\rm NS}=-\lambda\bm{\nabla}T-\bm{u}\cdot\bm{T}_{\rm NS},

(16)

where the first term is the Fourier heat flux and the second viscous heating. Note that the thermal conductivity $\lambda$ is tied to the relaxation time as,

\lambda=\tau p(c_{v}+r),

(17)

leading to fixed Prandtl number

{\rm Pr}=\frac{(c_{v}+r)p\tau}{(c_{v}+r)p\tau}=1,

(18)

and adiabatic exponent

\gamma=\frac{c_{v}+r}{c_{v}}.

(19)

In this study, we focus on the impact of various partitions of the energy (9) and thus use the simplest single relaxation time kinetic model for the sake of presentation. This certainly restricts the Prandtl number. Nevertheless, all energy partitions (9) discussed in the present study can be readily extended to a tunable Prandtl number, by introducing an intermediary state of relaxation via the quasi-equilibrium approach as discussed in [13]. These intermediary states are constructed via minimization of entropy under constraints defined by a set of select higher-order moments. The choice of additional constraints in the set as compared to the conserved moments, determines the fluxes with independent relaxation rates. Specific discussions on the quasi-equilibrium approach for independent control over bulk viscosity and Prandtl number can be found in [13, 14, 9, 10].

Now that the hydrodynamic limits of each realization has been clarified, let us go over minimum requirements for the discrete solver in each case; In discrete velocity methods such as LBM, the minimum number of degrees of freedom –i.e. discrete velocities– is determined by the number of moments of the equilibrium distribution function that need to be satisfied to recover the target hydrodynamic limit.

For instance, in the full energy model, the first reduced distribution function $f$ only needs to properly recover continuity and momentum balance equations to the Navier–Stokes order in the Chapman–Enskog expansion, which involves equilibrium moments up to order three, while the second distribution function $g$ requires accuracy of the equilibrium moments up to second order. This can be achieved either through a fourth-order quadrature, i.e. D $D$ Q $4^{D}$ lattice, or a third order quadrature, i.e. D $D$ Q $3^{D}$ lattice, with a correction term for the third-order moments [17, 18, 19]. For the internal energy model, the constraints on $f$ and $g$ remain the same, as the energy balance equation relies on a combination of the zeroth-order moments of both. The last realization, i.e. the non-translational splitting, defines energy as the sum of the second-order moment of $f$ and zeroth-order moment of $g$ , which brings in an additional constraint on the fourth-order moment of the $f$ distribution function,

\int_{\mathbb{R}^{D}}\bm{v}\otimes\bm{v}\otimes\bm{v}\otimes\bm{v}fd\bm{v}=\sum_{i=1}^{Q}\bm{c}_{i}\otimes\bm{c}_{i}\otimes\bm{c}_{i}\otimes\bm{c}_{i}f_{i},

and which increases the minimum lattice size to D $D$ Q $5^{D}$ . The present study will use fourth- and fifth-order quadrature-based lattices. Details of these lattices are given in Table 1.

Quadrature	$rT_{L}$	$c_{i\alpha}$	$w_{i\alpha}$	2-D lattice
D1Q3	1	$\{0,\pm\sqrt{3}\}$	$\{2/3,1/3\}$	D2Q9
D1Q4	1	$\{\pm\sqrt{3-\sqrt{6}},\pm\sqrt{3+\sqrt{6}}\}$	$\{(3+\sqrt{6})/12,(3-\sqrt{6})/12\}$	D2Q16
D1Q5	1	$\{0,\pm\sqrt{5-\sqrt{10}},\pm\sqrt{5+\sqrt{10}}\}$	$\{8/15,(7+2\sqrt{10})/60,(7-2\sqrt{10})/60\}$	D2Q25

Table 1: Details of quadratures used for simulations.

The above discussion establishes directly computational cost and memory footprint of each realization; This analysis has to be complemented with one further consideration which is the presence of source terms in Eq. (10b) involving space- and time-derivatives of $\bm{u}$ which brings in non-negligible computational cost, along with additional memory footprint due to the time-derivative. Finally, the presence of such non-local contributions, in the case of finite-differences approximations, as used for instance in [6], contributes to non-conservation issue which can become detrimental in high Mach number flows simulations. An overview of the properties of each splitting approaches is presented in Table 2.

Energy split	Eq. (9a)	Eq. (9b)	Eq. (9c)
Kinematic viscosity	Eq. (3)	Eq. (3)	Eq. (3)
Bulk viscosity	Eq. (14)	Eq. (14)	Eq. (14)
Thermal conductivity	Eq. (17)	Eq. (17)	Eq. (17)
Specific heat ratio	Eq. (19)	Eq. (19)	Eq. (19)
Minimum order of quadrature for NSF	four	four	five
Non-local source terms	No	Yes	No

Table 2: Overview of properties of different splitting strategies.

In the next section, we will probe points discussed here through numerical simulations.

III Numerical applications

In this section, simulation are conducted to verify theoretical results derived in the previous sections. First the numerical scheme is introduced and then simulation results are presented and discussed. For the sake of readability all variables will be presented in lattice units, i.e. normalized by grid and time-steps sizes, $\delta x$ and $\delta t$ ; Velocities are normalized by $\delta x/\delta t$ and temperatures $rT$ by $\delta x^{2}/\delta t^{2}$ .

III.1 Discrete velocity solver for highly compressible flows: PonD

As discussed in the introduction, given that final target of the present contribution is application to highly compressible flows, we use a numerical model developed in our group, [15], shown to be particularly well-adapted to such flows, i.e. the particles on Demand (PonD) method. To minimize errors in higher-order moments of the distribution function not supported by the lattice quadrature, in the PonD method [15], the reference frame is chosen via the local macroscopic quantities, i.e. fluid speed and temperature

rT^{\lambda}=rT,\ \bm{u}^{\lambda}=\bm{u},

(20)

and the particle velocity is locally defined as,

\bm{C}_{i}(rT^{\lambda},\bm{u}^{\lambda})=\sqrt{rT^{\lambda}}\bm{c}_{i}+\bm{u}^{\lambda}.

(21)

In so doing, equilibrium distribution functions are accurately evaluated and only depend on the local density

f_{i}^{eq}=\rho w_{i}

(22)

where $w_{i}$ are weights of the Gauss-Hermite quadrature [20]. Different from approaches such as [21, 22] where an adaptive reference-based Boltzmann equation is solved globally in the domain bringing in non-commutativity of the discrete velocities with the space and time-derivatives and resulting in additional terms in the Boltzmann equation involving derivatives of the reference frame velocity and temperature, in PonD at every point in space and time $\bm{x}$ and $t$ the Boltzmann equation is solved in the fixed reference frame of that point, here denoted as $\bar{\lambda}$ for the sake of readability, leading to the following local evolution equation:

D_{t}^{\bar{\lambda}}\mathcal{M}_{\lambda(\bm{x},t)}^{\bar{\lambda}}\{f_{i}^{\lambda(\bm{x},t)},g_{i}^{\lambda(\bm{x},t)}\}=\mathcal{M}_{\lambda(\bm{x},t)}^{\bar{\lambda}}\Omega_{\{f_{i}^{\lambda(\bm{x},t)},g_{i}^{\lambda(\bm{x},t)}\}},

(23)

where $\mathcal{M}_{\lambda(\bm{x},t)}^{\bar{\lambda}}$ is the reference frame transformation operator.

III.1.1 Frame transformation

In the previous section we have introduced the reference frame transformation, $\mathcal{M}_{\lambda(\bm{x},t)}^{\bar{\lambda}}$ , which allowing information to be transformed from one reference frame to other is an essential component of PonD. The main statement, allowing for such an operation is that discrete distribution functions can be equivalently written as the set of moment correctly matched by the quadrature and that these moment are frame-invariant, i.e.

\displaystyle\sum_{i=0}^{Q-1}f_{i}C_{ix}^{p}(rT^{\lambda},u_{x}^{\lambda})C_{iy}^{q}(rT^{\lambda},u_{y}^{\lambda})C_{iz}^{r}(rT^{\lambda},u_{z}^{\lambda})

\displaystyle=M^{\lambda}_{x^{p}y^{q}z^{r}},

(24)

This approach, also referred to as the moment matching method was developed and used in [15]. Here, in an effort to further control errors in higher-order moments not supported by the lattice during transformation, we make use of a reduced order Grad expansion [23], i.e. a regularized reconstruction of discrete populations [24, 25, 26]. The distribution functions in a reference frame $\bar{\lambda}$ can be computed from moment in a reference frame $\lambda$ as,

\displaystyle f_{i}^{\bar{\lambda}}=w_{i}\sum_{n=0}^{N}\frac{1}{n!rT_{L}}a^{\bar{\lambda}}_{(n)}:H_{i}^{(n)},

(25)

where $a^{\lambda}_{(n)}$ is the rank $n$ tensor of Hermite polynomials of the corresponding order in the reference frame $\lambda$ and $H_{i}^{(n)}$ the tensor of Hermite coefficients of the same order, $:$ the Frobenius inner product and $N$ the order of the Grad expansion. The Grad coefficients in reference frame $\bar{\lambda}$ can be computed from moments in the reference frame $\lambda$ as,

a_{0}^{\bar{\lambda}}=M^{\lambda}_{0},

(26)

a_{1}^{\bar{\lambda}}=\frac{1}{\sqrt{T^{\bar{\lambda}}/T_{L}}}\left(M^{\lambda}_{1}-M_{0}^{\lambda}u^{\bar{\lambda}}\right),

(27)

a_{2}^{\bar{\lambda}}=\frac{T_{L}}{T^{\bar{\lambda}}}\left(M^{\lambda}_{2}-M_{0}^{\lambda}rT^{\bar{\lambda}}\bm{I}-\sqrt{\frac{T^{\bar{\lambda}}}{T_{L}}}\overline{u^{\bar{\lambda}}\otimes a_{1}^{\bar{\lambda}}}-M_{0}^{\lambda}u^{\bar{\lambda}}\otimes u^{\bar{\lambda}}\right),

(28)

and

a_{3}^{\bar{\lambda}}=\frac{1}{{\left(T^{\bar{\lambda}}/T_{L}\right)}^{3/2}}\left(M^{\lambda}_{3}\right.\\ \left.-\frac{T^{\bar{\lambda}}}{T_{L}}\overline{u^{\bar{\lambda}}\otimes\left(M_{0}^{\lambda}T_{L}\bm{I}+a_{2}^{\bar{\lambda}}\right)}-T^{\bar{\lambda}}\sqrt{\frac{T^{\bar{\lambda}}}{T_{L}}}\overline{a_{1}^{\bar{\lambda}}\otimes\bm{I}}\right.\\ \left.-\sqrt{\frac{T^{\bar{\lambda}}}{T_{L}}}\overline{a_{1}^{\bar{\lambda}}\otimes u^{\bar{\lambda}}\otimes u^{\bar{\lambda}}}-M_{0}^{\lambda}u^{\bar{\lambda}}\otimes u^{\bar{\lambda}}\otimes u^{\bar{\lambda}}\right).

(29)

Here Eqs. (25) through (29) define the reference frame transformation operaor.

III.1.2 Streaming

While different space/time discretization approaches have been developed by our group in recent years, see for instance [27, 25], we will use the semi-Lagrangian discretization here. Propagation occurs along the co-moving particle velocities which are fully adaptive, and therefore no longer space filling,

\{f_{i},g_{i}\}^{\lambda}(\bm{x}+\bm{C}_{i}\delta t,t+\delta t)=\{f_{i}^{*},g_{i}^{*}\}^{\lambda}(\bm{x},t),

(30)

where $\{f_{i}^{*},g_{i}^{*}\}$ are post-collision distribution functions. As such, the streaming step is supplemented with a reconstruction step to get the discrete distribution functions on the grid points. Each of these populations exist in their local reference frames, which are distinct from one another and also distinct from the destination frame. We can write for a generalized $p$ -sized interpolation kernel $W(\bm{x})$ :

\displaystyle\{f_{i},g_{i}\}^{\bar{\lambda}}(\bm{x},t)=\sum_{s=0}^{p-1}W(\bm{x}-\bm{x}_{s})\mathcal{M}_{\lambda_{s}}^{\bar{\lambda}}\{f_{i},g_{i}\}^{\lambda_{s}}(\bm{x}_{s},t)

(31)

where $\bm{x}_{s}$ is the position of a node inside the interpolation kernel and $\lambda_{s}$ is the local reference frame at this node. In the context of the present work the reconstruction step is conducted through a 2nd-order Lagrangian interpolation stencil.

III.2 Probing hydrodynamic limits

III.2.1 Speed of sound: Effect of specific heats ratio

As a first step and to validate the dispersion properties of normal modes, we look at the temperature-dependence and effect of $c_{v}$ on the speed of sound. A freely travelling pressure front is simulated to that effect. A quasi-one-dimensional domain $L_{x}\times 1$ (with $L_{x}=800$ ) is separated into two parts with a pressure difference of $\Delta p=10^{-4}$ and a uniform initial temperature $T_{0}$ and velocity $\bm{u}_{0}=0$ . The speed of sound is computed by tracking the position of the shock front over time and compared with the theoretical value $c_{s}(T)=\sqrt{\gamma rT}$ . The simulations are performed with two different specific heat ratios, i.e. $\gamma=1.4$ and $1.8$ , in a wide range of temperature. From figure 1 we observe that all energy splitting strategies, as expected, can correctly capture the speed of sound at the considered range of temperatures. Note that temperatures are reported in non-dimensional form, i.e. $\theta=T/T_{L}$ and speed of sound is in lattice units.

Refer to caption — Figure 1: Speed of sound on the lattice plotted against different values of non-dimensional temperature. Red line indicates results for $\gamma=1.4$ and blue line $\gamma=1.8$ .

III.2.2 Measuring dissipation rates

Next we look into the dissipation rates of the three hydrodynamic modes, i.e. shear, normal and entropic. From the Chapman-Enskog analysis, the kinematic shear viscosity $\nu=\mu/\rho$ and bulk viscosity $\eta$ in all splittings are related to the relaxation coefficient $\tau$ as,


$\displaystyle\nu$	$\displaystyle=\frac{\mu}{\rho}=\tau rT,$	(32a)
$\displaystyle\frac{\eta}{\rho}$	$\displaystyle=\left(\frac{2}{D}-\frac{1}{c_{v}}\right)\tau rT,$	(32b)
$\displaystyle\alpha$	$\displaystyle=\frac{\lambda}{(c_{v}+r)\rho}=\tau rT.$	(32c)

To measure effective dissipation rates we conducted three set of simulations.

First, to measure the kinematic shear viscosity a plane shear wave is simulated. The wave is introduced via a small sinusoidal perturbation added to the initial velocity field. The initial conditions of the flows are

\rho=\rho_{0},T=T_{0},u_{x}=u_{0},u_{y}=A\sin\left(2\pi x/L_{x}\right).

(33)

where the initial density and temperature is $\left(\rho_{0},T_{0}\right)=\left(1.0,1.0\right)$ , the perturbation amplitude $A=1e-6$ and $u_{0}=0$ . The simulation domain is $L_{x}\times 1$ (with $L_{x}=1600$ ). The shear viscosity is measured by monitoring the time evolution of the maximum velocity and fitting an exponential function to it, i.e.

u_{x}^{\rm max}(t)\propto\exp{\left(-\frac{4\pi^{2}\nu}{L_{x}^{2}}t\right)}.

(34)

Note that convergence studies were conducted as preliminary runs and all results reported here correspond to converged simulations in time and space. The results from simulations for different Mach numbers as obtained using the different splitting strategies are displayed in Fig. (2).

The results clearly show the invariance of the effective kinematic viscosity with respect to the Mach number.

The bulk viscosity is measured through the decay rate of sound waves. To that end a small perturbation is added to the uniform initial density field [28, 19]. Note that the smallness of the perturbation ensures that the study is in the linear regime, excluding effects such as non-linear steepening. For a discussion on non-linear acoustics readers are referred to studies such as [29]. In the linear regime it is readily shown that wave dynamics are governed by the linear lossy wave equation, see [30]. The flow is initialized as

\rho=\rho_{0}+A\sin\left(2\pi x/L_{x}\right),T=T_{0},u_{x}=0,u_{y}=0,

(35)

where $\rho_{0}=1.0$ with initial amplitude $A=1e-6$ and the density perturbation is $\rho^{\prime}=\rho-\rho_{0}$ . The initial temperature is $T_{0}=1.0$ . The decay of energy $E(t)=u_{x}^{2}+u_{y}^{2}-u_{0}^{2}+c_{s}^{2}\rho^{\prime 2}$ over time is supposed to fit an exponential function depending upon an effective viscosity $\nu_{e}$ which is a combination of kinematic shear and bulk viscosity

\nu_{e}=\frac{4}{3}\nu+\frac{\eta}{\rho},

(36)

defined by [28] as

E(t)\propto\exp{\left(-\frac{4\pi^{2}\nu_{e}}{L_{x}^{2}}t\right)}.

(37)

The resolution is identical to the previous test case. Measured bulk viscosities for the different energy splits are shown in Fig. (3).

In agreement with the multi-scale analysis, all three schemes recover the predicted bulk viscosity and are invariant with respect to the Mach number.

To assess the thermal diffusivity $\alpha$ , a different type of perturbation is introduced into the initial state of the system,

\rho=\rho_{0}+Asin\left(2\pi x/L_{x}\right),T=\rho_{0}T_{0}/\rho,u_{x}=0,u_{y}=u_{0},

(38)

where the initial density and temperature are $\left(\rho_{0},T_{0}\right)=\left(1.0,1.0\right)$ , the perturbation amplitude $A=1e-6$ . The thermal diffusivity is measured through the time evolution of maximum temperature difference $T^{\prime}=T_{0}-T$ ,

T^{\prime}(t)\propto\exp{\left(-\frac{4\pi^{2}\alpha}{L_{x}^{2}}t\right)}.

(39)

Figure 4 plots the measured thermal diffusivity at different Mach numbers compared with the intended values $\alpha=10^{-3}$ .

III.2.3 Measuring viscous heating

To assess the proper recovery of the combined effect of viscous heating and thermal conductivity we next consider the thermal Couette flow. The configuration consists of a 2-D channel of height $L$ bounded by a stagnant wall (at the bottom) at temperature $T_{0}$ (set to $T_{0}=1$ here) and a moving wall (at the top), at velocity $u_{0}$ and temperature $T_{1}$ (set to $T_{1}=1.005$ here). The analytical solution for the temperature distribution in the channel is [31],

\frac{T-T_{0}}{T_{1}-T_{0}}=\frac{y}{L}+\frac{{\rm Pr}{\rm Ec}}{2}\frac{y}{L}\left(1-\frac{y}{L}\right),

(40)

where the Eckert number is defined as ${\rm Ec}=u_{0}^{2}/(c_{v}+r)\Delta T$ with $\Delta T=T_{1}-T_{0}$ . A first set of simulations was conducted at ${\rm Ec}=8$ using all splitting schemes, with (parameters). The results are displayed in Fig. 5.

Simulations with these set of parameters showed very good agreement with the analytical solution. To better show the effect of non-recovery of the fourth-order moment in the non-translational energy split a second set of simulation was conducted at ${\rm Ec}=72$ with stronger velocity and temperature gradients. The results are displayed in Fig. 6.

It is observed that the D2Q16 lattice combined with the non-translational split is the only scheme that leads to visible deviations from the reference analytical solution. This is to be expected as, per discussions in the previous sections, the fourth-order moment of the first-distribution function is not recovered on this lattice, and therefore leads to errors in the diffusive heat flux. Note, however that due to the locally-adaptive nature of PonD, these deviations are considerably reduced as compared to static reference frame realizations like the LBM. Finally, to further prove the source of these deviations, a simulation was also conducted using a higher-order lattice, i.e. D2Q25, and using the same grid-size. Results are displayed in Fig. 7, and we can see that restoring the fourth-order moment of the first-distribution function corrects these deviations.

III.2.4 2-D test-case: High Mach astrophysical jet

As a demonstration of applicability to high Mach numbers, we consider the case of an astrophysical jet of Mach $80$ , without radiative cooling [13]. This case is an example of actual gas flows revealed from images of the Hubble Space Telescope and therefore is of high scientific interest. Furthermore, given the rather large variation in temperature and pressure it is a challenging configuration to run using any numerical scheme. Following the configuration in (), the jet is initialized as a semi-circular region on the left boundary :

\{\rm{Ma},\rho,T\}=\{80,5,0.005\}.

The rest of the computational domain is at rest with

\{\rm{Ma},\rho,T\}=\{0,50,0.0005\}.

Outflow boundary conditions are used all around except on the left boundary where the prescribed conditions are imposed. The simulation was conducted with a resolution of $2000\times 1000$ with the diameter of the initial jet, $D_{jet}=100.$ The simulation was conducted using the total energy split and D2Q16 only, as based on previous discussions it was selected as the more efficient and accurate realization. Figure 8 shows the density, pressure and temperature field at lattice $t=240$ .

IV Conclusions and discussion

With the advent of more reliable discrete kinetic schemes for compressible flow simulations, development and tuning of efficient and accurate realizations has become a central topic in the literature. In the present contribution we aimed at discussing the different possible realizations in the context of double distribution function models with a focus on high speed flows. To that end multi-scale analyses of all schemes, considering internal degrees of freedom, were conducted. All energy splits were shown to recover consistent hydrodynamics in agreement with the single distribution function. Furthermore it was shown that one of the energy splits, i.e. non-translational internal, required higher-order quadratures to properly satisfy all equilibrium moments needed for the Navier-Stokes-Fourier equations. Simulations of the thermal Couette flow confirmed the issues with the fourth-order equilibrium moments, as for higher Eckert number the D2Q16 lattice led to non-negligible deviations which were not present for the D2Q25 lattice. The internal energy split on the other hand, while recovering proper hydrodynamics with fourth-order quadratures, brings in source terms that are non-local in both space and time. This impacts both computational load, memory requirements and also numerical properties. In our simulations, the added cost of the higher-order quadrature was approximately $\sim 60\%$ as compared to a fourth-order quadrature –note that this number will go up in 3-D simulations– while for the internal energy realization the source terms contributed to a $\sim 15\%$ increase in computation time as compared to the to total energy realization.

Acknowledgements.

This work was supported by European Research Council (ERC) Advanced Grant 834763-PonD. Computational resources at the Swiss National Super Computing Center CSCS were provided under the grant s1286.

References

Hosseini et al. [2024] S. A. Hosseini, P. Boivin, D. Thévenin, and I. Karlin, Lattice boltzmann methods for combustion applications, Progress in Energy and Combustion Science 102, 101140 (2024).
Philippi et al. [2006] P. C. Philippi, L. A. Hegele Jr, L. O. Dos Santos, and R. Surmas, From the continuous to the lattice boltzmann equation: The discretization problem and thermal models, Physical Review E 73, 056702 (2006).
Shan et al. [2006] X. Shan, X.-F. Yuan, and H. Chen, Kinetic theory representation of hydrodynamics: a way beyond the navier–stokes equation, Journal of Fluid Mechanics 550, 413 (2006).
Siebert et al. [2008] D. Siebert, L. Hegele Jr, and P. C. Philippi, Lattice boltzmann equation linear stability analysis: Thermal and athermal models, Physical Review E 77, 026707 (2008).
Frapolli et al. [2014] N. Frapolli, S. Chikatamarla, and I. Karlin, Multispeed entropic lattice boltzmann model for thermal flows, Physical Review E 90, 043306 (2014).
He et al. [1998] X. He, S. Chen, and G. D. Doolen, A novel thermal model for the lattice boltzmann method in incompressible limit, Journal of computational physics 146, 282 (1998).
Guo et al. [2007] Z. Guo, C. Zheng, B. Shi, and T. Zhao, Thermal lattice boltzmann equation for low mach number flows: Decoupling model, Physical Review E 75, 036704 (2007).
Karlin et al. [2013] I. Karlin, D. Sichau, and S. Chikatamarla, Consistent two-population lattice boltzmann model for thermal flows, Physical Review E 88, 063310 (2013).
Asinari and Karlin [2010] P. Asinari and I. V. Karlin, Quasiequilibrium lattice boltzmann models with tunable bulk viscosity for enhancing stability, Physical Review E 81, 016702 (2010).
Hosseini and Karlin [2023] S. Hosseini and I. Karlin, Lattice boltzmann for non-ideal fluids: Fundamentals and practice, Physics Reports 1030, 1 (2023).
Rykov [1975] V. Rykov, A model kinetic equation for a gas with rotational degrees of freedom, Fluid Dynamics 10, 959 (1975).
Kolluru et al. [2020] P. K. Kolluru, M. Atif, and S. Ansumali, Extended bgk model for diatomic gases, Journal of Computational Science 45, 101179 (2020).
Gorban and Karlin [1994] A. N. Gorban and I. V. Karlin, General approach to constructing models of the boltzmann equation, Physica A: Statistical Mechanics and its Applications 206, 401 (1994).
Ansumali et al. [2007] S. Ansumali, S. Arcidiacono, S. Chikatamarla, N. Prasianakis, A. Gorban, and I. Karlin, Quasi-equilibrium lattice boltzmann method, The European Physical Journal B 56, 135 (2007).
Dorschner et al. [2018] B. Dorschner, F. Bösch, and I. V. Karlin, Particles on demand for kinetic theory, Physical review letters 121, 130602 (2018).
Bhatnagar et al. [1954] P. L. Bhatnagar, E. P. Gross, and M. Krook, A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems, Physical review 94, 511 (1954).
Li et al. [2007] Q. Li, Y. He, Y. Wang, and W. Tao, Coupled double-distribution-function lattice boltzmann method for the compressible navier-stokes equations, Physical Review E 76, 056705 (2007).
Prasianakis and Karlin [2007] N. I. Prasianakis and I. V. Karlin, Lattice boltzmann method for thermal flow simulation on standard lattices, Physical Review E 76, 016702 (2007).
Hosseini et al. [2020] S. A. Hosseini, N. Darabiha, and D. Thévenin, Compressibility in lattice Boltzmann on standard stencils: effects of deviation from reference temperature, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 378, 20190399 (2020).
Frapolli et al. [2015] N. Frapolli, S. S. Chikatamarla, and I. V. Karlin, Entropic lattice boltzmann model for compressible flows, Physical Review E 92, 061301 (2015).
Kauf [2011] P. Kauf, Multi-scale approximation models for the Boltzmann equation (ETH Zurich, 2011).
Köllermeier [2013] J. Köllermeier, Hyperbolic approximation of kinetic equations using quadrature-based projection methods, Master’s thesis, RWTH Aachen University (2013).
Grad [1949] H. Grad, On the kinetic theory of rarefied gases, Communications on pure and applied mathematics 2, 331 (1949).
Zipunova et al. [2021] E. Zipunova, A. Perepelkina, A. Zakirov, and S. Khilkov, Regularization and the particles-on-demand method for the solution of the discrete Boltzmann equation, Journal of Computational Science 53, 101376 (2021).
Kallikounis et al. [2022] N. Kallikounis, B. Dorschner, and I. V. Karlin, Particles on demand for flows with strong discontinuities, Physical Review E 16, 015301 (2022).
Sawant et al. [2022] N. Sawant, B. Dorschner, and I. V. Karlin, Detonation modeling with the particles on demand method, AIP Advances 12, 075107 (2022).
Ji et al. [2024] Y. Ji, S. A. Hosseini, B. Dorschner, K. H. Luo, and I. V. Karlin, Eulerian discrete kinetic framework in comoving reference frame for hypersonic flows, Journal of Fluid Mechanics 983, A11 (2024).
Dellar [2001] P. J. Dellar, Bulk and shear viscosities in lattice boltzmann equations, Physical Review E 64, 031203 (2001).
Buick et al. [2000] J. Buick, C. Buckley, C. Greated, and J. Gilbert, Lattice Boltzmann BGK simulation of nonlinear sound waves: the development of a shock front, Journal of Physics A: Mathematical and General 33, 3917 (2000).
Lamb [1924] H. Lamb, Hydrodynamics (University Press, 1924).
Liepmann and Roshko [1957] H. W. Liepmann and A. Roshko, Elements of gasdynamics (Wiley, New York, 1957, 1957).

Appendix A Multi-scale analysis of double distribution function models

Let us consider the following general system of equations,

\mathcal{D}_{t}\{f,g\}=\frac{1}{\tau}\left(\{f^{\rm eq},g^{\rm eq}\}-\{f,g\}\right)+\{\mathcal{F},\mathcal{G}\},

(41)

where, $\mathcal{D}_{t}=\partial_{t}+\bm{v}\cdot\bm{\nabla}$ , $\mathcal{F}=0$ for all splitting approaches, and $\mathcal{G}$ is only non-zero for the internal energy splitting with,

\mathcal{G}=-f(\bm{v}-\bm{u})\cdot(\partial_{t}\bm{u}+\bm{v}\cdot\bm{\nabla}\bm{u}).

(42)

For the multi-scale analysis let us introduce the following parameters: characteristic flow velocity $\mathcal{U}$ , characteristic flow scale $\mathcal{L}$ , characteristic flow time $\mathcal{T}=\mathcal{L}/\mathcal{U}$ , characteristic density $\bar{\rho}$ , speed of sound of ideal gas $c_{s}=\sqrt{\gamma rT}$ . With these, the variables are reduced as follows (primes denote non-dimensional variables): time $t=\mathcal{T}t^{\prime}$ , space $\bm{x}=\mathcal{L}\bm{x}^{\prime}$ , flow velocity $\bm{u}=\mathcal{U}\bm{u}^{\prime}$ , particle velocity $\bm{v}=c_{s}\bm{v}^{\prime}$ , density $\rho$ = $\bar{\rho}\rho^{\prime}$ , distribution function $f=\bar{\rho}c_{s}^{-3}f^{\prime}$ . Furthermore, the following non-dimensional groups are introduced: Knudsen number ${\rm Kn}={\tau c_{s}}/{\mathcal{L}}$ and Mach number ${\rm Ma}={\mathcal{U}}/{c_{s}}$ . With this, the equations are rescaled as follows:

{\rm Kn}\,\mathcal{D}_{t}^{\prime}\{f^{\prime},g^{\prime}\}=\frac{1}{\tau^{\prime}}\left(\{f^{\rm eq^{\prime}},g^{\rm eq^{\prime}}\}-\{f^{\prime},g^{\prime}\}\right)+\{\mathcal{F}^{\prime},\mathcal{G}^{\prime}\},

(43)

Assuming ${\rm Kn}\sim\epsilon$ and dropping the primes for the sake of readability,

\epsilon\mathcal{D}_{t}\{f,g\}=\frac{1}{\tau}\left(\{f^{\rm eq},g^{\rm eq}\}-\{f,g\}\right)+\{\mathcal{F},\mathcal{G}\}.

(44)

Then introducing multi-scale expansions:

\{f,g\}=\{f^{(0)},g^{(0)}\}+\epsilon\{f^{(1)},g^{(1)}\}+\epsilon^{2}\{f^{(2)},g^{(2)}\}+O(\epsilon^{3}),

(45)

and,

\{\mathcal{F},\mathcal{G}\}=\epsilon\{\mathcal{F}^{(1)},\mathcal{G}^{(1)}\}+\epsilon^{2}\{\mathcal{F}^{(2)},\mathcal{G}^{(2)}\}+O(\epsilon^{3}),

(46)

the following equations are recovered at scales $\epsilon$ and $\epsilon^{2}$ :


$\displaystyle\epsilon$	$\displaystyle:\mathcal{D}_{t}^{(1)}\{f^{(0)},g^{(0)}\}=-\frac{1}{\tau}\{f^{(1)},g^{(1)}\}+\{\mathcal{F}^{(1)},\mathcal{G}^{(1)}\},$	(47a)
$\displaystyle\epsilon^{2}$	$\displaystyle:\partial_{t}^{(2)}\{f^{(0)},g^{(0)}\}+\mathcal{D}_{t}^{(1)}\{f^{(1)},g^{(1)}\}=-\frac{1}{\tau}\{f^{(2)},g^{(2)}\}$
	$\displaystyle+\{\mathcal{F}^{(2)},\mathcal{G}^{(2)}\},$	(47b)

with $\{f^{(0)},g^{(0)}\}=\{f^{\rm eq},g^{\rm eq}\}$ . Note that for this system, the solvability conditions are:


$\displaystyle\forall k>0,$	$\displaystyle\int_{\mathbb{R}^{D}}f^{(k)}d\bm{v}=0,$	(48a)
$\displaystyle\forall k>0,$	$\displaystyle\int_{\mathbb{R}^{D}}\bm{v}f^{(k)}d\bm{v}=0.$	(48b)

In addition, depending on the splitting, solvability conditions on energy are:

$\displaystyle\int_{\mathbb{R}^{D}}g^{(k)}d\bm{v}=0,\,$	$\displaystyle\forall k>0,$	(49)
$\displaystyle\int_{\mathbb{R}^{D}}g^{(k)}d\bm{v}=0,\,$	$\displaystyle\forall k>0,$	(50)
$\displaystyle\int_{\mathbb{R}^{D}}\left(g^{(k)}+\frac{\bm{v}^{2}}{2}f^{(k)}\right)d\bm{v}=0,\,$	$\displaystyle\forall k>0.$	(51)

Taking the moments, $\int\{f,\bm{v}f\}d\bm{v}$ , of the Chapman-Enskog-expanded equations at order $\epsilon$ :

	$\displaystyle\partial_{t}^{(1)}\rho+\bm{\nabla}\cdot\rho\bm{u}$	$\displaystyle=$	$\displaystyle 0,$		(52)
	$\displaystyle\partial_{t}^{(1)}\rho\bm{u}+\bm{\nabla}\cdot\rho\bm{u}\otimes\bm{u}+\bm{\nabla}\cdot p\bm{I}$	$\displaystyle=$	$\displaystyle 0.$		(53)

For the energy balance equation, taking these moments $\int_{\mathbb{R}^{D}}\{g,g+\frac{\bm{u}^{2}}{2}f,g+\frac{\bm{v}^{2}}{2}f\}d\bm{v}$ ,

\partial_{t}^{(1)}E+\bm{\nabla}\cdot E\bm{u}+\bm{\nabla}\cdot p\bm{u}=0.

(54)

Using the Euler level momentum balance and continuity equations, a balance equation for kinetic energy can be derived as,

\partial_{t}^{(1)}K+\bm{\nabla}\cdot K\bm{u}+\bm{u}\cdot\bm{\nabla}p=0,

(55)

which in turn can be used to derive a balance equation for internal energy,

\partial_{t}^{(1)}U+\bm{\nabla}\cdot U\bm{u}+p\bm{\nabla}\cdot\bm{u}=0,

(56)

and pressure,

\partial_{t}^{(1)}p+\bm{\nabla}\cdot p\bm{u}+\frac{r}{c_{v}}p\bm{\nabla}\cdot\bm{u}=0,

(57)

where

\frac{\partial U}{\partial T}=\rho c_{v}.

(58)

Going up one order in $\epsilon$ , at order $\epsilon^{2}$ the continuity equation is:

\partial_{t}^{(2)}\rho=0,

(59)

while for the momentum balance equation one has:

\partial_{t}^{(2)}\rho\bm{u}+\bm{\nabla}\left(\int_{\mathbb{R}^{D}}\bm{v}\otimes\bm{v}f^{(1)}\right)=0.

(60)

The second term can be further expanded using the first-order-in- $\epsilon$ equation as:

\int_{\mathbb{R}^{D}}\bm{v}\otimes\bm{v}f^{(1)}d\bm{v}=\\ -\tau\left[\partial_{t}^{(1)}\left(\rho\bm{u}\otimes\bm{u}+p\bm{I}\right)+\bm{\nabla}\cdot\rho\bm{u}\otimes\bm{u}\otimes\bm{u}\right.\\ \left.+\bm{\nabla}p\bm{u}+\bm{\nabla}p\bm{u}^{\dagger}+\bm{I}\bm{\nabla}\cdot p\bm{u}\right].

(61)

Next we expand the first term as:

\partial_{t}^{(1)}\left(\rho\bm{u}\otimes\bm{u}+p\bm{I}\right)=\bm{u}\otimes\partial_{t}^{(1)}\rho\bm{u}+\bm{u}\otimes\partial_{t}^{(1)}\rho\bm{u}^{\dagger}\\ -\bm{u}\otimes\bm{u}\partial_{t}^{(1)}\rho+\partial_{t}^{(1)}p\bm{I},

(62)

where we use,

\bm{u}\otimes\partial_{t}^{(1)}\rho\bm{u}=-\bm{u}\times\left[\bm{\nabla}\cdot\rho\bm{u}\otimes\bm{u}+\bm{\nabla}p\right],

(63)

and,

\bm{u}\otimes\bm{u}\partial_{t}^{(1)}\rho=-\bm{u}\otimes\bm{u}\bm{\nabla}\cdot\rho\bm{u},

(64)

to arrive at

\partial_{t}^{(1)}\left(\rho\bm{u}\otimes\bm{u}+p\bm{I}\right)=-\bm{\nabla}\cdot\rho\bm{u}\otimes\bm{u}\otimes\bm{u}-\bm{u}\otimes\bm{\nabla}p\\ -{\left(\bm{u}\otimes\bm{\nabla}p\right)}^{\dagger}+\partial_{t}^{(1)}p\bm{I},

(65)

where one can use Eq. (57) to get to,

\int\bm{v}\otimes\bm{v}f^{(1)}d\bm{v}=\\ -p\tau\left[\left(\bm{\nabla}\bm{u}+\bm{\nabla}\bm{u}^{\dagger}\right)-\frac{r}{c_{v}}\bm{\nabla}\cdot\bm{u}\bm{I}\right].

(66)

Plugging this final expression into the momentum balance equation at order $\epsilon^{2}$ ,

\partial_{t}^{(2)}\rho\bm{u}+\bm{\nabla}\cdot\bm{T}_{\rm NS}=0,

(67)

where

\bm{T}_{\rm NS}=\tau p\left[\bm{\nabla}\bm{u}+\bm{\nabla}\bm{u}^{\dagger}-\frac{2}{D}\bm{\nabla}\cdot\bm{u}\bm{I}\right]+\tau p\left(\frac{2}{D}-\frac{r}{c_{v}}\right)\bm{\nabla}\cdot\bm{u}\bm{I}.

(68)

Moving on to the energy balance equation at order $\epsilon^{2}$ , for the first split,

\partial_{t}^{(2)}E+\bm{\nabla}\cdot\left(\int_{\mathbb{R}^{D}}\bm{v}g^{(1)}d\bm{v}\right)=0,

(69)

where

\int_{\mathbb{R}^{D}}\bm{v}g^{(1)}d\bm{v}=-\tau\left[\partial_{t}^{(1)}\left(p+E\right)\bm{u}+\bm{\nabla}\cdot\left(\frac{p}{\rho}\left(E+p\right)\right.\right.\\ \left.\left.+p\bm{u}\times\bm{u}+\left(p+E\right)\bm{u}\times\bm{u}\right)\right].

(70)

To expand this equation, we multiply the pressure balance equation at Euler level by $\bm{u}$ and obtain the following balance equation using momentum and continuity balance,

\partial_{t}^{(1)}p\bm{u}+\bm{\nabla}\cdot p\bm{u}\otimes\bm{u}+\frac{p}{\rho}\bm{\nabla}p+\frac{r}{c_{v}}p\bm{u}\bm{\nabla}\cdot\bm{u}=0.

(71)

Using the same approach we can also derive the following additional balance equation,

\partial_{t}^{(1)}E\bm{u}+\bm{\nabla}\cdot\left(E+p\right)\bm{u}\otimes\bm{u}+\frac{E}{\rho}\bm{\nabla}p-p\bm{u}\cdot\bm{\nabla}\bm{u}=0.

(72)

Plugging these equations into Eq. (70) and after some algebra,

\int_{\mathbb{R}^{D}}\bm{v}g^{(1)}d\bm{v}=-\tau p\left(c_{v}+r\right)\bm{\nabla}T-\bm{u}\cdot\bm{T}_{\rm NS},

(73)

which leads to:

\partial_{t}^{(2)}E-\bm{\nabla}\cdot\left[\tau p\left(c_{v}+r\right)\bm{\nabla}T\right]+\bm{\nabla}\cdot\left(\bm{u}\cdot\bm{T}_{\rm NS}\right)=0.

(74)

For the second splitting approach integrating Eq. (47b),

\partial_{t}^{(2)}U+\bm{\nabla}\cdot\left[\int_{\mathbb{R}^{D}}\bm{v}g^{(1)}d\bm{v}\right]=\int_{\mathbb{R}^{D}}\mathcal{G}^{(2)}d\bm{v},

(75)

where using,

f^{(1)}=-\tau\left(\partial_{t}^{(1)}f^{(0)}+\bm{v}\cdot\bm{\nabla}f^{(0)}\right),

(76)

and,

g^{(1)}=-\tau\left(\partial_{t}^{(1)}g^{(0)}+\bm{v}\cdot\bm{\nabla}g^{(0)}-\mathcal{G}^{(1)}\right),

(77)

one arrives at,

\partial_{t}^{(2)}U-\bm{\nabla}\cdot\left[\tau p(c_{v}+r)\bm{\nabla}T\right]+\bm{T}_{\rm NS}:\bm{\nabla}\bm{u}=0.

(78)

For the third splitting approach,

\partial_{t}^{(2)}E+\bm{\nabla}\cdot\left[\int_{\mathbb{R}^{D}}\bm{v}\left(g^{(1)}+\frac{\bm{v}^{2}}{2}f^{(1)}\right)d\bm{v}\right]=0,

(79)

where using,

\int_{\mathbb{R}^{D}}\bm{v}\left(g^{(0)}+\frac{\bm{v}^{2}}{2}f^{(0)}\right)d\bm{v}=\left(E+p\right)\bm{u},

(80)

and,

\int_{\mathbb{R}^{D}}\bm{v}\otimes\bm{v}\left(g^{(0)}+\frac{\bm{v}^{2}}{2}f^{(0)}\right)d\bm{v}=\frac{p}{\rho}\left(E+p\right)+p\bm{u}\otimes\bm{u}\\ +\left(p+E\right)\bm{u}\otimes\bm{u},

(81)

the same equation as in Eq. (73) is recovered, which then leads to the energy balance of Eq. (74).