Compressible LBM

Every page in this chapter so far has assumed the flow is isothermal and weakly compressible: the temperature is uniform, the density barely moves, and the Mach number is held below the \(\mathrm{Ma} < 0.1\) threshold derived on the Chapman-Enskog page. That assumption covers the large majority of the incompressible flows Nassu targets, where the flow is far slower than the speed of sound and density changes are negligible. A few applications fall outside it: high-speed flows where compressibility is dynamically important, flows that contain shock waves, and non-isothermal flows where temperature differences alter the dynamics.

This page sets out the theory that extends the method to those regimes. It generalizes the equilibrium so that it carries a temperature, adds the moment correction the standard velocity sets need to stay accurate once the temperature varies, and supplements the mass and momentum balances with an energy equation and a shock-capturing term. It is an extension of the isothermal core rather than a replacement: the collision, regularization and LES machinery of the previous pages carry over, now with temperature threaded through them.

Compressible energy equation vs Boussinesq scalar

Reach for this fully compressible, energy-equation formulation only when compressibility is dynamically important: high-speed flows, shocks, or large temperature-driven density changes. Buoyancy-driven thermal flows where density barely moves (atmospheric stratification, urban heat island, indoor convection) are handled by the lighter route in the scalar-transport chapter, where temperature is carried as a transported passive scalar with a Boussinesq buoyancy force (see Thermal Extension (Boussinesq)).

Making the method suitable for compressible flow requires two changes to the isothermal core: a temperature-dependent equilibrium, and stabilising terms for the steep gradients that compressibility produces [1][2][3]. The first step is to let the equilibrium carry temperature variations, so the Maxwell-Boltzmann distribution is used in its more complete form [4]:

(1)\[ f^{\mathrm{eq}} = \rho\left[\frac{1}{2\pi c_{s}^{2}\left(\theta + 1\right)}\right]^{\frac{3}{2}}\mathrm{exp}\left[-\frac{\left(\xi_{\alpha}-u_{\alpha}\right)^{2}}{2 c_{s}^{2}\left(\theta + 1\right)}\right] \]

where \(\theta=(T/T_{0})-1\). When the Hermite polynomial expansion of equilibrium distribution function is performed for \(\theta \neq 0\), the coefficients \(\mathbf{a}^{(n),\mathrm{eq}}\) up to third-order are:

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

(3)\[ a^{(1),\mathrm{eq}}_{\alpha_{1}}=\rho u_{\alpha_{1}} \]

(4)\[ a^{(2),\mathrm{eq}}_{\alpha_{1}\alpha_{2}}=\rho u_{\alpha_{1}}u_{\alpha_{2}}+\rho c_{s}^{2}\theta\,\delta_{\alpha_{1}\alpha_{2}} \]

(5)\[ a^{(3),\mathrm{eq}}_{\alpha_{1}\alpha_{2}\alpha_{3}}=\rho u_{\alpha_{1}}u_{\alpha_{2}}u_{\alpha_{3}}+\rho c_{s}^{2}\theta\left[u\delta\right]_{\alpha_{1}\alpha_{2}\alpha_{3}} \]

where \(\left[u\delta\right]_{\alpha_{1}\alpha_{2}\alpha_{3}}=\left(u_{\alpha_{1}}\delta_{\alpha_{2}\alpha_{3}}+u_{\alpha_{2}}\delta_{\alpha_{1}\alpha_{3}}+u_{\alpha_{3}}\delta_{\alpha_{1}\alpha_{2}}\right)\). These are the same equilibrium moments used throughout the isothermal chapter (see moments collision), now carrying the temperature through the \(\rho c_{s}^{2}\theta\) terms; setting \(\theta = 0\) recovers the isothermal coefficients exactly.

Temperature as a moment of the populations

Unlike the Boussinesq scalar route, where temperature is a separately transported field, here the temperature is a moment of the same populations that carry mass and momentum, exactly as \(\rho\) and \(u_{\alpha}\) are. Its contracted second moment fixes it [5]:

(6)\[ \rho\left[D\left(\theta + 1\right) + u_{\alpha}u_{\alpha}\right] = \sum_{i} f_{i}\,c_{i,\alpha}c_{i,\alpha} \]

where \(D\) is the spatial dimension and \(u_{\alpha} = \rho^{-1}\sum_{i} f_{i}c_{i,\alpha}\) is the bare momentum velocity. Solving for the deviation gives the equivalent central-moment form

(7)\[ \theta = -1 + \frac{1}{D\rho}\sum_{i} f_{i}\left\lVert c_{i} - u\right\rVert^{2} \]

so \(\theta = 0\) whenever the populations sit at the reference temperature \(T_{0}\). He et al.[5] write the same relation in an absolute temperature \(T/T_{0} = \theta + 1\); the two differ only by this reference shift.

The D3Q27 velocity set recovers the isothermal third-order equilibrium moment exactly, which is why it is the preferred isothermal production set. The third-moment correction the standard sets carry measures the gap between the continuous Maxwell-Boltzmann third moment and the lattice one,

(8)\[ \Psi_{\alpha_{1}\alpha_{2}\alpha_{3}} = \Pi_{\alpha_{1}\alpha_{2}\alpha_{3}}^{\mathrm{eq,MB}} - \Pi_{\alpha_{1}\alpha_{2}\alpha_{3}}^{\mathrm{eq}} \]

where \(\Pi_{\alpha_{1}\alpha_{2}\alpha_{3}}^{\mathrm{eq}}=\sum f_{i}^{\mathrm{eq}}c_{i,\alpha_{1}}c_{i,\alpha_{2}}c_{i,\alpha_{3}}\), and \(\Pi_{\alpha_{1}\alpha_{2}\alpha_{3}}^{\mathrm{eq,MB}}=\int f^{\mathrm{eq}}c_{\alpha_{1}}c_{\alpha_{2}}c_{\alpha_{3}}d\mathbf{c}\).

On D3Q27 every mixed component of \(\Psi\) is zero, at any temperature: the temperature-dependent mixed third moments are recovered exactly, exactly as the isothermal ones are (D3Q19 misses them - that is the gap a correction repairs there). The only non-zero component on D3Q27 is the fully-diagonal

(9)\[ \Psi_{\alpha_{1}\alpha_{1}\alpha_{1}} = \rho u_{\alpha_{1}}\left(\theta+u_{\alpha_{1}}^{2}\right), \]

which is the weak-compressibility truncation of the third-order Hermite expansion: the \(u_{\alpha_{1}}^{3}\) piece is the same \(\mathcal{O}(\mathrm{Ma}^{3})\) error the isothermal scheme already tolerates (the diagonal lattice moment aliases through \(c_{\alpha}^{3}=c_{\alpha}\) on the \(\{-1,0,1\}\) axes), now with a thermal \(\theta u_{\alpha_{1}}\) addition. It is identical on D3Q19, so it is a universal Hermite-truncation residual, not a lattice-isotropy defect. This is pinned numerically in tests/mod/t_LBM/testLatticeMomentIsotropy.py.

Hence on D3Q27 the weakly-compressible thermal coupling does not apply a \(\Psi\) source term: the only residual it would remove is the diagonal weak-compressibility truncation that the isothermal scheme already accepts under \(\mathrm{Ma} < 0.1\), and the mixed parts are exact. The explicit correction below is retained for the fully compressible regime (high Mach, shocks) and for reduced lattices like D3Q19, where the mixed components do not vanish [1][2].

The correction term is decomposed in a mesoscopic \(\psi_{i}\) through a Hermite polynomial expansion:

(10)\[ \psi_{i} = -w_{i}\frac{H_{i,\alpha_{1}\alpha_{2}}^{(2)}}{2c_{s}^{4}}\frac{\partial}{\partial x_{\alpha_{3}}}\Psi_{\alpha_{1}\alpha_{2}\alpha_{3}} \]

Which, for the D3Q27 velocity results in:

(11)\[ \psi_{i} = \frac{w_{i}}{2c_{s}^{4}}\left\{ H_{i,xx}\frac{\partial}{\partial x}\left[\rho u_{x}\left(1-\frac{p}{\rho c_{s}^{2}} - u_{x}^{2}\right)\right]\right\} \]

where for a non-isothermal flow, the pressure is given by \(p = \rho c_{s}^{2}\left(\theta + 1\right)\) and the physical speed of sound by \(c_{sp}=\sqrt{p/\rho}\). The 0th and 1st moments of \(\psi_{i}\) are null, while the second moment is:

(12)\[ \sum_{i}\psi_{i}c_{i,\alpha_{1}}c_{i,\alpha_{2}} = - \frac{\partial}{\partial x_{\alpha_{3}}}\Psi_{\alpha_{1}\alpha_{2}\alpha_{3}} \]

The correction term is applied as a source term in LBE collision:

(13)\[ f_{i}^{*}=f_{i}^{\mathrm{eq}} + \left(1-\omega\right)\tilde{f}_{i}^{\mathrm{neq}} + \Delta t\left(1-\frac{\omega}{2}\right)F_i + \frac{\Delta t}{2}\psi_{i} \]

Where \(\tilde{f}_{i}^{\mathrm{neq}}\) is the regularized non-equilibrium population. The regularized population is calculated under the consideration that \(\Pi_{\alpha_{1}\alpha_{2}}^{\mathrm{neq}} = \sum \left(f_{i}^{\mathrm{neq}} + 0.5\Delta t \psi_{i}\right)c_{i,\alpha}c_{i,\beta}\).

For the macroscopics described collision, only \(\Pi_{\alpha_{1}\alpha_{2}}^{\mathrm{neq},*}\) will be affected by the addition of the correction term and will be calculated as:

(14)\[ \Pi_{\alpha_{1}\alpha_{2}}^{\mathrm{neq},*} = \left(1- \omega\right)\left(\Pi_{\alpha_{1}\alpha_{2}}^{\mathrm{neq}} - \frac{\Delta t}{2}\frac{\partial}{\partial x_{\alpha_{3}}}\Psi_{\alpha_{1}\alpha_{2}\alpha_{3}}\right) + \left(\Pi_{\alpha_{1}\alpha_{2}}^{\mathrm{eq}}- \Pi_{\alpha_{1}\alpha_{2}}^{\mathrm{eq}*}\right) + \left(1 - \frac{\omega}{2}\right)\left(F_{\alpha_{1}}u_{\alpha_{2}} + F_{\alpha_{2}}u_{\alpha_{1}}\right) - \frac{\Delta t}{2}\frac{\partial}{\partial x_{\alpha_{3}}}\Psi_{\alpha_{1}\alpha_{2}\alpha_{3}} \]

If the Chapman-Enskog expansion is performed taking into account the temperature, the viscosity is written in terms of mesoscopic variables as:

(15)\[ \eta = \rho c_{s}^{2}\Delta t\left(\frac{1}{\omega} - \frac{1}{2}\right)\left(\theta + 1\right) \]

Hence, the relaxation frequency must be tuned from a user-defined viscosity:

(16)\[ \omega = \left[\frac{1}{2}+\frac{\eta}{\rho c_{s}^{2}\Delta t\left(\theta + 1\right)}\right]^{-1} \]

In our Smagorinsky-LES formulation, the effective relaxation frequency is corrected to:

(17)\[ \omega^{*}=\frac{1}{2}\left(\omega_{0}^{-1}+\sqrt{\omega_{0}^{-2}+\frac{2(C_{\mathrm{S}}^{2}\sqrt{2Q_{\alpha\beta}Q_{\alpha\beta}})}{\Delta t \rho c_{s}^{4}\left(\theta + 1\right)^{2}}}\right) \]

where \(\omega_{0}\) is calculated from a reference viscosity \(\eta_{0}\).

Post-collision temperature update

In the macroscopic-storage collision the populations are not retained between steps; the macroscopics, \(\theta\) among them, are advanced directly through their post-collision values. Re-evaluating the temperature moment (6) on the post-collision populations and using the contracted second moment of the Guo force, \(\sum_{i} F_{i}c_{i,\alpha}c_{i,\alpha} = 2 F_{\alpha}u_{\alpha}\), gives the post-collision deviation

(18)\[ \theta^{*} = \theta + \frac{1}{D}\left(\left\lVert u\right\rVert^{2} - \left\lVert u^{*}\right\rVert^{2}\right) + \frac{2\Delta t}{\rho D}\left(1 - \frac{\omega}{2}\right)F_{\alpha}u_{\alpha} \]

where \(u_{\alpha} = \rho^{-1}\sum_{i} f_{i}c_{i,\alpha}\) is the bare momentum velocity and \(u^{*}_{\alpha} = \left(\rho u_{\alpha} + \Delta t F_{\alpha}\right)/\rho\) the post-collision velocity. The two contributions are the change in kinetic energy and the work done by the body force: together they are the compression work an isothermal collision does not need, and they keep the temperature consistent with the \(\rho c_{s}^{2}\theta\) pressure that (4) injects into the equilibrium. Omitting it leaves that thermal pressure unbalanced and a sustained run diverges. Equation (18) is the adiabatic parcel of the energy equation below; heat conduction and viscous dissipation are added on top.

The compression work uses the bare momentum velocity

Both \(\left\lVert u\right\rVert^{2}\) and the work term \(F_{\alpha}u_{\alpha}\) in (18) use the bare momentum velocity \(u_{\alpha} = \rho^{-1}\sum_{i} f_{i}c_{i,\alpha}\), the same velocity that closes the temperature moment (6), and not the Guo half-corrected velocity \(u_{\alpha} + \tfrac{\Delta t}{2}\rho^{-1}F_{\alpha}\) used elsewhere in the collision. The body force already enters (18) through \(u^{*}\) and through the explicit \(F_{\alpha}u_{\alpha}\) term; substituting the half-corrected velocity counts the force a third time and injects a spurious \(\mathcal{O}(\Delta t\,F_{\alpha}u_{\alpha})\) heating, of the same order as the physical work term.

The force-free compressible equilibrium and the diagnostic recovery (6) are those of He et al.[5]; the body force and the post-collision update (18) extend that scheme to Nassu’s Guo-forced collision and follow from (6) directly.

Validity and scale matching

Because the temperature is coupled directly into the equilibrium, the viscosity and the collision rather than carried as an independent passive scalar, the fluid and its energy equation are solved as a single coupled system, and their scales must be mutually consistent. This is what lets the method capture compressibility and thermo-acoustic effects the Boussinesq scalar route cannot, at the cost of extra constraints the user must respect:

• Small-deviation bound. The Hermite expansion is truncated at third order, so the equation of state \(p = \rho c_{s}^{2}\left(\theta + 1\right)\) is consistent only to leading order in \(\theta\); the model is trustworthy for \(\left|\theta\right| \ll 1\) (in practice \(\left|\theta\right| \lesssim 0.15\)).

• Local Mach number. The local sound speed is \(c_{s}\sqrt{\theta + 1}\), so a cold region (\(\theta < 0\)) lowers it and raises the effective Mach number. The \(\mathrm{Ma} < 0.1\) limit must be budgeted against the coldest expected \(\theta\): \(\left|u\right| < 0.1\,c_{s}\sqrt{1 + \theta_{\min}}\).

• Prandtl number. At ranks two and three the temperature relaxes on the same \(\omega\) as momentum, which pins the Prandtl number to order unity. An independent thermal diffusivity requires the conduction term of the energy equation below (a separate finite-difference parcel), not the collision alone.

• Buoyancy is not specified twice. The temperature already changes the pressure through the equation of state, so a direct-coupling run should drive buoyancy through that coupling (gravity acting on the temperature-induced density anomaly), and not additionally through a Boussinesq force, or the buoyant response is double-counted.

Energy Equation

The temperature-dependent equilibrium above closes the mass and momentum balances, but the temperature \(\theta\) it relies on must itself be advanced in time. That is the role of an energy equation, the third conservation law that an isothermal flow does not need. In order to solve for the temperature, the mass and momentum conservations are therefore supplemented by the entropy equation, written as [1][2]:

(19)\[ \frac{\partial s}{\partial t}+ u_{\alpha} \frac{\partial s}{\partial x_{\alpha}}=\frac{1}{\rho T}\frac{\partial}{\partial x_{\alpha}}\left(\lambda\frac{\partial T}{\partial x_{\alpha}}\right)+\frac{1}{\rho T}\sigma_{\alpha\beta}\frac{\partial u_{\alpha}}{\partial x_{\beta}} \]

where the entropy \(s=c_{v}\mathrm{ln}(p/\rho^{\gamma})\) is the entropy with \(c_{v}\) being the specific heat capacity at constant volume and \(\gamma\) the specific heat ratio. \(\lambda\) is the heat conductivity and \(\sigma_{\alpha\beta}\) is the stress tensor.

Shock Capturing Method

Once compressibility is dynamically important the flow can develop shock waves: near-discontinuities in pressure and density that a centred scheme cannot represent without spurious oscillations. The solution is to add a small, locally targeted artificial viscosity that switches on only where a shock is detected. To compute a flow containing shock waves, we adopt the shock sensor and associated viscosity used in the Jameson-Schmidt-Turkel (JST) scheme [6]:

(20)\[ \epsilon_{\alpha}=\kappa\left|\frac{p_{i-1}-2p_{i}+p_{i+1}}{p_{i-1}+2p_{i}+p_{i+1}}\right| \]

where \(i\) is the index along a cartesian grid line, \(\kappa\) is a free parameter that tunes the amount of artificial viscosity, and \(\epsilon_{\alpha}\) is the sensor evaluated along lattice direction \(\alpha\). Applying Equation (20) independently along \(x\), \(y\) and \(z\) yields a directional-sensor vector \(\boldsymbol{\epsilon}=\left(\epsilon_{x},\epsilon_{y},\epsilon_{z}\right)\), which is combined into a single effective relaxation frequency:

(21)\[ \omega_{e} = \left\{\frac{1}{\omega} + \left\lVert\boldsymbol{\epsilon}\right\rVert\right\}^{-1}, \qquad \left\lVert\boldsymbol{\epsilon}\right\rVert = \sqrt{\epsilon_{x}^{2}+\epsilon_{y}^{2}+\epsilon_{z}^{2}} \]

The directional sensors are combined through the Euclidean norm of \(\boldsymbol{\epsilon}\) rather than through their maximum \(\mathrm{max}\left[\epsilon_{x},\epsilon_{y},\epsilon_{z}\right]\), because a shock front can be misaligned with the lattice axes even though it tends to align with them. The norm accounts for all directions at once, while the maximum keeps only the strongest axial component.

---

[1] (1,2,3)

Yongliang Feng, Pierre Boivin, Jérôme Jacob, and Pierre Sagaut. Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows. Journal of Computational Physics, 394:82–99, 2019.

[2] (1,2,3)

S Guo, Yongliang Feng, Jérôme Jacob, Florian Renard, and Pierre Sagaut. An efficient lattice boltzmann method for compressible aerodynamics on d3q19 lattice. Journal of Computational Physics, 418:109570, 2020.

[3]

S Guo, Y Feng, and P Sagaut. Improved standard thermal lattice boltzmann model with hybrid recursive regularization for compressible laminar and turbulent flows. Physics of Fluids, 2020.

[4]

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.

[5] (1,2,3)

Bing He, Yingchun Chen, Weibing Feng, Qing Li, Anping Song, Yang Wang, Miao Zhang, and Wu Zhang. Compressible lattice Boltzmann method and applications. International Journal of Numerical Analysis and Modeling, 9(2):410–418, 2012.

[6]

Antony Jameson. Origins and further development of the jameson–schmidt–turkel scheme. AIAA Journal, 55(5):1487–1510, 2017.