(thermo_stability)=

# Balancing the dimensionless values and keeping the LBM stable

This is the practical centrepiece of the chapter. The {doc}`previous page <dimensionless>` showed that a buoyant flow is governed by $\mathrm{Pr}$, the geometry, and a driving group ($\mathrm{Ra}$ or $Q^*_H$), with the temperature ratio absorbed into the buoyancy velocity. This page turns that into the lattice knobs you actually set, and into the single constraint that decides whether a thermal run is stable: the buoyancy force per step must stay small. It mirrors the {ref}`units-page recipe <lattice_units>` for the isothermal case, adapted for the extra freedom buoyancy brings.

## The binding constraint: the force Mach number

An isothermal run is bounded by the velocity Mach number, $\mathrm{Ma} = U/c_s < 0.1$. A buoyant run inherits that bound, but the velocity is now *emergent*: the flow accelerates itself to the buoyancy velocity $U_b = \sqrt{G\,\beta\,\Delta T\, L}$ {eq}`td_buoyancy_velocity`. So the real constraint is on the buoyancy *force per step*, which must be small enough that the velocity it produces stays subsonic.

```{admonition} One knob too hard and the lattice goes supersonic
---
class: warning
---
The peak lattice velocity in a buoyant flow scales as $U_b \sim \sqrt{G_{\text{lat}}\,\beta\,\phi\,H_{\text{lat}}}$, where $G_{\text{lat}}$ is the lattice gravity, $\beta$ the expansion coefficient, $\phi$ the scalar amplitude and $H_{\text{lat}}$ the reference height in nodes. Demanding a physically large temperature ratio ($\beta\,\phi \sim 0.3$) at full lattice gravity ($G_{\text{lat}} = 1$) drives $U_b$ well above $c_s = 1/\sqrt{3}$: the run goes supersonic and the weakly-compressible approximation collapses. Keep $U_b$ in the same band as the isothermal lattice velocity, roughly $0.05$ to $0.15$.
```

## The lattice-gravity lever

The way out is the freedom the densimetric scaling hands you. Because the temperature ratio is absorbed into $U_b$, you do not have to shrink $\beta\,\phi$ to keep $U_b$ subsonic - you can instead lower the **lattice gravity** $G_{\text{lat}}$. The product $G_{\text{lat}}\,\beta\,\phi$ is what sets $U_b$, and matching the physical driving group ($Q^*_H$ or $\mathrm{Ra}$) fixes that product, not its individual factors. So:

- Set $\beta$ and the scalar amplitude $\phi$ to represent the *physical* temperature ratio ($\theta = \beta\,\phi \approx 0.3$ for a real fire), so the temperature field reconstructs correctly.
- Lower $G_{\text{lat}}$ (for example to $10^{-3}$) so that $U_b = \sqrt{G_{\text{lat}}\,\beta\,\phi\,H_{\text{lat}}}$ lands in the safe $0.05$-$0.15$ band.
- Check that the driving group still matches: $Q^*_H$ and $\mathrm{Ra}$ depend on $G\,\beta\,\Delta T$ together, and lowering $G_{\text{lat}}$ is compensated by the longer time the flow has (in lattice steps) to develop, so the matched non-dimensional flow is preserved.

```{admonition} Lattice gravity is a similarity lever, not a physical change
---
class: note
---
Lowering $G_{\text{lat}}$ does not make the simulated fire weaker. It rescales the lattice time step relative to the buoyant time scale so the same non-dimensional flow unfolds at a subsonic lattice velocity. It is the thermal analogue of choosing a small lattice velocity $U$ in the isothermal recipe: a unit choice that buys stability, undone exactly by the conversion factors when results are read back.
```

## The temperature ratio is emergent, not dialed

A recurring setup mistake is to treat the hot-layer temperature ratio as an input to be tuned to a target. It is not. With a heat *source* (rather than a fixed wall temperature), the steady temperature ratio emerges from the balance between the source and the buoyant venting; the knob you actually hold is the source strength.

```{admonition} Tune the source, read the emergent ratio
---
class: important
---
The GPU knob is the volumetric source rate (or the wall temperature, for fixed-temperature problems). Run the case, then *measure* the emergent effective ratio $\theta_{\text{eff}} = \beta\,(\phi_{\text{hot}} - \phi_{\text{amb}})$ off the field, and adjust the source until $\theta_{\text{eff}}$ reaches the physical target (for a compartment fire, $\theta_{\text{eff}} \approx 0.3$). You set $\beta$ and $G_{\text{lat}}$ once for the similarity match; the source rate is what you iterate, and $\theta_{\text{eff}}$ is an output you check, never a value you write directly.
```

## Dimensional reconstruction by similarity

Once a matched run has reached steady state, its dimensionless profiles (velocity normalised by $U_b$, temperature by $\Delta T$) are the physical profiles. Reconstruct dimensional fields from the run's *own* emergent ratio:

$$
\Delta T_{\text{phys}} = \theta_{\text{eff}}\, T_0, \qquad
U_{b,\text{phys}} = \sqrt{G_{\text{phys}}\, H_{\text{phys}}\, \theta_{\text{eff}}}
$$ (ts_reconstruct)

so a probed lattice velocity $u_{\text{norm}}$ (normalised by the lattice $U_b$) maps to $u_{\text{phys}} = u_{\text{norm}}\, U_{b,\text{phys}}$, and a normalised temperature $T_{\text{norm}} \in [0,1]$ maps to $T_{\text{phys}} = T_0 + \Delta T_{\text{phys}}\, T_{\text{norm}}$.

```{admonition} Temperature scales as the ratio, velocity as its square root
---
class: warning
---
Reconstruction respects the $\lambda$ versus $\sqrt{\lambda}$ split noted on the {ref}`dimensionless page <thermo_dimensionless>`: the temperature scale is $\theta_{\text{eff}}\,T_0$ (linear in the ratio) but the velocity scale is $\sqrt{G H\,\theta_{\text{eff}}}$ (square-root in the ratio). Multiplying every field by one common factor to "boost" a reduced-ratio run to the physical ratio is therefore incorrect - it would over-scale the velocity. Use {eq}`ts_reconstruct` with the two different exponents.
```

```{admonition} The Boussinesq ceiling
---
class: note
---
Reconstruction is only as good as the Boussinesq approximation underneath it. At $\theta_{\text{eff}} \approx 0.3$ the real density varies by tens of percent, which is the edge of the constant-density assumption {eq}`te_buoyancy`. The reconstructed fields inherit that modelling error; it is the intrinsic ceiling of a Boussinesq surrogate for a strongly heated flow, not a tuning artefact. A flow that needs more than this must move to the {doc}`variable-density low-Mach route <modeling_hierarchy>`, where the density is slaved to temperature through the equation of state rather than frozen.
```

## The balancing recipe

For a buoyant, heat-source-driven thermal run the workflow is:

1. Match $\mathrm{Pr}$ (set $\kappa = \nu/\mathrm{Pr}$) and the geometry; identify the physical driving group ($Q^*_H$ for a fire, $\mathrm{Ra}$ for fixed-temperature convection).
1. Set $\beta$ and the scalar amplitude to represent the physical temperature ratio $\theta \approx \beta\,\phi$.
1. Choose the lattice gravity $G_{\text{lat}}$ so the buoyancy velocity $U_b = \sqrt{G_{\text{lat}}\,\beta\,\phi\,H_{\text{lat}}}$ falls in the $0.05$-$0.15$ band (subsonic, weakly compressible).
1. Confirm the driving group matches the physical value at these choices (compute $Q^*_H$ or $\mathrm{Ra}$ in lattice units and compare).
1. Decide which groups to leave unmatched - typically the Reynolds/Grashof number, relying on Reynolds-independence of the gross exchange - and note it.
1. Run, measure the emergent $\theta_{\text{eff}}$, and iterate the source rate until it hits the physical target.
1. Reconstruct dimensional fields with {eq}`ts_reconstruct`, respecting the $\lambda$ / $\sqrt{\lambda}$ split.

```{admonition} The balance in one sentence
---
class: important
---
> Set $\beta$ to the physical temperature ratio, lower the lattice gravity until the buoyancy velocity is subsonic, match the driving group ($Q^*$ or $\mathrm{Ra}$) rather than the temperature ratio, tune the source so the emergent ratio is physical, and reconstruct with temperature linear and velocity square-root in that ratio.
```

This recipe is exactly what the {doc}`Steckler room-fire validation case </validation/buoyant_flows/03_steckler_room_fire/README>` follows: it matches $Q^*_H$ and $\mathrm{Pr}$, lowers the lattice gravity to keep the doorway jet subsonic, leaves the Reynolds number reduced, and reconstructs dimensional doorway profiles from the run's own emergent ratio.

(thermo_stability_lowmach)=

## The reduced-pressure mode in the variable-density closure

The {doc}`variable-density low-Mach route <modeling_hierarchy>` adds a stability concern the Boussinesq route does not have, because it changes what the fluid zeroth moment *means*. In that closure the streamed populations no longer carry density; they carry the **reduced hydrodynamic pressure** $\theta^h = p^h/(\rho c_s^2)$ (Taha et al. 2024, Eq. 33), and the density is closed separately from temperature through the equation of state $\rho = P/(rT)$. Temperature changes therefore drive a dilatation $\partial_\alpha u_\alpha = -\tfrac{1}{\rho}\,\mathrm{D}\rho/\mathrm{D}t$ that feeds the pressure field every step.

At the level of the continuous equations this coupling is benign. With gravity switched off it is the classical zero-Mach (Majda-Sethian) system: the temperature is parabolic and strictly dissipative, and the pressure acts as an **elliptic** Lagrange multiplier that enforces the divergence constraint instantaneously. That system has no growing high-wavenumber mode. The instability is introduced by the discretisation, not the model.

```{admonition} Why the bare scheme admits a checkerboard mode
---
class: warning
---
The LBM does not solve an elliptic pressure projection. It carries $\theta^h$ in the zeroth moment and relaxes it **acoustically**, at the finite lattice speed $c_s$. The explicit, one-step-lagged dilatation source is therefore injected into a finite-speed, un-projected pressure variable that has no damping at the grid scale. The result is an under-damped, highest-wavenumber ($k = \pi$, odd-even / checkerboard) acoustic mode: a smooth temperature bump seeds a smooth dilatation, the centred pressure-velocity coupling fails to damp the $k = \pi$ component, and it grows geometrically while the physical buoyant signal evolves underneath it. Shear viscosity does **not** reach this mode (it is a pressure checkerboard, not a shear mode), which is why raising $\tau$ does not cure it, and a continuous volumetric heat source makes it worse (more sustained dilatation pumps more acoustic energy, and $\mathrm{Ma}$ climbs).
```

### The bulk-viscosity stabiliser

The stabiliser damps exactly the mode the bare scheme fails to damp, and nothing else: it supplies the **bulk (volume) viscosity** the weakly-compressible lattice leaves out, applied inside the collision through `models.LBM.bulk_viscosity`. This is an independent relaxation rate $\omega_{\text{bulk}} \in (0,2)$ acting on the trace of the non-equilibrium stress, **off by default** so the isothermal and Boussinesq paths, and any low-Mach run that does not need it, are bit-identical. It reaches the $k = \pi$ reduced-pressure mode through the collision - the second-moment (stress) channel, which is non-zero at Nyquist.

The physical effect this realises is most transparent in its continuum analogue: a light Laplacian low-pass of the reduced pressure applied after its update,

$$
\theta^h \;\leftarrow\; \theta^h + \sigma\,\nabla^2\theta^h ,
\qquad 0 \le \sigma \le \frac{1}{2d},
$$ (ts_pressure_filter)

for spatial dimension $d$. In Fourier space such a filter multiplies a mode of wavenumber $k$ by $1 - 2\sigma\,(1 - \cos k)$ per axis: it is the identity at $k \to 0$ (the smooth buoyant signal passes through untouched) and it is most aggressive at $k = \pi$, where it scales the offending checkerboard component by $1 - 4\sigma$ per axis. The upper bound $\sigma \le 1/(2d)$ is the explicit-diffusion stability limit, the same constraint the energy Laplacian obeys; at $\sigma = 1/(2d)$ the $k = \pi$ mode is annihilated in one pass. The in-collision bulk viscosity delivers the same selective high-wavenumber absorption without that explicit-diffusion cap, because it is a relaxation with $\omega_{\text{bulk}} \in (0,2)$ rather than an explicit Laplacian step.

### What the stabiliser means physically

The low-pass model {eq}`ts_pressure_filter` makes the continuum meaning explicit. Adding $\sigma\,\nabla^2\theta^h$ to the pressure update is, term by term, a **diffusion of the hydrodynamic pressure**: over a step the reduced pressure evolves as $\partial_t \theta^h = \cdots + \sigma\,\nabla^2\theta^h$, a Laplacian relaxation with diffusivity $\sigma$.

A diffusion of pressure is precisely what **bulk (volume) viscosity** contributes to a real gas. The attenuation of a sound wave scales as $\delta\,k^2$, where the *sound diffusivity* $\delta = \tfrac{1}{\rho}\big(\tfrac{4}{3}\mu + \zeta\big)$ collects the shear and **bulk** viscosities $\mu$, $\zeta$. The weakly-compressible LBM is built with a near-zero bulk viscosity (the trace of the stress relaxes at the same rate as the shear), so it damps compression and dilatation waves far too weakly. The variable-density coupling then pours dilatation into that under-damped acoustic field every step, and the grid-scale sound has nowhere to go. The bulk-viscosity stabiliser **restores the missing acoustic absorption**: it raises the sound diffusivity at high wavenumber so the spurious compression mode is absorbed - the physical role $\zeta$ plays, supplied where the lattice left it out. The pressure low-pass {eq}`ts_pressure_filter` is the same effective bulk viscosity seen as a selective diffusion of the pressure field; the in-collision relaxation is the realisation that supplies it through the physical stress channel.

Seen another way, the effective bulk viscosity nudges the pressure toward its local neighbourhood mean, i.e. partially back toward the **elliptic** (instantaneous) pressure of the true low-Mach limit. The continuum low-Mach pressure communicates across the domain infinitely fast (the divergence constraint is global); the LBM replaces that with finite-speed acoustics; the absorption returns a little of that instantaneous, smoothed-out character to the field.

```{admonition} A pressure-side cousin of the regularized (RR/HRR) collision
---
class: note
---
This places the bulk-viscosity stabiliser in the same family as the **regularized** collision operators that stabilise LBM more generally. Recursive-regularised (RR) and hybrid (HRR) kernels cure spurious high-frequency and acoustic behaviour by projecting out, or selectively dissipating, the non-hydrodynamic **ghost moments** carried by the non-equilibrium populations - removing energy from modes that are numerical artefacts while leaving the resolved hydrodynamics untouched. The bulk viscosity does the same *kind* of thing, targeting the **trace** of the non-equilibrium stress, the channel that couples to the reduced pressure at high wavenumber. That is why it reaches the mode the shear-relaxation does not: the standard collision damps spurious content in the *deviatoric* stress, whereas the variable-density instability lives in the *pressure* field, so it needs a selective dissipation of the same spirit applied through the trace.
```

```{admonition} The stabiliser treats the symptom; the lattice corrections matter too
---
class: note
---
The bulk viscosity is a robust, route-agnostic stabiliser, but it is not the only lever. A pressure-carrying low-Mach LBM also relies on the lattice correction tensors of its reference scheme for isotropy and, in some formulations, for grid-scale dissipation. On the complete D3Q27 Hermite basis the D3Q19 isotropy patches are unnecessary and are dropped; on a reduced lattice (D3Q19) that shows the same mode, the first check is whether a dropped correction carried the missing dissipation, before reaching for a larger bulk viscosity. The dilatation per step $(\rho^{n} - \rho^{n-1})/\rho$ should also be kept well below the acoustic CFL: when a continuous heat source drives it too hard the velocity leaves the weakly-compressible band regardless of the stabiliser.
```

```{eval-rst}
.. footbibliography::
```