Recursive regularization step for high-order lattice Boltzmann methods

TL;DR

This paper introduces a recursive regularization step for high-order lattice Boltzmann methods, significantly improving their stability and accuracy in simulating complex fluid flows, especially under high Reynolds numbers.

Contribution

The paper develops a recursive computation of Hermite coefficients for LBM, reducing computational cost and enhancing stability and accuracy over standard regularization methods.

Findings

Enhanced stability range demonstrated in shear layer simulations at high Reynolds numbers.

Recursive regularization outperforms standard methods in thermal and compressible flow simulations.

Confirmed stability improvements through Sod shock tube tests with high-order lattices.

Abstract

A lattice Boltzmann method (LBM) with enhanced stability and accuracy is presented for various Hermite tensor-based lattice structures. The collision operator relies on a regularization step, which is here improved through a recursive computation of non-equilibrium Hermite polynomial coefficients. In addition to the reduced computational cost of this procedure with respect to the standard one, the recursive step allows to considerably enhance the stability and accuracy of the numerical scheme by properly filtering out second (and higher) order non-hydrodynamic contributions in under-resolved conditions. This is first shown in the isothermal case where the simulation of the doubly periodic shear layer is performed with a Reynolds number ranging from $10^4$ to $10^6$, and where a thorough analysis of the case at $Re=3\times 10^4$ is conducted. In the latter, results obtained using both regularization steps are compared against the BGK-LBM for standard (D2Q9) and high-order (D2V17 and D2V37) lattice structures, confirming the tremendous increase of stability range of the proposed approach. Further comparisons on thermal and fully compressible flows, using the general extension of this procedure, are then conducted through the numerical simulation of Sod shock tubes with the D2V37 lattice. They confirm the stability increase induced by the recursive approach as compared with the standard one.

Tables

Table 1. Table 1: Description of some common standard and high-order two-dimensional lattice structures (top), and their associated Hermite tensor basis (bottom). For each lattice structure, the convention E D , Q V superscript subscript 𝐸 𝐷 𝑄 𝑉 E_{D,Q}^{V} is adopted to summarize all the characteristics of interest, namely, the number of discrete velocities V 𝑉 V , the quadrature order Q 𝑄 Q , and the number of dimensions D 𝐷 D   Shan ( 2016 ) . Furthermore, p 𝑝 p stands for the number of discrete speeds of each velocity group, while their associated weights ω i subscript 𝜔 𝑖 \omega_{i} compose the right part of the table. The last row consists of the value of the normalization constant (the inverse of the lattice constant c s subscript 𝑐 𝑠 c_{s} ) needed for the on-grid property of the lattice structure ( ( ξ i , x , ξ i , y ) ∈ ℤ 2 subscript 𝜉 𝑖 𝑥 subscript 𝜉 𝑖 𝑦 superscript ℤ 2 (\xi_{i,x},\xi_{i,y})\in\mathbb{Z}^{2} ). Regarding the Hermite tensors, they are classified into two categories: those belonging to the basis (❍) and those which do not (✗). Lattice structures data are compiled from  Succi ( 2001 ); Philippi  et al. ( 2006 ) .

Group	${𝝃}_i$	$p$	$E_{2,5}^9$	$E_{2,7}^{17}$	$E_{2,9}^{37}$
$1$	$(0,0)$	$1$	$4/9$	$\left(575+193\sqrt{193}\right)/8100$	$0.23315066913235250228650$
$2$	$(1,0)$	$4$	$1/9$	$\left(3355-91\sqrt{193}\right)/18000$	$0.10730609154221900241246$
$3$	$(1,1)$	$4$	$1/9$	$\left(655+17\sqrt{193}\right)/27000$	$0.05766785988879488203006$
$4$	$(2,0)$	$4$	$1/36$		$0.01420821615845075026469$
$5$	$(2,1)$	$8$			$0.00535304900051377523273$
$6$	$(2,2)$	$4$		$\left(685-49\sqrt{193}\right)/54000$	$0.00101193759267357547541$
$7$	$(3,0)$	$4$		$\left(1445-101\sqrt{193}\right)/162000$	$0.00024530102775771734547$
$8$	$(3,1)$	$8$			$0.00028341425299419821740$
$\mathrm{\hspace{0.25em}\hspace{0.25em}1}/c_s$			$\sqrt{3}$	$\sqrt{5\left(25+\sqrt{193}\right)/72}$	$1.19697977039307435897239$

Table 2. Table 2: Description of some common standard and high-order three-dimensional lattice structures (top), and their associated Hermite tensor basis (bottom). For each lattice structure, the convention E D , Q V superscript subscript 𝐸 𝐷 𝑄 𝑉 E_{D,Q}^{V} is adopted to summarize all the characteristics of interest, namely, the number of discrete velocities V 𝑉 V , the quadrature order Q 𝑄 Q , and the number of dimensions D 𝐷 D   Shan ( 2016 ) . Furthermore, p 𝑝 p stands for the number of discrete speeds of each velocity group, while their associated weights ω i subscript 𝜔 𝑖 \omega_{i} compose the right part of the table. The last row consists of the value of the normalization constant (the inverse of the lattice constant c s subscript 𝑐 𝑠 c_{s} ) needed for the on-grid property of the lattice structure ( ( ξ i , x , ξ i , y , ξ i , z ) ∈ ℤ 3 subscript 𝜉 𝑖 𝑥 subscript 𝜉 𝑖 𝑦 subscript 𝜉 𝑖 𝑧 superscript ℤ 3 (\xi_{i,x},\xi_{i,y},\xi_{i,z})\in\mathbb{Z}^{3} ). Regarding the Hermite tensors, they are classified into two categories: those belonging to the basis (❍) and those which do not (✗). Lattice structures data are compiled from  Succi ( 2001 ); Shan ( 2016 ) .

Group	${𝝃}_i$	$p$	$E_{3,5}^{19}$	$E_{3,5}^{27}$	$E_{3,7}^{39}$	$E_{3,9}^{103}$
$1$	$(0,0,0)$	$1$	$1/3$	$8/27$	$1/12$	$0.032633351764471159466$
$2$	$(1,0,0)$	$6$	$1/18$	$2/27$	$1/12$	$0.097656833590334574221$
$3$	$(1,1,0)$	$12$	$1/36$	$1/54$
$4$	$(1,1,1)$	$8$		$1/216$	$1/27$	$0.028097750290257335627$
$5$	$(2,0,0)$	$6$			$2/135$	$0.001045259560430061466$
$6$	$(2,1,0)$	$24$				$0.005705329016894815990$
$7$	$(2,2,0)$	$12$			$1/432$	$0.000611939269829747839$
$8$	$(2,2,2)$	$8$				$0.000155964159374283722$
$9$	$(3,0,0)$	$6$			$1/1620$	$0.000284443251800055207$
$10$	$(3,1,1)$	$24$				$0.000130698375985191585$
$11$	$(3,3,3)$	$8$				$0.000001223194501323058$
$\mathrm{\hspace{0.17em}1}/c_s$			$\sqrt{3}$	$\sqrt{3}$	$\sqrt{3/2}$	$1.19697977039307435897239$

Equations

$${\partial }_{t}f+\mathbf{\xi }\cdot \nabla f={\Omega }_f,$$

$$\left\{\begin{array}[]{r @{\: = \:} l}\rho\>=\>&\displaystyle{\int f\,\mathrm{d}\bm{\xi}},\vspace{0.1cm}\\ \rho\bm{u}\>=\>&\displaystyle{\int f\bm{\xi}\,\mathrm{d}\bm{\xi}},\vspace{0.1cm}\\ 2\rho E\>=\>&\displaystyle{\int f\xi^{2}\,\mathrm{d}\bm{\xi}},\end{array}\right.$$

$$\int {\Omega }_f\Phi (\mathbf{\xi })\mathrm{d}\mathbf{\xi }=0,$$

$$f^{(eq)}=\frac{\rho }{(2\pi rT{)}^{D/2}}\exp \left(-\frac{c^2}{2rT}\right),$$

$${\Omega }_f^{BGK}=-\frac{1}{\tau }\left(f-f^{(eq)}\right),$$

$$f(\mathbf{x},\mathbf{\xi },t)=\omega (\mathbf{\xi })\sum _{n=0}^{\infty }\frac{1}{n!(r{T}_0{)}^n}{\mathbf{a}}^{(n)}(\mathbf{x},t):{\mathcal{H}}^{(n)}(\mathbf{\xi }),$$

$${\mathbf{a}}^{(n)}(\mathbf{x},t)=\int f(\mathbf{x},\mathbf{\xi },t){\mathcal{H}}^{(n)}(\mathbf{\xi })\mathrm{d}\mathbf{\xi },$$

$${\mathcal{H}}^{(n)}=\frac{(-r{T}_0{)}^n}{\omega (\mathbf{\xi })}{\nabla }_{\mathbf{\xi }}^n\omega (\mathbf{\xi })$$

$$\omega (\mathbf{\xi })=\frac{1}{(2\pi r{T}_0{)}^{D/2}}\exp \left(-\frac{{\xi }^2}{2r{T}_0}\right),$$

$$⟨g|h⟩\equiv \int \omega (\mathbf{\xi })g(\mathbf{\xi })h(\mathbf{\xi })\mathrm{d}\mathbf{\xi }.$$

$$\centering\left\{\begin{array}[]{r @{\: = \:} l}\bm{a}^{(0)}\>=\>&\rho,\vspace{0.1cm}\\ \bm{a}^{(1)}\>=\>&\rho\bm{u},\vspace{0.1cm}\\ \bm{a}^{(2)}\>=\>&\bm{\Pi}+\rho(\bm{u}^{2}-\bm{\delta}),\vspace{0.1cm}\\ \bm{a}^{(3)}\>=\>&\bm{Q}+\bm{ua}^{(2)}+(1-D)\rho\bm{u}^{3},\vspace{0.1cm}\\ \bm{a}^{(4)}\>=\>&\bm{R}-\bm{P\delta}+\bm{\delta}^{2},\end{array}\right.\@add@centering$$

$$\mathbf{\Pi }=\int f{\mathbf{c}}^2\mathrm{d}\mathbf{\xi },\mathbf{Q}=\int f{\mathbf{c}}^3\mathrm{d}\mathbf{\xi },\mathbf{R}=\int f{\mathbf{c}}^4\mathrm{d}\mathbf{\xi }.$$

$${\partial }_t\left({\mathbf{a}}^{(n)}\right)+\nabla \cdot \left({\mathbf{a}}^{(n+1)}\right)+r{T}_0\nabla {\mathbf{a}}^{(n-1)}={\Omega }^{BGK}.$$

$$\centering\left\{\begin{array}[]{r @{\: = \:} l}\bm{a}_{eq}^{(0)}\>=\>&\rho,\vspace{0.1cm}\\ \bm{a}_{eq}^{(1)}\>=\>&\rho\bm{u},\vspace{0.1cm}\\ \bm{a}_{eq}^{(2)}\>=\>&\rho\big{[}\bm{u}^{2}+rT_{0}\left(\theta-1\right)\bm{\delta}\big{]},\vspace{0.1cm}\\ \bm{a}_{eq}^{(3)}\>=\>&\rho\big{[}\bm{u}^{3}+rT_{0}\left(\theta-1\right)\bm{u\delta}\big{]},\vspace{0.1cm}\\ \bm{a}_{eq}^{(4)}\>=\>&\rho\big{[}\bm{u}^{4}+rT_{0}\left(\theta-1\right)\bm{u}^{2}\bm{\delta}+(rT_{0})^{2}\left(\theta-1\right)^{2}\bm{\delta}^{2}\big{]},\end{array}\right.\@add@centering$$

$${\partial }_t=ϵ{\partial }_{t_1}+{ϵ}^2{\partial }_{t_2},\nabla =ϵ{\nabla }_1.$$

$${\mathbf{a}}_0^{(n)}={\mathbf{a}}_{eq}^{(n)},$$

$${\partial }_{t_1}{\mathbf{a}}_0^{(n)}+{\nabla }_1\cdot \left({\mathbf{a}}_0^{(n+1)}\right)+r{T}_0{\nabla }_1{\mathbf{a}}_0^{(n-1)}=0,$$

$${\partial }_{t_2}{\mathbf{a}}_0^{(n)}=0.$$

$$\begin{array}[]{c}\partial_{t}{}(\rho)+\bm{\nabla}\cdot(\rho\bm{u})=0,\vspace{0.1cm}\\ \partial_{t}(\rho\bm{u})+\bm{\nabla}\cdot(\rho\bm{u}^{2}+p\bm{\delta})=0,\vspace{0.1cm}\\ \partial_{t}(\rho E)+\bm{\nabla}\cdot[(\rho E+p)\bm{u}]=0,\end{array}$$

$$\begin{array}[]{l @{,\:} r}f=f^{(0)}+f^{(1)}\hfil,\>&f^{(0)}\gg f^{(1)}\sim\mathcal{O}(\epsilon),\vspace{0.1cm}\\ \bm{a}^{(n)}=\bm{a}_{0}^{(n)}+\bm{a}_{1}^{(n)}\hfil,\>&\bm{a}_{0}^{(n)}\gg\bm{a}_{1}^{(n)}\sim\mathcal{O}(\epsilon),\end{array}$$

$${\partial }_{t_1}{\mathbf{a}}_0^{(n)}+{\nabla }_1\cdot \left({\mathbf{a}}_0^{(n+1)}\right)+r{T}_0{\nabla }_1{\mathbf{a}}_0^{(n-1)}=-\frac{{\mathbf{a}}_1^{(n)}}{\tau },$$

$${\partial }_{t_2}{\mathbf{a}}_0^{(n)}+{\partial }_{t_1}{\mathbf{a}}_1^{(n)}+{\nabla }_1\cdot \left({\mathbf{a}}_1^{(n+1)}\right)+r{T}_0{\nabla }_1{\mathbf{a}}_1^{(n-1)}=0.$$

$$\begin{array}[]{c}\partial_{t}(\rho)+\bm{\nabla}\cdot(\rho\bm{u})=0,\vspace{0.1cm}\\ \partial_{t}(\rho\bm{u})+\bm{\nabla}\cdot(\rho\bm{u}^{2}+p\bm{\delta})=-\bm{\nabla}\cdot\bm{a}_{1}^{(2)},\vspace{0.1cm}\\ \partial_{t}(\rho E)+\bm{\nabla}\cdot[(\rho E+p)\bm{u}]=-\dfrac{1}{2}\mathrm{Tr}\left(\bm{\nabla}\cdot\bm{a}_{1}^{(3)}\right),\end{array}$$

$$\bm{a}_{1}^{(0)}=\bm{a}_{1}^{(1)}=\mathrm{Tr}\big{(}\bm{a}_{1}^{(2)}\big{)}=\bm{0}.$$

$$\begin{array}[]{l @{\: = \:} l @{\: = \:} l @{\: \equiv\:} l}\bm{\Pi}^{(1)}\hfil\>=\>&\displaystyle{\int f^{(1)}\bm{c}^{2}\mathrm{d}\bm{\xi}}\hfil\>=\>&\displaystyle{\int f^{(1)}\bm{\xi}^{2}\mathrm{d}\bm{\xi}}+0\hfil\>\equiv\>&\bm{a}_{1}^{(2)},\vspace{0.1cm}\\ \bm{Q}^{(1)}\hfil\>=\>&\displaystyle{\int f^{(1)}\bm{c}^{3}\mathrm{d}\bm{\xi}}\hfil\>=\>&\displaystyle{\int f^{(1)}\bm{\xi}^{3}\mathrm{d}\bm{\xi}}+0\hfil\>\equiv\>&\bm{a}_{1}^{(3)},\end{array}$$

$$\left\{\begin{array}[]{r @{\: = \:} l}\bm{\Pi}^{(1)}\>=\>&-\tau p\>\left[\bm{S}-\left(\frac{2}{D}\bm{\nabla}\cdot\bm{u}\right)\bm{\delta}\right],\vspace{0.1cm}\\ \bm{q}^{(1)}\>=\>&-\tau pc_{p}\bm{\nabla}T+\bm{u}\cdot\bm{\Pi}^{(1)},\end{array}\right.$$

$$\begin{array}[]{c}\partial_{t}(\rho)+\bm{\nabla}\cdot(\rho\bm{u})=0,\vspace{0.1cm}\\ \partial_{t}(\rho\bm{u})+\bm{\nabla}\cdot(\rho\bm{u}^{2})=\bm{\nabla}\cdot(\bm{\sigma}),\vspace{0.1cm}\\ \partial_{t}(\rho E)+\bm{\nabla}\cdot[\rho E\bm{u}]=\bm{\nabla}\cdot(\lambda^{\prime}\bm{\nabla}T)+\bm{\nabla}\cdot(\bm{\sigma}\cdot\bm{u}),\end{array}$$

$${\Omega }^{MRT}=-\sum _{n=0}^{\infty }\frac{1}{{\tau }_n}\frac{1}{n!(r{T}_0{)}^n}{\mathbf{a}}^{(n)}:{\mathcal{H}}^{(n)}.$$

$$f^N=\omega \sum _{n=0}^N\frac{1}{n!(r{T}_0{)}^n}{\mathbf{a}}^{(n)}:{\mathcal{H}}^{(n)}.$$

$${\mathbf{a}}^{(M)}=\int {\mathcal{H}}^{(M)}f\mathrm{d}\mathbf{\xi },$$

Full text

Regularization, Stability, High-order LBM, Compressible, Shock.

I INTRODUCTION

Over the past two decades, the Lattice Boltzmann Method (LBM) has emerged as an interesting candidate for Computational Fluid Dynamics (CFD) and beyond. Despite a first restriction to isothermal and weakly compressible flows, its range of applicability in both physics and engineering has grown in such a way that it is now possible to simulate very complex phenomena including turbulence Yu et al. (2005); Sagaut (2010), combustion Ponce Dawson et al. (1993); Filippova and Hänel (2000); Lin et al. (2016), multiphase interactions Shan and Chen (1993, 1994); Chen et al. (2014); Mazloomi et al. (2015), hemodynamics Ye et al. (2016), magnetohydrodynamics Dellar (2002, 2005, 2009, 2010), relativistic flows Succi et al. (2014); Mohseni et al. (2015) and even quantum systems Succi (2002); Lapitski and Dellar (2011).

From the numerical point of view, the LBM requires a strong coupling between discretizations of the velocity and the physical spaces. This is usually done using a cartesian grid coupled with an octree-based refinement technique Lagrava (2012). Combining these numerical tools with the particulate nature and the local dynamics of the LBM, the simulation of complex phenomena around realistic geometries is greatly eased Sengissen et al. (2015); Casalino et al. (2014); Ribeiro et al. (2017). Moreover, the method is simple to implement, induces a very low computational cost per degree of freedom, and presents a compact stencil, all of which contribute to its intrinsic advantage for parallel computations Schornbaum and Rüde (2016). All these key points make the LBM of great interest for both academic and industry groups.

The LBM derives from the Boltzmann Equation (BE), the milestone of the kinetic theory of gases Boltzmann (1872). The BE describes the balance between the transport and the collision of packets of particles through the evolution of the Velocity Distribution Function (VDF). The latter, usually written $f(\bm{x},\bm{\xi},t)$, can be seen as the probability density of finding fictive particles at a given point $(\bm{x},t)$ and with a given speed $\bm{\xi}$. Hydrodynamic variables (such as density $\rho$, momentum $\rho\bm{u}$ and total energy $\rho E$) are then recovered through the computation of the average, over the velocity space, of their mesoscopic counterparts. The Lattice Boltzmann Equation (LBE) results from a velocity discretization of the BE, meaning that the degrees of freedom allowed to the transport phenomenon are restricted to a finite set of velocities (directions and norms are fixed).

In order to recover the proper macroscopic set of equations from the LBE, the lattice of discrete speeds $\bm{\xi}_{i}$ and their associated equilibrium states $f^{eq}_{i}$ must be chosen accordingly. To do so, three main approaches are proposed in the literature.

The first suggestion is based on the expansion of the equilibrium state $f^{eq}_{i}$ into a polynomial, which is performed through a Taylor expansion, assuming that the Mach number is small. The polynomial coefficients are defined a posteriori by requiring that the continuum limit of the LBM recovers the Navier-Stokes-Fourier set of equations using a predefined lattice of discrete velocities Qian et al. (1992). This way of building a LBM does not allow to easily link lattices with their associated equilibrium states in a systematic way, especially when the complexity of the macroscopic behavior of interest increases.

The second approach is called the Entropic LBM (ELBM) and relies on ensuring the validity of the H-theorem using a discretized velocity space. This is done through a careful evaluation of both the equilibrium state $f^{eq}_{i}$ and the relaxation time $\tau$ Karlin et al. (2013); Frapolli et al. (2014, 2016), leading to an unconditionally stable LBM. Despite a broad range of applicability, the ELBM needs to solve a minimization problem for any grid point and at any time step, hence highly increasing the computational cost per degree of freedom. Such a drawback was recently overcome through the use of an approximated analytic solution of the minimization problem Karlin et al. (2014), but there seems to be some theoretical work remaining to properly define the validity range of this approximation Mattila et al. (2015).

Finally, the last idea originates from Grad’s work Grad (1949a) and was reintroduced into the LBM framework more recently Shan and He (1998); Shan et al. (2006); Philippi et al. (2006). It relies on the projection of the BE onto the Hilbert space spanned by Hermite polynomials. Such a procedure allows to build a systematic link between the kinetic theory, fluid mechanics and beyond. Indeed, once the macroscopic behavior of interest is chosen, the associated Hermite polynomial basis is then fixed. Typically, a second-order approximation is mandatory for weakly compressible and isothermal flows, a third-order one for isothermal flows, a fourth-order for thermal and fully compressible flows (Navier-Stokes-Fourier) and so on. The final step is to build the associated lattice of speeds by solving a linear system of equations. This system is constrained in order to keep the orthogonality properties of the previously chosen Hermite polynomial basis. This way of proceeding allows to straightforwardly build high-order LBMs Philippi et al. (2006).

The present work is based on this Hermite basis framework and focuses on the regularization of the pre-collision VDFs, which consists in filtering out non-hydrodynamic contributions of the streaming step to stabilize the LBM. Originally, this Projection-based Regularization (PR) procedure was designed for the simulation of particulate suspension Ladd and Verberg (2001). In fact, the stabilization property of the PR was only recently understood. Latt & Chopard Latt and Chopard (2005, 2006) emphasized its importance for the simulation of flows at high Reynolds number, while Chen and coworkers Chen et al. (2006) focused on flows at high Knudsen number. It is now one of the standard stabilization procedures available for both academic studies Zhang et al. (2006); Niu et al. (2007); Montessori et al. (2015) and engineering computations Sengissen et al. (2015); Casalino et al. (2014); Ribeiro et al. (2017). Even more recently, a Recursive Regularized (RR) collision operator was introduced by Malaspinas Malaspinas (2015) for isothermal LBMs only. The main improvement provided by the RR approach concerned its spectral properties and thus its enhanced stability range, regarding both Reynolds and Mach numbers, for a relatively low additional computational cost. Basically, the RR and PR processes differ in the way off-equilibrium Hermite coefficients (denoted $\bm{a}_{1}$ hereafter) are computed.

This paper studies the impact of non-equilibrium Hermite coefficients on the stability range of standard and high-order LBMs. Furthermore, a thorough mathematical derivation of the RR procedure extension to thermal and fully compressible flows is provided and numerically validated. The rest of the paper is organized as follows. In Sec. II, the recovery of the Navier-Stokes set of equations through the projection of the BE onto the Hermite polynomial basis and the concept behind both regularization steps are reminded. Sec. III is dedicated to numerics, and starts with the space/time discretization of the LBE (Sec. III.1). In Sec. III.2, the discrepancies between the stability range of the PR and the RR versions of the regularized collision operator are highlighted for the simulation of isothermal flows using standard and high-order lattices. In Sec. III.3, the extension of the RR step to thermal and fully compressible high-order LBMs is presented and illustrated on two configurations of the Sod shock tube. For the sake of completeness, appendixes on the mathematical proof of the Hermite coefficients recursive formulas (App. A & B), implementation details (App. C) and Hermite tensor basis for standard and high-order lattice structures (App. D) are finally provided.

II THEORETICAL BACKGROUND

In this paper, all quantities are defined in the $D$-dimensional Cartesian space $\mathbb{R}^{D}$ and by definition, for any vector $\bm{v}$, $v^{2}=\bm{v}\cdot\bm{v}$ is the square of the norm of $\bm{v}$. Properties concerning tensor products and Hermite polynomials in $\mathbb{R}^{D}$are based on Grad’s and Shan’s works Grad (1949b); Shan et al. (2006).

II.1 Boltzmann equation

In kinetic theory, gases are modeled by the VDF $f(\bm{x},\bm{\xi},t)$ describing the probability density of finding a fictive particle at position $\bm{x}$, time $t$ and with a mesoscopic velocity $\bm{\xi}$. When no external accelerations are considered, this VDF evolves through time and space in accordance to the force-free form of the BE:

[TABLE]

where $\cdot$ denotes the scalar product over $\mathbb{R}^{D}$, $\bm{\nabla}$ is the gradient operator associated to the physical space and $\Omega_{f}$ is the collision operator. The macroscopic quantities of interest (density $\rho$, momentum $\rho\bm{u}$, and total energy $\rho E$) are recovered averaging their mesoscopic counterparts over the velocity space:

[TABLE]

with integrals computed over $\mathbb{R}^{D}$. Hereafter, integration bounds will be omitted for the sake of clarity.

Regarding the collision model $\Omega_{f}$, it must satisfy the conservation of mass, momentum and total energy:

[TABLE]

with $\Phi(\bm{\xi})=(1,\bm{\xi},\xi^{2}/2)$. This collision process induces a relaxation of the VDF to the local thermodynamic equilibrium

[TABLE]

where $\bm{c}=\bm{\xi}-\bm{u}$, $r$ is the gas constant and $T$ the thermodynamic temperature. Throughout this paper, only linearized collision operators, such as the single relaxation time collision term of Bhatnagar-Gross-Krook (BGK) Bhatnagar et al. (1954)

[TABLE]

will be considered.

II.2 Projection onto the Hermite polynomial basis

In the present context, solutions of Eq. (1) are sought in the form of Hermite polynomials Shan and He (1998),

[TABLE]

where $:$ stands for the full contraction of indexes, $T_{0}$ is a reference temperature, $\bm{a}^{(n)}$ is the (tensor of) coefficient(s) related to the Hermite tensor $\bm{\mathcal{H}}^{(n)}$, both being $n-$rank tensors, and $\omega(\bm{\xi})$ is the weight function. They are defined as follows:

[TABLE]

where

[TABLE]

with

[TABLE]

$\bm{\nabla}^{n}_{\bm{\xi}}$ being the $n$-th derivative with respect to the velocity space. Eq. (6) can be seen as the decomposition of $f$ onto an orthogonal polynomial basis, since Hermite tensors are orthogonal with respect to the following scalar product:

[TABLE]

Thus $\bm{a}^{(n)}$ can simply be obtained as a projection of $f$ onto this orthogonal basis $\bm{a}^{(n)}=\langle\bm{\mathcal{H}}^{(n)}|f/\omega\rangle$. For the sake of clarity, the list of function variables will be omitted throughout the rest of the paper, except in Sec III.1 where the space/time discretization of the LBE is presented.

By construction, $\bm{a}^{(n)}$ can be linked to the familiar hydrodynamic moments Shan et al. (2006); Malaspinas (2009):

[TABLE]

where $\bm{\delta}$ is the identity matrix, $\bm{\delta}^{2}$ is the fourth-order identity tensor, and

[TABLE]

In order to create a systematic link between the BE and its macroscopic counterpart, Eq. (1) is projected onto the Hermite tensor basis using the projection operator $\langle\cdot|\cdot\rangle$ defined in Eq. (10):

[TABLE]

Here the BGK approximation (5) is adopted for the computation of the collision term, i.e., $\Omega^{BGK}=-(\bm{a}^{(n)}-\bm{a}_{eq}^{(n)})/\tau$, where $\bm{a}_{eq}^{(n)}=\langle\bm{\mathcal{H}}^{(n)}|f^{(eq)}/\omega\rangle$. The computation of these coefficients is straightforward since the expression of $f^{(eq)}$ is known. Up to the fourth order, one obtains

[TABLE]

with $\theta=T/T_{0}$, $T_{0}$ being a reference temperature.

II.3 Chapman-Enskog expansion

After expressing the BE in the Hermite tensor basis and defining the related coefficients $\bm{a}^{(n)}$, the macroscopic behavior linked to Eq. (13) is recovered by applying a separation of scales through the Chapman-Enskog expansion Chapman and Cowling (1970). For this purpose, the time derivative is expanded in powers of the Knudsen number $\epsilon$ as in Chapman and Cowling (1970):

[TABLE]

First, let us assume that $f$ and $\bm{a}^{(n)}$ are at equilibrium. Injecting Eq. (15) into Eq. (13), this separation of scales leads to

$\epsilon^{0}$:

[TABLE]

$\epsilon^{1}$:

[TABLE]

$\epsilon^{2}$:

[TABLE]

Adding all contributions, the macroscopic equations corresponding to the conservation of mass, momentum and total energy are finally recovered, using Eq. (14) for $n=0$, $n=1$ and $n=2$ respectively,

[TABLE]

where $p=\rho rT$ is the thermodynamic pressure. Hence, assuming $f$ and $\bm{a}^{(n)}$ are at equilibrium leads to Euler’s set of equations.

In a second step, the contributions of order $\mathcal{O}(\epsilon)$ are taken into account in the definition of $f$ and $\bm{a}^{(n)}$

[TABLE]

assuming the continuum limit $\epsilon\ll 1$. Applying the same multiscale analysis leads to

$\epsilon^{1}$:

[TABLE]

$\epsilon^{2}$:

[TABLE]

Therefore the set of macroscopic equations (19) becomes

[TABLE]

since the collision model must satisfy Eq. (3) or equivalently,

[TABLE]

Now, only $a_{1,\alpha\beta}^{(2)}$ and $a_{1,\alpha\beta\beta}^{(3)}$ remain to be computed. They can either be computed thanks to Eq. (13) or noticing that,

[TABLE]

where $\bm{\Pi}^{(1)}$ and $\bm{q}^{(1)}=\mathrm{Tr}(\mathbf{Q}^{(1)})/2$ are the standard second- and third-order off-equilibrium moments at the Navier-Stokes level Malaspinas (2009),

[TABLE]

with $\bm{S}=\bm{\nabla u}+(\bm{\nabla u})^{T}$, $\cdot^{T}$ standing for the transpose operator. $c_{p}=(1+D/2)r$ is the heat capacity at constant pressure.

Finally, injecting Eq. (26) in Eq. (23) allows to link the Hermite formulation of the BE (13) to the following macroscopic set of equations, namely, Navier-Stokes-Fourier equations,

[TABLE]

where $\bm{\sigma}=\overline{\bm{\Pi}}-p\bm{\delta}$ is the stress tensor, $\overline{\bm{\Pi}}=\mu\left[\bm{S}-\left(\frac{2}{D}\bm{\nabla}\cdot\bm{u}\right)\bm{\delta}\right]$ is the traceless viscous stress tensor, $\mu=\tau p$ and $\lambda^{\prime}=\tau pc_{p}$ being the dynamic viscosity and the thermal conductivity coefficients.

Before proceeding any further, several remarks can be pointed out.

First, only coefficients $\bm{a}_{0}^{(n)}$ ($\bm{a}_{1}^{(n)}$) up to the fourth order (third order) are needed to recover the Navier-Stokes set of equations. And more generally, the Chapman-Enskog expansion of the BE at order $k$ needs coefficients of the Hermite expansion to include terms at order $n+k$. This is mandatory to properly recover the hydrodynamic behavior of the BE Shan et al. (2006).

Second, a strong assumption is made by using the BGK collision operator, which allows only one parameter to represent the physical behavior of the fluid: the single relaxation time $\tau$. This choice induces a coupling between the momentum and the energy relaxation processes since the Prandtl number is fixed at $Pr=\mu c_{p}/\lambda^{\prime}=1$. To overcome this deficiency, a more sophisticated collision operator may be employed, such as the general Multi-Relaxation Time (MRT) collision term expressed in the Hermite tensor basis Shan and Chen (2007). The latter reads

[TABLE]

Here choosing $\tau_{2}=\mu/p$ and $\tau_{3}=\lambda/pc_{p}$ allows to tune the Prandtl number which now equals $Pr=\tau_{2}/\tau_{3}$. A double distribution function LBM could also be used to overcome this difficulty (see for example Guo et al. (2007)).

Third, as depicted in Eq. (27), the computed viscous stress tensor is traceless, which means that only fluids without bulk viscosity can be simulated. This is another consequence of the BGK approximation, and again this can be avoided by using a more sophisticated collision operator belonging to the Hermite polynomial expansion framework Krüger et al. (2017).

Last, the specific heat ratio has a fixed value $\gamma=(D+2)/D$. To overcome this issue, one can employ a second distribution function to take into account the energy evolution linked to internal degrees of freedom (rotational and vibrational) of molecules Rykov (1975); Nie et al. (2008a); Frapolli et al. (2016).

In the rest of the paper, only single-VDF-based LBMs coupled with the single-relaxation-time approximation are considered.

II.4 Truncation of the VDF

As briefly discussed in Sec. II.3, a finite expansion of the VDF in Hermite tensors, up to the fourth order, is sufficient to recover the macroscopic behavior of Navier-Stokes’ equations from the BE. Here particular attention is paid to error terms arising, at the macroscopic level, from a wrong truncation of the VDF. From now on, let $f^{N}$ denote the truncation of $f$, to the order $N$, of the Hermite polynomial development

[TABLE]

The key element is to ensure that this truncation still allows the proper conservation of the macroscopic moment of the VDF. Introducing $\bm{a}^{(M)}$ as the $M$th-order moment of $f$:

[TABLE]

orthogonality properties of the Hermite polynomials lead to Shan et al. (2006)

[TABLE]

This means that $f$ and $f^{N}$ share the same moments up to the order $M$, if the truncation order $N$ is at least equal to the highest moment $M$ that we need to conserve.

As mentioned in Sec. II.3, moments of $f^{(0)}$ up to the fourth order are sufficient to recover the proper macroscopic behavior of the BE. Consequently, truncation errors appear at the macroscopic level when the equilibrium VDF is truncated to an order lower than $N=4$. These error terms can be evaluated quite easily using Eq. (13). For $N=3$, the only term that cannot be used for the computation of $\bm{a}_{1}^{(3)}$ is $\bm{\nabla}\cdot\bm{a}_{0}^{(4)}$. Hence, the error term introduced by a third-order truncation modifies the heat flux $\bm{q}=\mathrm{Tr}(\mathbf{Q})/2$:

[TABLE]

For $N=2$, too many error terms are introduced in the energy equation and the concept of temperature should not be used. This is why this truncation is restricted to isothermal, or more precisely athermal LBMs ($\theta=1$). Furthermore, the computation of $\bm{a}_{1}^{(2)}$ is also affected by the truncation. As for the heat flux when $N=3$, the deviation also comes from the divergence term $\bm{\nabla}\cdot(\bm{a}_{0}^{(3)})$ which cannot be used anymore for the computation of $\bm{a}_{1}^{(2)}$. The viscous stress tensor is thus modified:

[TABLE]

The famous $\mathcal{O}(Ma^{3})$ error term is then recovered from the above procedure, explaining why second-order truncation of the VDF are restricted to the simulation of weakly compressible flows Qian et al. (1992).

In the most general case, the truncation criterium (31) becomes Shan et al. (2006)

[TABLE]

if $N+k\geq M$. Here, we recover the fact that the Navier-Stokes macroscopic behavior is achieved if $N=4$ at equilibrium ($k=0$), while $N=3$ is sufficient for the non-equilibrium VDF $f^{(1),N}$. From the macroscopic point of view, there is no reason why $f^{(0),N}$ and $f^{(1),N}$ should be truncated at the same order $N$. But it will be shown hereafter (Sec. II.5) that keeping the same order $N$ for both $f^{(0),N}$ and $f^{(1),N}$ is possible thanks to a recursive computation of these coefficients.

II.5 Regularized collision operator

Both truncations $f^{(0),N}$ and $f^{(1),N}$ will now be considered. The reason lies in the definition of the regularization step Latt and Chopard (2005, 2006); Chen et al. (2006). This collision operator is based on the regularization of the non-equilibrium part of the pre-collision distribution functions, and aims at filtering out non-hydrodynamic sources. To do so, $f^{N}$ is reconstructed before the collision step discarding $\mathcal{O}(\epsilon^{k})$ contributions ($k\geq 2$):

[TABLE]

with

[TABLE]

and where the sum begins at $n=2$ due to mass and momentum conservation (24). Post-collision VDFs are then defined as follows,

[TABLE]

Here coefficients $\bm{a}_{1}^{(n)}$ are the only missing information required to reconstruct $f^{(1),N}$. Originally in Latt and Chopard (2005, 2006), this stabilization technique was used for the simulation of isothermal and weakly compressible flows for which only second-order terms were kept. Off-equilibrium coefficients $\bm{a}_{1}^{(2)}$ were then computed projecting the VDFs onto the second-order Hermite polynomials

[TABLE]

assuming that $f^{N}-f^{(0),N}\approx f^{(1),N}$. But doing so, only contributions belonging to the Hilbert space that were not taken into account are filtered, while some non-hydrodynamic contributions are still hidden in $\bm{a}_{1}^{(2)}$. In the rest of the paper, this method will be referred as to the projection-based regularization (PR) process.

Later, Malaspinas Malaspinas (2015) proposed a complete regularization procedure based on recursive properties of the off-equilibrium coefficients, allowing an enhancement of accuracy and stability compared to BGK and standard MRT models. Regarding the recursive property of this stabilization procedure, it comes from the fact that non-equilibrium coefficients $\bm{a}_{1}^{(n)}$ are computed using a recursive formula flowing from the Chapman-Enskog expansion. This step provides a proper way to filter out non-hydrodynamic sources.

It must be understood that this recursive regularization (RR) approach was introduced in the context of isothermal and weakly compressible flows for standard LBMs, while in the present paper the recursive formula is extended to the thermal and fully compressible case (see App. A for its derivation) and applied to high-order LBMs. The latter reads

[TABLE]

where $\bm{a}_{1}^{(2)}$ and $\bm{a}_{1}^{(3)}$ can either be computed thanks to the projection of $f^{(1)}$ onto the Hermite polynomial basis or using finite differences. In the isothermal case ($\theta=1$), Malaspinas’ recursive relation Malaspinas (2015) is recovered for $n\geq 3$:

[TABLE]

Again, $\bm{a}_{1}^{(2)}$ coefficients can be computed by projection or using finite differences.

Furthermore, a correct evaluation of $\bm{a}_{1}^{(2)}$ requires a proper evaluation of $\bm{a}_{0}^{(3)}$ since

[TABLE]

Thus, in the isothermal case, the equilibrium VDF should theoretically be developed up to the third order. Similarly, in the thermal case, correct calculations of $\bm{a}_{1}^{(2)}$ and $\bm{a}_{1}^{(3)}$ require a development of $f^{(0)}$ up to the fourth order, which is no more binding than the condition to recover the thermal and fully compressible Navier-Stokes-Fourier equations.

II.6 Discretization of the velocity space

In order to numerically solve the BE, a discretization of the velocity space is necessary. It consists in keeping only a discrete set of $V$ velocities $\bm{\xi}_{i},\ i\in\llbracket 1,V\rrbracket$, ensuring the preservation of $f_{i}$’s moments, from the continuum velocity space to the discrete one. To do so, a Gauss-Hermite quadrature is applied Shan et al. (2006); Nie et al. (2008b):

[TABLE]

Here $\bm{\mathcal{H}}_{i}^{(M)}=\bm{\mathcal{H}}^{(M)}(\bm{\xi}_{i})$, $f_{i}^{(0),N}=\frac{\omega_{i}}{\omega(\bm{\xi}_{i})}f^{(0),N}(\bm{\xi}_{i})$, $\omega_{i}$’s are the gaussian weights of the quadrature, and $Q$ is the order of accuracy of the quadrature. The resolution of this quadrature problem aims at finding lattices, i.e., discrete weights $\omega_{i}$’s associated with discrete velocities $\bm{\xi_{i}}$’s, allowing the conservation of Hermite polynomials orthogonality properties up to a requested order $M$. Then, the $M$th-order moment can be exactly recovered if the equilibrium function is truncated to the order $N\geq M$, and if the quadrature order is greater than $2N-1$. Some common two- and three-dimensional lattices, satisfying Eq. (42) for several $N$, are detailed in App. D.

In the rest of the paper, the regularization procedure will be assessed on three lattices: the common D2Q9 lattice and the high-order D2V17 and D2V37 lattices. The D2V37 lattice, for which $Q=9$, allows the conservation of moments up to $M=4$. Therefore, this lattice is able to reproduce the macroscopic behavior of the fully compressible set of Navier-Stokes-Fourier equations. On the contrary, the D2V17 lattice flows from a seventh-order quadrature, which means that the equilibrium VDF can be truncated up to the third order only, leading to truncation error terms in the energy equation (II.4). This lattice should then be restricted to isothermal flows, even though the simulation of flows with weak temperature fluctuations could also be valid. Eventually, the case of the D2Q9 lattice, with a quadrature order of $Q=5$, allows developments up to the second order only, which leads to the famous compressibility error in the momentum equation (II.4).

II.7 Regularization step: Projection vs Recursivity

In the general case of a lattice including $V$ discrete speeds, the dimension of the associated Hilbert space is also $V$. Assuming the quadrature order of this lattice structure is $Q$, then Hermite polynomials up to the order $N=(Q-1)/2$ will be orthogonal to each other Nie et al. (2008b) and may form part of a basis $\mathcal{B}=\mathcal{B^{\mathcal{H}}}\cup\mathcal{B^{\overline{\mathcal{H}}}}$, where

[TABLE]

is a sub-set of $\mathcal{B}$ entirely composed of Hermite polynomials, while elements of $\mathcal{B^{\overline{\mathcal{H}}}}$ are linearly independent of each others and may not be Hermite polynomials. Using this decomposition, the polynomial coefficients of $f^{(1)}$ can also be recast into two sub-sets:

[TABLE]

where coefficients $\bm{b}_{1}^{(n)}$ ($n>N$) are supposedly related to non-hydrodynamic behaviors.

The purpose of the PR approach is to keep only Hermite polynomial coefficients, since the mathematical expression of $\bm{b}_{1}^{(n)}$ ($n>N$) is usually unknown. After the PR procedure, remaining polynomial coefficients of $f^{(1),N}_{i,PR}$ are

[TABLE]

where contributions from $\mathcal{B}^{\overline{\mathcal{H}}}$ have been completely filtered out. Nevertheless, spurious sources coming from the approximation $f_{i}^{(1)}\approx(f_{i}-f_{i}^{(0)})$ may still be hidden in $\bm{a}_{1,PR}^{(n)}$. Hence, this approach reduces the order of the polynomial development, and filters out spurious contributions originating from presumably non-Hermite polynomials. In the particular case of $\mathcal{B}^{\overline{\mathcal{H}}}$ being empty, this regularization step reduces to the standard BGK collision model if $f_{i}^{(1)}$ is projected onto the complete basis $\mathcal{B}$. We will see later that $\mathcal{B}_{D2Q9}$ and $\mathcal{B}_{D3Q27}$ belong to this particular case.

Regarding now the RR approach, it further filters out high-order contributions, left by the PR approach, recomputing most coefficients $\bm{a}_{1}^{(n)}$ ($n\leq N$) by a Chapman-Enskog expansion, and without assuming $f^{(1)}_{i}\approx(f_{i}-f^{(0)}_{i})$. For $f^{(1),N}_{i,RR}$ the remaining coefficients are then

[TABLE]

where $\bm{a}_{1,RR}^{(n)}=\bm{a}_{1,PR}^{(n)}$ for $n\leq 2$ ($n\leq 3$) in the isothermal (thermal) case, whereas high-order coefficients are recomputed using Eqs. (40) and (II.5). Eventually, while working on a complete basis, this approach is the only one leading to the expected filtering behavior.

For the velocity sets of interest, the following observations can be made. The D2V37 (D2V17) is built ensuring that Hermite polynomials orthogonality properties are preserved, up to $N=4$ ($N=3$), during the velocity space discretization Philippi et al. (2006). Thus it is known for sure that

[TABLE]

and

[TABLE]

whereas the true form of $\mathcal{B}_{D2V37}^{\overline{\mathcal{H}}}$ and $\mathcal{B}_{D2V17}^{\overline{\mathcal{H}}}$ are unknown. Hence, both PR and RR approaches first discard coefficients related to $\mathcal{B}_{D2V37}^{\overline{\mathcal{H}}}$ and $\mathcal{B}_{D2V17}^{\overline{\mathcal{H}}}$. Then, $\bm{a}_{1}^{(n)}$ ($n\leq N$) are computed either by the PR or by the RR approach. In the particular case of the D2Q9 lattice, Hermite polynomials up to $N=2$ were first considered for the polynomial expansion Qian et al. (1992), leading to

[TABLE]

Nevertheless, the development of the VDF can be extended including some third- and fourth-order terms, which also satisfy the orthogonality property conservation (42):

[TABLE]

and $\mathcal{B}^{\overline{\mathcal{H}},Complete}_{D2Q9}=\O$. In other words, the complete basis $\mathcal{B}_{D2Q9}$ can be derived thanks to the tensor properties of the D2Q9 lattice, which allows the orthogonality conservation of every second-order Hermite polynomial per direction. The same applies to the D3Q27 lattice.

Malaspinas has recently shown that developing $f_{i}^{(0),N}$ onto these complete basis helps reducing the $\mathcal{O}(Ma^{3})$ error term in the Navier Stokes’ equations by removing all non-diagonal terms Malaspinas (2015). The same development was also done for $f_{i}^{(1),N}$ during the regularization process. And even if it was not explicitly specified in Malaspinas (2015), one should notice that doing a full projection of the non-equilibrium VDF onto $\mathcal{B}_{D2Q9}$, in order to compute the 9 off-equilibrium coefficients, has no impact at all on the off-equilibrium VDF. As an analogy, it would be pointless to project a 2D vector onto the orthonormal basis of the 2D physical space using the Euclidean scalar product. Hence a complete PR approach would be absolutely useless, since it would not filter out any physical information.

The rest of the paper will aim at showing some interesting properties of the RR procedure:

It enhances numerical stability compared to the PR procedure for both standard and high-order LBMs, 2. 2.

It is the only way to filter out non-hydrodynamic spurious sources without discarding any Hermite coefficients, 3. 3.

The extension of the RR collision model (II.5) helps improving numerical stability for the simulation of thermal and fully compressible flows.

III NUMERICS & VALIDATION

In this section, the principle of the space/time discretization of the LBE is first recalled. Then, improvements regarding numerical stability and accuracy, induced by the use of the RR collision model, are confirmed for both standard and high-order LBMs. For the sake of clarity, the truncation notation $\cdot^{N}$ is dropped.

III.1 Key points on the space/time discretization

Let us start from the force-free lattice-Boltzmann equation with a general collision operator,

[TABLE]

This equation is a first-order partial differential equation. The LHS term is linear and corresponds to the convection of the discrete VDFs $f_{i}$’s at constant speed $\xi_{i}$, while the RHS term is non-linear and translates the rate of change of $f_{i}$’s induced by collisions. Since these two contributions behave differently from the mathematical point of view, two different time integrations are used: (a) the method of characteristics allows to exactly integrate the convection term between $t$ and $t+\Delta t$, whereas (b) the trapezoidal rule ensures a second-order accuracy in time for the collision term integration. These techniques lead to

[TABLE]

where the space discretization error $\mathcal{O}(\Delta x^{2})$ comes from (a): $\Delta\bm{x}=\bm{\xi_{i}}\Delta t$. This space/time discretization results in an implicit formulation since $\Omega_{i}\left(t+\Delta t\right)$ depends on $f\left(\bm{x}+\bm{\xi}_{i}\Delta t,t+\Delta t\right)$. Nevertheless, a change of variables compliant with the conservation of mass, momentum and total energy allows to get around it Dellar (2001)

[TABLE]

and leads to an explicit numerical scheme,

[TABLE]

To be fully consistent with the new set of VDFs $\overline{f_{i}}$’s, the collision operator needs to be slightly modified. For the BGK collision model,

[TABLE]

with $\overline{\tau}=\tau/\Delta t\,+\,1/2$. Dropping the overline notation, and adopting the shorthand notation $\Delta x=\Delta t=1$, we end up with the famous second-order accurate and explicit numerical scheme

[TABLE]

All results presented hereafter will be based on the coupling between this numerical scheme and several discretizations of the velocity space, namely, the D2Q9, the D2V17 and the D2V37 velocity sets. These lattice structures were built using an on-grid condition Philippi et al. (2006),

[TABLE]

allowing to achieve the full potential of the method of characteristics. The renormalization constant (or lattice constant) $c_{s}$ can further be identified as the isothermal lattice speed of sound, i.e., $c_{s}=\sqrt{rT_{0}}\Delta t/\Delta x$ Viggen (2014). Its value is recalled for each lattice structure in App. D. In the rest of the paper, all necessary quantities used to define each numerical case will be expressed in lattice and dimensionless units.

III.2 Isothermal LBM ($\theta=1$)

In a first time, the case of isothermal flows is considered. Lattice structures with increasing complexity (D2Q9, D2V17 and D2V37) are successively used to point out the capability of the RR procedure to deal with both standard and high-order LBMs in a straightforward way. All velocity set definitions are recalled in App. D.

The double shear layer is a well-known test case which allows to quantify the stability of numerical schemes as a first step Minion and Brown (1997). This flow is composed of two longitudinal shear layers, located at $y=L/4$ and $y=3L/4$, in a 2D doubly periodic domain with $(x,y)\in[0,L]^{2}$. Furthermore, a transverse perturbation is superimposed to the flow. This leads to the roll-up of the shear layers, and the generation of two counter-rotating vortices by the Kelvin-Helmholtz instability mechanism. In addition, any numerically induced disturbances may lead to the formation of further spurious vortices, or in the worst case, make the simulation reach its stability threshold. This is why this test case is an excellent candidate to evaluate the stability of numerical schemes. The initial state is defined by,

[TABLE]

and

[TABLE]

where $(x^{*},y^{*})=(x/L,y/L)$, and $u_{0}$ is the characteristic speed. $k$ is related to the width of the shear layers while $\delta$ controls the amplitude of the transverse perturbation. Here the case of thin shear layers, where $(k,\delta)=(80,0.05)$, is considered. The Reynolds number, which is the ratio between convective and diffusive phenomena, is first fixed to a moderate value of $Re=u_{0}L/\nu=3\times 10^{4}$ (cf Fig. 1). This is sufficient to reach the stability limit of the BGK collision model for standard LBMs when an under-resolved mesh ($L=128$) is considered Bösch et al. (2015).

The stability range of the proposed model is then evaluated varying the Reynolds number from $10^{4}$ to $10^{6}$. Simulations using the D2Q9, the D2V17 and the D2V37 lattices are performed with a free-stream Mach number $M_{0}=u_{0}/c_{s}=0.2,$ $0.35$ and $0.57$ respectively. These values are chosen in order to: (i) properly distinguish the impact of each collision operator on the numerical stability of each LBM, and (ii) reduce the impact of the $\mathcal{O}(Ma^{3})$ error encountered when the D2Q9 lattice is employed. The initialization step is achieved using the approximation $f_{i}\approx f_{i}^{(0)}+f_{i}^{(1)}$, with $f_{i}^{(1)}$ computed using analytic formulas for the velocity gradients. This allows to reduce spurious oscillations at the beginning of the simulation as demonstrated in Mattila et al. (2015). Finally, it should be noted that for all simulations below, $f_{i}^{(0)}$ is expanded to the maximal authorized order, i.e., $N=3$ for the D2V17, and $N=4$ for both the D2Q9 and the D2V37, using their associated Hermite tensor basis (see App. C & D for more details on this last point).

Extensive results concerning the case at $Re=3\cdot 10^{4}$ are compiled in Figs. 2, 3 and 4 for the D2Q9, the D2V17 and the D2V37 lattices respectively. All models are compared to a reference solution obtained using the BGK collision model with $L=2048$. Several partial conclusions can be drawn from these results.

Firstly, correct kinetic-energy-related evolutions are recovered for all computations, even using a relatively coarse mesh ($L=128$). Regarding mean and standard deviation of the enstrophy, a convergence study has been conducted, leading to the choice of a centered and fourth-order-accurate finite-difference scheme for the gradient evaluation. Nevertheless small errors are still observed, when the Mach number is increased, even for $L=256$. Thus recovering the proper evolution of the mean and the standard deviation of the enstrophy is more difficult than obtaining the correct evolution of the same quantities in the case of the kinetic energy. The study of the enstrophy evolution should then be preferred for the accuracy evaluation of LBMs.

Secondly, the new RR approach is more stable than the PR one, at least for the present lattices. As explained in Sec. II.7, this originates from a better computation of non-equilibrium Hermite coefficients $\bm{a}_{1}^{(n)}$ from the kinetic theory point of view. This point is further highlighted in Fig. 5, where the maximal achievable Mach number $M_{0}^{Max}$ allowing a stable simulation, up to $t/t_{c}=2$, is plotted for Reynolds numbers ranging from $10^{4}$ to $10^{6}$, in the under-resolved configuration $L=128$.

The stability criterium that has been chosen is based on the mean kinetic energy $\langle u^{2}\rangle$: a computation is considered to remain stable if $\langle u^{2}(t\leq 2t_{c})\rangle<\langle u^{2}(t=0)\rangle$. For a proper comparison, the standard PR and the most stable RR versions associated to each lattice structure are compared. Results obtained in the particular case of $Re=3\times 10^{4}$ seem to be extendable to a wide range of Reynolds number.

Thirdly, the PR 4 and the BGK collision models give exactly the same results in the particular case of the D2Q9. This confirms what was anticipated in Sec. II.7, i.e., when the VDF is expanded over the complete Hermite basis then the PR 4 and the BGK reduce to the same collision operator. Therefore the only way to reach the full potential of the D2Q9 lattice is to use the RR procedure.

Lastly, the computational overhead for our non-optimized implementation is $1.5\leq t_{RR}/t_{BGK}\leq 2$, the lower limit is for the D2Q9 and the upper for the D2V37. This is far from being excessive considering the tremendous improvements obtained in terms of numerical stability. Regarding the regularization step itself, the RR approach is always faster (about $20\%$) than the standard PR one.

III.3 Fully compressible LBM

To further validate the general RR process, the simulation of thermal and fully compressible flows is now considered. This is done through the numerical computation of the famous Sod shock tube Sod (1978), using the most compact lattice structure allowing to ensure the preservation of the orthogonality of all fourth-order Hermite tensors: the D2V37 velocity set Philippi et al. (2006); Shan (2016). For this purpose, Hermite coefficients $\bm{a}_{0}^{(n)}$ and $\bm{a}_{1}^{(n)}$ are computed thanks to our proposed extensions (74) and (II.5).

This 1D Riemann problem consists in a closed tube divided into two regions by a thin membrane. Each region is filled with the same gas but has different thermodynamic properties (density $\rho$, temperature $T$, pressure $P$ and velocity $u$). At the initialization, the breakdown of the membrane induces a strong acceleration of the flow, from the high-pressure side to the low-pressure one, whose purpose is to equalize the pressure inside the tube. This leads to the generation and the propagation of three characteristic waves: (1) the compression of the gas creates a shock wave which propagates towards the low-pressure side, (2) the expansion of the gas towards the high-pressure side induces the propagation of the expansion or rarefaction wave, and (3) the separation between the two waves, namely, the contact discontinuity. The latter can be seen as a fictitious diaphragm traveling at a constant speed towards the low-pressure side.

In this paper, two different configurations are studied. They share the same pressure ratio but differ when it comes to their temperature or density ones:

[TABLE]

[TABLE]

where subscripts $L$ and $R$ stand for the left and the right states respectively. To avoid the contribution of the boundary conditions, the computation takes place in a periodic domain of length $2L_{x}$ centered around the location $x=L_{x}/2$ where the discontinuity between the two states belongs. The simulation domain is then spatially discretized using $L_{x}=400$ grid points. Such a coarse mesh will allow to further highlight the numerical stability issues encountered computing discontinuities.

Results obtained using the PR and the RR processes, at order $N=4$, are plotted along $[0,L_{x}]$ in Fig. 6. Even though both models are able to properly reproduce the generation and the propagation of all the characteristic waves of this 1D Riemann problem, the PR procedure introduces a coupling between high-order and Navier-Stokes physics in the form of standing waves, whereas the RR procedure completely filters them out even if small over/undershoots still remain. The latter can be attenuated using either a finer grid or a shock-capturing technique, such as a shock sensor Bogey et al. (2009).

It must be noted that without regularization steps, the standard LBM encountered severe numerical stability issues for the present configurations, and could not be stabilized even using extremely fine meshes (more than $10000$ points in the longitudinal direction). This further points out that the RR collision operator allows to get more stable solutions without degrading the accuracy of the numerical scheme.

IV CONCLUSION

Despite a wide range of validity of the standard lattice Boltzmann method (LBM), the simulation of certain flows remains a tedious task: (1) weakly compressible flows at high Reynolds numbers, and (2) fully compressible flows including discontinuities such as shockwaves. In this context, the present work focuses on the extension of the recursive regularization step to high-order LBMs. New LBMs, with increased stability range, are then obtained filtering out the non-hydrodynamic contributions induced by the streaming step. This technique relies on two points: (a) the computation of non-equilibrium coefficients $\bm{a}_{1}$ through the Chapman-Enskog expansion, eased thanks to (b) the recursive properties of Hermite polynomial basis. This procedure was originally derived in the particular case of isothermal and weakly compressible flows simulations, with standard lattice structures such as the D2Q9 and the D3Q27. Here we further validate it using high-order LBMs (the D2V17 and the D2V37) for the simulation of both isothermal flows and thermal fully compressible flows. The latter was possible thanks to our derivation of general recursive formulas for the computation of Hermite coefficients $\bm{a}_{0}$ and $\bm{a}_{1}$. Strong improvements in terms of numerical stability are confirmed for both kind of simulations, and at a relatively low computational overhead.

Flows including external accelerations, boundary conditions, and LBMs with two sets of populations are currently under investigation. Theoretical derivations and numerical validations, concerning the extension of the RR procedure to deal with these topics, will be presented in the near future.

Acknowledgements.

C. C. wants to thank D. Ricot & O. Malaspinas for fruitful conversations about the lattice Boltzmann theoretical background, and regularization steps. This work was supported by the French funded project CLIMB in the framework of the Programme d’Investissement et d’Avenir - Calcul Intensif et Simulation Numérique.

Appendix A Proof of the recursive formula

The aim of this section is to prove the following recursive relation:

[TABLE]

where some mathematical notations are used for the sake of clarity: (1) $\beta_{i}$ is used when the index $\alpha_{i}$ is omitted, e.g., $a_{0,\beta_{i}}^{(n)}\equiv a_{0,\alpha_{1}..\alpha_{i-1}\alpha_{i+1}..\alpha_{n}}^{(n)}$, and (2) if $\alpha_{i}$ and $\alpha_{j}$ are omitted then $\beta_{ij}$ is used.

The following relations are needed to prove (II.5):

• •

The Hermite coefficient based LBE Shan et al. (2006)

[TABLE]

where Einstein summation notation is used on subscript $\gamma$,

• •

The recursive formula for Hermite coefficients at equilibrium (App. B)

[TABLE]

• •

Euler’s equations

[TABLE]

Let us now move on to the proof itself.

STEP 1: Boltzmann equation associated to Hermite coefficients (63) at order $(n-1)$.

[TABLE]

where the derivation by parts rule is used.

STEP 2: Remove all time derivatives $\partial_{t}$.

First, we use Eq. (64),

[TABLE]

for $n\geq 3$. Then, thanks to (63) one can replace $\partial_{t}a_{0,{\alpha_{1}..\alpha_{n}}}^{(n)}$ by

[TABLE]

and

[TABLE]

In the same spirit as for Chapman-Enskog expansion, Euler’s equation of momentum (Eq. (66)) is used to replace the time derivative $\partial_{t}u_{\alpha_{n}}$. In addition to that, removing the implicit summation on $\gamma$ leads to

[TABLE]

so that

[TABLE]

Let us now focus on $(iv)$ and $(v)$:

[TABLE]

[TABLE]

Once again, one can replace $\partial_{t}a_{0,{\beta_{i}}}^{(n-2)}$ using (63),

[TABLE]

Moreover, since $\partial_{\alpha_{n}}p=\rho c_{s}^{2}\partial_{\alpha_{n}}\theta+c_{s}^{2}\theta\partial_{\alpha_{n}}\rho$ we can rearrange several terms as follows,

[TABLE]

By using (64), one gets

[TABLE]

The last time derivative $\partial_{t}\theta$ is removed thanks to (67). Further simplifications, using the implicit summation on $\gamma$ in the second term, lead to

[TABLE]

Eventually, regrouping all terms leads to

[TABLE]

STEP 3: Dispose of $\partial_{\alpha_{n}}a_{0,{\alpha_{1}..\alpha_{n-1}}}^{(n-1)}$ via the definition of the Maxwell-Boltzmann equilibrium state $f^{(0)}$.

By the chain rule,

[TABLE]

and since $f^{(0)}=\frac{\rho}{(2\pi c_{s}^{2}\theta)^{D/2}}\exp\left[-\frac{(\bm{\xi-u})^{2}}{2c_{s}^{2}\theta}\right]$, the following derivatives can be obtained:

[TABLE]

Hence

[TABLE]

and by integration by parts, one can show that

[TABLE]

Then, removing the implicit summation on $\gamma$ leads to

[TABLE]

Injecting this expression in (A) gives

[TABLE]

Finally, knowing that

[TABLE]

leads to (A), which is valid for $n\geq 4$.

❏

Particularity of the isothermal case

The same kind of recursive relation can be obtained in the isothermal case, after noticing that this assumption has two important consequences in the previous derivation:

• •

the temperature is constant ($\theta=1$),

• •

the term $-\dfrac{2}{D}\delta_{\alpha_{i}\alpha_{n}}\partial_{\gamma}u_{\gamma}$ disappears in (71), leading to a non-zero bulk viscosity $\mu_{b}=\dfrac{2}{D}\mu$ Dellar (2001).

Thus equation (71) becomes

[TABLE]

and knowing that

[TABLE]

Malaspinas’ recursive formula Malaspinas (2015) can be recovered:

[TABLE]

Appendix B Recursive formula over equilibrium coefficients

The aim of this section is to prove the following relation

[TABLE]

which simplifies in

[TABLE]

for the isothermal case ($\theta=1$). To this end, the Rodrigues relation over Hermite polynomials will be used:

[TABLE]

Using this relation, one gets, for $n\geq 2$:

[TABLE]

where $\bm{c}=\bm{\xi}-\bm{u}$. By integration by parts, the first term gives

[TABLE]

Finally, using the fact that

[TABLE]

one gets

[TABLE]

so that

[TABLE]

❏

Appendix C Implementation details

This appendix summarizes all the necessary material to properly implement the recursive regularization procedure, in both 2D and 3D, for some common and high-order lattice structures. As a reminder, the purpose of this step is to rebuilt the discrete VDF $f_{i}$ knowing the macroscopic properties of the flow (density $\rho$, velocity $\bm{u}$ and temperature $\theta$) and their gradients:

[TABLE]

Hermite tensors:

The definition of Hermite tensors up to the fourth order (sufficient for the recovery of Navier-Stokes’ equations) is

[TABLE]

where the following mathematical notations $\xi_{i,\alpha\beta}\equiv\xi_{i,\alpha}\xi_{i,\beta}$, $\xi_{i,\alpha\beta\gamma}\equiv\xi_{i,\alpha}\xi_{i,\beta}\xi_{i,\gamma}$, $\xi_{i,\alpha\beta\gamma\delta}\equiv\xi_{i,\alpha}\xi_{i,\beta}\xi_{i,\gamma}\xi_{i,\delta}$ is used for the sake of compacity.

In 2D, $(\alpha,\beta,\gamma,\delta)\in\{x,y\}^{4}$ and thus only two different values can be chosen for each index. This simplifies the definition of Hermite polynomials into:

[TABLE]

In 3D, $(\alpha,\beta,\gamma,\delta)\in\{x,y,z\}^{4}$ and therefore one more value can be chosen for each index. This adds the following Hermite polynomials:

[TABLE]

Due to symmetry properties of Hermite tensors ($\mathcal{H}^{(1)}_{i,y}$, $\mathcal{H}^{(1)}_{i,z}$, $\mathcal{H}^{(2)}_{i,yy}$, $\mathcal{H}^{(2)}_{i,zz}$,$\mathcal{H}^{(2)}_{i,xz}$,$\mathcal{H}^{(2)}_{i,yz}$, $\mathcal{H}^{(3)}_{i,yyy}$, $\mathcal{H}^{(3)}_{i,zzz}$, $\mathcal{H}^{(3)}_{i,xxz}$, $\mathcal{H}^{(3)}_{i,yyx}$, $\mathcal{H}^{(3)}_{i,yyz}$, $\mathcal{H}^{(4)}_{i,yyyy}$, $\mathcal{H}^{(4)}_{i,zzzz}$, $\mathcal{H}^{(4)}_{i,xxxz}$, $\mathcal{H}^{(4)}_{i,yyyx}$, $\mathcal{H}^{(4)}_{i,yyyz}$, $\mathcal{H}^{(4)}_{i,xxzz}$, $\mathcal{H}^{(4)}_{i,yyzz}$, $\mathcal{H}^{(4)}_{i,yyxz}$, $\mathcal{H}^{(4)}_{i,zzxy}$) can be obtained changing $x$ by $y$ or $z$ in the above formulas. By definition, the very same property is also available for $a_{0,\bm{\alpha}}^{(n)}$ and $a_{1,\bm{\alpha}}^{(n)}$.

Hermite coefficients at equilibrium:

The recursive definition of Hermite coefficients at equilibrium is first recalled (74). Up to the fourth order, they read

[TABLE]

This becomes in 2D,

[TABLE]

And the 3D extension adds the following coefficients,

[TABLE]

First-order off-equilbrium Hermite coefficients:

The recursive definition of first-order (with respect to the Knudsen number) off-equilibrium Hermite coefficients (Eq. (II.5)), up to the fourth order (with respect to the VDF moments), reads

[TABLE]

with $S_{\alpha\beta}=\partial_{\alpha}u_{\beta}+\partial_{\beta}u_{\alpha}$, and $\mu=\rho c_{s}^{2}\theta\tau$ is the dynamic viscosity.

For the 2D case, this simplifies into,

[TABLE]

And for the 3D extension, more coefficients need to be taken into account,

[TABLE]

Here attention must be paid to the way $a_{1,\alpha\beta}^{(2)}$ and $a_{1,\alpha\beta\gamma}^{(3)}$ are computed. One can use finite differences to evaluate the velocity and temperature gradients, but it will degrade the accuracy of the algorithm. Instead, we propose the following idea for LBM allowing the preservation of moments of the VDF up to the fourth order ($D2V37$ and $D3Q103$):

Compute $a_{1,\alpha\beta}^{(2)}$ using the definition of the second-order moment of the VDF,

[TABLE]

This is justified by the fact that the error introduced by the approximation $f_{i}=\sum_{n}f_{i}^{(n)}\approx f_{i}^{(0)}+f_{i}^{(1)}$, with $f_{i}^{(n)}~{}\sim\mathcal{O}\left(\epsilon^{n}\right)$, should be small enough when the LBM allows to conserve moments of the VDF up to the fourth order. 2. 2.

Compute $a_{1,\alpha\beta\gamma}^{(3)}$ using the definition of the third-order moment of the VDF,

[TABLE] 3. 3.

Average the contributions of all $a_{1,\alpha\beta\gamma}^{(3)}$ to the temperature gradients,

[TABLE] 4. 4.

Reconstruct $a_{1,\alpha\beta\gamma}^{(3)}$ using Eq. (84) thanks to Eqs. (86) & (89).

Regularized collision operator:

The regularization procedure of pre-collision VDFs can be interpreted as a particular collision step where

[TABLE]

with $A^{(n)}_{\bm{\alpha}}$ being the number of times each Hermite tensors appears in the expansion. As an example, $\mathcal{H}^{(3)}_{i,xxy}$ appears three times in the above development since $\mathcal{H}^{(3)}_{i,xxy}=\mathcal{H}^{(3)}_{i,xyx}=\mathcal{H}^{(3)}_{i,yxx}$. This property is taken into account imposing $A^{(3)}_{xxy}=3$. And more generally,

[TABLE]

It flows from Pascal’s triangle and pyramid rules (also called binomial and trinomial expansions) using $n_{x}$, $n_{y}$ and $n_{z}$ the number of occurrences of $x$, $y$ and $z$ in $\bm{\alpha}=(\alpha_{1},..,\alpha_{n})$ respectively. It should be noted that Eq. (90) is the most general form of the regularized collision step, which allows to decouple the relaxation process specific to each Hermite coefficients. A similar formula was previously given by Shan & Chen Shan and Chen (2007) where all Hermite coefficients belonging to the same expansion order $n$ were related by the same relaxation time, i.e., $\tau_{n,\bm{\alpha}}\longrightarrow\tau_{n}$. The latter formulation allows to preserve isotropy properties of the LBM and is thus preferred.

Appendix D Lattice structures and Related Hermite tensors

Here the link between lattice structures and Hermite tensors is emphasized. To properly choose which Hermite polynomials should be taken into account in the expansion of $f_{i}$’s, the preservation of the orthogonality property of these polynomials, with respect to the weighted scalar product, is considered Philippi et al. (2006). As a reminder, this condition is as follows

[TABLE]

with $c_{s}^{2}=rT_{0}$. Results concerning some standard and high-order lattice structures are reported in Tabs. 1 & 2, for the 2D & 3D cases respectively. The convention from Shan et al. (2006) was used to describe the properties of each lattice structure: $E_{D,Q}^{V}$ where $D$ is the number of physical dimensions, $Q$ is the degree of precision of the quadrature, and $V$ is the number of discrete velocities. Furthermore, all velocities obtained by cyclic permutations and/or reflections with respect to each axis are omitted for the sake of clarity.

To conclude this appendix, expressions of the VDFs for the most complex lattice structures described herein ($E_{2,9}^{37}$ and $E_{3,9}^{103}$) are given by:

[TABLE]

and for the 3D extension,

[TABLE]

In the case of $f^{(1)}_{i}$ the very same terms, as for the equilibrium part $f^{(0)}_{i}$, have to be taken into account.