Comprehensive comparison of collision models in the lattice Boltzmann framework: Theoretical investigations

TL;DR

This paper provides a unified mathematical framework for various collision models in the lattice Boltzmann method, comparing their mathematical properties and impacts on simulation accuracy across multiple lattice configurations.

Contribution

It introduces a common formalism for diverse collision models, analyzes their relationships, and evaluates the effects of different polynomial bases and regularization on LBMs.

Findings

Unified formalism for collision models in LBM

Impact of polynomial basis and regularization on accuracy

Links established between different collision models

Abstract

Over the last decades, several types of collision models have been proposed to extend the validity domain of the lattice Boltzmann method (LBM), each of them being introduced in its own formalism. The present article proposes a formalism that describes all these methods within a common mathematical framework, and in this way allows us to draw direct links between them. Here, the focus is put on single and multirelaxation time collision models in either their raw moment, central moment, cumulant or regularized form. In parallel with that, several bases (non orthogonal, orthogonal, Hermite) are considered for the polynomial expansion of populations. General relationships between moments are first derived to understand how moment spaces are related to each other. In addition, a review of collision models further sheds light on collision models that can be rewritten in a linear matrix form. More quantitative mathematical studies are then carried out by comparing explicit expressions for the post collision populations. Thanks to this, it is possible to deduce the impact of both the polynomial basis (raw, Hermite, central, central Hermite, cumulant) and the inclusion of regularization steps on isothermal LBMs. Extensive results are provided for the D1Q3, D2Q9, and D3Q27 lattices, the latter being further extended to the D3Q19 velocity discretization. Links with the most common two and multirelaxation time collision models are also provided for the sake of completeness. The present work ends by emphasizing the importance of an accurate representation of the equilibrium state, independently of the choice of moment space. As an addition to the theoretical purpose of the present article, general instructions are provided to help the reader with the implementation of the most complicated collision models.

Tables

Table 1. Table 1: Can X recover the behavior of Y by using a particular set of relaxation frequencies? Yes (✓), Partially ( ∼ similar-to \bm{\sim} ), No (❏). Collision models considered in this summary are based on: raw moments (RM), Hermite moments (HM), central moments (CM), central Hermite moments (CHM), both projection based (PR) and recursive regularization (RR) steps, cumulants (K). MRT and TRT models corresponds to collision models expressed within an orthogonal basis derived from the Gram-Schmidt orthogonalization procedure. RR ∗ superscript RR \text{RR}^{*} is the multirelaxation time version of the RR collision model. Only links that do not impact the resulting physics ( ω ν ≠ 1 subscript 𝜔 𝜈 1 \omega_{\nu}\neq 1 and ω ν b ≠ 1 subscript 𝜔 subscript 𝜈 𝑏 1 \omega_{\nu_{b}}\neq 1 ) are reported here. Furthermore, a collision model X will be considered to partially recover the behavior of another one Y, if the condition to enforce Δ ​ f ∗ = 0 Δ superscript 𝑓 0 \Delta f^{*}=0 implies that at least of one of ω n Y superscript subscript 𝜔 𝑛 𝑌 \omega_{n}^{Y} cannot be chosen freely anymore, with the most constraining case being the reduction to a single relaxation time approach. Eventually, all reported results are valid for both the D2Q9 and the D3Q27 lattices. They further remain valid for the D3Q19 lattice if and only if one compares collision models relying on the same equilibrium state.

	BGK	MRT	TRT	RM	HM	CM	CHM	PR	RR	${\text{RR}}^{\ast }$	K
BGK	–	$\sim$	$\sim$	$\sim$	$\sim$	$\sim$	$\sim$	❏	❏	❏	❏
MRT	✓	–	✓	$\sim$	$\sim$	$\sim$	$\sim$	❏	❏	❏	❏
TRT	✓	$\sim$	–	$\sim$	$\sim$	$\sim$	$\sim$	❏	❏	❏	❏
RM	✓	$\sim$	✓	–	$\sim$	$\sim$	$\sim$	❏	❏	❏	❏
HM	✓	$\sim$	✓	$\sim$	–	$\sim$	$\sim$	✓	❏	❏	❏
CM	✓	$\sim$	$\sim$	$\sim$	$\sim$	–	$\sim$	❏	❏	❏	❏
CHM	✓	$\sim$	$\sim$	$\sim$	$\sim$	$\sim$	–	❏	✓	❏	❏
PR	❏	❏	❏	❏	$\sim$	❏	❏	–	❏	❏	❏
RR	❏	❏	❏	❏	❏	❏	$\sim$	❏	–	$\sim$	❏
${\text{RR}}^{\ast }$	❏	❏	❏	❏	❏	❏	❏	❏	✓	–	❏
K	❏	❏	❏	❏	❏	❏	❏	❏	❏	❏	–

Table 2. Table 2: Description of standard lattice structures of interest. Weights w i subscript 𝑤 𝑖 w_{i} related to velocity groups compose the right part of the table. Here, the cyclic permutation is implied, which means for example that ( 1 , 0 , 0 ) 1 0 0 (1,0,0) stands for all six possibilities ( ± 1 , 0 , 0 ) plus-or-minus 1 0 0 (\pm 1,0,0) , ( 0 , ± 1 , 0 ) 0 plus-or-minus 1 0 (0,\pm 1,0) and ( 0 , 0 , ± 1 ) 0 0 plus-or-minus 1 (0,0,\pm 1) in the 3D case, while only ( ± 1 , 0 ) plus-or-minus 1 0 (\pm 1,0) , ( 0 , ± 1 ) 0 plus-or-minus 1 (0,\pm 1) and ( ± 1 ) plus-or-minus 1 (\pm 1) are implied in the 2D and 1D cases respectively. All velocity discretizations share the same lattice constant c s = 1 / 3 subscript 𝑐 𝑠 1 3 c_{s}=1/\sqrt{3} .

Group	${𝝃}_i$	D1Q3	D2Q9	D3Q19	D3Q27
$1$	$(0,0,0)$	$2/3$	$4/9$	$1/3$	$8/27$
$2$	$(1,0,0)$	$1/6$	$1/9$	$1/18$	$2/27$
$3$	$(1,1,0)$		$1/36$	$1/36$	$1/54$
$4$	$(1,1,1)$				$1/216$

Equations

$$f_{i}=\sum_{n}c_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}i,}n}P_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}i,}n},$$

$$\mathcal{M}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X})=\displaystyle{\int}\exp(\xi{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X})f(\xi)\>\mathrm{d}\xi,$$

$$\Phi_{M}=\exp(\xi{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X})=\displaystyle{\sum_{n}}\dfrac{(\xi{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X})^{n}}{n!},$$

$$\mathcal{M}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X})=\displaystyle{\sum_{n}}\left(\displaystyle{\int}\xi^{n}f(\xi)\>\mathrm{d}\xi\right)\dfrac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{n}}{n!}=\displaystyle{\sum_{n}}M_{n}\dfrac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{n}}{n!},$$

$$M_n={\mathcal{M}}^{(n)}(0),$$

$$\mathcal{H}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X})=\displaystyle{\int}\exp(\xi{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}-c_{s}^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{2}/2)f(\xi)\>\mathrm{d}\xi,$$

$$H_0(\xi )$$

$$H_1(\xi )$$

$$H_2(\xi )$$

$$H_3(\xi )$$

$$H_4(\xi )$$

$$H_5(\xi )$$

$$H_6(\xi )$$

$$\Phi_{H}=\exp(\xi{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}-c_{s}^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{2}/2)=\Phi_{M}\exp(-c_{s}^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{2}/2).$$

$$\dfrac{\mathrm{d}^{(k)}}{\mathrm{d}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{k}}e^{-c_{s}^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{2}/2}|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}=0}=\left\{\begin{array}[]{c l}1&\text{if }k=0\\ (-c_{s}^{2})^{k/2}(k-1)!!&\text{if }k\text{ is even}\\ 0&\text{otherwise}\end{array}\right.$$

$$(k-1)!!=(k-1)\times (k-3)\times ...\times 3\times 1.$$

$$\mathcal{H}^{(n)}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X})=\sum_{k=0}^{n}\binom{n}{k}\bigg{(}\dfrac{\mathrm{d}^{(k)}}{\mathrm{d}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{k}}e^{-c_{s}^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{2}/2}\bigg{)}\mathcal{M}^{(n-k)}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}),$$

$$A_n={\mathcal{H}}^{(n)}(0).$$

$$\mathcal{M}^{(n)}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X})=\sum_{k=0}^{n}\binom{n}{k}\bigg{(}\dfrac{\mathrm{d}^{(k)}}{\mathrm{d}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{k}}e^{c_{s}^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{2}/2}\bigg{)}\mathcal{H}^{(n-k)}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}),$$

$$\dfrac{\mathrm{d}^{(k)}}{\mathrm{d}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{k}}e^{c_{s}^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{2}/2}|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}=0}=\left\{\begin{array}[]{c l}1&\text{if }k=0\\ (c_{s}^{2})^{k/2}(k-1)!!&\text{if }k\text{ is even}\\ 0&\text{otherwise}\end{array}\right.$$

$$\displaystyle\widetilde{\mathcal{M}}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X})$$

$$\displaystyle=\displaystyle{\sum_{n}}\left[\displaystyle{\int}(\xi-M_{1})^{n}f(\xi)\>\mathrm{d}\xi\right]\dfrac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{n}}{n!}$$

$$\displaystyle=\displaystyle{\sum_{n}}\widetilde{M}_{n}\dfrac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{n}}{n!}.$$

$$\widetilde{\mathcal{M}}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X})=\exp(-M_{1}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X})\mathcal{M}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}).$$

$$\widetilde{\mathcal{M}}^{(n)}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X})=\sum_{k=0}^{n}\binom{n}{k}(-M_{1})^{n-k}\mathcal{M}^{(k)}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}).$$

$${\tilde{M}}_n=\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)(-M_1{)}^{n-k}{M}_k.$$

$$\mathcal{M}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X})=\exp(M_{1}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X})\widetilde{\mathcal{M}}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}).$$

$$M_n=\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)(M_1{)}^{n-k}{\tilde{M}}_k.$$

$${\tilde{A}}_n=\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)(-A_1{)}^{n-k}{A}_k,$$

$$A_n=\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)(A_1{)}^{n-k}{\tilde{A}}_k.$$

$$\mathcal{K}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X})=\ln[\mathcal{M}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X})].$$

$$\displaystyle\displaystyle{\sum_{n}}K_{n}\dfrac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{n}}{n!}$$

$$\displaystyle=\ln(M_{0})+\ln\bigg{(}1+\displaystyle{\sum_{p\geq 1}}\dfrac{M_{p}}{M_{0}}\dfrac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{p}}{p!}\bigg{)}$$

$$\displaystyle=\ln(M_{0})+\displaystyle{\sum_{q\geq 1}}\dfrac{(-1)^{q-1}}{q}\bigg{(}\displaystyle{\sum_{p\geq 1}}\dfrac{M_{p}}{M_{0}}\dfrac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{p}}{p!}\bigg{)}^{q}$$

$$\displaystyle=\displaystyle{\sum_{q\geq 1}}\dfrac{(-1)^{q-1}}{q}\bigg{(}\displaystyle{\sum_{p\geq 1}}M_{p}\dfrac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{p}}{p!}\bigg{)}^{q},$$

$$\ln(1+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}Y})=\displaystyle{\sum_{q\geq 1}}\dfrac{(-1)^{q-1}}{q}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}Y}^{q},$$

Full text

Lattice Boltzmann method, BGK, MRT, TRT, Regularization, Hermite moment, Raw moment, Central moment, Central Hermite moment, Cumulant.

I INTRODUCTION

During the past three decades, the lattice Boltzmann method (LBM) was proven to be of particular interest to the field of Computational Fluid Dynamics. This started showing its efficiency for the simulation of isothermal and weakly compressible flows, where it is now considered to be a particularly tough challenger for more conventional fluid solvers based on the solving of Navier-Stokes-Fourier equations Manoha and Caruelle (2015); Rumsey et al. (2018). In the mean time, the validity domain of LBMs was extended to a very large set of phenomena which are not restricted anymore to Fluid Mechanics Succi (2018). Despite these successes, severe stability issues were also encountered when simulating high Reynolds number flows with the BGK-LBM Ricot et al. (2009); Marié and Gloerfelt (2017). A great number of collision models have been proposed to circumvent this issue. These models are either based on static or dynamic single (multiple) relaxation time(s).

Most common LBMs based on dynamic relaxation times can be divided into two categories: (1) Large Eddy simulation (LES) based LBMs, and (2) Entropic LBMs (ELBMs). LES-LBMs rely on the well known fact that underresolved turbulent scales must be accounted for through the use of a subgrid scale model Sagaut (2006, 2010). This model mainly consists in mimicking the behavior of small scales through the dissipation of structures at the grid cut off size. This can be done modifying the relaxation time to take into account the additional eddy viscosity which usually scales as the strain rate tensor.

Regarding ELBMs, these models ensure the H-theorem to be valid after the velocity discretization of the Boltzmann equation Karlin et al. (1998); Boghosian et al. (2001). This is done solving a minimization problem at each grid point and time step. This leads to a variable relaxation time that locally self-adjusts to the flow. Hence a non constant dynamic viscosity is obtained with ELBMs, especially when underresolved mesh grids are used for the simulation of high Reynolds (turbulent) flows. In that sense, ELBMs seem very similar to LES-LBMs. In fact, Malaspinas et al. proved the behavior of the ELBM was sharing similarities with the standard Smagorinsky subgrid scale model Sagaut (2006) since its dynamic viscosity also scales as the strain rate tensor, and its value tends towards zero when the resolution of the mesh grid is increased Malaspinas et al. (2008). Nevertheless, solving a minimization problem at each grid node and for each time step induces a nonnegligible extra CPU cost, which is higher than for subgrid scale models, while sharing an equivalent accuracy Hamdi et al. (2018).

Hence, one may prefer to rely on KBC models to improve the numerical stability of LBMs at a low CPU cost. These collision models are based on an approximation to the minimization problem, and use an analytic formula for the computation of the dynamic relaxation time Karlin et al. (2014); Bösch et al. (2015, 2015); Dorschner et al. (2016, 2018); Flint and Vahala (2018); Bösch et al. (2018). They can further decouple the relaxation of shear modes from acoustic and ghost modes, hence, freeing themselves from the generation of spurious vortices induced by the use of a varying shear viscosity in underresolved conditions. Nevertheless, it is worth noting that some theoretical work remains to be done in order to properly define the validity range of the approximated minimization problem Mattila et al. (2015). Eventually, one can further find in the literature stabilization techniques based on other kinds of dissipation control Brownlee et al. (2006, 2007, 2008); Gorban and Packwood (2014).

All the aforementioned collision models, despite showing a great stability gain as compared to the BGK operator, cannot be easily studied from a mathematical viewpoint due to the dynamic computation of the relaxation time.

Regarding LBMs based on static relaxation times, d’Humières originally proposed to apply the collision step within the moment space to improve the numerical stability of the resulting LBM through the increase of the number of free parameters. This ended up in a multirelaxation time (MRT) collision model for the D3Q19 velocity discretization d’Humières (1992). With this approach, several relaxation times controlled the relaxation process of moments toward their respective equilibrium state. This model was further popularized by the 2D formulation (D2Q9) of Lallemand and Luo Lallemand and Luo (2000). In their work, stability gains were obtained through fine tuning of these free parameters via a linear stability analysis. Nevertheless, finding optimal values for all relaxation times rapidly becomes a tedious task, especially in the 3D case. This is why two relaxation time (TRT) models were then proposed by Ginzburg et al. Ginzburg et al. (2008a, b). In their most simple and naive form, these models can be related to both Lallemand’s and Luo’s work in 2D, and d’Humière et al. 3D extension d’Humières (2002), where a particular distinction is made between the relaxation times of even ($\tau^{+}$) and odd moments ($\tau^{-}$). While the former is related to the dynamic viscosity, the second one is a free parameter. Using very thorough mathematical derivations, general relationships between $\tau^{+}$ and $\tau^{-}$ were proposed through the so called ‘magic number’ $\Lambda$ to improve both the accuracy and the stability of these TRT-LBMs d’Humières and Ginzburg (2009); Ginzburg et al. (2010); Kuzmin et al. (2011). Both types of collision model were shown to be of particular interest for the simulation of various types of phenomena Hao and Cheng (2009); Xu and Sagaut (2011); Wang et al. (2013); Servan-Camas and Tsai (2009); Kuzmin et al. (2011); Silva and Semiao (2017).

To further improve the numerical stability of MRT-LBMs, several authors proposed to perform the collision step in the moment space associated to the comoving reference frame Geier et al. (2006, 2009); Geller et al. (2013); Lycett-Brown and Luo (2014); Dubois et al. (2015a); Ning et al. (2016); Fei and Luo (2017); De Rosis (2017a); Fei et al. (2018); De Rosis and Luo (2019); De Rosis et al. (2019). These models were originally named as cascaded LBM by Geier et al. Geier et al. (2006) due to the construction of high order moments from lower order ones in a cascade like way. Despite an improved stability for high Reynolds flow simulations, this model was still suffering from Galilean invariance defects that were partially corrected later by the same authors Geier et al. (2009). Nonetheless, the complex implementation of the cascaded process made it difficult to extend this collision model to both a broader range of physics, and to other kinds of lattices. This deficiency was overcome, for example, by Lycett-Brown and Luo who simplified the derivation of the cascaded mechanism, and applied it to the simulation of multiphase flows Lycett-Brown and Luo (2014, 2016). In the mean time, Dubois et al. properly reformulated the cascaded model in terms of collisions performed in the central moment (CM) space Dubois et al. (2015a). The latter being well known from the point of view of Statistics, more solid mathematical foundations were then provided to these models. The same authors also quantified the stability gain obtained with this kind of collision models through a linear stability analysis Dubois et al. (2015b). In their work, they further confirmed the superiority, in terms of linear stability, of the tensor product basis as compared to the more standard ones used in common MRT models. All of these led to the extension of CM-LBMs to very different fields of research such as shallow water and magnetohydrodynamics among others De Rosis (2017b, a); De Rosis et al. (2018).

The last version of these MRT-LBMs is based on a collision step occurring this time in the cumulant space. By definition, cumulants are another kind of statistical quantities, such as raw and central moments, that allow the description of both continuous and discrete probability distribution functions Cook (1951). They were originally used for the simulation of gas dynamics by Seeger et al. Seeger and Hoffmann (2000, 2002); Seeger (2003). Their method was based on solving equations obtained by taking cumulants of the Boltzmann equation. In the LBM context, Geier et al. did confirm the stability improvement induced by the collision model based on cumulants instead of raw moments Geier et al. (2015). Nonetheless, it was done in an extremely reduced configuration where the diffusive scaling was adopted, thus reducing the validity of the comparison to flows governed by the incompressible Navier-Stokes equations Junk et al. (2005); Krüger et al. (2017). Hence, no information regarding the generation and propagation of acoustic waves can be deduced from either their work or from those that followed Geier et al. (2017a, b); Geier and Pasquali (2018). More surprisingly, no comparison with central moments was provided for the most complicated case, i.e, the study of the flow around a sphere. As a consequence, it is still not clear if cumulant methods are (or not) always more stable than CM-LBMs at the time of writing.

In parallel to MRT-like collision operators, regularization steps were proposed to increase the stability of the BGK-LBM in the low viscosity limit. This kind of LBMs was originally proposed by Ladd et al. to reduce the memory consumption Ladd and Verberg (2001). Nevertheless, it was rapidly shown to be able to improve the numerical stability by filtering out nonhydrodynamic modes (Latt and Chopard Latt and Chopard (2005, 2006)), and ensuring the rotational invariance of the numerical scheme (Chen et al. Chen et al. (2006)). This model was further extended to high order velocity discretizations by Zhang et al. Zhang et al. (2006). A recursive version of these models was recently derived for standard Malaspinas (2015), and high order LBMs Coreixas et al. (2017). This recursive approach was proven to increase the numerical stability of LBMs as compared to the original regularized collision model Coreixas (2018). It was further shown to be competitive against state of the art Navier-Stokes and LBM solvers for the simulation of jet noise Brogi et al. (2017). This confirmed its capability for high fidelity simulations, including computational aeroacoustics of moderate to high Reynods number flows. In the mean time, Mattila et al. worked on a regularization step operating in the comoving reference frame Mattila et al. (2017). They noticed that by neglecting particular central Hermite coefficients, they were able to recover the recursive formulation of Malaspinas. These coefficients being related to diffusive phenoma, they were supposed to be negligible for high Reynolds number flows. In the end, all the aforementioned regularized models were also employed to simulate a large panel of distinct phenomena Niu et al. (2007); Viggen (2013); Wang et al. (2015); Coreixas et al. (2017); Brogi et al. (2017); Ba et al. (2018); Montessori et al. (2018).

Despite a wide variety of collision models, only very few comprehensive and consistent comparative studies have been conducted in the past. Most of them are restricted to less than four collision models, or are only based on numerical test cases, and contain only very few sound comments on their respective theoretical discrepancies. The former issue most likely comes from the fact that each type of collision model is usually introduced in its own framework. Consequently, it is complicated to express all the collision models within the very same framework in order to eventually compare them to each other from a general viewpoint. Regarding the lack of theoretical comparisons, it can be explained by the nonnegligible mathematical background required to properly derive links between collision models. Nonetheless, rigorous comparisons relying on different tools are also available in the literature. As an example, several authors used linear stability analysis to quantify the numerical stability of collision models Lallemand and Luo (2000); Siebert et al. (2008); Adhikari and Succi (2008); Marié et al. (2009); Ricot et al. (2009); Kuzmin et al. (2011); Xu and Sagaut (2011); Dubois et al. (2015b); Hosseini et al. (2017); Chávez-Modena et al. (2018); Coreixas (2018); Wissocq et al. (2019). In addition, it is worth noting that studies dedicated to the derivation of links between several kinds of LBMs from a (more or less) theoretical viewpoint can also be found in the literature. Among them, one may refer to Refs. Asinari (2008); Ginzburg et al. (2008a); Gorban and Packwood (2014); Geier et al. (2015); Dubois et al. (2015a); Coreixas et al. (2017) and therein references.

Knowing all of this, the present work proposes a formalism that describes all static methods within a common mathematical framework, and in this way allows us to draw direct links between them. This is done considering the LBM in its general form, i.e, with no restriction to the incompressible regime. In addition, it is intended to answer several questions that arise from the literature. For instance, one can wonder if all collision models can be rewritten in a linear matrix form as originally proposed by Higuera et al. Higuera and Jimenez (1989), and further discussed in the context of LBMs by several authors Benzi et al. (1992); d’Humières (1992). If this assumption were to be true, then it would be possible to easily switch from one collision model to another through (1) matrix products and (2) adjustments of relaxation frequencies. One can also wonder if the equilibrium state does or not depend on the Galilean invariance properties of the collision model, and consequently, if the resulting macroscopic behavior of the LBM also depends on the collision model Geier et al. (2006, 2009, 2015, 2017a). As a consequence, the present work can be considered as the first (theoretical) stage of a large project that intends to compare collisions models from both a theoretical and a numerical viewpoint.

The rest of the paper is organized as follows. Recalls on Statistics are first provided in Sec. II. Thanks to them, the way all kinds of statistical quantities are linked is determined by relating their moment generating function to each other. All collision models are then thoroughly reviewed, and potential links between the models that have been previously pointed out in the literature are recalled (Sec. III). In this way, a first picture of all the families of collision models is drawn. In the following sections (IV, V and VI), relationships derived between moments are used to carry out, in a straightforward way, the comparative study of post collision populations for lattices of increasing complexity (D1Q3, D2Q9, D3Q27 and D3Q19). All algebraic manipulations are detailed for the most simple lattices, hence allowing the reader to readily extend the present approach to any kind of LBMs. This comparative study ends with a reflection on macroscopic equations resulting from the LBM. General conclusions regarding the main discrepancies and similarities between collision models are finally made in Sec. VII.

For the sake of completeness, several appendices are also provided to help the interested reader properly understand all concepts addressed in the present work. Lattice structures and partial Bell polynomials are described in Apps. A and B. Relationships between one-, two- and three-dimensional statistical quantities are gathered in Apps. C, D and E respectively. Linear transformations leading to the matrix form of collision models are compiled in App. F for the D2Q9 lattice. Explanations regarding the way to build them are also provided, which allows a straightforward calculation of these matrices for any type of lattice. The D3Q27 formulations of all types of populations are recalled in App. G. It is also explained how to easily derive their D3Q19 counterparts. Eventually, the most general form of equilibrium states is presented in App. H.

To help the reader with the implementation of the most complicated collision models, general instructions are finally provided as Supplemental material \bibnote[SupMat]See Supplemental Material at [COREIXAS_PRE_2019_Theo_SupMat.pdf] for instructions regarding the implementation of most collision models considered in the present work.. More specifically, pseudo codes, describing the computation of each type of moment, as well as, equilibrium and post collision populations, are supplied for both the D2Q9 \bibnote[SupMatQ9]See Supplemental Material for instructions regarding the implementation of D2Q9 collision models at [Comprehensive_Collision_Models_D2Q9.txt]. and the D3Q27 \bibnote[SupMatQ27]See Supplemental Material for instructions regarding the implementation of D3Q27 collision models at [Comprehensive_Collision_Models_D3Q27.txt]. formulations.

II Theoretical background on Statistics

II.1 Motivations

As it will be shown in Sec. IV.1, the choice of the moment space directly impacts the derivation of LBMs through their populations

[TABLE]

where c_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}i,}n} is a coefficient related to the moment P_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}i,}n}, both yet to be defined. To properly quantify differences between LBMs, it is then essential to understand how all kinds of moments are derived, and how they are linked to each other. In the present work, it is proposed to relate families of moments through their sole common feature, namely, their generating function. This method further simplifies their calculation in any $D$-dimension, $D$ being the number of physical dimensions.

II.2 Moment generating function

To properly understand how relationships between raw, Hermite, central, central Hermite moments and cumulants are derived, it is necessary to introduce the notion of (raw) moment generating function (MGF) $\mathcal{M}$ of a probability distribution function $f$,

[TABLE]

By expanding the generating function of monomials in power series,

[TABLE]

the MGF can be directly linked to raw moments $M_{n}$ of $f$ as follows,

[TABLE]

where $n!=n\times(n-1)\times...\times 2\times 1$. Hence, (raw) moments of order $n$ are obtained through the partial derivatives of the MGF (2) evaluated at {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}=0,

[TABLE]

$\mathcal{M}^{(n)}$ being the $n$th derivative of $\mathcal{M}$ with respect to {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}. In the rest of the paper, $M_{0}=1$ will be assumed. This corresponds to the assumption of normalized moments, as already used by several authors in the lattice Boltzmann community Karlin and Asinari (2010); Lycett-Brown and Luo (2014); Bösch et al. (2015).

II.3 Hermite moment generating function

Hermite moments are defined through their generating function,

[TABLE]

where \Phi_{H}=\exp(\xi{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}-c_{s}^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{2}/2) is the generating function of Hermite polynomials as they are usually defined in the lattice Boltzmann framework, i.e, taking into account the lattice constant $c_{s}$ Coreixas et al. (2017). As a reminder, Hermite polynomials up to $n=6$ are

[TABLE]

Since (Hermite) polynomials are linear combinations of monomials, then relationships between Hermite and raw moments will be the same as those between their polynomial counterparts. To derive these relationships, let us start noticing that both polynomial generating functions are linked as follows,

[TABLE]

Since these polynomials are obtained by deriving $n$ times $\Phi_{H}$ ($\Phi_{M}$) about {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}=0, the only expression that needs to be evaluated is the $n$th derivative of \exp(-c_{s}^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{2}/2). While its general formulation is quite complicated, it can easily be shown by induction that

[TABLE]

with ‘!!’ standing for the double factorial, i.e.,

[TABLE]

Hence,

[TABLE]

with $\binom{n}{k}=n!/k!(n-k)!$ being the binomial coefficient, and where Hermite moments are defined as

[TABLE]

The inversion formula is simply obtained changing $-c_{s}^{2}$ to $+c_{s}^{2}$ in Eq. (7),

[TABLE]

where now,

[TABLE]

II.4 Central moment generating functions

Regarding central moments $\widetilde{M}_{n}$, their generating function reads as

[TABLE]

A simple way to link raw and central moments is to express the central moment generating function as

[TABLE]

Using the binomial formula, the $n$th derivative of \widetilde{\mathcal{M}}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}) then becomes

[TABLE]

Eventually, evaluating the $n$th derivative about {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}=0 leads to

[TABLE]

The inversion formula is simply obtained noticing that

[TABLE]

Following the very same steps as before, one ends up with

[TABLE]

Since Hermite polynomials are linear combinations of monomials, the above derivation can be applied to link raw and central moments in the Hermite expansion framework,

[TABLE]

where $\widetilde{A}_{n}$ are central Hermite moments, and $A_{1}=M_{1}$ since $H_{1}(\xi)=\xi$. Finally, the inversion formula is

[TABLE]

II.5 Cumulant generating function

When it comes to cumulants, their generating function $\mathcal{K}$ is defined as Cook (1951)

[TABLE]

Its power series expansion allows us to link cumulants and raw moments as follows,

[TABLE]

where the series expansion of the logarithmic part was computed using

[TABLE]

with {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}Y}=\sum_{p\geq 1}M_{p}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}^{p}/p!.

To determine the inversion formula, the starting point is

[TABLE]

Taking its power series expansion, the following expression is derived,

[TABLE]

In their current formulation, it is quite difficult to use Eqs. (15) and (17). This is why two ways of simplifying these formulas, as originally proposed by Kendall Kendall (1940), are presented hereafter:

‘Handmade’ expansion of the double sum: It consists on expanding the right hand side of these formulas assuming truncated sums. Relationships are then obtained identifying polynomials coefficients of $t^{n}$ with those of $t^{p+q}$, where $p+q=n$ (uniqueness of the power series expansion). In other words, by considering $p\leq n$ and $q\leq n$, one just needs to collect all coefficients of $t^{n}$. Even if this method is rather general, it rapidly leads to lengthy calculations even for small values of $n$.

Partial Bell polynomials: They are used in Combinatorial Mathematics to study set partitions Bell (1934); Withers (1994). These polynomials are defined as follows,

[TABLE]

where $B_{n,k}$ are partial Bell polynomials. They can be individually computed thanks to the recursive formula

[TABLE]

with

[TABLE]

Expressions up to $n=6$, and further details concerning these polynomials, are given in App. B. Injecting Eq. (18) into Eq. (17), the $n$th raw moment is simply expressed as

[TABLE]

Likewise, the $n$th cumulant is obtained via the inversion formula of partial Bell polynomials Kendall (1940),

[TABLE]

To obtain formulas between cumulants and other statistical quantities, one simply needs to switch between the quantity of interest to its raw moment counterpart, and then use the above formulas. Relationships between all kinds of moments, up to $n=6$, are detailed in App. C.

II.6 Multivariate extensions

In the $D$-dimensional case, the multivariate generating moment function is now defined as

[TABLE]

where \bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}}=({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}_{1},...,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}_{D}) and $\bm{\xi}=(\xi_{1},...,\xi_{D})$ are now vectors. As in the univariate case, moments are obtained through successive differentiations of \mathcal{M}(\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}}) about \bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}X}}=\bm{0}. For the most complicated case ($D=3$), this leads to

[TABLE]

$M_{pqr}$ is a raw moment of order $n=p+q+r$. $\mathcal{M}^{(p,q,r)}$ is the $n$th derivative of $\mathcal{M}$, where $p$, $q$, $r$ are the number of differentiations with respect to $\xi_{x}$, $\xi_{y}$ and $\xi_{z}$.

While rederiving all relationships between statistical quantities using Eq. (21) is relatively straightforward, the computation of multivariate cumulants rapidly becomes lengthy. It is proposed hereafter to rely on a rather elegant way to extend the derivation of cumulant formulas to the $D$-dimensional case. The most general way relies on advanced knowledge of set partitioning. Here, it is preferred to use a simpler, yet rigorous, method based on a particular kind of differential operator. This approach was originally proposed by Kendall Kendall (1940), and explained in more details by Cook Cook (1951). Let us start from

[TABLE]

Each component of $M_{n}$ or $K_{n}$ depends on $\xi_{x}^{n}$ in such a way that Eq. (23) can be rewritten as

[TABLE]

Applying the operator $\Delta_{y}^{x}\equiv\xi_{y}\partial/\partial\xi_{x}$ to the variable of $M(\xi_{x}^{3})$ leads to Kendall (1940); Cook (1951),

[TABLE]

When it comes to the RHS terms of Eq. (24), the following formulas are obtained,

[TABLE]

Using Eqs. (25)-(28), $M_{21}$ is eventually linked to cumulants through

[TABLE]

Applying $\Delta_{y}^{x}$ a second time leads to

[TABLE]

while applying it a third time further gives the expression

[TABLE]

The formula for $M_{30}$ is simply obtained noticing that $\xi_{x}^{3}=\xi_{x}^{3}\xi_{y}^{0}$, and then taking [math] as second subscript for all terms in Eq. (23),

[TABLE]

As a way to confirm the validity of the above third order bivariate moments, one just has to sum all subscripts for each moment. If the bivariate formula is correct, then one should obtain the univariate formula for $M_{3}$ as a result.

To derive trivariate ($D=3$) formulas, two operators of the same kind are required, namely, $\Delta_{z}^{x}$ and $\Delta_{z}^{y}$. Using them, it is straightforward to obtain trivariate expressions, such as,

[TABLE]

and

[TABLE]

Once again, one can check the validity of the above formulas summing all subscripts for each moment, and comparing the result with the univariate case, i.e, $M_{3}$. While Kendall’s diffential operators might appear rather empirical, it is important to know that it relies on rigorous mathematical derivations originating from the multivariate version of the Taylor series expansion of generating functions Kendall (1940); Cook (1951). Using this method, only univariate formulas up to $n=4$ ($n=6$) are needed to derive bivariate (trivariate) statistical quantities of interest for the D2Q9 (D3Q19, D3Q27) lattice. Univariate formulas are summarized in App. C, whereas expressions needed for the implementation of the D2Q9-, D3Q19- and D3Q27-LBMs are compiled in Apps. D and E.

II.7 Application to the normal distribution

In Kinetic Theory, the normal (or Gaussian) distribution plays a major role in the modeling of the collision term. Indeed, the relaxation process is usually based on the BGK assumption which says that collisions make distribution functions tend towards their equilibrium state Bhatnagar et al. (1954). The latter follows a normal distribution function and is usually named as the Maxwell-Boltzmann equilibrium (distribution) function. It reads as

[TABLE]

where $D$ is the number of physical dimensions, $\rho$ being the density, $\bm{u}$ the macroscopic velocity, $\bm{\xi}$ the mesoscopic velocity, $\theta=T/T_{0}$ the reduced temperature, and $c_{s}$ the isothermal (or Newtonian) sound speed. In the present work, only isothermal LBMs are considered, which implies that $T=T_{0}$ so that $\theta=1$.

Injecting Eq. (35) in the definition of raw, central, Hermite, central Hermite moments and cumulants, one can derive their equilibrium counterparts. Another way to obtain all equilibrium quantities is to start with equilibrium raw moments in either the univariate or the multivariate case, and then use corresponding formulas (Apps. C-E) to compute other families of moments.

Starting with one-dimensional equilibrium raw moments, one obtains up to $n=4$

[TABLE]

Related central moments are then

[TABLE]

while equilibrium Hermite moments read as,

[TABLE]

and equilibrium central Hermite moments are expressed as,

[TABLE]

Eventually, equilibrium cumulants are

[TABLE]

Central moments are usually preferred to their counterparts in the reference frame at rest because, for the former, the collision is applied in a moment space that is not impacted by Galilean invariance issues Geier et al. (2006, 2009); Ning et al. (2016). Nonetheless, lattice dependent defects still remain in $\widetilde{M}^{eq}_{2p}$ ($p>0$) through terms proportional to $c_{s}$ that are nonnegligible in the weakly compressible limit, since $c_{s}\gg u$ for the latter. Thus one may wonder if, in the isothermal case, it would not be wiser to apply the collision in the central Hermite moment space instead. In addition, one of the reason behind the use of cumulants is that they allow us to easily quantify deviations from a normal distribution since $K^{eq}_{n}=0$ for $n\geq 3$ McCullagh (2018). But once again, it seems more logical to rely on central Hermite moments since $\widetilde{A}^{eq}_{n}=0$ for $n\geq 1$ in the isothermal case. Nevertheless, one should keep in mind that $\widetilde{A}^{eq}_{n}\neq 0$ for $n\geq 2$ in the thermal case. Hence, cumulants may be a preferable choice for the latter case.

Multivariate formulations of all kinds of equilibrium moments are simply obtained thanks to the isotropy of the multivariate normal distribution. Hence, multivariate equilibrium moments are built through products of univariate ones.

Henceforth, only equilibrium moments computed via the Maxwell-Boltzmann distribution function will be considered. But it is also worth noting that depending of the physics of interest, other equilibrium distribution functions might be used, e.g., those derived from the principle of maximum entropy Kogan (1965); Pressé et al. (2013). Interestingly, the latter methodology was used in the lattice Boltzmann community for the derivation of exact Ansumali et al. (2003); Asinari and Karlin (2009) and approximated Frapolli et al. (2015); Atif et al. (2018) equilibrium distribution functions.

III Review of collision models

The LBM is a very specific numerical scheme employed to solve the lattice Boltzmann equation – a set of Boltzmann equations resulting from the discretization of the velocity space.

III.1 The lattice Boltzmann method

This method describes the space and time evolution of the velocity distribution function (VDF) $f_{i}(\bm{x},\bm{\xi}_{i},t)$, also referred to as ‘population’ in the present work. Roughly speaking, this quantity can be assimilated to the number of fictive particles at a point $(\bm{x},t)$ and characterized by a given velocity $\bm{\xi}_{i}$. Its space and time evolution is obtained following two (dimensionless) successive steps

[TABLE]

where $f^{*}_{i}$ are post collision VDFs. Eq. (36) corresponds to the streaming step that propagates $f_{i}$ to their neighboring nodes following the direction $\bm{\xi}_{i}$. Eq. (37) is the collision step that locally takes into account the rate of change of $f_{i}$ due to collisions, where $\Omega$ is the collision model in its general form, $f_{i}^{eq}$ is the equilibrium VDF, and $f_{i}^{neq}$ is the deviation from the equilibrium state. Altogether, they form the famous ‘Collide and Stream’ algorithm which shows both CPU time efficiency and accuracy Krüger et al. (2017).

Hereafter, the discussion will be restricted to the D2Q9 lattice. This is done for the sake of simplicity, and because it is sufficient to explain the differences between collision models of interest in the present work.

III.2 BGK collision model

In the particular case of the BGK collision operator,

[TABLE]

with $\omega_{\nu}=1/(\tau+1/2)$ the collision frequency, $\tau$ being the relaxation time, $\nu$ the kinematic viscosity, and $f^{eq,2}_{i}$ is the standard second order equilibrium VDF Qian et al. (1992); Shan et al. (2006). Using this collision model, severe stability issues have been reported in the zero viscosity limit Succi (2018). While one can improve the numerical stability of the LBM changing either the space or the time discretizations of the streaming and collision steps Lee et al. (2006); Surmas et al. (2009); Guo et al. (2011, 2013), it is usually preferred to derive new collision models introducing new forms for either $\Omega$ d’Humières (1992); Lallemand and Luo (2000); Ginzburg et al. (2008a); Geier et al. (2006); Dubois et al. (2015a); Geier et al. (2015) or $f_{i}^{neq}$ Latt and Chopard (2005, 2006); Malaspinas (2015); Coreixas et al. (2017). Doing so, both the efficiency and the accuracy can be preserved in most of the cases, and the stability of the resulting LBM is drastically improved.

III.3 Raw moment space

To circumvent the stability issue, it has first been proposed to apply the collision step in the moment space through a MRT collision model Lallemand and Luo (2000),

[TABLE]

where $\bm{S}$ is the collision matrix which may have a diagonal form depending on the choice of moments. $\bm{M}$ and $\bm{M}^{-1}$ are commonly defined as orthogonal matrices allowing us to move from the 9-dimensional velocity space

[TABLE]

to the 9-dimensional orthogonal moment space

[TABLE]

and vice versa. The polynomials basis related to this moment space is

[TABLE]

and it is built using the Gram-Schmidt orthogonalization procedure based on a unweighted scalar product.

In the present case, the collision matrix is defined as Lallemand and Luo (2000)

[TABLE]

Hence, the linear matrix form of $\bm{\Omega}$ allows the proper control of the relaxation process of each moment. This approach can drastically increase the numerical stability of the MRT-LBM, but this will depend on the choice of both the moment space, $\bm{M}$ and $\bm{M}^{-1}$ Dubois et al. (2015b), and the relaxation frequencies composing the collision matrix $\bm{S}$ Lallemand and Luo (2000). This last point becomes a huge drawback in either the 3D case or when dealing with high order LBMs, since the number of parameters that need to be fine tuned rapidly becomes very large.

To reduce the number of free parameters in a general way, the TRT-LBM was proposed Ginzburg et al. (2008a, b). In the context of the orthogonal basis $\mathcal{B}_{\mathrm{MRT}}$ (38), it relies on the following decomposition into a symmetric (even moments), and an antisymmetric (odd moments) sub-bases Krüger et al. (2017)

[TABLE]

with

[TABLE]

and

[TABLE]

Hence, the relaxation of even and odd moments is decoupled using two relaxation frequencies. In the above context, this leads to

[TABLE]

with $\omega^{+}=\omega_{\nu}$, $(1/\omega^{-}-1/2)=\Lambda/(1/\omega^{+}-1/2)$, and $\Lambda$ the ‘magic parameter’ that controls both accuracy and stability of the TRT-LBM d’Humières and Ginzburg (2009); Ginzburg et al. (2010); Kuzmin et al. (2011); Silva and Semiao (2017).

It is worth noting that, by relying on its original set of relaxation frequencies, the orthogonal MRT model is not well suited for aeroacoustic simulations. Indeed, the standard value of the relaxation frequency related to acoustic modes ($\omega_{e}$) is always far from $2$, hence leading to huge values of the bulk viscosity Lallemand and Luo (2000); d’Humières (2002). This further increases the stability of the corresponding LBM by overdissipating acoustic waves Marié et al. (2009); Xu and Sagaut (2011); Gendre et al. (2017). But from the theoretical viewpoint, the major defect of the above models is that they originally relied on polynomial bases including high order terms that are not compliant with the velocity discretization. Taking the example of $q_{x}$, its associated vector is

[TABLE]

and it reduces to

[TABLE]

due to the aliasing defects $\xi_{i,x}^{3}=\xi_{i,x}$ and $\xi_{i,y}^{3}=\xi_{i,y}$ of the D2Q9 lattice, where $\xi_{i,x},\xi_{i,y}\in\{0,\pm 1\}$. The same issue is also encountered with $M_{q_{y}}$ and $M_{\varepsilon}$. Unfortunately, this aliasing problem remains present for D3Q19 and D3Q27 lattices where it might be even more dominant. Nevertheless, both defects can easily be circumvented by using: (1) the correct value of the bulk viscosity to compute $\omega_{e}$, and (2) a polynomial basis that is not affected by aliasing issues (see for example Eq. (39)).

It is also worth noting that these models were shown to recover the behavior of the BGK collision operator based on $f^{eq,2}_{i}$ if one neglects higher order velocity dependent terms. Knowing that a fourth order equilibrium state $f^{eq,4}_{i}$ increases the linear stability of the LBM Coreixas (2018); Wissocq et al. (2019), it would be interesting to determine if these orthogonal MRT and TRT models recover the behavior of the corresponding BGK-LBM. A few words regarding this possible link can be found in Sec. V.3.

In the rest of the paper, LBMs based on the orthogonal basis (38), and its 3D extensions d’Humières (2002); Suga et al. (2015); Fakhari et al. (2017), will be denoted as orthogonal MRT and TRT models.

To find polynomial bases compliant with the order of accuracy of the velocity discretization, it was proposed to build $D$-dimensional bases following tensor product rules of 1D bases Karlin and Asinari (2010). Since the D2Q9 lattice is a second order velocity discretization Shan et al. (2006), these rules lead to the construction of a basis only composed of polynomials of the form $\xi_{i,x}^{n}\xi_{i,y}^{m}$ with $n\leq 2$ and $m\leq 2$. The most natural basis compliant with the D2Q9 lattice is then

[TABLE]

where the subscript TP stands for the tensor product formalism. In fact, any basis composed of linear combinations of $\mathcal{B}_{\mathrm{TP}}$ elements would also be correct. Hence, its Hermite counterpart (HTP) could also be used Malaspinas (2015); Coreixas et al. (2017),

[TABLE]

If one wants to further decouple shear and bulk viscosities, one simply needs to change $\xi_{i,x}^{2}$ and $\xi_{i,y}^{2}$ ($H_{i,20}$ and $H_{i,02}$) to $\xi_{i,x}^{2}+\xi_{i,y}^{2}$ and $\xi_{i,x}^{2}-\xi_{i,y}^{2}$ ($H_{i,20}+H_{i,02}$ and $H_{i,20}-H_{i,02}$) as proposed in Refs. Karlin et al. (2011a, 2014); Lycett-Brown and Luo (2014).

This TP approach paved the way for a systematic derivation of new MRT-LBMs as introduced hereafter. It is finally worth noting that other kinds of TRT models were also proposed in different contexts Ansumali et al. (2007); Adhikari and Succi (2008), but they will not be studied in the present work since they do not introduce new concepts as compared to those already contained within the orthogonal TRT model.

III.4 Central moment space

Despite a nonnegligible increase in stability of both the MRT- and TRT-LBMs, issues were still encountered when simulating high Reynolds number flows. It was then suggested by several authors to improve the numerical stability of MRT-LBMs by applying the collision step in the comoving reference frame Geier et al. (2006, 2009); Geller et al. (2013); Lycett-Brown and Luo (2014); Dubois et al. (2015a); Ning et al. (2016); Fei and Luo (2017); De Rosis (2017a); Fei et al. (2018); De Rosis and Luo (2019); De Rosis et al. (2019).

These models rely on two modifications of the previous MRT-LBMs. First, the collision is applied in the central moment (CM) space based on either the TP Geier et al. (2006, 2009); Lycett-Brown and Luo (2014); Fei and Luo (2017); De Rosis (2017a); Fei et al. (2018); De Rosis and Luo (2019); De Rosis et al. (2019) or the orthogonal MRT formalism Dubois et al. (2015a, b), i.e,

[TABLE]

where the tilde stands for the velocity shift of each vector of the basis: $\widetilde{\xi}_{i,x}=\xi_{i,x}-u_{x}$ and $\widetilde{\xi}_{i,y}=\xi_{i,y}-u_{y}$. Second, the equilibrium state is also expanded following the TP formalism, leading to the inclusion of third and fourth order terms in the definition of $f_{i}^{eq}$. It is important to note that while these terms were added in an a posteriori manner in the first (cascaded) model Geier et al. (2006, 2009), the TP formalism offers an a priori and systematic way to include high order terms without exceeding the order of accuracy of the lattice. By extending the previous reasoning on equilibrium moments to the equilibrium population itself, one can derive new equilibrium states that includes high-order Hermite polynomials. This improves the Galilean invariance of the resulting LBM in a systematic and a priori way, using either the CM formalism Huang et al. (2018); De Rosis and Luo (2019); De Rosis et al. (2019) or any other one. This point will further be investigated in Sec. VI.4.

The CM-LBM can be summarized by

[TABLE]

where $\bm{N}$ and $\bm{N^{-1}}$ are the velocity dependent matrices used to move from the raw to the central moment space, and vice versa Dubois et al. (2015a); Fei and Luo (2017). In the particular case of a flow at rest ($u_{x}=u_{y}=0$), these matrices reduce to the identity matrix, and the orthogonal MRT (or TP) collision model is recovered. $f^{eq,\mathrm{ext}}$ is the extended equilibrium state that includes up to fourth order (six order) terms for the D2Q9 (D3Q27) velocity discretization Huang et al. (2018); De Rosis and Luo (2019); De Rosis et al. (2019). Both of these modifications lead to a drastic stability gain, and especially when the TP formalism is adopted for the moment space Dubois et al. (2015b). This is why only this formalism will be considered in the present work.

Eventually, it is important to note that despite a further stability increase, these models were mostly validated imposing acoustically related moments (those including the trace of second order moments) to their equilibrium value in the single phase and isothermal case Geier et al. (2006, 2009); Geller et al. (2013); Ning et al. (2016). This translates in a severe overdissipation of acoustic waves. As a consequence, it is of paramount importance to properly decouple the stability increase induced by the CM-LBM itself, from the one induced by the bulk viscosity. A few results about this problem can be found in Ref. Dubois et al. (2015b), but more in depth investigations are still required to properly decouple these two stabilization mechanisms. For example, it would be interesting to quantify the impact of the moment space, and of relaxation frequencies, on both academic and realistic configurations that would include acoustically related phenomena.

III.5 Regularization steps and Hermite moments

In parallel of the derivation of CM-LBMs, another type of collision models have been proposed to increase the numerical stability of the BGK-LBM without accounting for collisions in the moment space. Instead, these regularized collision models aim at filtering out nonhydrodynamic contributions that appear during the streaming step Latt and Chopard (2005, 2006); Zhang et al. (2006); Niu et al. (2007); Coreixas et al. (2017); Coreixas (2018). This is done projecting the nonequilibrium VDF on Hermite polynomials up to a given order $N$. Hence, this regularized collision model can be summarized as

[TABLE]

where ‘:’ is the tensor index contraction, $w_{i}$ are the quadrature weights, $c_{s}$ is the lattice constant, $\bm{a}_{n}$ are Hermite coefficients, and the subscript PRN stands for the projection based regularization at order $N$, with $N=2$ in the original model Latt and Chopard (2005, 2006), whereas $N\geq 3$ for high order LBMs Zhang et al. (2006); Niu et al. (2007); Coreixas et al. (2017); Coreixas (2018). In the original PR framework, Hermite nonequilibrium coefficients read as

[TABLE]

The above description corresponds to the way the regularization step was first introduced. In fact, it can be reinterpreted in terms of a MRT collision model based on $\mathcal{B}_{\mathrm{HTP}}$ (40) and the corresponding extended equilibrium state Coreixas (2018)

[TABLE]

and where, $\bm{S_{\mathrm{PR2}}}=\mathrm{diag}(0,0,0,\omega_{\nu},\omega_{\nu},\omega_{\nu},\omega_{3},\omega_{3},\omega_{4})$, with $\omega_{3}=\omega_{4}=1$. This is the sole set of relaxation frequencies that will be considered for this collision model in the present work. It is interesting to note that a more general MRT-LBM based on the Hermite moment space was also introduced in Refs. Shan and Chen (2007); Chen et al. (2014).

Despite its great success in various fields Niu et al. (2007); Wang et al. (2015); Ba et al. (2018); Montessori et al. (2018), it was recently proven that this model does not filter out all nonhydrodynamic contributions Malaspinas (2015); Coreixas et al. (2017). The reason lies in the fact that part of them are still present in $\bm{a}^{neq}_{n}$ via $f_{i}-f^{eq,4}_{i}$. Indeed, the latter only reduces to the first order (Navier-Stokes) nonequilibrium VDF $f^{(1)}_{i}$ in the continuum limit, which might no be valid anymore when the mesh becomes very coarse, i.e, nonnegligible Knudsen number based on the grid cell size. To counter this defect, its was proposed to impose the correct nonequilibrium part through a recursive computation of $\bm{a}^{(1)}_{n}$ using formulas derived from the Chapman-Enskog expansion Malaspinas (2015); Coreixas et al. (2017); Brogi et al. (2017); Coreixas (2018). This recursive regularized (RR) collision model is defined as

[TABLE]

with $\bm{S}=\omega_{\nu}\bm{I}$ in the original model, $\bm{I}$ being the identity matrix. The recursive formulas for the computation of $\bm{a}^{(1)}_{n}$ were derived in Refs. Malaspinas (2015); Coreixas et al. (2017). They are recalled in the case of the D2Q9 lattice in Eq. (82), whereas formulas for the D3Q27 are summarized in Eq. (89). The single relaxation time (SRT) PR collision model was shown to recover the behavior of the BGK operator when it is applied to the complete Hermite polynomial basis of the D2Q9 lattice, i.e, $\omega_{4}=\omega_{3}=\omega_{\nu}$ Coreixas et al. (2017). When it comes to the RR procedure, no direct link with the BGK operator is known at the time of writing.

Before moving to LBMs based on cumulants, it is worth noting that extensions of the PR model to other moment spaces have been proposed in the literature. A first extension was proposed in the TRT context Silva and Semiao (2017). Imposing $\Lambda_{q}^{-}\Lambda_{\varepsilon}^{+}=1/4$, it was noticed that nonhydrodynamic contributions could be filtered out in the same spirit as for the original PR model. Unfortunately, this model suffers the same problems as the orthogonal TRT, since it also relies on the same polynomial basis $\mathcal{B}_{\mathrm{TRT}}$. Another extension was proposed in the framework of collisions occurring in the comoving reference frame. Originally, the purpose of such an extension was to improve the Galilean invariance of LBMs relying on the central Hermite polynomial expansion of VDFs Chen et al. (2015). This formalism was further used to increase the stability of the PR model for the simulation of high Reynolds number flows Mattila et al. (2017). Discarding several diffusive terms in an adhoc manner, the authors recovered recursive formulas obtained from the Chapman-Enskog expansion in the context of the isothermal RR collision model. This suggests that there might be an underlying link between the CHM and the RR frameworks. More information will be given in Sec. V.4 regarding this last point.

III.6 Cumulant space

The most recent improvement of the MRT-LBM is based on a collision step occurring in the cumulant space Geier et al. (2015, 2017a, 2017b). The main advantage of cumulants, as compared to standard raw and central moments, is their ability to quantify the deviation of a probability density function with respect to a Gaussian distribution McCullagh (2018), in our case, the normalized equilibrium VDF $f^{eq}/\rho$. This was confirmed in Sec. II.7 where the only nonzero equilibrium cumulants were the first, and the second ones:

[TABLE]

By definition, cumulants belong to a type of statistical quantities quite different from raw and central moments. Nevertheless, they directly flow from (the logarithm of) the moment generating function as explained in Sec. II. Besides, they also share some similarities with central moments in the sense that their second and third order terms are strictly equivalent, as already precised in Refs. Kendall (1940); Cook (1951), and further illustrated in Apps. C-E. This means that without doing any calculations, it is known for sure that CM-LBM and cumulant-LBM (K-LBM) are identical for the D1Q3 lattice.

Due to the nonlinear relationships between cumulants and the other moment like quantities for $n\geq 4$, it is not possible to express the collision model in a linear matrix form. Nonetheless, one obtains a very simple algorithm using formulas derived to link central moments with respect to the cumulant one, and vice versa (Apps. C-E). According to Refs. Geier et al. (2015, 2017a, 2017b), the K-LBM cannot recover the behavior of the BGK-LBM. This point will be further studied in Sec. V.5.

III.7 Partial conclusions

Almost all collision models can be rewritten in a linear matrix form,

[TABLE]

with the exception of the K-LBM. Hence, it is already known for sure that this LBM belongs to a completely different family of LBMs. For collision models that satisfy Eq. (42), another distinction can be made between those relying on the approximation of the continuum limit $f^{neq}_{i}=f_{i}-f^{eq}_{i}$, and the RR collision model which imposes $f^{neq}_{i}=f^{(1)}_{i}$ instead. Regarding models operating in the comoving reference frame, their linear matrix form depends on the local velocity. Hence it seems complicated, at first sight, to find relationships with other collision models expressed in the reference frame at rest. Consequently, four different groups seem to emerge from the above review: (1) RM, HM and PR, (2) CM and CHM, (3) RR and (4) cumulant based collision models. Nevertheless, the possible relationship between RR- and CHM-LBMs remains to be addressed.

The rest of the paper will mainly be devoted to the validation of the above assumptions. To do so, it will be confirmed if the form of the equilibrium state is independent of the moment space or not. In the meantime, possible sets of relaxation frequencies allowing to link different types of collision model will be sought.

For the sake of completeness, the derivation of all linear transformations matrices $\bm{M}$ and $\bm{M^{-1}}$ are detailed for the D2Q9 lattice in App. F. General explanations regarding their construction are also provided, which makes possible their derivation for both the D3Q19 and D3Q27 lattices in a straightforward manner.

IV Theoretical comparison using the D1Q3 lattice

The purpose of this section is to compare, in the one-dimensional case, the most common collision models encountered in the LBM framework through their equilibrium, pre and post collision populations. To do so, the derivation of raw, Hermite, central, central Hermite LBMs is first recalled before moving to the comparative study.

IV.1 Discrete moments and populations

Let us consider the one-dimensional discretization of the velocity space using the following three discrete velocities $(\xi_{0},\xi_{\pm 1})=(0,\pm 1)$, namely, the D1Q3 lattice Krüger et al. (2017). This model has three degrees of freedom meaning it allows, at best, to preserve the first three moments of the VDF during the velocity discretization, i.e, $\forall n\in\{0,1,2\}$,

[TABLE]

where $P_{n}$ is a normalized moment of order $n$, yet to be defined. $\mathcal{P}_{n}$, and $\mathcal{P}_{i,n}$ (its discrete counterpart), also depend on the framework on which the collision model relies. Four different moment spaces will be considered hereafter. Hence, $\mathcal{P}_{i,n}$ can be defined using the

Raw moment (RM) space

[TABLE] 2. 2.

Hermite moment (HM) space

[TABLE] 3. 3.

Central moment (CM) space

[TABLE] 4. 4.

Central Hermite moment (CHM) space

[TABLE]

where $H_{i,n}$ is the discrete version of the Hermite polynomial of degree $n$, as introduced in Eq. (5). The tilde symbol stands for the central moment version of the polynomial basis considered for the collision process.

Once the moment space has been chosen, Eq. (43) is enforced adopting the correct form of $f_{+1}$, $f_{0}$ and $f_{-1}$. Assuming the collision takes place in the RM space, the following system needs to be solved

[TABLE]

which gives

[TABLE]

with $\sigma=\pm 1$. Corresponding equilibrium VDFs are then obtained computing equilibrium moments $M_{n}^{eq}$ that satisfy the moment conservation rule (43). By using the isothermal moments of the continuous Maxwell equilibrium VDF (35), one obtains

[TABLE]

which leads to

[TABLE]

$c_{s}$ is the isothermal (or Newtonian) speed of sound, and it corresponds to the lattice constant when the lattice Boltzmann unit system is adopted Viggen (2014).

Instead of performing collisions in the moment space using the linear matrix form of the collision operator, it is preferred here to work directly with VDFs. Post collision VDFs then read as

[TABLE]

with $M_{n}^{*}$ the post collision moment defined as

[TABLE]

where $\omega_{n}=1/(\tau_{n}+1/2)$ is the relaxation frequency associated to $M_{n}$, and $\tau_{n}$ the corresponding relaxation time. In the case of the D1Q3 lattice, only mass and momentum are collision invariants, thus one can freely choose the value of their corresponding relaxation frequency. Most of the time, one imposes $\omega_{0}=\omega_{1}=0$ in the absence of external forces Ginzburg et al. (2008a).

Injecting equilibrium moments (49) into the definition of post collision moments (51), post collision VDFs (50) now read as

[TABLE]

where $\omega_{\nu}=\omega_{2}$ is the relaxation time related to the kinematic viscosity $\nu$ through $1/\omega_{\nu}=\nu/c_{s}^{2}+1/2$.

IV.2 Impact of the moment space (D1Q3)

By solving Eq. (43) in different moment frameworks, various forms of post collision VDFs are obtained. In the HM framework,

[TABLE]

since $A_{0}^{*}=A_{0}^{eq}=1$, $A_{1}^{*}=A_{1}^{eq}=u_{x}$, $A_{2}^{eq}=u^{2}_{x}$. Using central moments, these VDFs are then

[TABLE]

where $\widetilde{M}_{0}^{*}=\widetilde{M}_{0}^{eq}=1$, $\widetilde{M}_{1}^{*}=\widetilde{M}_{1}^{eq}=0$, $\widetilde{M}_{2}^{eq}=c_{s}^{2}$. Eventually, the CHM framework leads to

[TABLE]

using $\widetilde{A}_{0}^{*}=\widetilde{A}_{0}^{eq}=1$, $\widetilde{A}_{1}^{*}=\widetilde{A}_{1}^{eq}=0$, $\widetilde{A}_{2}^{eq}=0$.

Regarding equilibrium moments, it is interesting to note that one can simply derive them applying two simple rules to Eq. (49). First, one can switch from equilibrium raw or central moments to their Hermite counterparts by neglecting $c_{s}$-dependent terms. The reason behind this is that Hermite polynomials reduce to monomials when terms proportional to $c_{s}$ are discarded. Second, the change from equilibrium raw (or Hermite) moments to their central (Hermite) versions is done discarding velocity dependent terms. This is explained by the fact that in the comoving reference frame, equilibrium moments are velocity independent. Consequently, all isothermal equilibrium moments, but the zeroth, are null in the CHM framework.

At first sight, all moment spaces seem to lead to different expressions for VDFs. Nonetheless, links can be drawn between all the above populations. Starting from post collision VDFs obtained within the RM framework (50),

[TABLE]

and

[TABLE]

Thus, the definition of post collision VDFs originating from the HM framework (54) are recovered from those of the RM approach. Moving now to the CM framework (55),

[TABLE]

and

[TABLE]

which proves that post collision VDFs originating from the CHM framework (56) are also recovered.

The question remains regarding the validity of this result when the reference frame is changed. Starting with,

[TABLE]

and using Eq. (11), one obtains $M_{1}^{*}=\widetilde{M}_{1}^{*}+u_{x}$ and $M_{2}^{*}=\widetilde{M}_{2}^{*}+u_{x}^{2}$. This leads to

[TABLE]

and it results in the equivalence between raw and central moment frameworks for the D1Q3 lattice.

In summary, raw, Hermite, central and central Hermite frameworks recovers the very same behavior (BGK) when the D1Q3 lattice is employed. Imposing $\omega_{\nu}=1$, one can further confirm that all models share the very same equilibrium state.

It is important to note that the K-LBM reduces to the CM-LBM as far as the D1Q3 lattice is employed since only moments up to the second order are considered for this velocity discretization. To quantify the possible benefit of the cumulant collision model in 1D, one must consider LBMs derived from high order velocity discretizations, such as the (zero-one-three) D1Q5 or the D1Q7 lattices Chikatamarla and Karlin (2009); Karlin and Asinari (2010); Frapolli (2017).

The last point that need to be checked is the way regularized collision operators are related to the above collision models. To answer this question, it is necessary to link populations derived through the Gauss-Hermite quadrature Grad (1949); Shan and He (1998); Shan et al. (2006) with those obtained via the conservation of Hermite moments (54).

IV.3 Gauss-Hermite quadrature and regularized collision models

Another way to derive populations associated to a given velocity discretization is based on the Gauss-Hermite (GH) quadrature Grad (1949); Shan and He (1998); Shan et al. (2006). This mathematical tool allows the exact preservation of certain properties (moments up to a certain order $N$) during the discretization of the velocity space. This method defines, in a unique and systematic way, the discrete populations $f_{i}$ as

[TABLE]

where the quadrature weights $w_{i}$ and the lattice constant $c_{s}$ can be obtained, for example, ensuring the preservation of Hermite polynomial orthogonality properties up to the order $N$ Philippi et al. (2006). For the D1Q3 lattice, the above expression reduces to,

[TABLE]

with $w_{0}=2/3$, $w_{\sigma}=1/6$, $c_{s}=1/\sqrt{3}$ and $a^{eq}_{n}=\rho u_{x}^{n}$ in the isothermal case Malaspinas (2015). Thus,

[TABLE]

using $a^{*}_{2}=\sum_{i}f^{*}_{i}H_{i,2}=\rho A^{*}_{2}$. Furthermore,

[TABLE]

Since both PR and RR collision models rely on this Gauss-Hermite quadrature, and that RR and PR approaches are equivalent for moments of order $n\leq 2$ Coreixas et al. (2017), it can be concluded that both regularization steps recover the behavior of all aforementioned collision models.

All in all, every collision model considered in the present work reduces to the BGK collision model, as far as the D1Q3 lattice is employed. The question is now to determine if the present conclusion remains valid in both 2D and 3D cases. To answer it, the $D$-dimensional extension based on the tensor product formulation will be adopted for both monomial Karlin and Asinari (2010), and Hermite polynomial bases Coreixas (2018). This step will allow us to properly understand which moments should be included in the derivation of populations, in the case of tensor product based LBMs such as the D2Q9 or the D3Q27 lattices (Secs. V and  VI respectively).

V Theoretical comparisons using the D2Q9 lattice

V.1 $D$-dimensional extension

The $D$-dimensional extension is done following the tensor product rules proposed in Refs. Karlin and Asinari (2010); Karlin et al. (2011b). In the most general (3D) case, they read as

Lattice

[TABLE] 2. 2.

Population

[TABLE]

with $\phi_{i}=f_{i}/\rho$, $\phi_{j}=f_{j}/\rho$, $\phi_{k}=f_{k}/\rho$, 3. 3.

Weight

[TABLE]

If not otherwise stated, the $3$-tuple $(i,j,k)$ (or $2$-tuple $(i,j)$) will be used instead of the conventional single index $i$ to describe a discrete velocity and its associated population.

Before moving to the derivation of post collision VDFs in the 2D case, it is important to note that by following the above rules, it can be shown (App. H) that all equilibrium states derived following these rules recover the exact same form whatever the polynomial basis used, and in any $D$-dimension. This important result flows from the fact that they already share the same equilibrium state in the 1D case, as demonstrated in Sec. IV.

Thanks to tensor product rules (67) and (68), populations evolving in the RM space now read as

[TABLE]

with $(\sigma,\lambda)\in\{-1,+1\}^{2}$. This leads to

[TABLE]

where the isotropy of the VDF itself is enforced up to the second order (in each direction) imposing Karlin and Asinari (2010)

[TABLE]

with $(p,q)\in\{0,1,2\}^{2}$.

Hence, post collision populations for the RM framework read as

[TABLE]

where post collision raw moments are defined as

[TABLE]

with $\nu_{b}$ being the bulk viscosity, and $\nu$ the kinematic viscosity. $\omega_{3}$ and $\omega_{4}$ are the relaxation frequencies associated to third and fourth order moments respectively. Furthermore, corresponding equilibrium moments are

[TABLE]

From this, the collision matrix cannot be diagonal if one requires $\nu_{b}\neq\nu$ Krüger et al. (2017). A simple remedy consists in replacing $(\xi_{i,x}^{2},\xi_{i,y}^{2})$ in $\mathcal{B}_{\mathrm{TP}}$ by $(\xi_{i,x}^{2}+\xi_{i,y}^{2},\xi_{i,x}^{2}-\xi_{i,y}^{2})$ as proposed in Ref. Bösch et al. (2015).

From a more general viewpoint, attention must be paid to the polynomial basis that is chosen since it may introduce extra computations as compared to the original one, and this overhead can become nonnegligible in the 3D case. As an example, this was investigated in Ref. Fei et al. (2018) by comparing the CPU times required by several CM-LBMs to simulate the very same configuration.

From now on, let us determine if previous conclusions regarding the impact of the polynomial basis in 1D are still valid in the 2D case.

V.2 Impact of the moment space (D2Q9)

Regarding Hermite, central and central Hermite moment spaces, their post collision populations are obtained following the very same tensor product rules as before (Eqs. (67) and (68)). In the HM framework, one obtains

[TABLE]

with

[TABLE]

While formulas expressed in RM and HM frameworks are rather simple, their counterparts in the comoving reference frame rapidly become lengthy. It is even worse in the 3D case, as detailed in App. G. Populations based on the Gauss-Hermite quadrature do recover the very same expression as those expressed within the HM framework. Hence they will only be considered in Sec. V.4 to further study possible links between RR- and CHM-LBMs.

For the D2Q9 lattice, HM equilibrium moments necessary for the collision step are defined as,

[TABLE]

while CM formulas read as

[TABLE]

Finally, expressions for the CHM framework are

[TABLE]

Now that all required formulas have been provided, let us start with the comparison between post collision VDFs obtained with both the RM and the HM frameworks. Considering the population at rest in the RM framework,

[TABLE]

so that

[TABLE]

It flows from these comparisons that it is not possible anymore to freely switch between raw and Hermite moments. This is explained by the fact that $M^{*}_{22}$ is not only linked to $A^{*}_{22}$ but also to $A^{*}_{20}+A^{*}_{02}$, while these two kinds of moments are relaxed using two different frequencies, namely, $\omega_{4}$ and $\omega_{\nu_{b}}$. The equivalency between these two approaches is then lost when $\omega_{4}\neq\omega_{\nu_{b}}$. One can further suppose that discrepancies between these approaches will tend to zero when phenomena linked to compressibility effects are negligible, i.e, $A^{neq}_{20}+A^{neq}_{02}\approx 0$. These deviations cannot be observed through the comparison of pre collision or equilibrium VDFs, hence the use of post collision VDFs for the comparison of different frameworks.

Since relationships between raw and Hermite moments are the same as those between central and central Hermite moments, it is known for sure that CM and CHM frameworks only merge for the same reasons as before.

Let us continue with the comparison between RM and CM frameworks. A similar computation shows that

[TABLE]

with

[TABLE]

As expected, even more deviations are present for other populations since they contain third order RM terms ($M_{21}$ and $M_{12}$) that will impact the relaxation of second order CM terms ($\widetilde{M}_{20}$ and $\widetilde{M}_{02}$). It is then clear that RM and CM frameworks only share the same behavior when $\omega_{4}=\omega_{3}=\omega_{\nu_{b}}=\omega_{\nu}$, and a fortiori, when the behavior of the BGK collision is recovered by both models. On the contrary, if a MRT approach is adopted then these deviations can only be neglected in the case of a flow at rest ($u_{x}=u_{y}=0$), for which both frameworks merge as already mentioned in Sec. III.4. Eventually, the case where $\widetilde{M}_{20}^{neq}+\widetilde{M}_{02}^{neq}$, $\widetilde{M}_{20}^{neq}-\widetilde{M}_{02}^{neq}$ and $\widetilde{M}_{11}^{neq}$ are simultaneously negligible is very restrictive, since it would imply that both compressibility and shear effects are also negligible. To the best of the authors’ knowledge, this is very unlikely to happen. As a consequence, this possibility will not be considered in the rest of the paper.

Once again, since relationships between Hermite and central Hermite moments are the same as those between raw and central moments, one can affirm that above conclusions can be extended to HM and CHM frameworks:

[TABLE]

In summary, the choice of the moment space used for the D2Q9-LBM does have an impact on the resulting post collision VDFs, when either a MRT approach is considered, or the flow is not at rest. More precisely, the deviation between collision models expressed in either the reference frame at rest or the comoving reference frame are not constant, and they depend on the local flow velocity for MRT approaches. These discrepancies between models naturally emerge from the 2D representation, while it was not present with the D1Q3 lattice.

V.3 Orthogonal MRT and TRT models

These two collision models are based on the expansion of populations over a polynomial basis built through the Gram-Schmidt orthogonalization procedure (38). The use of this basis leads to the following formulas for VDFs

[TABLE]

where the superscript LL stands for Lallemand’s and Luo’s MRT model introduced in Ref. Lallemand and Luo (2000). Interestingly, it is possible to show that the equilibrium form of Eq. (77) equals $f_{i}^{eq,4}$ if one does not neglect high order velocity terms in the equilibrium moments

[TABLE]

By fixing all relaxation frequencies to $\omega_{\nu}$, the orthogonal LL-MRT recovers the behavior of the BGK-LBM based on $f_{i}^{eq,4}$. By extension, the orthogonal TRT-LBM also recovers its behavior fixing $\omega^{+}=\omega^{-}=\omega_{\nu}$. Thus, it is quite surprising that both models usually neglect $\mathcal{O}(u^{3})$ and $\mathcal{O}(u^{4})$ terms while the use of $f_{i}^{eq,4}$ might improve the linear stability of the resulting LBM, as already proven for the BGK operator and for both regularization steps in Ref. Coreixas (2018).

The question is now to determine if these collision models are further related to other models based on different moment spaces. Rewriting LL moments within the HM framework,

[TABLE]

Thus the only issue that will be encountered with post collision VDFs will originates once again from the relaxation of the fourth order term ${\varepsilon}^{*}$. Indeed,

[TABLE]

with $\Delta f^{*,\mathrm{HM}}_{\varepsilon}=3(\omega_{\varepsilon}-\omega_{\nu_{b}})\rho(A_{20}^{neq}+A_{02}^{neq})$ and $\omega_{\varepsilon}=\omega_{4}$ so that,

[TABLE]

As for RM and HM frameworks, the LL-MRT-LBM reduces to the HM-MRT-LBM if the collision step of the fourth order and the bulk moments follow the very same relaxation rate, i.e, $\omega_{\varepsilon}=\omega_{\nu_{b}}$. Once again, the same result is obtained when compressibility effects are negligible ($A_{20}^{neq}+A_{02}^{neq}\approx 0$). One can derive the very same constraints for the LL-MRT-LBM to recover the behavior of the RM-LBM. In the same spirit, one can further come to the conclusion that LL and CM frameworks merge if a SRT collision model is employed, since it reduces both approaches to the BGK-LBM. In the case of a MRT collision model, the LL-LBM reduces to the CM-LBM in the zero velocity limit, as it was also shown in the context of the RM collision model (74). In the particular case of the orthogonal TRT approach, the relaxation of all odd (or even) moments are linked together. Hence, the condition $\omega_{\varepsilon}=\omega_{\nu_{b}}$ is always satisfied for this collision model. In addition, the behavior of the BGK collision model is also recovered when $\omega^{+}=\omega^{-}=\omega_{\nu}$.

V.4 Regularized collision models

The interesting thing about the Gauss-Hermite quadrature is that it remains valid whatever the number of dimensions or the lattice considered Grad (1949); Shan and He (1998); Shan et al. (2006). Nevertheless, it was not used (in its original form) to define which high order components should be included in the definition of populations. In the context of the D2Q9 lattice, the original basis used for the polynomial expansion was

[TABLE]

One can notice that $\mathcal{B}_{\mathrm{Q9}}$ is not a ‘basis’ from the mathematical viewpoint since it does not contain nine elements. In order to improve the description of populations, it was proposed by several authors to include three more elements Dellar (2002, 2014); Malaspinas (2015); Coreixas et al. (2017). While the first approaches added these elements through a Gram-Schmidt orthogonalization procedure, the last two relied on Hermite polynomials that were orthogonal with respect to those belonging to $\mathcal{B}_{\mathrm{Q9}}$. Nonetheless, they all end up with the same complete (and true) basis

[TABLE]

that corresponds to $\mathcal{B}_{\mathrm{HTP}}$ (40), as explained in Ref. Adhikari and Succi (2008). This basis is the starting point to compare all kinds of regularization steps against each other. Let us first consider projection based regularized (PR) models in both the reference frame at rest Latt and Chopard (2005, 2006); Chen et al. (2006)

[TABLE]

and in the comoving reference frame where central Hermite coefficients are defined as Chen et al. (2015); Mattila et al. (2017)

[TABLE]

so that Hermite coefficients are now computed using relationships between Hermite and central Hermite moments (App. D). This leads to

[TABLE]

where $\widetilde{a}_{00}=\rho$ and $\widetilde{a}_{10}=\widetilde{a}_{01}=0$ have been assumed. Restricting ourselves to the isothermal case, equilibrium coefficients are then computed as

[TABLE]

while nonequilibrium coefficients are

[TABLE]

In Ref. Mattila et al. (2017), the authors suggested to discard part of nonequilibrium (diffusive) terms $\widetilde{a}^{neq}$ based on the fact that these terms should be negligible for high Reynolds number flows. Doing so, they noticed that by imposing $\widetilde{a}^{neq}_{22}=\widetilde{a}^{neq}_{21}=\widetilde{a}^{neq}_{12}=0$, they were able to recover the recursive formulas of the RR approach Malaspinas (2015),

[TABLE]

where $a^{neq}_{20}$, $a^{neq}_{02}$ and $a^{neq}_{11}$ are computed using Eq. (41). While this result was obtained in an adhoc manner, it can be shown that it originates from particular values of relaxation times in the CHM framework. Indeed, post collision central Hermite coefficients are expressed as

[TABLE]

where $\overline{\tau}_{pq}=1/2+\tau_{pq}$ is the discrete relaxation time, with $(p,q)\in\{0,1,2\}^{2}$ for the D2Q9 lattice. Considering $\overline{\tau}_{21}=\overline{\tau}_{12}=\overline{\tau}_{22}=1$, corresponding coefficients are fixed to their equilibrium values. In the CHM framework, all equilibrium moments are zero with the exception of $\widetilde{a}^{eq}_{00}$ (see Sec. II.7). Hence, RR formulas (82) are indeed recovered when Eqs. (81) and (83) are employed to compute $a^{*}_{pq}$ with this particular set of relaxation times. The SRT-RR-LBM then corresponds to a PR approach in the comoving reference frame where high order ($p+q\geq 3$) nonequilibrium contributions are filtered out.

Furthermore, both the MRT-RR-LBM (with $\overline{\tau}_{\nu_{b}}^{\mathrm{RR}}=1$) and the cascaded based LBM share the same behavior when $\overline{\tau}_{\nu_{b}}^{\mathrm{CM}}=\overline{\tau}_{21}^{\mathrm{CM}}=\overline{\tau}_{12}^{\mathrm{CM}}=\overline{\tau}_{22}^{\mathrm{CM}}=1$. This result is pertinent if and only if the reader has no interest in acoustically related phenomena.

Eventually, since it is mandatory to use several relaxation times in the CHM framework to recover $f_{i}^{*,\mathrm{RR}}$, then the RR collision model cannot recover the behavior of the BGK collision operator. Hence the MRT-RR procedure does belong to another category of collision models.

V.5 Cumulant space

Here, the discrepancies between collision steps occurring in the CM and the cumulant spaces are investigated. As demonstrated in Apps. D and E, cumulants are not equivalent anymore to central moments when high order terms ($n\geq 4$) have to be taken into account in the expansion of populations. In the context of the D2Q9 lattice, it has been shown that one should include the fourth order term $\widetilde{M}_{22}$. If one follows the definition of post collision cumulants, as they are defined in Refs Geier et al. (2015, 2017a, 2017b), one obtains

[TABLE]

Post collision populations then read as

[TABLE]

with $U_{x}=1-u_{x}^{2}$, $\sigma_{x}=\sigma+u_{x}$, $\sigma_{2x}=\sigma+2u_{x}$, $U_{y}=1-u_{y}^{2}$, $\lambda_{y}=\lambda+u_{y}$ and $\lambda_{2y}=\lambda+2u_{y}$. From this, one can see that deviations from their CM counterparts come from the fact that

[TABLE]

while

[TABLE]

where for the sake of simplicity the relaxation frequency of $K_{11}$, $K_{20}$ and $K_{02}$ was taken as $\omega_{\nu}$. The latter assumption amounts to impose $\omega_{\nu_{b}}=\omega_{\nu}$ Krüger et al. (2017).

Eventually, post collision populations can be rewritten as

[TABLE]

with

[TABLE]

From this, one needs to impose $\omega_{4}=\omega_{\nu}=1$ to make the deviation $\Delta f^{*,\mathrm{K}}_{22}$ disappear. In that particular case, the kinematic viscosity cannot be chosen freely. Hence, it is not possible to recover the behavior of the BGK collision operator using cumulants, as already assumed in Refs. Geier et al. (2015, 2017a, 2017b). As a consequence, the K-LBM belongs to another category of LBMs, different from those presented in Secs. V.2-V.4.

Before moving to the investigation of 3D models, the interested reader may refer to the Supplemental material \bibnotemark[SupMatQ9] for more information regarding the way to implement D2Q9-LBMs in either their RM, HM, CM, CHM, K or RR formulation. Thanks to them, one can further derive instructions for orthogonal models (LL and TRT) by switching from HM to LL moments using Eq. (78).

VI Further investigations

This section is dedicated to several studies. First, collision models will be compared in the 3D case considering both the D3Q27 and the D3Q19 lattices. Second, the impact of collision models on the resulting macroscopic equations will be investigated.

VI.1 D3Q27-LBMs

The most straightforward way to build LBMs based on the D3Q27 lattice is to consider it as a tensor product of three D1Q3 lattices (one for each spatial direction). Doing so, it is possible to rely on rules defined in Eqs. (67)-(69). This approach then leads to populations that are expanded over tensor product like polynomial bases. For raw moments,

[TABLE]

using the shorthand notation $\xi_{pqr}=\xi_{i,x}^{p}\xi_{i,y}^{q}\xi_{i,z}^{r}$. For other moment spaces, one just needs to change from monomials ($\xi_{pqr}$) to either Hermite polynomials ($H_{i,pqr}$), central monomials ($\widetilde{\xi}_{pqr}$) or central Hermite polynomials ($\widetilde{H}_{i,pqr}$). One can notice that by using the tensor product rules a true basis is obtained, i.e, it contains the same number of elements as the number of velocities in the lattice. In addition, one can prove that once again population obtained through the Gauss-hermite quadrature do share the same expression as those obtained via the tensor product rules.

Knowing the basis corresponding to each framework, one can easily derive post collision populations for the D3Q27 lattice (App. G). The comparison of collision models follow the same steps as those presented in both the 1D and the 2D cases, with the exception that now calculations are more complex. Taking the example of the RM framework, post collision populations are

[TABLE]

with $(\sigma,\lambda,\delta)\in\{-1,+1\}^{3}$. Before moving to the comparison itself, one must define the number of relaxation frequencies required for the collision step. To do so, one simply needs to categorize all types of moments employed in Eq. (86). As an example, one can suppose that $M_{210}$ and all its cyclic permutations ($M_{201},M_{120},M_{102},M_{021}$ and $M_{012}$) should be relaxed to their equilibrium value via the very same relaxation frequency, the latter being different for the relaxation of $M_{111}$. Consequently, one possible set of relaxation frequencies would be

[TABLE]

which is a particular case of the collision matrix adopted in Ref. Fei et al. (2018). Here, $\omega_{\nu_{b}}=\omega_{\nu}$ is assumed, and this corresponds to $s_{+}=s_{2}$ and $s_{-}=0$ in their framework. This is done for the sake of simplicity, and more importantly, without significantly impacting the conclusions that are drawn below.

Hereafter, the comparative study will be simplified by starting from $f^{*,\mathrm{RM}}_{i}$ (86), and using relationships between each kind of statistical quantities (App. E) to derive deviations between each type of collision models. If one considers the HM framework, deviations with the RM approach come from fourth and higher order moments. Considering the following moments

[TABLE]

is then sufficient to determine all discrepancies between both families of moments. As previously discussed for the D2Q9 lattice, if one discards flows where compressibility and shear phenomena are negligible, then formulas on fourth order moments (88a and 88b) lead to the constraints $\omega_{5}=\omega_{1}$ and $\omega_{6}=\omega_{2}$ to enforce the equivalence between RM and HM frameworks. Expressions for fifth and six order moments (88c and 88d) further implies that $\omega_{7}=\omega_{3}$ and $\omega_{8}=\omega_{5}=\omega_{1}$, since it is quite difficult to derive further constraints based on the physical meaning of these moments, at least in a straightforward manner. Considering that $\omega_{1}$ and $\omega_{2}$ are related to the dissipation of acoustic and shear waves, only $\omega_{3}$ and $\omega_{4}$ are free parameters. As for the D2Q9 lattice, both RM and HM frameworks do recover the behavior of the BGK collision operator when only one relaxation frequency is used. The very same conclusions can be drawn regarding CM and CHM frameworks.

For RM and CM frameworks, the relaxation of each raw moment $M_{pqr}$ impacts the relaxation of all lower order moments $\widetilde{M}_{p^{\prime}q^{\prime}r^{\prime}}$ such as $p^{\prime}\leq p$, $q^{\prime}\leq q$ and $r^{\prime}\leq r$. Hence, the configuration for which both approaches merge corresponds to the use of a single relaxation frequency (BGK). One can further show that if a MRT approach is considered, then both frameworks recover the same behavior if and only if the flow is at rest. Once again, one comes to the same conclusions for HM and CHM frameworks.

Orthogonal MRT and TRT collision models are based on the Gram-Schmidt orthogonalization procedure Suga et al. (2015); Fakhari et al. (2017). The latter does not strongly impact the relationships between families of moments since it relies on linear transformations. Hence, it can be shown that above conclusions remain valid for orthogonal MRT collision models. Once again, orthogonal TRT-LBMs recover the behavior of the HM-LBMs since the relaxations of even moments are controlled by the same relaxation frequency. If one further imposes $\omega^{+}=\omega^{-}=\omega_{\nu}$ then the orthogonal TRT approach mimics the behavior of the BGK operator, as it was already the case for the D2Q9 lattice.

When it comes to the PR collision model, it leads to different results when different moment spaces are considered, since it relies on the equilibriation of higher than second order moments ($\omega_{3}=\omega_{4}=\omega_{5}=\omega_{6}=\omega_{7}=\omega_{8}=1$). The question is then to determine if it is still linked to the RR approach when central Hermite moments are employed for the relaxation process. Using relationships between central Hermite and Hermite moments (Sec. E), one can show that the PR collision model applied in the CHM framework recovers the behavior of the SRT-RR-LBM. Indeed, by following the same steps as in Sec. V.4, one obtains

[TABLE]

where other formulas are derived through cyclic permutations. The above relationships do match recursive formulas flowing from the Chapman-Enskog expansion at the Navier-Stokes-Fourier level Malaspinas (2015); Coreixas et al. (2017). Nevertheless, the equivalence is once again lost when several relaxation times are employed for the RR collision model.

The last comparison is related to cumulants and central moments. In the 2D case, it was shown that the deviation was only impacting the sole fourth order term, while now fourth, fifth and six order terms are impacted. In addition, the deviation was directly proportional to the square of the relaxation frequency, whereas it is now proportional to up to the cube of the relaxation frequency. Due to these nonlinear deviations, it is not possible to link CM- and K-LBMs.

In summary, no further links can be derived between collision models for the D3Q27 lattice as compared to the 2D case. This was to be expected from the way the D3Q27 lattice is built using tensor product rules (67 and 68). Nevertheless, these rules do not hold anymore in the case of the D3Q19 velocity discretization. Consequently, one may wonder if the aforementioned conclusions are still valid for this particular LBM.

VI.2 D3Q19-LBMs

To derive populations for the D3Q19 lattice, it is proposed to start from the D3Q27-LBM and then to discard particular terms that are not compliant with the velocity discretization using a ‘pruning method’ Karlin and Asinari (2010); Karlin et al. (2011b). The latter is of particular interest since it also allows us to determine which moments should be included in the expansion of populations. Hereafter, another way to choose the number of moments that are necessary for the D3Q19 lattice is presented.

It has been known for a long time that the projection of either the D3Q15 or the D3Q19 lattices, onto the 2D velocity space, leads to the D2Q9 lattice Krüger et al. (2017). Nonetheless, only the D3Q19 velocity discretization contains all the discrete velocities of the D2Q9 lattice in each of its planes $(x,y)$, $(x,z)$ and $(y,z)$. Hence all moments encountered in the definition of populations in the 2D case should also be accounted for in the present case. In other words, $\mathcal{B}^{\mathrm{Q19}}$ should, at least, contain all the 3D versions of monomials encountered in Eq. (39). By 3D version, it is meant that, for example, $\xi_{120}$ should be included since $\xi_{12}$ belongs to the 2D basis. Starting from Eq. (39) and considering all cyclic permutations, one ends up with the following polynomial basis

[TABLE]

expressed in the RM framework, and where the shorthand notation $\xi_{pqr}=\xi_{i,x}^{p}\xi_{i,y}^{q}\xi_{i,z}^{r}$ was used for the sake of clarity. This simple yet rigorous reasoning leads to a polynomial basis composed of nineteen elements. Thus one can simply move from D3Q27-LBMs to their D3Q19-LBMs discarding $\xi_{111}$, $\xi_{211}$, $\xi_{121}$, $\xi_{112}$, $\xi_{221}$, $\xi_{212}$, $\xi_{122}$, $\xi_{222}$ and then using relationships detailed in App. E. Doing so, it is possible to extend several collision models to the D3Q19-LBM even if they were originally developed within the framework of the D3Q27 velocity discretization, such as both RR- and K-LBMs (see App. G). In addition, this basis was recently used to derive a more efficient CM-LBM based on the D3Q19 lattice Fei et al. (2018). Nonetheless, it should be noted that, at the time of writing, the accuracy of this D3Q19 formulation as not yet been compared to its D3Q27 counterpart.

Using the above methodology, one can see that $\mathcal{B}^{\mathrm{Q19}}$ still contains fourth order terms. Consequently, all links that have been drawn for the D3Q27 lattice remain valid for the D3Q19 velocity discretization. One just needs to be careful and only compare collision models based on the same equilibrium state (see App. H for both choices obtained with the D3Q19 lattice).

General instructions for the coding of 3D models, as well as pseudo codes of the D3Q27 formulations, are provided in Supplemental materials \bibnotemark[SupMat] and \bibnotemark[SupMatQ27] respectively.

VI.3 Partial conclusions

Tab. 1 summarizes all derived links between collision models. More precisely, it details all collision models whose behavior can be entirely or partially recovered using another collision model with a particular set of collision frequencies. This is done considering only cases where the resulting physics is not impacted by the choice of the relaxation frequencies. Taking the example of the D2Q9 lattice, if one applies the collision step within the CM framework then one can entirely recover the results of the BGK collision model using only one relaxation frequency $\omega^{\mathrm{CM}}_{\nu}$. Besides, the CM-LBM can partially recover the behavior of the CHM-LBM if one imposes that $\omega^{\mathrm{CM}}_{\nu_{b}}=\omega^{\mathrm{CM}}_{4}=\omega^{\mathrm{CHM}}_{4}=\omega^{\mathrm{CHM}}_{\nu_{b}}$. In other words, this implies that the bulk viscosity cannot be chosen independently of the relaxation coefficient of fourth order moments. The number of free parameters is then reduced, hence the partial recovery of the CHM-LBM behavior by the CM-LBM. The most constraining case is the reduction to a single relaxation time approach, where the behaviors of RM, HM, orthogonal MRT and TRT models are also recovered. Eventually, it would be possible to get the same results as those obtained with the (SRT) RR approach by imposing $\omega^{\mathrm{CM}}_{\nu_{b}}=\omega^{\mathrm{RR}}_{\nu_{b}}=1$, but this would lead to an extreme overdissipation of acoustic waves Marié et al. (2009); Xu and Sagaut (2011); Gendre et al. (2017). Thus, this link is not considered in Tab. 1.

VI.4 Macroscopic behavior

To conclude this comparative study, let us have a look at macroscopic equations flowing from all the different collision models.

VI.4.1 Motivations

Previous investigations suggest that the choice of moment space has a major impact on the macroscopic behavior of LBMs, and more specifically on their Galilean invariance properties (e.g., Refs. Geier et al. (2006, 2009); Premnath and Banerjee (2011); Geier et al. (2015); Ning et al. (2016); Geier et al. (2017a); Geier and Pasquali (2018) among others). Knowing that these collision models naturally rely on an extended equilibrium state, one might wonder if the improved macroscopic behavior of these models come from either (1) high-order velocity terms of their equilibrium state, (2) the moment space, or (3) both of them. With this idea in mind, the hydrodynamic limit of LBMs will be studied in a general manner, meaning that the asymptotic study will be conducted before the (space-time) numerical discretization of the collision model. This is explained by the fact that one would have to properly distinguish errors resulting from every possible discretization technique (finite difference, finite volume, finite element, discontinuous Galerkin, etc) in order to objectively quantify the impact of each collision model on the resulting macroscopic behavior, and this is out of the scope of the present work.

In the following, a brief review on the origin of errors encountered in the context of LBMs is first proposed, whereas the investigation of the macroscopic behavior of LBMs is conducted in the second part of this section.

VI.4.2 Velocity and space-time discretization errors

The LBM relies on two types of discretization, namely, the velocity and the space-time discretizations. To correctly make the distinction between errors that emerge from both of them, let us recall the two main steps that are required for the design of any numerical scheme, and a fortiori for the derivation of LBMs.

The first step consists in selecting the mathematical model corresponding to the desired level of approximation to reality Hirsch (2007). In the following, and for the sake of simplicity, let us start from the force-free Boltzmann equation (BE),

[TABLE]

where $\mathcal{C}$ is the general form of the collision term. The physics governed by Eq. (90) goes far beyond the validity of standard macroscopic equations of interest, namely, the compressible Navier-Stokes-Fourier equations Huang (1987). As a consequence, this mathematical model is too detailed for the level of approximation to reality that is required in the present framework. To simplify it, a physical discretization is applied to the BE. It consists in drastically reducing the number of possible velocities $\xi$ for the propagation of populations $f$. This is done in such a way that the physics of interest is not lost during the overall process. One then ends up with the discrete velocity Boltzmann equation, also known as lattice Boltzmann equation (LBE),

[TABLE]

with $V$ being the number of discrete velocities required to recover the macroscopic behavior of interest Shan et al. (2006).

Once this set of partial differential equations has been chosen, the second step consists in choosing a type of space and time discretization in order to be able to numerically solve it Hirsch (2007). Among the wide panel of numerical discretizations available in the literature, the most commonly used approach is the ‘Collide and Stream’ algorithm (Eqs. (37) and (36)).

To determine the hydrodynamic limit of the LBM, one can either start from the LBE or from its numerical discretization – the LBM itself. Hence at least two types of asymptotic study are possible. The first methodology is the most commonly used in the lattice Boltzmann community, as it naturally flows from Statistical Physics, and it is based on the Chapman-Enskog expansion of the LBE Chapman and Cowling (1970); Huang (1987). The second approach is more generally used for the evaluation of numerical errors introduced during the space and time discretization of a given set of equations. It is based on the Taylor expansion of these discretized equations, and it allows the user to rigorously evaluate the order of accuracy of the numerical scheme of interest Hirsch (2007). In the lattice Boltzmann framework, the Taylor expansion is applied to the LBM itself by considering either a diffusive Junk et al. (2005) or an acoustic Dubois (2008) scaling. While the former leads to the derivation of error terms with respect to the incompressible Navier-Stokes equations, the latter derives them in the context of their weakly compressible formulation.

Since the Chapman-Enskog expansion is conducted before the space-time discretization of the LBE, results obtained from this asymptotic study remain valid whatever the numerical discretization considered (finite difference, finite volume, finite element, discontinuous Galerkin, etc). On the contrary, the Taylor expansion leads to the derivation of error terms that flow from both velocity and space-time discretizations. Hence, it is not possible to determine in a straightforward manner which discretization is related to error terms obtained with this expansion. For all of these reasons, the Chapman-Enskog expansion of the LBE seems to be the best approach to draw general conclusions about the hydrodynamic limit of LBMs. Results obtained hereafter will then flow from this asymptotic study. For the sake of completeness, it is also worth noting that the Chapman-Enskog could also be applied to the LBM itself, in order to further derive numerical errors introduced by the space-time discretization Otomo et al. (2017); Gendre et al. (2017).

VI.4.3 Macroscopic equations

In the present context, only second order velocity discretizations (D1Q3, D2Q9, D3Q19 and D3Q27) are considered. Let us start with their second order equilibrium state Qian et al. (1992); Shan et al. (2006),

[TABLE]

The corresponding (isothermal) macroscopic equations recovered through the Chapman-Enskog expansion Chapman and Cowling (1970) are

[TABLE]

where $\bm{\Pi}^{\prime}=\mu\left[\bm{\nabla}\bm{u}+(\bm{\nabla}\bm{u})^{T}\right]+\mathcal{O}(\mathrm{Ma}^{3})$ is the viscous stress tensor with the isothermal and weakly compressible limitation. The superscript $T$ stands for the transpose operator and $\mu=\tau p$ is the dynamic viscosity. For the D1Q3 lattice, the compressibility error reads as

[TABLE]

In the case of the D2Q9 lattice, each component of $\bm{\Pi}$ is affected as follows

[TABLE]

When it comes to both the D3Q19 and the D3Q27 lattice, even more error terms are obtained

[TABLE]

It is important to note that the general form of the above error terms is independent of the collision model framework. Only the truncation order of the equilibrium state $f_{i}^{eq}$ does have an impact on the compressibility error terms, since Krüger et al. (2017); Coreixas (2018)

[TABLE]

It is not possible to get rid of third order terms proportional to $u_{\alpha}^{3}$ ($\alpha=x,y,z$), due to the aliasing defect $\xi_{i,\alpha}^{3}=\xi_{i,\alpha}$ that is specific to second order velocity discretizations. Nonetheless, by taking into account third order moments compliant with these velocity discretizations, one can discard the influence of all non diagonal error terms as already pointed out in several studies Dellar (2002); Házi and Kávrán (2006); Dellar (2014); Malaspinas (2015); Coreixas et al. (2017); Huang et al. (2018); De Rosis and Luo (2019); De Rosis et al. (2019). This result is valid for both D2Q9 and D3Q27 lattices, as long as, third order velocity dependent terms are not discarded in the definition of their equilibrium state. For the D3Q19 lattice, error terms proportional to $\partial_{z}(u_{x}u_{y}u_{z})$ are still present since $\xi_{i,111}$ is not taken into account in $\mathcal{B}_{\mathrm{Q19}}$. The interested reader may refer to App. H where full forms of equilibrium states are compiled.

Regarding the impact of the moment space used for the collision process, it is important to understand that all moment spaces, even the cumulant one when fourth- and higher-order cumulants are discarded, can be related to each other through linear transformation matrices (see App. F). In this context, the only difference between the macroscopic behavior of LBMs comes from their original equilibrium state. As compared to standard collision models (BGK, MRT, PR, etc) which originally relied on a second-order equilibrium state, the CM-, CHM-, K- and RR-LBMs were derived using extended ones that include high-order velocity terms. It is then only natural that the latter models reduce the number of velocity dependent error terms present in the viscous stress tensor. Nonetheless, the present reasoning also suggests that by keeping high-order velocity terms, even for standard collision models, then one would also improve the macroscopic behavior of the resulting LBM. This was discussed by several authors Dellar (2002); Házi and Kávrán (2006); Dellar (2014); Malaspinas (2015); Coreixas et al. (2017); Huang et al. (2018); De Rosis and Luo (2019); De Rosis et al. (2019), and further confirmed through the linear stability analysis of several LBMs Coreixas (2018).

Consequently, only the equilibrium state, which depends on the lattice (and not on the collision model), does impact the resulting macroscopic equations recovered by the isothermal and weakly compressible LBMs. For the sake of fairness, however, it is also important to understand that by changing the collision model then one also modifies the numerical properties of the LBM. This is explained by the fact that the only error introduced by the numerical discretization comes from the collision term Dellar (2013); Krüger et al. (2017). Hence, one can show that the numerical behavior of LBMs is drastically impacted by the choice of both the moment space and the relaxation parameters d’Humières and Ginzburg (2009); Geier et al. (2015, 2017a); Geier and Pasquali (2018). Nevertheless, these errors should not be attributed to a physical problem, such as the Galilean invariance issue, but rather to a purely numerical defect, in order not to mislead the reader on this particularly complex topic.

VII CONCLUSION

The BGK-LBM has been really successful during the last decades due to both its efficiency and accuracy. Nevertheless, this single relaxation time (SRT) collision model is known to suffer from severe stability issues during the simulation of high Reynolds number flows. To circumvent this deficiency numerous collision models has been developed. They can be classified according to the number of relaxation times they rely on, and via the moment space used for the derivation of populations. Most of the time, all these collision models are presented in their own framework with only very few suggestions regarding their possible link with already known models. While authors meanly concentrate on the comparison of collision models through a list of numerical test cases, they seldomly explain observed discrepancies from a theoretical viewpoint.

In this context, the present work focused on the understanding of fundamental differences between the most common collision models. To do so, an extensive search of links between these collision operators was performed in the case of standard LBMs of increasing complexity (D1Q3, D2Q9, D3Q19 and D3Q27 lattices). After drawing relationships between all moment spaces (raw (RM), central (CM), Hermite (HM), central Hermite (CHM), cumulants (K)) through the use of their corresponding moment generating function, a thorough review of collision models was conducted. The latter drew a first picture of known links between collision models and further showed that all these collision models can be rewritten in a linear matrix form with the exception of the K-LBM.

In a general way, it was also demonstrated that all collision models recover the very same second, fourth and six order equilibrium states when the D1Q3, the D2Q9 and the D3Q27 lattices are respectively employed for the velocity discretization of the Boltzmann equation. Nonetheless, the use of the D3Q19 lattice led to two different equilibrium states because of its non compliance with tensor product rules. Interestingly, the use of these extended equilibrium states might improve the linear stability of the corresponding LBMs, as recently demonstrated for both BGK and regularized D2Q9-LBMs Coreixas (2018).

Regarding the mathematical comparison of collision models, while it was not possible to find discrepancies between them in the one-dimensional case, deviations started appearing in the two-dimensional case. Using the D2Q9 lattice and with the assumption of a SRT operator, all models recovered the behavior of the BGK-LBM but both RR- and K-LBMs. Using several relaxation times, a partial equivalency between raw and Hermite frameworks was obtained in both the reference frame at rest (RM and HM) and in the comoving reference frame (CM and CHM), when the collision of bulk and fourth order moments were sharing the same relaxation frequency ($\omega_{\nu_{b}}=\omega_{4}$). Still using a multirelaxation time (MRT) approach, it was confirmed that raw and central frameworks can only share the same behavior if the simulated flow is at rest, i.e, imposing $u_{x}=u_{y}=0$.

In addition, the projection based regularization (PR) step was recasted in a collision step occurring in the HM framework, where nonhydrodynamic contributions are actually filtered out imposing $\omega^{\mathrm{HM}}_{3}=\omega^{\mathrm{HM}}_{4}=1$. Most importantly, the SRT-RR-LBM was reinterpreted as an extension of the PR approach to the comoving reference frame (CHM space) where third and fourth order contributions were also discarded ($\omega^{\mathrm{CHM}}_{3}=\omega^{\mathrm{CHM}}_{4}=1$). Imposing $\omega_{\nu_{b}}^{\mathrm{CM}}=\omega_{3}^{\mathrm{CM}}=\omega_{4}^{\mathrm{CM}}=1$, the cascaded-LBM was further demonstrated to correspond to the MRT-RR-LBM where $\omega_{\nu_{b}}^{\mathrm{RR}}=1$. The latter condition is viable if and only if one is not interested in the simulation of acoustically related phenomena. This might explain why the cascaded collision model was usually employed with these collision frequencies in order to simulate incompressible flows. In conclusion, the most general formulation of the MRT-RR-LBM was shown to belong to a completely different kind of collision models.

Discrepancies between orthogonal (MRT and TRT) and nonorthogonal approaches (RM and HM frameworks) were also highlighted. For the orthogonal MRT, they were shown to originate from the orthogonalization procedure, which is based on the construction of an orthogonal basis through linear combinations of monomials, eventually leading to spurious entanglements between moments of different orders. From this perspective, the orthogonalization procedure seemed to induce more issues than it solved. The orthogonal version of the TRT-LBM was shown to be able to cancel out these spurious entanglements.

A deviation between the CM- and the K-LBMs appeared due to the inclusion of the fourth order central moment in the expansion of post collision populations. This discrepancy originates from the nonlinear relaxation of cumulants related to this central moment, eventually leading to a radically different behavior as compared to all aforementioned collision models.

The above results were further confirmed for the D3Q27 lattice where the number of deviations increased due to the presence of more high order moments. Guidelines concerning the extension of all collision models to the D3Q19 lattice were also provided. In particular, these explanations make possible the derivation of both RR- and K-LBMs in a straightforward manner (App. G). The comparative study then led to the very same conclusions as before due to the presence of fourth order contributions in the expansion of populations. Nevertheless, equilibrium states obtained through Gauss-Hermite quadrature and tensor product rules were shown to not be equivalent anymore, meaning that only collision models based on the same equilibrium state should be compared to each other.

It was finally shown that the reduction of error terms in the macroscopic equations recovered by the LBM – those originating from the limited accuracy of standard lattices and not from the numerical discretization – was only due to the form of the equilibrium state and not to moment space employed for the collision process.

For the sake of completeness, three Supplemental materials are also provided to help the reader with the coding of collision models considered in the present work. While the general methodology is discussed in \bibnotemark[SupMat], pseudo codes dedicated to the coding of the D2Q9 and the D3Q27 formulations are provided in \bibnotemark[SupMatQ9] and \bibnotemark[SupMatQ27] respectively.

Regarding future works, it is planned to further compare collision models through the evaluation of their linear stability domain. This study will include both results on eigenvalues Coreixas (2018) and eigenvectors Wissocq et al. (2019) of corresponding linearized LBMs. Not only will it allow the proper classification of collision models based on their stability domain, but it will also lead to a better understanding of the stabilization mechanism specific to each model. These linear stability analysis should also help to: (1) better quantify errors introduced by the numerical discretization of the collision term, and (2) find optimal values of collision frequencies to further improve the stability of LBMs in the low-viscosity regime. Collision models offering the best trade off between stability and accuracy will then be compared to dynamic models through both academic and realistic configurations.

In parallel to this, it is also planned to continue the work initiated in Ref. Coreixas (2018). This should lead to further extensions of all aforementioned collision models to high order LBMs, for which the simulation of compressible flows induces even more severe stability issues.

Acknowledgements.

Fruitful discussions with G. Wissocq, T. Astoul, F. Renard, D. Ricot, F. Gendre and O. Malaspinas are gratefully acknowledged. The authors would also like to thank the reviewers for their pertinent remarks and insightful questions. This work was supported by the Swiss National Science Fund SNF through the Sinergia project 170930, ‘A 3D Cell-Based Simulation Framework for Morphogenetic Problems’. All mathematical derivations were conducted using ‘Mathematica’ software.

Appendix A Lattices

Henceforth, the main characteristics of standard velocity discretizations of interest are compiled. Fig. 1 contains the representation of each velocity discretization, whereas Tab. 2 gives further information regarding their velocity sets $\bm{\xi}_{i}$, their associated Gauss-Hermite weights $w_{i}$ and their lattice constant $c_{s}$.

Appendix B Bell polynomials

These polynomials are used in combinatorial mathematics to study set partitions Bell (1934); Withers (1994). The Bell polynomial of degree $n$ is defined as

[TABLE]

where $B_{n,k}$ are partial Bell polynomials, and $x_{1},...,x_{n}$ are $n$ variables. Each $B_{n,k}$ corresponds to the partitioning of a set composed of $n$ elements into $k$ non-empty subsets.

As an example, let us consider the set $\mathcal{S}=\{a,b,c\}$, and let us find all its possible partitionings (described via $B_{3}$). They are three ways to partition $\mathcal{S}$:

One subset composed of three elements

[TABLE] 2. 2.

Two subsets composed of one and two elements

[TABLE] 3. 3.

Three subsets composed of one element

[TABLE]

Through $x_{1}$, $x_{2}$ and $x_{3}$, partial Bell polynomials $B_{3,1}$, $B_{3,2}$ and $B_{3,3}$ allow the mathematical description partitionings $\mathcal{S}_{1}$, $\mathcal{S}_{2}$ and $\mathcal{S}_{3}$ respectively. Assuming $x_{p}$ indicates the presence of a subset composed of $p$ elements, then

[TABLE]

since only one subset, composed of three elements, is included in $\mathcal{S}_{1}$. In addition,

[TABLE]

because $\mathcal{S}_{2}$ is composed of (three) subsets divided into two blocks of length two and one respectively. Eventually,

[TABLE]

since $\mathcal{S}_{3}$ contains three singletons. Hence, the corresponding Bell polynomial is

[TABLE]

Applying the very same reasoning for $1\leq n\leq 6$, Bell polynomials read as

[TABLE]

Defining $y_{n}={\sum_{k=1}^{n}B_{n,k}(x_{1},...,x_{n-k+1})},$ the inversion formula, allowing to express $x_{n}$ with respect to $y_{1},...,y_{n}$, is Bell (1934); Withers (1994)

[TABLE]

Finally, one obtains relationships between raw moments and cumulants (Eqs. 100 and  101) replacing $(x_{n},y_{n})$ by $(M_{n},K_{n})$ in the above formulas.

Appendix C Univariate formulas

In this appendix, formulas used to link all statistical quantities of interest between each other are compiled in the univariate case.

Relationships between raw and Hermite moments flows from formulas expressing Hermite polynomials with respect to monomials (5). Up to $n=6$, they read as

[TABLE]

and

[TABLE]

Continuing with relationships between raw and central moments, they are simply obtained using recursive formulas introduced in Eq. (11),

[TABLE]

and Eq. (10),

[TABLE]

It is important to note that replacing $M$ and $\widetilde{M}$ by their counterparts in the Hermite framework ($A$ and $\widetilde{A}$) the above formulas remain valid. This is explained by the fact that when the change from the reference frame at rest to the comoving one is done then a simple shift of $\exp(-M_{1}t)$ is applied to the raw (or Hermite) moment generating function.

Considering now relationships between raw moments and cumulants, formulas derived thanks to Bell polynomials (App. B) lead to

[TABLE]

and

[TABLE]

Regarding formulas between central moments and cumulants, one just need to move to the comoving reference frame. In fact, this amounts to neglecting terms proportional to $u_{x}$ in the above formulas Kendall (1940); Cook (1951). Hence

[TABLE]

and

[TABLE]

Appendix D Bivariate formulas

Henceforth, bivariate formulations of relationships between statistical quantities of interest are presented.

Starting with raw and Hermite moments, the most straightforward way to compute their bivariate relationships is to take advantage of the orthogonality properties of Hermite tensors. Hence, they are simply computed as Coreixas (2018),

[TABLE]

and

[TABLE]

One must be careful regarding the above method. Indeed, one must wait until the very end before replacing zeroth and first order moments by their values. As an example, $M_{11}\neq u_{x}u_{y}(=M_{10}M_{01})$ since $M_{11}^{neq}=M_{11}-M^{eq}_{11}$ is related to viscous phenomena and is usually non zero.

The simplest way to extend relationships between raw and central moments to the bivariate case is to rely on the corresponding binomial formula,

[TABLE]

Inversion formulas are then obtained discarding the minus sign,

[TABLE]

In the case of the D2Q9 lattice ($p\leq 2$ and $q\leq 2$), this leads to

[TABLE]

and

[TABLE]

Due to the fact that there is a linear transformation allowing us to move from monomials to Hermite polynomials, these relationships remain valid for Hermite and central Hermite moments. Hence, one simply needs to replace $M$ and $\widetilde{M}$ by $A$ and $\widetilde{A}$ in the above formulas to recover relationships between Hermite and central Hermite moments.

To link raw moments and cumulants in the bivariate case, differential operators introduced in Sec. II.6, and proposed by Kendall Kendall (1940); Cook (1951), are employed. Corresponding relationships are then

[TABLE]

and

[TABLE]

For central moments, one can again discard terms proportional to $u_{x}$ and $u_{y}$. Doing so, one obtains

[TABLE]

and

[TABLE]

where it is clear that deviations from central moments start appearing in fourth order cumulants.

From the implementation point of view, the conversion from central moments to cumulants is clearly simpler than the one using raw moments. Apart from that, there is no fundamental reason for the use of central moments instead of raw ones. Furthermore, formulas for post collision populations are more complex in the CM framework (73b) than in the RM one (71). In the 2D case, it is then not clear if the use of central moments should be preferred or not.

Appendix E Trivariate formulas

Hereafter, the methodology introduced in App. D is further extended to the trivariate case. Starting with raw and Hermite moments, orthogonality properties of Hermite polynomials are, once again, used to relate these two families of moments. For the D3Q27 lattice, corresponding formulas are

[TABLE]

For raw and central moments,

[TABLE]

Inversion formulas are once again obtained discarding the minus sign,

[TABLE]

Thus, relationships obtained for the D3Q27 lattice are

[TABLE]

and

[TABLE]

When it comes to raw moments and cumulants, Kendall’s differential operators are used to derive trivariate formulas from the univariate case Kendall (1940); Cook (1951). Up to $K_{222}$ and $M_{222}$, this leads to

[TABLE]

and

[TABLE]

Relationships between central moments and cumulants are finally obtained discarding velocity dependent terms in the above formulas,

[TABLE]

Appendix F Linear transformation matrices

This appendix compiles all linear transformations allowing us to write the collision step in a matrix form in the context of the D2Q9 velocity discretization, i.e,

[TABLE]

where {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{f}^{*}}, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{f}^{eq}} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{f}^{neq}} are vectors composed of all nine post collision, equilibrium and nonequilibrium populations respectively. Starting with the set of orthogonal polynomials (38), it is explained how to compute $\bm{M}$ and $\bm{M^{-1}}$ from a general point of view Lallemand and Luo (2000). This allows a straightforward extension to any kind of lattice, in either 2D or 3D. These linear transformation matrices (LTMs) are then given for RM, HM, CM and CHM frameworks. For the last two, the concept of shifting matrices is used to move from the reference frame at rest to the comoving one Fei and Luo (2017). Eventually, these matrices are provided in the context of the Gauss-Hermite quadrature on which regularized collision models are based. Hereafter, the single index $i$ is used to describe each discrete velocity $\bm{\xi}_{i}$ and its related population $f_{i}$.

F.1 Construction of LTMs

Let us start with the orthogonal basis proposed by Lallemand and Luo Lallemand and Luo (2000),

[TABLE]

where each of its component are defined as,

[TABLE]

Due to the aliasing defect of the D2Q9 lattice ($\xi_{i,x}^{3}=\xi_{i,x}$ and $\xi_{i,y}^{3}=\xi_{i,y}$), $M_{q_{x}}$, $M_{q_{y}}$ and $M_{q_{\varepsilon}}$ can further be simplified into

[TABLE]

On can see that the ordering of polynomials proposed here is different from the one introduced in the original paper. The present ordering is based on a progressive increase of the degree of each polynomial. Doing so, the comparison with other frameworks, such as cascaded or regularized collision models, is simplified.

In addition, the standard ordering of discrete velocities is considered. In other words, velocities are gathered within several groups depending on their norm, and they are then read ‘counterclockwise’. This leads to

[TABLE]

The ordering of both sets of discrete velocities and moments is the first source of mistakes when deriving LTMs. It is then of uttermost importance to check them before moving to the construction of $\bm{M}$ and $\bm{M^{-1}}$. In the present work, chosen orderings are described in Eqs. (118) and (38). While the LTM to change from populations to moments is derived replacing each $\bm{\xi}_{i}$ by its value in the chosen set of polynomials, its inverse is simply derived using standard linear algebra libraries. This eventually leads to

[TABLE]

and

[TABLE]

where the subscript LL stands for the orthogonal basis proposed by Lallemand and Luo. The corresponding collision matrix is diagonal, and it reads as $\bm{S_{\mathrm{LL}}}=\mathrm{diag}(0,0,0,\omega_{e},\omega_{\nu},\omega_{\nu},\omega_{q},\omega_{q},\omega_{\varepsilon})$. The collision frequency $\omega_{\nu}$ controls the relaxation of shear related phenomena, whereas $\omega_{e}$ is related to the attenuation of acoustic waves. $\omega_{q}$ and $\omega_{\varepsilon}$ further control the dissipation of third and fourth order moments respectively. The latter are usually considered as free parameters that can be tuned to increase the stability of the LBM without impacting the resulting physics.

F.2 Application to RM and HM frameworks

In the RM framework, monomials compose the polynomials basis

[TABLE]

which is now nonorthogonal. Following the same steps as before, new LTMs are obtained. They read as

[TABLE]

In addition, the collision matrix can have two forms depending on the prerequisites of the reader. If one is not interested in acoustically related phenomena, one can impose $\omega_{\nu_{b}}=\omega_{\nu}$ which leads to $\bm{S_{RM}}=\rm{diag}(0,0,0,\omega_{\nu},\omega_{\nu},\omega_{\nu},\omega_{3},\omega_{3},\omega_{4})$. Otherwise, the collision matrix is only block-diagonal Krüger et al. (2017); Fei et al. (2018)

[TABLE]

with

[TABLE]

For the HM framework, one can either use the above method to derive LTMs, or rely on relationships obtained for raw and Hermite moments (Eqs. (104) and (105)). The second method was originally used to build matrices allowing the shift from raw to central moments, and vice versa Fei and Luo (2017). Hereafter, it will be employed to derive matrices allowing to link raw and Hermite moments. To do so, one simply needs to rewrite relationships between both kinds of moment in the form of a row-column product. As an example,

[TABLE]

and

[TABLE]

Eventually,

[TABLE]

where $\bm{B_{\mathrm{RM}}^{\mathrm{HM}}}$ and $\bm{B_{\mathrm{HM}}^{\mathrm{RM}}}=\bm{(B_{\mathrm{RM}}^{\mathrm{HM}})^{-1}}$ allow us to move from raw to Hermite moments, and vice versa. The lower triangular form of these matrices translates the dependency of high order Hermite moments with respects to lower order raw moments and conversely.

Post collision populations are finally obtained through

[TABLE]

with $\bm{M_{\mathrm{HM}}}=\bm{B_{\mathrm{RM}}^{\mathrm{HM}}M_{\mathrm{RM}}}$ and $\bm{M_{\mathrm{HM}}^{-1}}=\bm{(B_{\mathrm{RM}}^{\mathrm{HM}}M_{\mathrm{RM}})^{-1}}=\bm{M_{\mathrm{RM}}^{-1}B_{\mathrm{HM}}^{\mathrm{RM}}}$. Once again, $\bm{S_{\mathrm{HM}}}$ is diagonal if one assumes $\omega_{\nu_{b}}=\omega_{\nu}$, and block-diagonal otherwise.

F.3 Change of reference frame

To derive LTMs corresponding to the CM framework, one starts from matrices relating populations to raw moments (121). Then, one further uses shifting matrices to move from raw to central moments. These shifting matrices are obtained thanks to relationships compiled in Eqs. (106) and (107). One must be careful and use binomial formulas where zeroth and first order moments are not replaced by their values, i.e, $M_{00}=\widetilde{M}_{00}=1$, $M_{10}=u_{x}$, $\widetilde{M}_{10}=0$, etc. As an example,

[TABLE]

and

[TABLE]

By rewriting zeroth through third order moment in the same spirit, one ends up with the following shifting matrices

[TABLE]

[TABLE]

LTMs for the CM framework then read as

[TABLE]

and

[TABLE]

To derive the LTMs for the CHM framework, one simply needs to start from $\bm{M_{\mathrm{CM}}}$ and $\bm{M_{\mathrm{CM}}^{-1}}$ and then to switch from CMs to CHMs. The latter transformation can be done using

[TABLE]

since relationships between RMs and HMs are the same as those between their counterparts in the comoving reference frame. Hence,

[TABLE]

and

[TABLE]

Appendix G D3Q27 and D3Q19 formulations of populations

This appendix is dedicated to the derivation of 3D populations for raw, Hermite, central and central Hermite moments. They are obtained enforcing the isotropy of $f_{i}$ up to the second order (in each direction),

[TABLE]

with $(p,q,r)\in\{0,1,2\}^{3}$. This isotropy condition is valid for all kinds of moments considered hereafter, with the exception of post collision cumulants due to their nonlinear relationship with central moments.

Assuming that $(\sigma,\lambda,\delta)\in\{\pm 1\}^{3}$, populations evolving in the raw moment space read as

[TABLE]

whereas in the case of Hermite moments they read as

[TABLE]

with $C=1-c_{s}^{2}$. It is interesting to note that $f_{i}^{\mathrm{HM}}$ reduce to $f_{i}^{\mathrm{RM}}$ imposing $c_{s}=0$ in the above formulas. This is due to the fact that $H_{i,pqr}=\xi_{i,x}^{p}\xi_{i,y}^{q}\xi_{i,z}^{r}$ when terms proportional to $c_{s}$ are discarded in the definition of Hermite polynomials.

[TABLE]

with $U_{x}=1-u_{x}^{2}$, $U_{y}=1-u_{y}^{2}$, $U_{z}=1-u_{z}^{2}$, $\sigma_{x}=\sigma+u_{x}$, $\sigma_{2x}=\sigma+2u_{x}$, $\lambda_{y}=\lambda+u_{y}$, $\lambda_{2y}=\lambda+2u_{y}$, $\delta_{z}=\delta+u_{z}$, $\delta_{2z}=\delta+2u_{z}$. Discarding velocity dependent terms in $f_{i}^{\mathrm{CM}}$ allows us to recover the definitions of $f_{i}^{\mathrm{RM}}$. This is another simple way that permits to easily check the validity of the above formulas.

[TABLE]

with

[TABLE]

Here, if one neglects $c_{s}$-dependent terms, the definition of $f_{i}^{\mathrm{CM}}$ is recovered. On the contrary, discarding velocity dependent terms leads to the definition of $f_{i}^{\mathrm{HM}}$. Eventually, neglecting both kinds of terms allows us to recover the definition of $f_{i}^{\mathrm{RM}}$. Regarding regularized LBMs, their populations are exactly the same as those expressed within the HM framework. For the K-LBM, one simply needs to replace central moments of order $n\geq 4$ in $f_{i}^{\mathrm{CM}}$ using corresponding formulas detailed in Eq. (117).

In the particular case of the D3Q19 velocity discretization, the easiest way to derive populations is to start from their D3Q27 counterparts and then to neglect moments that are not compliant with this lattice. One further needs to choose to either rely on the Gauss-Hermite formalism or not since this choice directly impacts the form of the equilibrium state (see App. H and more specifically Eqs. (159)-(165) and (166)-(172)). Once this is done, one simply needs to use relationships derived in App. E to switch from one framework to another one. Taking the example of $f_{i}^{\mathrm{CM}}$ (126), one needs to discard $M_{111}$, $M_{211}$, $M_{121}$, $M_{112}$, $M_{221}$, $M_{212}$, $M_{112}$, $M_{222}$ as explained in Sec. VI.2. This leads to

[TABLE]

One can notice that populations corresponding to discrete velocities $(\pm 1,\pm 1,\pm 1)$ are null. This confirms the validity of the proposed choice of raw moments. Replacing raw moments by their Hermite counterpart further results in

[TABLE]

with $C_{2}=1-2c_{s}^{2}$. For the CM framework, populations read as

[TABLE]

where $U_{xy}=1-u_{x}^{2}-u_{y}^{2}$, $U_{xz}=1-u_{x}^{2}-u_{z}^{2}$ and $U_{yz}=1-u_{y}^{2}-u_{z}^{2}$. To derive populations corresponding to the K-LBM, one simply needs to replace $\widetilde{M}_{220}$, $\widetilde{M}_{202}$ and $\widetilde{M}_{022}$ by their cumulant counterparts using Eq. (117). Eventually, populations defined in the CHM framework are

[TABLE]

with

[TABLE]

If one chooses to start from the Gauss-Hermite formalism, then populations read as

[TABLE]

with $n=p+q+r$. It is preferred here to rely on $\bm{A}_{pqr}$ instead of $\bm{a}_{pqr}=\rho\bm{A}_{pqr}$ to ease the change of framework through formulas compiled in App. E. From this, it is clear that the D3Q19 formulation of both regularization steps is based on the following populations

[TABLE]

where both the lattice constant and the weights have been replaced by their values (Tab. 2). Depending on the way $A_{pqr}^{neq}$ will be computed, Eq. (137) will either lead to the PR ($A_{pqr}^{neq}=A_{pqr}-A_{pqr}^{eq}$) or to the RR ($A_{pqr}^{neq}=A_{pqr}^{(1)}$ computed via Eq. (89)) collision models.

One can see that $f^{\mathrm{HM,Q19}}_{i}$ and $f^{\mathrm{GH,Q19}}_{i}$ lead to two families of LBMs that are based on two different kinds of equilibrium states (see App. H). The origin of this mismatch comes from (1) the weights which differ between the D3Q19 and the D3Q27 lattices, and (2) the difference between high-order contributions that are discarded. The former point is crystal clear when one looks at the values of the weights (see Tab. 2). Regarding the latter point, it flows from the fact that when neglecting $M_{222}$ in the unweighted formalism, then it affects, for example, $A_{222}$ and all other even Hermite moments ($A_{220}$, $A_{200}$, $A_{000}$, and their cyclic permutations). On the contrary, only the contribution of $A_{222}$ is discarded in the Gauss-Hermite formalism. The same reasoning applies to all moments that are neglected when moving from the D3Q27 formulation to the D3Q19 one.

What is interesting about these two families of LBMs is that one can derive new collision models by simply switching from $A_{pqr}$ to the statistical quantity of interest. As an example, one could derive a new D3Q19-K-LBM based on either $f^{\mathrm{HM,Q19}}_{i}$ or $f^{\mathrm{GH,Q19}}_{i}$ by first moving from $A_{pqr}$ to $M_{pqr}$, then to $\widetilde{M}_{pqr}$, and finally to $K_{pqr}$ using successively Eqs. (112), (113) and (117).

In the end, one may wonder what formalism should be preferred when it comes to the D3Q19 lattice. In fact, populations obtained in the unweighted formalism (131)-(134) can be shown to be solutions of the problem (43). This leads to the recovery of the proper set of $19$ moments, and as a consequence, these models have the correct macroscopic behavior. On the contrary, the present Gauss-Hermite formulation (137) leads to spurious coupling between hydrodynamic and high-order moments, and a fortiori, to an improper macroscopic behavior. Indeed,

[TABLE]

meaning the viscous stress tensor will depend on fourth-order nonequilibrium moments. Further couplings can be observed among third-order moments,

[TABLE]

and also fourth-order moments,

[TABLE]

These coupling defects are usually not encountered since, in most cases, the Hermite polynomial expansion is done up to the second-order for the D3Q19 lattice Shan et al. (2006); Krüger et al. (2017). Nevertheless, a few authors have also proposed to include third-order terms in its equilibrium state to improve the macroscopic behavior of the resulting LBM Házi and Kávrán (2006); Dellar (2014); Feng et al. (2016).

There are several possibilities to deal with these spurious couplings. First, one can impose third- and fourth-order terms to their equilibrium value, which leads to a PR collision model based on an extended equilibrium state. One then ends up with the correct viscous stress tensor. The second option consists in discarding fourth-order moments, and using orthogonal third-order Hermite polynomials Jacob et al. (2018),

[TABLE]

Equilibriating fourth-order moments also works here. The last way to obtain a stable numerical scheme is to rely on a fully orthogonal basis. By orthogonalizing fourth-order Hermite polynomials with second-order ones using the weighted scalar product, the following formulas are derived

[TABLE]

with the superscript $4o2$ standing for ‘order $4$ orthogonal with respect to order $2$’. By further orthogonalizing them between each other, one obtains

[TABLE]

where the superscript $\mathrm{FO}$ means ‘Fully Orthogonal’. Eventually, by combining Eqs. (141) and (143) with the second order polynomial basis, the new D3Q19-GH formulation of the populations reads

[TABLE]

where $A_{pqr}^{\mathrm{FO}}$ are Hermite polynomial coefficients corresponding to the fully orthogonal set of Hermite polynomials $H_{i,pqr}^{\mathrm{FO}}$. Among the three above propositions, this last model was identified as the most (linearly) stable one when coupled with the BGK collision operator (all relaxation frequencies equal to $\omega_{\nu}$).

Appendix H Equilibrium distribution functions

This appendix is dedicated to the equilibrium state recovered by all LBMs considered in the present work. Starting from the one-dimensional case, it was demonstrated that all LBMs recovered the following formulas (Sec. IV)

[TABLE]

for the D3Q27 velocity discretization, where $C=1-c_{s}^{2}$. For all above formulas, it is considered that $(\sigma,\lambda,\delta)=\{\pm 1\}^{3}$.

Knowing that all models recover the same equilibrium state in the one-dimensional case, the above construction then ensures that all LBMs also share the exact same equilibrium state for both two- and three-dimensional cases. In particular, this is true for models based on orthogonal bases (MRT and TRT), the regularized LBMs and the K-LBM. For regularized LBMs, it is explained by the fact that in the one-dimensional case, either solving Eq. (73) or using the Gauss-Hermite quadrature formulation (63) leads to the very same formalism. Tensor product rules compiled in Eqs. (67)-(69) allow us to extend the previous result to both two- and three-dimensional cases. For cumulants, deviations with respect to the CM formalism can only come from equilibrium cumulants of order $n\geq 4$, but they are null by definition if $n\geq 3$.

One must be careful regarding the equilibrium state derived in the case of the D3Q19 lattice. Indeed, if one starts from the RM framework and discards equilibrium raw moments via instructions given in Sec. VI.2, then one obtains the following type of equilibrium state

[TABLE]

where the lattice constant has been replaced with its value $c_{s}^{\mathrm{Q19}}=1/\sqrt{3}$. Unfortunately, these equilibrium VDFs do not match anymore those obtained through the Gauss-Hermite quadrature,

[TABLE]

where both the lattice constant and weights have been replaced with their value (Tab. 2) to ease the comparison with $f^{eq,\mathrm{RM}}_{i}$.

From this, one obtains two different kinds of equilibrium states whose difference scales as $\mathcal{O}(u^{2})$. This means that these equilibrium VDFs should lead to different results even for low Mach number flows. One can further expect that they will also result in different linear stability domains for the D3Q19 lattice, as it was already confirmed for the D2Q9 lattice Coreixas (2018).