Recursive Regularized-BGK (RR-BGK)

Note

The following development is more involved and may require familiarity with some specific topics to be fully understood. The Hermite-polynomial machinery it uses was motivated on the Lattice Boltzmann Equation page.

The idea of regularization is to keep, in the non-equilibrium populations, only the moments that carry real physics and discard the ghost content. Hermite polynomials are the natural tool: each Hermite order corresponds to one moment order, so truncating the Hermite expansion at order \(N\) keeps exactly the moments up to that order and nothing else. The “recursive” in RR-BGK refers to the fact that the higher-order coefficients can be built from the lower ones by a recurrence, which makes the third-order reconstruction cheap. This is the operator introduced by Malaspinas[1] and Mattila et al.[2], which Nassu uses by default for LES.

The functional form of the equilibrium distribution function \(f_i^\mathrm{eq}\) is taken as a mesoscopic velocity expansion in Hermite polynomials up to order \(N\), as suggested by Malaspinas[1] and Mattila et al.[2]:

(1)\[ f_{i}^{\mathrm{eq}} = w_{i}\sum_{n=0}^{N}\frac{1}{n!}\mathbf{H}_{i}^{(n)}:\mathbf{a}^{(n),\mathrm{eq}} \]

in which “:” stands for a full index contraction, and \(\mathbf{H}_{i}^{(n)}\) is the three-dimensional Hermite polynomial associated to the mesoscopic velocity vector \(\xi_{i,\alpha}=c_{i,\alpha}/c_{s}\). The Hermite polynomials up to order \(N\) are given by:

(2)\[ H^{(0)}\left(\mathbf{\xi}_{i}\right) = 1, \]

(3)\[ H_{i,\alpha}^{(1)}\left(\mathbf{\xi}_{i}\right) = \xi_{i,\alpha}, \]

(4)\[ H_{i,\alpha_{1}\alpha_{2}}^{(2)}\left(\mathbf{\xi}_{i}\right) = \xi_{i,\alpha_{1}}\xi_{i,\alpha_{2}}-\delta_{\alpha_{1}\alpha_{2}}, \]

(5)\[ H_{i,\alpha_{1}\alpha_{2}\alpha_{3}}^{(3)}\left(\mathbf{\xi}_{i}\right) = \xi_{i,\alpha_{1}}\xi_{i,\alpha_{2}}\xi_{i,\alpha_{3}} -\left(\delta_{\alpha_{1}\alpha_{2}}\xi_{i,\alpha_{3}} + \delta_{\alpha_{1}\alpha_{3}}\xi_{i,\alpha_{2}} + \delta_{\alpha_{2}\alpha_{3}}\xi_{i,\alpha_{1}}\right), \]

(6)\[ H_{i,\alpha_{0}...\alpha_{n}}^{(n+1)}\left(\mathbf{ \xi_{i}}\right)=\xi_{i,\alpha_{0}}H_{i,\alpha_{1}...\alpha_{n}}^{(n)}\left(\mathbf{ \xi_{i}}\right)-\sum_{i=1}^{n}\delta_{\alpha_{0}\alpha_{i}}H_{i,\alpha_{1}...\alpha_{i-1}\alpha_{i+1}...\alpha_{n}}^{(n-1)}\left(\mathbf{\xi}_{i}\right) \]

Due to orthogonality of Hermite polynomials, the value of the coefficients \(\mathbf{a}^{(n),\mathrm{eq}}\) can be found with:

(7)\[ \mathbf{a}^{(n),\mathrm{eq}}=\sum_{i=0}^{q-1}\mathbf{H}_{i}^{(n)}f_{i}^{\mathrm{eq}}, \]

where \(q\) is the number of velocity directions of the velocity set adopted. Knowing that the moments of equilibrium populations are:

(8)\[ \sum_{i=0}^{q-1}f_{i}^{\mathrm{eq}}=\rho, \]

(9)\[ \sum_{i=0}^{q-1}f_{i}^{\mathrm{eq}}\xi_{i,\alpha}=\frac{1}{c_{s}}\rho u_{\alpha}, \]

(10)\[ \sum_{i=0}^{q-1}f_{i}^{\mathrm{eq}}\xi_{i,\alpha}\xi_{i,\beta}=\frac{1}{c_{s}^{2}}\rho u_{\alpha}u_{\beta} + \rho \delta_{\alpha\beta}. \]

The values of \(\mathbf{a}^{(n),\mathrm{eq}}\) can be written recursively as:

(11)\[ a^{(0),\mathrm{eq}}=\rho, \]

(12)\[ a^{(n),\mathrm{eq}}_{\alpha_{1}...\alpha_{n}}=\frac{1}{c_{s}}u_{\alpha_{n}}a^{(n-1),\mathrm{eq}}_{\alpha_{1}...\alpha_{n-1}}. \]

We adopt an expansion up to third order to represent the equilibrium populations, which gives:

(13)\[ f_{i}^{\mathrm{eq}}=w_{i}\left[\rho + \frac{\mathbf{c}_{i}\cdot \left(\rho\mathbf{u}\right)}{c_{s}^{2}} + \frac{1}{2}\mathbf{H}_{i}^{(2)}:\mathbf{a}^{(2),\mathrm{eq}} + \frac{1}{6}\mathbf{H}_{i}^{(3)}:\mathbf{a}^{(3),\mathrm{eq}}\right] \]

or:

(14)\[ f_{i}^{\mathrm{eq}}=w_{i}\left\{\rho + \frac{c_{i,\alpha} \left(\rho u_{\alpha}\right)}{c_{s}^{2}} + \frac{1}{2}\frac{\left(c_{i,\alpha}c_{i,\beta}-c_{s}^{2}\delta_{\alpha\beta}\right)\left(\rho u_{\alpha}u_{\beta}\right)}{c_{s}^{4}} + \frac{1}{6}\frac{\left[c_{i,\alpha}c_{i,\beta}c_{i,\gamma}-c_{s}^{2}\left(\delta_{\alpha\beta}c_{i,\gamma} + \delta_{\alpha\gamma}c_{i,\beta} + \delta_{\beta\gamma}c_{i,\alpha}\right)\right]\left[\rho u_{\alpha}u_{\beta}u_{\gamma}\right]}{c_{s}^{6}}\right\} \]

Following the same approach, the non-equilibrium populations are written as:

(15)\[ f_{i}^{\mathrm{neq}} = w_{i}\sum_{n=0}^{N}\frac{1}{n!}\mathbf{H}_{i}^{(n)}:\mathbf{a}^{(n),\mathrm{neq}} \]

Being the first moments of non-equilibrium populations:

(16)\[ \sum_{i=0}^{q-1}f_{i}^{\mathrm{neq}}=0, \]

(17)\[ \sum_{i=0}^{q-1}f_{i}^{\mathrm{neq}}\xi_{i,\alpha}=-\frac{\Delta t}{2c_{s}}F_{\alpha}, \]

(18)\[ \sum_{i=0}^{q-1}f_{i}^{\mathrm{neq}}\xi_{i,\alpha}\xi_{i,\beta}=\frac{1}{c_{s}^{2}}\Pi_{\alpha\beta}^{\mathrm{neq}}, \]

where \(\Pi_{\alpha\beta}^{\mathrm{neq}}=\sum_{i=0}^{q-1}f_{i}^{\mathrm{neq}}c_{i,\alpha}c_{i,\beta}\). From those moments, the coefficients \(\mathbf{a}^{(n),\mathrm{neq}}\) can be written recursively as:

(19)\[ a^{(0),\mathrm{neq}}=0, \]

(20)\[ a_{\alpha_{1}}^{(1),\mathrm{neq}}=-\frac{\Delta t}{2c_{s}}F_{\alpha_{1}}, \]

(21)\[ a_{\alpha_{1}\alpha_{2}}^{(2),\mathrm{neq}}=\frac{1}{c_{s}^{2}}\sum_{i}f_{i}^{\mathrm{neq}}c_{i,\alpha_{1}}c_{i,\alpha_{2}}, \]

(22)\[ a_{\alpha_{1}\alpha_{2}\alpha_{3}}^{(3),\mathrm{neq}} = \frac{1}{c_{s}}u_{\alpha_{3}}a_{\alpha_{1}\alpha_{2}}^{(2),\mathrm{neq}}+ \frac{1}{c_{s}}\left[ u_{\alpha_{1}}a_{\alpha_{2}\alpha_{3}}^{(2),\mathrm{neq}}+ u_{\alpha_{2}}a_{\alpha_{1}\alpha_{3}}^{(2),\mathrm{neq}} \right] \]

(23)\[ a_{\alpha_{1}...\alpha_{n}}^{(n),\mathrm{neq}}=\frac{1}{c_{s}}\sum_{k=1}^{n} u_{\alpha_{k}}\, a_{\alpha_{1}...\widehat{\alpha_{k}}...\alpha_{n}}^{(n-1),\mathrm{neq}} \]

where the sum runs symmetrically over all \(n\) indices: each index \(\alpha_{k}\) in turn carries the \(u\) factor while the remaining \(n-1\) indices (the hat \(\widehat{\alpha_{k}}\) denotes the omitted index) index the lower-order coefficient \(a^{(n-1),\mathrm{neq}}\). There is no privileged index; the operation is a symmetric sum over all indices, not a cyclic permutation.

Note

In Nassu the recursion is only ever evaluated up to third order: the Hermite basis is truncated at 3rd order, so Eq. (22) is the highest-order non-equilibrium coefficient actually built. Eq. (23) is given only to show the general structure.

The non-equilibrium populations truncated up to a third order Hermite polynomial expansion will be given by:

(24)\[ f_{i}^{\mathrm{neq}}=w_{i}\left[\mathbf{H}_{i}^{(1)}:\mathbf{a}^{(1),\mathrm{neq}} + \frac{1}{2}\mathbf{H}_{i}^{(2)}:\mathbf{a}^{(2),\mathrm{neq}} + \frac{1}{6}\mathbf{H}_{i}^{(3)}:\mathbf{a}^{(3),\mathrm{neq}} \right] \]

or:

(25)\[ f_{i}^{\mathrm{neq}}=w_{i}\left\{\frac{c_{i,\alpha}}{c_{s}}\left[-\frac{\Delta t}{2 c_{s}}F_{\alpha}\right] + \frac{1}{2}\left[\frac{c_{i,\alpha}c_{i,\beta}-c_{s}^{2}\delta_{\alpha\beta}}{c_{s}^{2}}\right]\left[\frac{1}{c_{s}^{2}}\Pi_{\alpha\beta}^{\mathrm{neq}}\right] + \frac{1}{6}\left[\frac{c_{i,\alpha}c_{i,\beta}c_{i,\gamma} - c_{s}^{2}\left(\delta_{\alpha\beta}c_{i,\gamma} + \delta_{\alpha\gamma}c_{i,\beta} + \delta_{\beta\gamma}c_{i,\alpha}\right)}{c_{s}^{3}}\right]\left[\frac{u_{\gamma}\Pi_{\alpha\beta}^{\mathrm{neq}} + \left(u_{\alpha}\Pi_{\beta\gamma}^{\mathrm{neq}} + u_{\beta}\Pi_{\alpha\gamma}^{\mathrm{neq}}\right)}{c_s^{3}}\right] \right\} \]

Important

Through this projection in an Hermite polynomial basis, the stability of LBM is greatly increased.

The discrete form of lattice Boltzmann equation is then written as:

(26)\[ f_i(\boldsymbol{x} + \boldsymbol{c}_i \Delta t, t+\Delta t) = f_i(\boldsymbol{x}, t)^{\mathrm{eq}} + \left(1-\omega\right)f_i(\boldsymbol{x}, t)^{\mathrm{neq}} + \Delta t\left(1-\frac{\omega}{2}\right)F_i(\boldsymbol{x}, t) \]

which describes the flow evolution.

Note

A simpler formulation of this development is the regularized-BGK (R-BGK) collision operator, in which the truncation of Hermite polynomials is performed only up to second order. This operator is also available in the present solver.

Hybrid Recursive Regularized-BGK (HRR-BGK)

A hybrid variant of the operator (HRR-BGK) exists in the literature [3]. It blends the mesoscopic rate-of-strain \(S_{\alpha\beta}\) with a finite-difference estimate \(S_{\alpha\beta}^{\mathrm{FD}}\) computed from the velocity field, through a parameter \(\sigma\):

(27)\[ S_{\alpha\beta}^{\mathrm{HRR}}=\sigma S_{\alpha\beta} + \left(1-\sigma\right)S_{\alpha\beta}^{\mathrm{FD}} \]

The intent is to add a small, controlled amount of dissipation from the smoother finite-difference estimate.

Note

HRR-BGK is not part of Nassu’s production path. The default and validated production operator is RR-BGK; HRR-BGK is available in the code but has not been thoroughly tested or validated, and is documented here only for completeness.

---

[1] (1,2)

Orestis Malaspinas. Increasing stability and accuracy of the lattice boltzmann scheme: recursivity and regularization. 2015. URL: https://arxiv.org/abs/1505.06900, doi:10.48550/ARXIV.1505.06900.

[2] (1,2)

Keijo K. Mattila, Paulo C. Philippi, and Luiz A. Hegele. High-order regularization in lattice-boltzmann equations. Physics of Fluids, 29(4):046103, 2017. URL: https://doi.org/10.1063/1.4981227, arXiv:https://doi.org/10.1063/1.4981227, doi:10.1063/1.4981227.

[3]

Jérôme Jacob, Orestis Malaspinas, and Pierre Sagaut. A new hybrid recursive regularised bhatnagar–gross–krook collision model for lattice boltzmann method-based large eddy simulation. Journal of Turbulence, 19(11-12):1051–1076, 2018. doi:10.1080/14685248.2018.1540879.