Three-dimensional central-moments-based lattice Boltzmann method with external forcing: A consistent, concise and universal formulation

TL;DR

This paper introduces a unified, mathematically rigorous formulation of the three-dimensional central-moments-based lattice Boltzmann method with external forcing, enhancing its consistency, simplicity, and applicability to complex physics.

Contribution

It combines recent systematic derivations of Galilean invariant central moments with external forcing into a compact, universal lattice Boltzmann framework.

Findings

Provides a consistent derivation of Galilean invariant forcing terms

Develops a simple, compact algorithm for 3D CM-LBM with external forcing

Enhances the applicability of CM-LBM to high-Reynolds number flows

Abstract

The cascaded or central-moments-based lattice Boltzmann method (CM-LBM) is a robust alternative to the more conventional BGK-LBM for the simulation of high-Reynolds number flows. Unfortunately, its original formulation makes its extension to a broader range of physics quite difficult. To tackle this issue, a recent work [A. De Rosis, Phys. Rev. E 95, 013310 (2017)] proposed a more generic way to derive concise and efficient three-dimensional CM-LBMs. Knowing the original model also relies on central moments that are derived in an adhoc manner, i.e., by mimicking those of the Maxwell-Boltzmann distribution to ensure their Galilean invariance a posteriori, a very recent effort [A. De Rosis and K. H. Luo, Phys. Rev. E 99, 013301 (2019)] was proposed to further generalize their derivation. The latter has shown that one could derive Galilean invariant CMs in a systematic and a priori manner by taking into account high-order Hermite polynomials in the derivation of the discrete equilibrium state. Combining these two approaches, a compact and mathematically sound formulation of the CM-LBM with external forcing is proposed. More specifically, the proposed formalism fully takes advantage of the D3Q27 discretization by relying on the corresponding set of 27 Hermite polynomials (up to the sixth order) for the derivation of both the discrete equilibrium state and the forcing term. The present methodology is more consistent than previous approaches, as it properly explains how to derive Galilean invariant CMs of the forcing term in an a priori manner. Furthermore, while keeping the numerical properties of the original CM-LBM, the present work leads to a compact and simple algorithm, representing a universal methodology based on CMs and external forcing within the lattice Boltzmann framework.

Tables

Table 1. Table 1: Orszag-Tang vortex: findings from different LB schemes and reference spectral values from  Dellar ( 2002 ) in terms of peak value of the current, j max subscript 𝑗 max j_{\mathrm{max}} , and vorticity, ζ max subscript 𝜁 max \zeta_{\mathrm{max}} , at two representative time instants.

	$t$	Dellar (2002)	$M_1$	$M_2$	$M_3$	$M_4$	$M_5$	$M_6$	$M_7$
$j_{\max }$	0.5	18.24	18.26	18.26	18.26	18.26	18.26	18.27	18.26
	1	46.59	46.66	46.66	46.66	46.66	46.66	46.69	46.66
${\zeta }_{\max }$	0.5	6.758	6.756	6.756	6.756	6.756	6.756	6.755	6.756
	1	14.20	14.10	14.10	14.11	14.10	14.10	14.10	14.10

Table 2. Table 2: Rayleigh-Taylor instability: time evolution of the position of the spike and the bubble of the interface at salient time instants. Present results are compared those from (i) a MRT LB study  Saito  et al. ( 2017 ) , (ii) a BGK LB model for multiphase flows  He  et al. ( 1999 ) , (iii) a phase-field MRT LB scheme  Wang  et al. ( 2016 ) , and (iv) a solution of the coupled Navier-Stokes-Cahn-Hilliard equations  Lee and Kim ( 2013 ) .

	Spike					Bubble
$t/t_0$	Present	Saito et al. (2017)	He et al. (1999)	Wang et al. (2016)	Lee and Kim (2013)	Present	Saito et al. (2017)	He et al. (1999)	Wang et al. (2016)	Lee and Kim (2013)
0.0	1.897	1.895	1.887	1.888	1.904	2.100	2.095	2.092	2.096	2.101
0.5	1.845	1.864	1.839	1.860	1.869	2.101	2.131	2.113	2.129	2.131
1.0	1.753	1.763	1.744	1.755	1.776	2.202	2.230	2.229	2.228	2.224
1.5	1.591	1.587	1.555	1.569	1.618	2.345	2.377	2.372	2.364	2.372
2.0	1.378	1.357	1.312	1.325	1.396	2.516	2.535	2.545	2.524	2.535
2.5	1.121	1.085	1.022	1.037	1.149	2.682	2.695	2.693	2.672	2.688
3.0	0.791	0.788	0.712	0.740	0.863	2.834	2.847	2.846	2.824	2.856
3.5	0.537	0.481	0.390	0.419	0.572	2.997	2.999	3.009	2.984	3.009
4.0	0.233	0.160	0.060	0.090	0.271	3.184	3.179	3.178	3.164	3.181

Equations

$$|c_{ix}⟩$$

$$|c_{iy}⟩$$

$$|c_{iz}⟩$$

$$f_i(\mathbf{x}+{\mathbf{c}}_i,t+1)=f_i(\mathbf{x},t)$$

$$+\left(1-\frac{\omega }{2}\right){\mathcal{F}}_i(\mathbf{x},t),$$

$$f_i^{\star }(\mathbf{x},t)=f_i(\mathbf{x},t)$$

$$+\left(1-\frac{\omega }{2}\right){\mathcal{F}}_i(\mathbf{x},t),$$

$$f_i(\mathbf{x}+{\mathbf{c}}_i,t+1)=f_i^{\star }(\mathbf{x},t),$$

$$\rho$$

$$\rho \mathbf{u}$$

$$f_i^{eq}$$

$$\displaystyle\quad+\frac{1}{2c_{s}^{6}}\bigg{[}\mathcal{H}_{ixxy}^{(3)}u_{x}^{2}u_{y}+\mathcal{H}_{ixxz}^{(3)}u_{x}^{2}u_{z}+\mathcal{H}_{ixyy}^{(3)}u_{x}u_{y}^{2}+\mathcal{H}_{ixzz}^{(3)}u_{x}u_{z}^{2}+\mathcal{H}_{iyzz}^{(3)}u_{y}u_{z}^{2}+\mathcal{H}_{iyyz}^{(3)}u_{y}^{2}u_{z}+2\mathcal{H}_{ixyz}^{(3)}u_{x}u_{y}u_{z}\bigg{]}$$

$$\displaystyle\quad+\frac{1}{4c_{s}^{8}}\bigg{[}\mathcal{H}_{ixxyy}^{(4)}u_{x}^{2}u_{y}^{2}+\mathcal{H}_{ixxzz}^{(4)}u_{x}^{2}u_{z}^{2}+\mathcal{H}_{iyyzz}^{(4)}u_{y}^{2}u_{z}^{2}+2\bigg{(}\mathcal{H}_{ixyzz}^{(4)}u_{x}u_{y}u_{z}^{2}+\mathcal{H}_{ixyyz}^{(4)}u_{x}u_{y}^{2}u_{z}+\mathcal{H}_{ixxyz}^{(4)}u_{x}^{2}u_{y}u_{z}\bigg{)}\bigg{]}$$

$$\displaystyle\quad+\frac{1}{4c_{s}^{10}}\bigg{[}\mathcal{H}_{ixxyzz}^{(5)}u_{x}^{2}u_{y}u_{z}^{2}+\mathcal{H}_{ixxyyz}^{(5)}u_{x}^{2}u_{y}^{2}u_{z}+\mathcal{H}_{ixyyzz}^{(5)}u_{x}u_{y}^{2}u_{z}^{2}\bigg{]}$$

$$\displaystyle\quad+\frac{1}{8c_{s}^{12}}\mathcal{H}_{ixxyyzz}^{(6)}u_{x}^{2}u_{y}^{2}u_{z}^{2}\bigg{\}},$$

$$\mathcal{F}$$

$$\displaystyle=-\bm{F}\cdot\nabla_{\bm{c}}\bigg{(}w(\bm{c})\sum_{n=0}^{\infty}\frac{1}{n!}\bm{a}_{eq}^{(n)}\cdot\mathcal{H}^{(n)}\bigg{)}$$

$$\displaystyle=-\bm{F}\cdot\sum_{n=0}^{\infty}\frac{1}{n!}\bm{a}_{eq}^{(n)}\cdot\nabla_{\bm{c}}\bigg{(}w(\bm{c})\mathcal{H}^{(n)}\bigg{)}$$

$$=\omega (\mathbf{c})\sum _{n=1}^{\infty }\frac{1}{n!}\rho \left[\mathbf{F}{\mathbf{u}}^{(n-1)}\right]\cdot {\mathcal{H}}^{(n)}$$

$$[\mathbf{F}\mathbf{u}\mathbf{u}]=\mathbf{F}\mathbf{u}\mathbf{u}+\mathbf{u}\mathbf{F}\mathbf{u}+\mathbf{u}\mathbf{u}\mathbf{F}.$$

$${\mathcal{F}}_i=$$

$$\displaystyle+2\mathcal{H}_{iyz}^{(2)}(u_{y}F_{z}+u_{z}F_{y})\bigg{]}+\frac{1}{2c_{s}^{6}}\bigg{[}\mathcal{H}_{ixxy}^{(3)}(2u_{x}u_{y}F_{x}+u_{x}^{2}F_{y})+\mathcal{H}_{ixxz}^{(3)}(2u_{x}u_{z}F_{x}+u_{x}^{2}F_{z})+\mathcal{H}_{ixyy}^{(3)}(2u_{x}u_{y}F_{y}+u_{y}^{2}F_{x})$$

$$+{\mathcal{H}}_{ixzz}^{(3)}(2u_x{u}_z{F}_z+u_z^2{F}_x)+{\mathcal{H}}_{iyzz}^{(3)}(2u_y{u}_z{F}_z+u_z^2{F}_y)+{\mathcal{H}}_{iyyz}^{(3)}(2u_y{u}_z{F}_y+u_y^2{F}_z)+2{\mathcal{H}}_{ixyz}^{(3)}(u_x{u}_y{F}_z+u_x{F}_y{u}_z$$

$$\displaystyle+F_{x}u_{y}u_{z})\bigg{]}+\frac{1}{4c_{s}^{8}}\bigg{[}\mathcal{H}_{ixxyy}^{(4)}(2u_{x}^{2}u_{y}F_{y}+2u_{x}u_{y}^{2}F_{x})+\mathcal{H}_{ixxzz}^{(4)}(2u_{x}^{2}u_{z}F_{z}+2u_{x}u_{z}^{2}F_{x})+\mathcal{H}_{iyyzz}^{(4)}(2u_{y}^{2}u_{z}F_{z}+2u_{y}u_{z}^{2}F_{y})$$

$$+2{\mathcal{H}}_{ixyzz}^{(4)}(2u_x{u}_y{u}_z{F}_z+u_x{F}_y{u}_z^2+F_x{u}_y{u}_z^2)+2{\mathcal{H}}_{ixyyz}^{(4)}(2u_x{u}_y{u}_z{F}_y+u_x{u}_y^2{F}_z+F_x{u}_y^2{u}_z)$$

$$\displaystyle+2\mathcal{H}_{ixxyz}^{(4)}(2u_{x}u_{y}u_{z}F_{x}+u_{y}u_{x}^{2}F_{z}+u_{x}^{2}u_{z}F_{y})\bigg{]}+\frac{1}{4c_{s}^{10}}\bigg{[}\mathcal{H}_{ixxyyz}^{(5)}(2u_{x}u_{y}^{2}u_{z}F_{x}+2u_{x}^{2}u_{y}u_{z}F_{y}+u_{x}^{2}u_{y}^{2}F_{z})$$

$$\displaystyle+\mathcal{H}_{ixxyzz}^{(5)}(2u_{x}u_{y}u_{z}^{2}F_{x}+u_{x}^{2}u_{z}^{2}F_{y}+2u_{x}^{2}u_{y}u_{z}F_{z})+\mathcal{H}_{ixyyzz}^{(5)}(u_{y}^{2}u_{z}^{2}F_{x}+2u_{x}u_{y}u_{z}^{2}F_{y}+2u_{x}u_{y}^{2}u_{z}F_{z})\bigg{]}$$

$$\displaystyle+\frac{1}{8c_{s}^{12}}\mathcal{H}_{ixxyyzz}^{(6)}(2u_{x}u_{y}^{2}u_{z}^{2}F_{x}+2u_{x}^{2}u_{y}u_{z}^{2}F_{y}+2u_{x}^{2}u_{y}^{2}u_{z}F_{z})\bigg{\}},$$

$$|{\bar{c}}_{ix}⟩$$

$$|{\bar{c}}_{iy}⟩$$

$$|{\bar{c}}_{iz}⟩$$

$$|{\bar{T}}_0⟩$$

$$|{\bar{T}}_1⟩$$

$$|{\bar{T}}_2⟩$$

$$|{\bar{T}}_3⟩$$

$$|{\bar{T}}_4⟩$$

$$|{\bar{T}}_5⟩$$

$$|{\bar{T}}_6⟩$$

$$|{\bar{T}}_7⟩$$

$$|{\bar{T}}_8⟩$$

$$|{\bar{T}}_9⟩$$

$$|{\bar{T}}_{10}⟩$$

$$|{\bar{T}}_{11}⟩$$

$$|{\bar{T}}_{12}⟩$$

$$|{\bar{T}}_{13}⟩$$

$$|{\bar{T}}_{14}⟩$$

$$|{\bar{T}}_{15}⟩$$

$$|{\bar{T}}_{16}⟩$$

$$|{\bar{T}}_{17}⟩$$

Full text

††preprint: APS/123-QED

I Introduction

Nowadays, the lattice Boltzmann method (LBM) is a consolidated approach to simulate fluid flows Benzi et al. (1992); Succi (2001, 2015, 2016); Krüger et al. (2016); Succi (2018). The basic idea is to model the fluid through populations (or distributions) of fictitious particles moving along the links of a fixed Cartesian lattice. The space and time evolution of the distributions is predicted by solving the lattice Boltzmann equation that involves a collide-and-stream process, where the collision stage retains the flow physics. Then, a proper treatment of the collision process is instrumental to obtain accurate predictions of the fluid dynamics.

The impressive popularity of the LBM mainly stems from the intrinsic simplicity of the BGK collision operator Bhatnagar et al. (1954). In short, all the populations are forced to relax (with a common unique rate) towards a discrete equilibrium state derived by applying a Gauss-Hermite quadrature to the continuous Maxwellian distribution Shan et al. (2006). Unfortunately, it is well known to be prone to numerical instability in the low-viscosity regime, thus becoming unsuitable for the prediction of turbulent flows. A possible alternative is represented by the multiple-relaxation-time (MRT) LBM, that suggests to perform the collision in a space of raw (or absolute) moments d’Humières (2002). While second-order ones are related to the flow physics and relax with a frequency directly linked to the fluid kinematic viscosity, higher-order moments are related to non-hydrodynamic modes Latt and Chopard (2006), and their relaxation rates can be considered as free parameters. The MRT formulation allows to directly tune the relaxation frequencies related to these high-order contributions, eventually leading to an improved stability of the resulting LBM. In addition, the original MRT-LBM further increases the dissipation of acoustic waves, hence reducing their impact on the surrounding fluid dynamics Marié et al. (2009). Usually, both BGK and MRT LBMs relax to equilibrium states derived through a second-order truncated Taylor expansion in the local Mach number of the above-mentioned continuous Maxwellian distribution. As a consequence, these are Galilean invariant only up to the second order. Unfortunately, the lack of Galilean invariance at higher orders leads to numerical instability in the low viscous regime Nie et al. (2008); Coreixas (2018).

To cope with this problem, Geier et al. Geier et al. (2006) proposed LBMs based on the concept of central moments (CMs), the latter being obtained by shifting the lattice directions by the local fluid velocity. In the actual practice, it is more convenient to operate in terms of raw moments, thus requiring a transformation between the two types of moments. This is achieved by using binomial formulas, that generate central moments of a certain order as a function of lower order raw moments only. Consequently, the collision in the CM space shows a pyramidal hierarchical topology, where the post-collision state of CMs is constructed starting from the lowest order, and then proceeding in ascending sequence, hence the name “cascaded” lattice Boltzmann method (CLBM). Undoubtedly, the CLBM drastically outperforms both BGK and MRT in terms of stability Geier et al. (2006, 2007); Geier (2008); Asinari (2008); Geier et al. (2009, 2015); De Rosis and Lévêque (2016); Geier et al. (2017a, b); Geier and Pasquali (2018); Fei and Luo (2017); Fei et al. (2018); Shah et al. (2017); Kumar et al. (2017); Sharma et al. (2017); Fei and Luo (2018a, b); Safari et al. (2018); Premnath and Banerjee (2009, 2011); Ning et al. (2016); Hajabdollahi and Premnath (2017). The latter point can partially be explained by the positive hyperviscosity naturally introduced by CM-LBMs, and that can further be adjusted, through fine tuning of relaxation times related to high-order CMs, in order to find the best compromise between stability and accuracy Chávez-Modena et al. (2018).

An important step in the design of the CLBM is the match between continuous and discrete central moments, that enforces the Galilean invariance at all orders. In other words, CMs are forced to relax to the equilibrium state of the continuous Maxwellian distribution, rather than to its discrete counterpart. Consequently, this methodology assumes that equilibrium CMs are unchanged by the velocity discretization of the Boltzmann equation, which might not be true depending on both equilibrium CMs of interest and on the considered lattice of discrete velocities.

More recently, we approached central moments from a different viewpoint De Rosis (2016, 2017a, 2017b, 2017c). Given a certain lattice, our methodology consists of building a transformation matrix allowing us to move from the space of populations to the one of central moments and vice versa. The resultant algorithm loses the above-mentioned pyramidal cascaded structure and, as a consequence, it can be interpreted as a “non-cascaded” way to apply the collision step in CM space De Rosis (2017d). In addition, this methodology is based on CMs computed from the polynomial expansion of populations instead of those derived from the continuous ones. Thus, one does not need anymore to choose which continuous CMs should be kept for the collision step. Consequently, the proposed approach allows the derivation of CM-LBMs for any lattice of discrete velocities in a straightforward manner. This was thoroughly demonstrated by successfully recovering different sets of governing equations with this approach, hence allowing the simulation of a rich variety of physics problems such as shallow waters De Rosis (2017e), magnetohydrodynamic De Rosis et al. (2018), and multicomponent flows Saito et al. (2018) among others.

Nonetheless, it is important to understand that the above attempt to derive a CMs-based scheme in the D3Q27 space was relying on the classical second-order truncated equilibrium state De Rosis (2017a). Even if the stability of the algorithm was greatly improved, it involved a huge computational cost Fei et al. (2018). This is mainly due to the use of an improper (or incomplete) equilibrium state, which eventually leads to a non-negligible number of non-zero velocity-dependent equilibrium CMs. Indeed, as pointed out by Malaspinas Malaspinas (2015) and Coreixas et al. Coreixas et al. (2017); Coreixas (2018), the full potential of the D3Q27 velocity space can only be achieved by employing the correct (or complete) equilibrium populations accounting for the whole set of 27 Hermite polynomials. By including these high-order polynomials, the methodology outlined in De Rosis (2017a) was proven to lead to Galilean invariant CMs that were used for the derivation of a force-free D3Q27-CM-LBM De Rosis and Luo (2019). This formulation was further shown to recover the behavior of the original cascaded LBM by only relying on the correct set of Hermite polynomials. Since this family of polynomials is tightly linked to the velocity discretization of the Boltzmann equation Grad (1949); Shan and He (1998); Shan et al. (2006), the correct set of Hermite polynomials is known in an a priori way, hence ensuring its consistency for any kind of lattices.

When it comes to external or internal forces, it is well known that they are ubiquitous in many fluid systems. Indeed, gravity, Coriolis or Lorentz forces are just few examples of the fields that a fluid might undergo. It is then of paramount importance to properly take into account these forces in the lattice Boltzmann framework. To meet this need, a seminal contribution has been proposed by Guo et al. Guo et al. (2002) in 2002, where the authors analyzed the discrete lattice effects on the forcing term in the BGK-LBM. Later, Guo $\&$ Zheng Guo and Zheng (2008) extended the approach to the MRT-LBM. Concerning the adoption of CMs, Premnath $\&$ Banerjee derived forcing schemes based on the cascaded approach for two-and three-dimensional simulations Premnath and Banerjee (2009, 2011). In particular, they approximated the forcing term according to the formula in the pioneering work by He et al. He et al. (1998). In 2017, this forcing term has been adopted by Fei $\&$ Luo Fei and Luo (2017) within a framework where equilibrium CMs were computed from the continuous Maxwellian distributions. The same force treatment has been adopted in De Rosis (2017b), where equilibrium CMs were derived from the classical second-order-truncated expansion. More recently, it has been demonstrated that the forcing term can be written as a Hermite polynomial expansion up to the fourth order in the D2Q9 space Huang et al. (2018), where the well-known formula by Guo et al. Guo et al. (2002) is recovered if Hermite polynomials of order higher than two are neglected.

In this paper, we derive a three dimensional CMs-based model with external forcing in an a priori manner. Following the idea originally proposed by Shan et al. Shan et al. (2006), and further adopted for the D2Q9 space by Huang et al. Huang et al. (2018), the forcing term is written as an expansion of Hermite polynomials. By adopting the very same set of Hermite polynomials for both the equilibrium state and the forcing term, it is thoroughly shown that corresponding CMs become Galilean invariant, i.e., all remaining error terms depending on the mean flow vanish. As a consequence, the number of non-zero CMs of the forcing term drops to nine, which drastically simplify its implementation.

The rest of the paper is organized as follows. The adopted methodology is first devised in Sec. II, and the impact of high-order Hermite polynomials on CMs of the forcing term is thoroughly quantified. In Sec. III, several numerical experiments of increasing complexity highlight the consistency, accuracy and universality of the proposed model. A comparative study with six other forcing schemes is also provided to further evaluate the numerical properties of the present approach. Some conclusions are finally drawn in Section IV. For the sake of completeness, details regarding Hermite polynomials, the D2Q9-CM-LBM with external forcing, the lattice Boltzmann method for the magnetic field, and the CMs-based color-gradient method are reported in Apps. A, B, C, and  D, respectively.

II Lattice Boltzmann schemes with higher-order Hermite polynomials

Let us consider an Eulerian basis $\bm{x}=[x,y,z]$ and the D3Q27 space Succi (2001). The lattice Boltzmann equation predicts the space and time evolution of the particle distribution functions $\displaystyle|f_{i}\rangle=\left[f_{0},\,\ldots,\,f_{i},\,\ldots,\,f_{26}\right]^{\top}$ that collide and stream along the generic link $i=0\ldots 26$ with the discrete veloctities $\displaystyle\bm{c}_{i}=[|c_{ix}\rangle,\,|c_{iy}\rangle,\,|c_{iz}\rangle]$ defined as

[TABLE]

Let us employ the symbols $|\bullet\rangle$ and $\top$ to denote a column vector and the transpose operator, respectively. Within the BGK approximation Bhatnagar et al. (1954), the numerical discretization of the lattice Boltzmann equation with external forcing reads as follows:

[TABLE]

where the lattice unit system was assumed. As usual, this numerical scheme can be divided into two parts, i.e., collision:

[TABLE]

and streaming:

[TABLE]

where the superscript $\star$ denotes post-collision quantities here and henceforth. If not otherwize stated, the dependence on the space and the time will be tacitly assumed in the rest of the paper. The term $\mathcal{F}_{i}$ accounts for external body forces $\bm{F}=[F_{x},\,F_{y},\,F_{z}]$ and it will be discussed later. The fluid density $\rho$ and velocity $\bm{u}=[u_{x},\,u_{y},\,u_{z}]$ are computed as

[TABLE]

respectively. The second term of the right-hand side of Eq. (3) is the BGK collision operator, that forces all the populations to relax to a discrete local equilibrium $f_{i}^{eq}$ with the same rate $\omega=\left(\nu/c_{s}^{2}+1/2\right)^{-1}$, $\nu$ being the fluid kinematic viscosity. Following the works by Malaspinas Malaspinas (2015) and Coreixas Coreixas et al. (2017); Coreixas (2018), the full potential of the D3Q27 lattice can only be achieved by correctly expanding the equilibrium distribution onto the complete basis composed of 27 Hermite polynomials

[TABLE]

where the weights are $w_{0}=8/27$, $w_{1\ldots 6}=2/27$, $w_{7\ldots 18}=1/54$, $w_{19\ldots 26}=1/216$ and $c_{s}=1/\sqrt{3}$ is the lattice sound speed. $\mathcal{H}_{i}^{(n)}$ denotes the $n$th-order Hermite polynomial tensor. Coefficients before these tensors are $(n_{x}!n_{y}!n_{z}!c_{s}^{2n})^{-1}$ where $n_{x}$, $n_{y}$ and $n_{z}$ are the number of occurrences of $x$, $y$ and $z$ respectively. Full expressions of the Hermite polynomials are compiled in App. A.

Before going any further, it is important to understand that in the present context the sole purpose of the extended equilibrium state (7) is to derive Galilean invariant equilibrium and forcing CMs. Hence, only these CMs are required for the implementation of the D3Q27-CM-LBM with or without external forcing. In addition, one can note that Eq. (7) recovers the classical second-order truncated equilibrium if $\mathcal{H}_{i}^{(3)}$, $\mathcal{H}_{i}^{(4)}$, $\mathcal{H}_{i}^{(5)}$ and $\mathcal{H}_{i}^{(6)}$ are neglected.

Regarding the forcing term, in the continuous Boltzmann equation it is expressed as $\mathcal{F}=\displaystyle-\bm{F}\cdot\bm{\nabla}_{\bm{c}}f$, where $\nabla_{\bm{c}}$ is the derivative with respect to the mesoscopic velocity $\bm{c}$ Shan et al. (2006). Based on the Chapman-Enskog expansion, the Navier-Stokes-Fourier equations can be recovered if the force term is approximated as $\displaystyle-\bm{F}\cdot\bm{\nabla}_{\bm{c}}f^{eq}$. Then, it is possible to write an Hermite expansion of this forcing term based on the expansion of the equilibrium state Shan et al. (2006)

[TABLE]

where both Rodrigues’ formula and $\bm{a}_{eq}^{(n-1)}=\rho\bm{u}^{(n-1)}$ have been used Malaspinas (2015). $\bm{u}^{(n-1)}$ is the velocity tensor of rank $(n-1)$. The square bracket $[\bullet]$ stands for cyclic permutations, e.g., for $n=3$,

[TABLE]

By adopting the extended equilibrium state related to the D3Q27 lattice (7), the expansion of the complete forcing term reads as

[TABLE]

where $1/c_{s}^{2n}$ terms appear due to the use of discrete Hermite polynomials $\mathcal{H}_{i}^{(n)}$. Once again, the full form (9) is only used to derive Galilean invariant CMs of the forcing term. As a consequence, only these CMs are required for the implementation of the D3Q27-CM-LBM with external forcing, while the extended equilibrium state only serves a theoretical purpose in the present context. Regarding the forcing term itself, one can notice that the popular formula by Guo et al. Guo et al. (2002) is recovered if terms proportional to $\mathcal{H}_{i}^{(3)}$, $\mathcal{H}_{i}^{(4)}$, $\mathcal{H}_{i}^{(5)}$ and $\mathcal{H}_{i}^{(6)}$ are neglected.

Now, let us focus on central moments. The pillar to design any CMs-based collision operator is to shift the lattice directions by the local fluid velocity Geier et al. (2006). Therefore, it is possible to define $\displaystyle\bar{\bm{c}}_{i}=[|\bar{c}_{ix}\rangle,\,|\bar{c}_{iy}\rangle,\,|\bar{c}_{iz}\rangle]$, where

[TABLE]

One possibility is to adopt the basis proposed in Ref. De Rosis (2017a), which is $\bar{\mathcal{T}}=\left[\bar{T}_{0},\,\ldots,\,\bar{T}_{i},\,\ldots,\,\bar{T}_{26}\right]$ with

[TABLE]

Let us first recall important results concerning the force-free CM-LBM. Pre-collision, equilibrium and post-collision CMs are defined as

[TABLE]

respectively. The first two quantities are evaluated by applying the transformation matrix $\bar{\mathcal{T}}$ to the corresponding distributions, i.e,

[TABLE]

where $|f_{i}^{eq}\rangle=[f_{0}^{eq},\,\ldots f_{i}^{eq},\,\ldots f_{26}^{eq}]^{\top}$. By adopting a Hermite polynomial expansion of the equilibrium state up to $\mathcal{H}_{i}^{(6)}$ (Eq. (7)), equilibrium CMs read as follows

[TABLE]

with $k_{1\ldots 8}^{eq}=k_{10\ldots 16}^{eq}=k_{19\ldots 25}^{eq}=0$. Notably, only five equilibrium CMs assume non-zero values, and they have the same form as those of the continuous Maxwellian. In other words, Galilean invariant equilibrium CMs are obtained when the full set of Hermite polynomials is considered, as originally demonstrated in Ref. De Rosis and Luo (2019).

The post-collision CMs can be written as

[TABLE]

or

[TABLE]

where $\mathbf{I}$ is the $27\times 27$ unit tensor and

[TABLE]

is the $27\times 27$ relaxation matrix adopted in the present work. The latter allows us to increase the numerical stability through (1) the overdissipation of acoustic waves d’Humières (2002), and (2) the equilibriation of high-order moments Latt and Chopard (2006). The post-collision state of CMs then becomes

[TABLE]

where

[TABLE]

and $\displaystyle k_{1\ldots 3}^{\star}=k_{10\ldots 16}^{\star}=k_{19\ldots 25}^{\star}=0$. Eventually, the post-collision populations are reconstructed as

[TABLE]

with $|f_{i}^{\star}\rangle=[f_{0}^{\star},\,\ldots f_{i}^{\star},\,\ldots f_{26}^{\star}]^{\top}$, before being streamed (see Eq. (4)). One can notice that the computation of only five pre-collision CMs and ten post-collision CMs is required in the present work. Hence, not only the two previous choices (equilibrium (7) and relaxation matrix (17)) lead to a robust and accurate CM-LBM De Rosis and Luo (2019), but they also lead to a very concise formulation. The CPU time required for the computation of the remaining CMs can further be reduced by adopting the fast CM transform Geier et al. (2015).

In presence of an external force, the post-collision CMs can be rewritten as:

[TABLE]

or

[TABLE]

where the CMs of the forcing term are

[TABLE]

Now, let us proceed by assuming the presence of Hermite polynomials of progressively higher order. If only $\mathcal{H}_{i}^{(1)}$ and $\mathcal{H}_{i}^{(2)}$ are considered, Eq. (8) leads to the following CMs for the forcing term:

[TABLE]

with $R_{0}=R_{4\ldots 9}=0$. One can notice that most of these CMs are not Galilean invariant since they contain velocity-dependent terms. By further taking into account third order Hermite polynomials, terms proportional to the square of the velocity vanish, and the corresponding CMs become

[TABLE]

with now $R_{0}=R_{4\ldots 9}=R_{13\ldots 16}=0$. If $\mathcal{H}_{i}^{(4)}$ polynomials are further considered in the definition of the forcing term (8), we obtain

[TABLE]

with $R_{0}=R_{4\ldots 9}=R_{13\ldots 22}=0$. By adding $\mathcal{H}_{i}^{(5)}$, CMs are

[TABLE]

Finally by matching the Hermite polynomial expansion of the forcing term to the one of the extended equilibrium state (7), the following expressions are derived:

[TABLE]

with $R_{0}=R_{4\ldots 9}=R_{13\ldots 22}=R_{26}=0$, and where all the velocity-dependent terms have vanished, eventually leading to a fully Galilean invariant forcing term.

Before going any further, let us point out some interesting properties of the present approach. First, the same expressions of $|\mathcal{R}_{i}\rangle$ in Eqs. (28) can be achieved when CMs are computed from the continuous Maxwellian distribution Fei et al. (2018). This is consistent with a recent paper by De Rosis $\&$ Luo De Rosis and Luo (2019), where it is shown that the CMs of the discrete equilibrium distribution reduce to those of the continuous counterpart when the full set of Hermite polynomials is considered. Hence the proposed approach consistently lead to Galilean invariant forcing terms in accordance with the velocity discretization of interest. This methodology can then be extended to any kind of lattices in a straightforward and a priori manner, i.e., without having to rely on the CMs of the Maxwellian equilibrium state. Second, as an alternative to the present Hermite expansion of the forcing term, the latter could have been expressed by the formula $\displaystyle\bm{F}\cdot(\bm{c}_{i}-\bm{u})f_{i}^{eq}/\left(\rho c_{s}^{2}\right)$ in He et al. He et al. (1998). However, we find that, in this case, the CMs of the force show a dependence on the third power of the velocity. This is due to the fact that the forcing term $\displaystyle-\bm{F}\cdot\bm{\nabla}_{\bm{c}_{i}}f^{eq}$ is only approximated as $\displaystyle-\bm{F}\cdot\frac{\bm{c}_{i}-\bm{u}}{c_{s}^{2}}f^{eq}$ in the pioneering work by He et al. He et al. (1998). This is consistent with the analysis presented in Huang et al. (2018) for the D2Q9 space, and such a problem is avoided in the present work thanks to the extended Hermite polynomial expansion.

Eventually, post-collision CMs accounting for sixth-order Hermite polynomials in both equilibrium populations and forcing term are:

[TABLE]

with $k_{13\ldots 16}^{\star}=k_{19\ldots 22}^{\star}=0$. Once again, only five pre-collision CMs need to be computed in the present work. Nevertheless, the inclusion of the forcing term now requires the computation of nine more post-collision CMs, and this brings the number of non-zero post-collision CMs to nineteen instead of ten for the force-free case (18).

To conclude, a time iteration of the present algorithm can be summarized as follows:

Computation of macroscopic variables (6); 2. 2.

Evaluation of pre-collision CMs $k_{4\ldots 8}$ (II); 3. 3.

Calculation of post-collision CMs (29); 4. 4.

Reconstruction of post-collision populations (20); 5. 5.

Streaming (4).

For the sake of completeness, a script is provided to allow the reader to perform all the symbolic manipulations required for the implementation of the present approach 111See Supplemental Material at [D3Q27_CentralMomentsWithForcingScheme.m] for performing all the symbolic manipulations to easily compute all formulas required for the implementation (the transformation matrix, the post-collision CMs, etc)..

III Numerical tests

In this Section, the accuracy, stability and universality of the present approach are discussed using several numerical testcases. Unless differently specified, the density field is initialized as $\rho(\bm{x},t=0)=\rho_{0}$, with $\rho_{0}=1$, and the fluid is initially at rest, i.e., $\bm{u}(\bm{x},t=0)=0$. Velocity boundary conditions are assigned through the regularized technique Latt et al. (2008). Six test cases of increasing complexity are investigated: (1) four-rolls mill, (2) Hartmann flow, (3) two-, (4) three-dimensional Orszag-Tang vortex, (5) static droplet, and (6) the Rayleigh-Taylor instability.

The first three tests are two-dimensional (i.e, only one point is used for the spatial discretization in the $z$ direction), while the last one is a fully three-dimensional case. In the first scenario, the effect of a force driving the flow is simulated. In the second, third and fourth ones, an electrically conductive fluid is considered and the proposed treatment of the forcing scheme is applied to the Lorentz force. A self-consistent internal force is considered in order to simulate a two-components flow in the fifth scenario. Finally, the flexibility of the present approach to deal with more sophisticated LB models is shown through simulations of the Rayleigh-Taylor instability by means of the color-gradient method.

In all the tests, a diffusive scaling is adopted for the computation of the time step. In addition, all quantities are expressed in lattice units hereafter.

To further evaluate the robustness and the accuracy of the present method (referred to as $\mathrm{M_{1}}$), a comparative study with six other forcing schemes is conducted:

• –

$\mathrm{M_{2}}$: a recent effort by Fei $\&$ Luo Fei and Luo (2017), where equilibrium CMs are derived from the continuous Maxwellian distribution (equivalent to the six-order truncated expansion of the discrete equilibrium state), and CMs of the forcing terms are computed from the approximation by He et al. He et al. (1998) ;

• –

$\mathrm{M_{3}}$: the methodology discussed in De Rosis (2017b), where only a second-order truncated expansion is used for the discrete equilibrium populations, while the forcing scheme is the same as for $\mathrm{M_{2}}$ (i.e., He et al. He et al. (1998));

• –

$\mathrm{M_{4}}$: a model using sixth-order truncated expansion for the discrete equilibrium populations and a second-order force treatment (i.e., Guo et al. Guo et al. (2002));

• –

$\mathrm{M_{5}}$: a model using the second-order truncated expansion for both the discrete equilibrium populations and the forcing term (i.e., Guo et al. Guo et al. (2002));

• –

$\mathrm{M_{6}}$: the cascaded approach proposed by Premnath $\&$ Banerjee Premnath and Banerjee (2009);

• –

$\mathrm{M_{7}}$: the exact difference method by Kupershtokh $\&$ Medvedev Kupershtokh and Medvedev (2006).

Comparisons between methods $\mathrm{M_{1}}$, $\mathrm{M_{4}}$ and $\mathrm{M_{5}}$ are proposed to better understand the contribution of each high-order modification (equilibrium or forcing term) in the context of the Hermite polynomial expansion. $\mathrm{M_{2}}$ and $\mathrm{M_{4}}$, as well as $\mathrm{M_{3}}$ and $\mathrm{M_{5}}$, further highlight the impact of the forcing scheme since they are based on the same equilibrium CMs. On the contrary, comparing $\mathrm{M_{2}}$ and $\mathrm{M_{3}}$ (or $\mathrm{M_{4}}$ and $\mathrm{M_{5}}$) gives more information about the influence of equilibrium CMs. Eventually, both $\mathrm{M_{6}}$ and $\mathrm{M_{7}}$ were also considered for the present comparative study since the former was derived in the context of CMs, whereas the latter is a common forcing scheme used in many applications of the lattice Boltzmann framework.

III.1 Two-dimensional four-rolls mill

Let us consider a square periodic box of size $[x,y]\in[0:N]^{2}$, $N$ being the number of points in each spatial direction. This test simulates the effect of four rotating cylinders that are replaced by a constant body force driving the flow. It reads as

[TABLE]

This test is known as four-rolls mill Malaspinas et al. (2010), and admits a steady state. When the latter is reached, the following velocity field is recovered

[TABLE]

where $u_{0}=10^{-3}$ and $\psi=2\pi/N$. A contour plot of the velocity magnitude is depicted in FIG. 1, exhibiting the pattern that is typical of this kind of flow.

In the present study, a Reynolds number $\displaystyle\mathrm{Re}=u_{0}N/\nu$ of 100 is considered. To assert the accuracy of the present numerical scheme, a convergence analysis is carried out by varying $N$ as $N=\left[8,\,16,\,32,\,64,\,128,\,256\right]$. Let us collect the velocity field from Eq. (31) and the one from our numerical experiments in the vectors $\mathbf{r}_{\mathrm{an}}$ and $\mathbf{r}_{\mathrm{num}}$, respectively. By doing so, it is then possible to compute the $L_{2}$-norm of the relative error $\varepsilon$ between present findings and the analytical solution as

[TABLE]

As shown in FIG. 2, the expected second-order accuracy property of the present method ($\mathrm{M_{1}}$) is successfully recovered, as an excellent convergence rate equal to 1.998 is found.

The approach $\mathrm{M_{6}}$ manifests a poor performance for the coarset mesh grid ($N=8$), while other models show the same accuracy. The latter result can simply be explained by the fact that discrepancies between these models are proportional to powers of the velocity.

By increasing the mean velocity from $u_{0}=10^{-3}$ to $10^{-2}$, these discrepancies can be exerted, as illustrated in FIG. 3.

Generally speaking, the proposed approach shows a desirable convergence rate $\sim 2$ for each value of $N$, while the other methodologies require finer grid meshes ($N\geq 64$) to recover the expected behavior. By looking at the results more closely, it is clear that $\mathrm{M_{6}}$ generates large errors for the lowest values of $N$. The exact difference method $\mathrm{M_{7}}$ deviates from the ideal second-order of convergence when $N=16$. Non-negligible errors are experienced for $N=32$ by adopting $\mathrm{M_{2}}-\mathrm{M_{5}}$. Surprisingly, $\mathrm{M_{3}}$ and $\mathrm{M_{5}}$ show the lowest mismatch while they are only second-order approaches regarding both the equilibrium state and the forcing scheme. In fact, this can be explained by the fact that $\mathrm{M_{2}}$ and $\mathrm{M_{4}}$ are inconsistent models since they do not share the same order for both the equilibrium and the forcing terms. Consequently, one can easily show that high-order CMs of these terms contain velocity-dependent terms that violate the principle of Galilean invariance.

III.2 Two-dimensional Hartmann flow

Let us now consider the flow of an electrically conductive fluid of magnetic resistivity $\eta=\nu$, and which is subjected to a magnetic field $\bm{b}$. The interest reader shall refer to App. C for further details about the computation of the magnetic field $\bm{b}$ in the lattice Boltzmann framework. Let us assume a rectangular domain of length $L_{x}$ and height $L_{y}$. Initial conditions consist of

[TABLE]

We reproduce the so-called Hartmann flow, that it is analogous to the Poiseuille flow, but in the former case: (1) the fluid is assumed to be conductive, and (2) a constant uniform vertical magnetic field ($b_{y0}$) is enforced at the bottom and top walls, where the no-slip condition is enforced too. At the other sides, the magnetic field is assumed to be periodic. Magnetic boundary conditions are assigned according to Dellar (2013). The analytical solution of the Hartmann flow reads as follows:

[TABLE]

where $y^{\prime}=2y-L_{y}$ and the Hartmann number is defined as $\mathrm{Ha}=b_{y0}L_{y}/\sqrt{4\rho_{0}\nu\eta}$ Dellar (2002). Moreover, Eq. (33) is prescribed at the inlet section, while an outflow boundary condition is applied at the opposite side of the domain as $\bm{n}^{\ddagger}\cdot\nabla\bm{u}=0$, $\bm{n}^{\ddagger}$ being the unit vector normal to the boundary. The presence of the magnetic field gives rise to the Lorentz force $\bm{F}=\bm{j}\times\bm{b}$, $\bm{j}$ being the electric current that is computed directly from the populations Pattison et al. (2008).

In FIG. 4, results obtained with the present approach are compared to the analytical solution (33) for four values of the Hartmann number, i.e., $\mathrm{Ha}=1,\,3,\,10,\,20$. From them, it is clear that $\mathrm{M_{1}}$ is able to correctly recover the physics of the Hartmann flow when the grid mesh is fine enough. The impact of the grid mesh, as well as, the convergence order of the present method is further studied by evaluating the relative error (32) for different values of $L_{y}\in[5:1025]$, and using the analytical solution (33). In FIG. 5, the relative error is plotted against the number of points required to discretize the vertical dimension. Interestingly, a poor convergence is experienced for low values of $L_{y}$, especially for high values of $\mathrm{Ha}$. This behavior can be explained by the presence of progressively thinner Hartmann layers, the latter requiring a larger number of grid points to be successfully and accurately captured.

This test is further considered to compare the behavior of the proposed approach with other forcing schemes. Before moving to the results, one should note that the Hartmann flow is a purely one-dimensional problem (i.e., $u_{y}=u_{z}=0$). Hence, only terms properly to $u_{x}$ remain in Eqs. (7, 9). As a consequence, it is expected that $\mathrm{M}_{1}$ and $\mathrm{M}_{4}$ reduce to $\mathrm{M}_{5}$, whereas $\mathrm{M}_{2}$ should recover the behavior of $\mathrm{M}_{3}$. Indeed, for any combination of $\mathrm{Ha}$ and $L_{y}$, it has been confirmed that high-order approaches lead exactly to the same results, in terms of relative errors and velocity profiles, as their low-order versions. In the following, $\mathrm{M_{1}}$, $\mathrm{M_{3}}$, $\mathrm{M_{6}}$ and $\mathrm{M_{7}}$ are then compared against each other for the most difficult configuration ($\mathrm{Ha}=20$). In addition, the maximal number of points in the vertical direction is increased ($L_{y}^{\mathrm{max}}=4097$) in order to exert discrepancies between all the models.

Results compiled in FIG. 6 show that the exact difference method exhibits the best performance by keeping the second-order accuracy even for the finest grid mesh. On the contrary, a deterioration of the convergence rate is experienced by all the other methods. More precisely, $\mathrm{M}_{1}$ and $\mathrm{M}_{3}$ show non-negligible errors for $L_{y}\geq 1025$. Eventually, the cascaded model $\mathrm{M}_{6}$ starts deviating from the expected second-order convergence for $L_{y}\geq 129$, which confirms the poor behavior of this forcing scheme, as already pointed out in the four-mill test case.

III.3 Two-dimensional Orszag-Tang vortex

The next test case consists in the simulation of the Orszag-Tang vortex problem Dellar (2002); Orszag and Tang (1979). Its simulation represents a hard task since smaller and smaller structures appear in the domain as the Reynolds number grows. This test case is then a perfect candidate to assess (1) the accuracy of the present method to capture the existence of fine flow features, as well as (2) its ability to deal with high local gradients that arise during the simulation De Rosis et al. (2018).

In what follows, a square periodic box of length $2\pi$ is considered with a Reynolds number of $\mathrm{Re}=1600$, and initial conditions are set to

[TABLE]

In FIG. 7, the mean kinetic energy $\langle u^{2}\rangle$ (normalized by its initial value $\langle u_{0}^{2}\rangle$) is reported for different values of the Mach number, i.e., $\mathrm{Ma}=u_{0}/c_{s}=[0.07,\,0.14,\,0.28]$ with $c_{s}=1/\sqrt{3}$. Moreover, it is set $b_{0}=u_{0}$. Experiments are carried out by adopting Hermite polynomials up to the second and sixth orders for both equilibrium populations and forcing term, corresponding to the models $M_{1}$ and $M_{5}$, respectively. Here, the stability properties of the proposed algorithm clearly manifest. In fact, the two strategies are both stable for $\mathrm{Ma}=0.07$ and 0.14. By further increasing the Mach number, $M_{5}$ becomes unstable whereas the present approach, $M_{1}$, leads to the correct behavior.

We also use this scenario to evaluate the accuracy of different LB approaches. In TABLE 1, the peak values of the current $j_{\mathrm{max}}$ and vorticity $\zeta_{\mathrm{max}}$ are compared to reference high-resolution spectral values in Dellar (2002). All the models show very good accuracy, with the relative errors lower than $0.22\%$. To further evaluate their robustness, the Mach number is increased up to $\mathrm{Ma}=0.28$ for all the involved forcing schemes. From this, it is found that only $\mathrm{M_{5}}$ manifests stability issues.

Now that the good numerical properties of our approach have been rigorously demonstrated in the 2D case, its validation will further be carried out considering 3D flows. Due to the non-negligible CPU time required to run these simulations, $\mathrm{M_{1}}$ will only be compared to the second-order formulations $\mathrm{M_{4}}$ and $\mathrm{M_{5}}$, with the exception of the Rayleigh-Taylor instability test case, where data are already available in the literature.

III.4 Three-dimensional Orszag-Tang vortex

This test involves an Orszag-Tang vortex that develops in a three-dimensional cubic periodic box of length $2\pi$ with the initial conditions defined according to Mininni et al. (2006), i.e.,

[TABLE]

with $u_{0}=b_{0}=0.0203$. Similarly to the two-dimensional case, here the flow is characterized by the presence of intermittent small-scale structures. However, the topology of the reconnection of the magnetic field is considerably more complicated in three dimensions Greene (1988). Firstly, the accuracy of the method is assessed by running a simulation at $\mathrm{Re=100}$ with three different grid sizes, i.e., $N=\left[32,\,64,\,128\right]$.

In FIG. 8, the time evolution of the peak value of the current magnitude ($j_{\mathrm{max}}$) is depicted for the different grid dimensions. As $N$ increases, the numerical solution gets progressively closer to the reference one obtained from a high-resolution pseudo-spectral simulation De Rosis (2017a). This is further quantified by defining the vectors $\mathbf{r}_{\mathrm{ref}}$ and $\mathbf{r}_{\mathrm{num}}$, gathering findings from the reference pseudo-spectral analysis and our results, respectively. Then, we introduce $\varepsilon_{\infty}$ as the $L_{\infty}$-norm of the relative discrepancy, i.e.,

[TABLE]

that is of order $10^{-2}$ for the finest grid, thus highlighting the accuracy of the method. Then, in order to evaluate the stability of the algorithm, higher values of the Reynolds number are considered, i.e., $\mathrm{Re}=[570,\,1040,\,3040,\,5600]$. In FIG. (9), the time evolution of $j_{\mathrm{max}}$ is reported for all these Reynolds numbers. From this, it is clear that the present approach ($\mathrm{M_{1}}$) remains stable even for the highest value of the Reynolds number. To further highlight the robustness of $\mathrm{M_{1}}$, the case $\mathrm{Re}=570$ (the one that is expected to be the least prone to the onset of instability) is also simulated using this time the BGK collision operator Pattison et al. (2008). One can notice that, for the latter collision model, the run blows-up at $t\sim 0.43$. Conversely, the adoption of CMs (either the sixth- or the second-order one, i.e. $\mathrm{M_{1}}$, $\mathrm{M_{4}}$ and $\mathrm{M_{5}}$) overcome this problem, allowing us to simulate a larger time span for the whole desired set of $\mathrm{Re}$. This is consistent with the observations in De Rosis et al. (2018) regarding two-dimensional magnetohydrodynamic flows simulated by a CMs-based D2Q9 model. Nonetheless, by increasing the Mach number, $\mathrm{M_{5}}$ becomes unstable as it was already the case in Sec. III.3. It should be noted that this paper delivers the first tangible solution of the three-dimensional high-Reynolds Orszag-Tang vortex within the framework of the lattice Boltzmann method.

Finally, a contour plot of the current magnitude is sketched in FIG. (10) at different time instants for the configuration $\mathrm{Re}=570$ with $N=128$. One can appreciate that the present method is able to capture the existence of fine flow features, as well as the presence of strong gradients.

III.5 Static droplet

In the following, the behavior of the present model is further investigated in the context of a multiphase flow simulated through the popular Shan-Chen model Shan and Chen (1993). Specifically, the latter introduces an interaction force

[TABLE]

that mimics the molecular interactions leading to phase segregation. $\psi=1-e^{-\rho}$ is the so-called pseudo-potential and $G$ is a parameter controlling the strength of the interaction.

Let us consider a periodic domain consisting of $200$ lattice points in each direction. A droplet of radius $R$ and density equal to $\rho_{0}$ is placed in the center of the domain, while the density is set to $\rho_{0}/8$ elsewhere. The parameter $G$ is set equal to -6. The pseudo-potential LBM is known to be prone to numerical instability due to the presence of spurious velocity currents. In FIG. 11, the time evolution of the mean kinetic energy is reported for numerical experiments carried out by using the models $M_{1}$, $M_{4}$ and $M_{5}$. Here, the beneficial effect of using sixth-order Hermite polynomials for the forcing term are clearly visible. Indeed, the run corresponding to the poorest model (i.e., $M_{5}$) blows-up after $\sim 120$ time iterations. Conversely, the adoption of sixth-order Hermite polynomials allows us to simulate a considerably larger time span, with the kinetic energy that tends to vanish after initial large-amplitude oscillations. This behavior is slightly impacted by the form of the forcing term, as both $M_{1}$ and $M_{4}$ lead to stable simulations. A small mismatch between the two forcing schemes is found in the early stage of the simulation, where larger velocity magnitude are encountered. They emphasize the impact of higher-order Hermite polynomials on the forcing term, and suggest a slightly better behavior of $M_{1}$, as compared to $M_{4}$, since the former leads to a lower maximal kinetic energy.

III.6 Rayleigh-Taylor instability

To conclude the numerical investigation of the present method, the simulation of the Rayleigh-Taylor instability is carried out. The fundamental mechanism, upon which this kind of instability relies, can be summarized as follows. Let us consider a heavy fluid on top of a lighter one, both being immiscible fluids that are separated by an interface. Assuming that they are subject to a gravitational field, their interface will not be able to stay in its equilibrium state. As a consequence, plumes (spikes) will start forming and flowing downward (upward).

This Rayleigh-Taylor instability mechanism is simulated hereafter by adopting the color-gradient method Latva-Kokko and Rothman (2005); Reis and Phillips (2007); Leclaire et al. (2012). Recently, Saito et al. Saito et al. (2018) have proposed a central-moments-based color-gradient model by using third- and second-order terms for the equilibrium and forcing, respectively. Here, the approach is recasted within the present framework, hence adopting sixth-order Hermite polynomials for both the equilibrium and the forcing term. This problem, as well as its solution procedure by the color-gradient method, clearly shows the universality of the present approach. The CM-based color-gradient method is briefly outlined in App. D.

Let us consider a three-dimensional domain of size $W\times 4W\times W$, with $W=64$. A fluid of density $\rho_{h}$ is placed over a lighter one of density $\rho_{l}=1$. The initial flow field is perturbed as

[TABLE]

No-slip walls are enforced at the top and bottom sections, while periodic boundaries are assumed at the other sides of the domain. A gravitational body force is considered as

[TABLE]

with $\bm{g}=(0,\,-g,\,0)$, and $g$ chosen so that $\sqrt{gW}=0.04$ Latva-Kokko and Rothman (2005). The problem is governed by two dimensionless parameters: the Reynolds number

[TABLE]

and the Atwood number

[TABLE]

that are set here to $\mathrm{Re}=512$ and $\mathrm{At}=0.5$, respectively. In addition, the characteristic time used for the time evolution of the Rayleigh-Taylor instability is defined as

[TABLE]

In the following, let us denote the bottom and top points of the interface as spike and bubble, respectively. In FIG. 12, the evolution of the interface between the two fluids is sketched at salient time instants. They show that the edge of the spike rolls up at $t/t_{0}=3$. This is consistent with another LB study Saito et al. (2017), where the Rayleigh-Taylor instability has been investigated by means of the multiple-relaxation-time kernel.

More quantitative results, regarding the time evolution of the position of the two reference points, are reported in TABLE 2, and depicted in FIG. 13. Present findings are further compared to several models to assess its accuracy: (1) a MRT LB study Saito et al. (2017), (2) a BGK LB model for multiphase flows He et al. (1999), (3) a phase-field MRT LB scheme Wang et al. (2016) and (4) a solution of the coupled Navier-Stokes-Cahn-Hilliard equations Lee and Kim (2013). From this, it is clear that the present method shows a pretty good agreement with data from the literature. This confirms the good numerical properties of the proposed approach, as well as, its universality.

IV Conclusions

In the present work, a systematic way to derive CM-LBMs with external forcing has been proposed. It naturally flows from a previous work based on the D3Q27 velocity discretization De Rosis and Luo (2019), where Galilean invariant CMs were derived from the discrete equilibrium state instead of the continuous Maxwell Boltzmann distribution. Here, by relying on the very same set of 27 Hermite polynomials, Galilean invariant CMs of the forcing term are also obtained but without any assumption on these CMs. By further equilibrating acoustic and high-order CMs, the present method only requires the computation of five pre-collision CMs. Both points eventually lead to a very concise algorithm.

Interestingly, it is also shown that the present method recovers the behavior of the cascaded LBM Fei et al. (2018). But differently from the latter, the proposed methodology allows the derivation of forcing terms in an a priori way thanks to the tight link existing between the Hermite polynomial expansion and the lattice of discrete velocities. As a consequence, its extension to any kind of velocity discretizations can be done in a systematic and straightforward manner.

The simulation of the four-rolls mill, the Hartmann flow, the two- and three-dimensional Orszag-Tang vortex, a static droplet and the Rayleigh-Taylor instability were conducted to evaluate both the numerical behavior as well as the universality of the present method. By further comparing the latter with six other forcing schemes on several of these testcases, its good accuracy and robustness properties were properly demonstrated.

All in all, the consistency with the cascaded LBM, the compactness of the algorithm as well as its excellent numerical properties lead the present scheme to be a very good candidate to perform multiphysics simulations with (or without) the presence of external forces. Regarding future works, the extension to both high-order lattices and to other kinds of collision space (raw moment, Hermite moment, central Hermite moment and cumulant Coreixas et al. (2019)) is currently under investigation. Corresponding results shall be presented in the near future.

Acknowledgments

This article is based upon work from COST Action MP1305, supported by COST (European Cooperation in Science and Technology). The authors would like to express their gratitude to the reviewers for their relevant remarks and insightful questions. A.D.R. is grateful to Dr. Shimpei Saito for very valuable hints and discussions related to the setup of the simulations of the Rayleigh-Taylor instability.

Appendix A Hermite polynomials

Discrete Hermite polynomials are defined as Malaspinas (2015); Coreixas et al. (2017)

[TABLE]

with the corresponding Gaussian weight function

[TABLE]

where $D$ is the number of physical dimensions and $c_{s}$ is the lattice constant. In the 3D case, discrete Hermite tensors can be computed using tensor products of their 1D formulation,

[TABLE]

where $w^{\prime}$ is the one-dimensional version of the Gaussian weight (43), and $(nx,ny,nz)$ are the number of occurrences of indexes $(x,y,z)$. Since the D3Q27 lattice is a tensor product of three D1Q3 lattices, then only 1D Hermite polynomials of degree up to $n=2$ can be used for the construction of their 3D counterparts. This means that the degree of each Hermite polynomial will be at most two per direction. Indeed,

[TABLE]

Appendix B Two-dimensional model

Let us consider an Eulerian basis $\bm{x}=[x,y]$. Here, we derive the two-dimensional LB model by means of the D2Q9 velocity discretization, where lattice directions $\displaystyle\bm{c}_{i}=[|c_{ix}\rangle,\,|c_{iy}\rangle]$ are

[TABLE]

with $i=0\dots 8$. Let us define the particle distribution functions $\displaystyle|f_{i}\rangle=[f_{0},\,\ldots f_{i},\,\ldots f_{8}]^{\top}$ and the velocity vector $\bm{u}=[u_{x},u_{y}]$. In two dimensions, the equilibrium distributions $|f_{i}^{eq}\rangle=[f_{0}^{eq},\,\ldots f_{i}^{eq},\,\ldots f_{8}^{eq}]^{\top}$ can be expanded onto a basis of Hermite polynomials $\mathcal{H}^{(n)}$ up to the fourth order Malaspinas (2015); Coreixas et al. (2017); Coreixas (2018), i.e

[TABLE]

where $w_{0}=4/9$, $w_{1\ldots 4}=1/9$, $w_{5\ldots 8}=1/36$ and $c_{s}=1/\sqrt{3}$ is the lattice sound speed  Succi (2001). Again, the development of a central-moments-based collision operator begins by shifting the lattice directions by the local fluid velocity, i.e., $\displaystyle\bar{\bm{c}}_{i}=[|\bar{c}_{ix}\rangle,\,|\bar{c}_{iy}\rangle]$, where

[TABLE]

Let us choose the basis $\bar{\mathcal{T}}=\left[\bar{T}_{0},\,\ldots,\,\bar{T}_{i},\,\ldots,\,\bar{T}_{8}\right]$ which components are

[TABLE]

Let us collect pre-collision, equilibrium and post-collision CMs as

[TABLE]

respectively. The first two quantities are evaluated by applying the transformation matrix $\bar{\mathcal{T}}$ to the corresponding distribution, that is

[TABLE]

By adopting equilibrium populations with the full set of Hermite polynomials (see Eq. (46)), equilibrium CMs are

[TABLE]

and $k_{1,2,4,5,6,7}^{eq}=0$. In analogy with the 3D case, it should be noted that (1) the equilibrium state is fully Galilean invariant since it does not show any dependence on the fluid velocity, and (2) the discrete equilibrium CMs have the same form of the continuous counterpart when the full set of Hermite polynomials is considered.

Again, the post-collision CMs can be written as

[TABLE]

or

[TABLE]

where $\mathbf{I}$ is the ($9\times 9$) unit tensor, $\displaystyle\bm{\Lambda}=\mathrm{diag}\left[1,\,1,\,1,\,1,\,\omega,\,\omega,\,1,\,1,\,1\right]$ is a ($9\times 9$) relaxation matrix.

Now, let us define the two-dimensional forcing term $\mathcal{F}_{i}$. Specifically, we employ the expression adopted by Huang et al. Huang et al. (2018):

[TABLE]

where the square bracket in Hermite coefficient denotes permutations. Notice that the popular formula by Guo et al. Guo et al. (2002) is recovered if $\mathcal{H}_{i}^{(3)}$ and $\mathcal{H}_{i}^{(4)}$ are neglected. The CMs of the forcing term can be computed as

[TABLE]

and read as follows:

[TABLE]

and $R_{0,3,4,5,8}=0$. The resultant expressions of the post-collision CMs are:

[TABLE]

Again, one can appreciate that the post-collision state in terms of CMs show simple expressions. Then, the post-collision populations are reconstructed as usual:

[TABLE]

with $|f_{i}^{\star}\rangle=[f_{0}^{\star},\,\ldots f_{i}^{\star},\,\ldots f_{8}^{\star}]^{\top}$, that are eventually streamed.

Appendix C Lattice Boltzmann method for the evolution of the magnetic field

The evolution of the magnetic field $\bm{b}$ is predicted by a second LBM based on the D3Q7 velocity space, that is $\forall l\in\llbracket 0,6\rrbracket$,

[TABLE]

with $\omega_{m}$ the relaxation frequency that is related to the magnetic resistivity as $\displaystyle\eta=c_{s}^{2}\left(\frac{1}{\omega_{m}}-\frac{1}{2}\right)$, where $c_{s}=1/2$ for the D3Q7 lattice. For this particular lattice, discrete velocities $\displaystyle\bm{\xi}_{l}=[|\xi_{lx}\rangle,\,|\xi_{ly}\rangle,\,|\xi_{lz}\rangle]$ are defined as

[TABLE]

Vector-valued populations $\bm{h}_{l}$ are necessary due to the anti-symmetry of the electric tensor Dellar (2002). Their equilibrium state is defined as Pattison et al. (2008)

[TABLE]

where $w_{0}=1/4$, $w_{1\ldots 6}=1/8$ and $(\alpha,\beta)\in\{x,y,z\}^{2}$. Regarding magnetic boundary conditions, we follow the approach in Dellar (2013) that is built on the fact that the magnetic field can be computed as

[TABLE]

Specifically, let us assume that one wants to enforce a certain magnetic field $\bm{b}^{\dagger}=[b_{x}^{\dagger},\,b_{y}^{\dagger},\,b_{z}^{\dagger}]$ to the bottom section of the considered domain. In this case, the only population to be assigned is $\bm{h}_{2}$, that is evaluated as

[TABLE]

Appendix D Central-moments-based color-gradient method

Let us consider two immiscible fluids, namely, red and blue. The evolution of populations $f_{i}^{k}$ is

[TABLE]

where $k=r$ for the red fluid, and $k=b$ for the blue one. Moreover, it is possible to define the total distribution functions as $\displaystyle f_{i}=f_{i}^{r}+f_{i}^{b}$. The collision process $\Omega_{i}^{k}$ is composed of three sub-stages:

[TABLE]

where $\left(\Omega_{i}^{k}\right)^{(1)}$, $\left(\Omega_{i}^{k}\right)^{(2)}$ and $\left(\Omega_{i}^{k}\right)^{(3)}$ are the single-phase, perturbation and recoloring operators, respectively. Within the CM-based framework Saito et al. (2018), the single-phase collision operator can be written as

[TABLE]

with $|\mathcal{F}_{i}\rangle$, $\bar{\mathcal{T}}^{\top}$ and $\Lambda$ that obey Eqs. (9), (II) and (17), respectively. Macroscopic variables are given by

[TABLE]

where $\rho_{k}$ is the density of the fluid $k$, $\rho$ is the total mass density, $\bm{u}$ is the total momentum and $\bm{F}$ is a body force. Distributions relax to an enhanced local equilibrium Leclaire et al. (2013); Saito et al. (2017) that, by adopting sixth-order Hermite polynomials, read as follows:

[TABLE]

Here, quantities $\varphi_{i}$ and $\Phi_{i}$ are

[TABLE]

and

[TABLE]

where $\otimes$ is the tensor (outer) product, and : stands for the contraction of index (Frobenius inner product). The tensor $\mathbf{G}$ is defined as

[TABLE]

where gradients are computed by a second-order isotropic central scheme Liu et al. (2012). $\bar{\nu}$ is the kinematic viscosity interpolated by Leclaire et al. (2012)

[TABLE]

To distinguish the two components, the order parameter $\phi$ is introduced, that is

[TABLE]

The values $\phi=1,-1$, and [math] correspond to a purely red fluid, a purely blue fluid, and the interface, respectively. To obtain a stable interface, the density ratio between the fluids must be taken into account as follows to obtain a stable interface Grunau et al. (1993):

[TABLE]

where the superscript “0” indicates the initial value of the density at the beginning of the simulation Leclaire et al. (2017). The pressure of the fluid is given as an isothermal equation of state:

[TABLE]

for the D3Q27 lattice, where $c_{s}^{k}$ is the speed of sound of the fluid $k$ Leclaire et al. (2012), $\bar{\alpha}$ is interpolated by

[TABLE]

with $\alpha_{b}=8/27$ and $c_{s}^{b}=1/\sqrt{3}$ Saito et al. (2018).

Following Brackbill et al. (1992); Reis and Phillips (2007); Liu et al. (2012), the interfacial tension is modeled by the so-called perturbation operator:

[TABLE]

where

[TABLE]

Eventually, the following recoloring operator is applied to promote phase segregation and maintain the interface:

[TABLE]

where $\gamma=0.7$ Leclaire et al. (2012); Saito et al. (2017, 2018), $\theta_{i}$ is the angle between $\nabla\phi$ and $\bm{c}_{i}$, which is defined by

[TABLE]