(rr_bgk)= # 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 {ref}`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 {footcite:t}`Malaspinas2015-1048` and {footcite:t}`Matilla2017-046`, 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 {footcite:t}`Malaspinas2015-1048` and {footcite:t}`Matilla2017-046`: $$ f_{i}^{\mathrm{eq}} = w_{i}\sum_{n=0}^{N}\frac{1}{n!}\mathbf{H}_{i}^{(n)}:\mathbf{a}^{(n),\mathrm{eq}} $$ (lb_equilibrium-hermite-expansion) 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: $$ H^{(0)}\left(\mathbf{\xi}_{i}\right) = 1, $$ (lb_H0) $$ H_{i,\alpha}^{(1)}\left(\mathbf{\xi}_{i}\right) = \xi_{i,\alpha}, $$ (lb_H1) $$ 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}}, $$ (lb_H2) $$ 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), $$ (lb_H3) $$ 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) $$ (lb_Hn) Due to orthogonality of Hermite polynomials, the value of the coefficients $\mathbf{a}^{(n),\mathrm{eq}}$ can be found with: $$ \mathbf{a}^{(n),\mathrm{eq}}=\sum_{i=0}^{q-1}\mathbf{H}_{i}^{(n)}f_{i}^{\mathrm{eq}}, $$ (lb_a-equilibrium) where $q$ is the number of velocity directions of the velocity set adopted. Knowing that the moments of equilibrium populations are: $$ \sum_{i=0}^{q-1}f_{i}^{\mathrm{eq}}=\rho, $$ (lb_mom0-equilibrium) $$ \sum_{i=0}^{q-1}f_{i}^{\mathrm{eq}}\xi_{i,\alpha}=\frac{1}{c_{s}}\rho u_{\alpha}, $$ (lb_mom1-equilibrium) $$ \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}. $$ (lb_mom2-equilibrium) The values of $\mathbf{a}^{(n),\mathrm{eq}}$ can be written recursively as: $$ a^{(0),\mathrm{eq}}=\rho, $$ (lb_a-eq0) $$ 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}}. $$ (lb_a-eqn) We adopt an expansion up to third order to represent the equilibrium populations, which gives: $$ 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] $$ (lb_equilibrium-pop) or: $$ 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\} $$ (lb_equilibrium-pop-exp) Following the same approach, the non-equilibrium populations are written as: $$ f_{i}^{\mathrm{neq}} = w_{i}\sum_{n=0}^{N}\frac{1}{n!}\mathbf{H}_{i}^{(n)}:\mathbf{a}^{(n),\mathrm{neq}} $$ (lb_non-equilibrium-hermite-expansion) Being the first moments of non-equilibrium populations: $$ \sum_{i=0}^{q-1}f_{i}^{\mathrm{neq}}=0, $$ (lb_mom0-non-equilibrium) $$ \sum_{i=0}^{q-1}f_{i}^{\mathrm{neq}}\xi_{i,\alpha}=-\frac{\Delta t}{2c_{s}}F_{\alpha}, $$ (lb_mom1-non-equilibrium) $$ \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}}, $$ (lb_mom2-non-equilibrium) 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: $$ a^{(0),\mathrm{neq}}=0, $$ (lb_a-neq0) $$ a_{\alpha_{1}}^{(1),\mathrm{neq}}=-\frac{\Delta t}{2c_{s}}F_{\alpha_{1}}, $$ (lb_a-neq1) $$ 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}}, $$ (lb_a-neq2) $$ 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] $$ (lb_a-neq3) $$ 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}} $$ (lb_a-neqn) 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. {math:numref}`lb_a-neq3` is the highest-order non-equilibrium coefficient actually built. Eq. {math:numref}`lb_a-neqn` 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: $$ 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] $$ (lb_non-equilibrium-pop) or: $$ 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\} $$ (lb_non-equilibrium-pop-exp) ```{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: $$ 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) $$ (lb_evolution-lbe) 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 {footcite:p}`jacob2018new`. 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$: $$ S_{\alpha\beta}^{\mathrm{HRR}}=\sigma S_{\alpha\beta} + \left(1-\sigma\right)S_{\alpha\beta}^{\mathrm{FD}} $$ (hrr-bgk) 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. ``` ```{eval-rst} .. footbibliography:: ```