Thermally Compressible Flow

Overview

The standard LBM BGK solver is designed for isothermal incompressible flow. For applications where the fluid density varies due to temperature or pressure changes — such as natural convection, gas-phase combustion, or packed-bed reactor flow — a thermally compressible extension is available.

Enable it with Do_ThermalComp = true in the input file. The thermal compressible (TC) variant uses an additional set of methods in FlowSolverLBM and requires coupling to the Temperature and Composition modules.

---

Thermodynamic pressure

The key concept in the TC extension is the separation of pressure into a thermodynamic and hydrodynamic part:

$$\begin{array}{}\text{(1)}& p_{\text{total}}=p_{\text{th}}+p_{\text{hydro}}\end{array}$$

• $p_{\text{th}}$ — spatially uniform thermodynamic pressure (set globally)
• $p_{\text{hydro}}$ — local hydrodynamic pressure fluctuation solved by LBM

For an ideal gas mixture:

$$\begin{array}{}\text{(2)}& p_{\text{th}}=\rho \frac{R}{M_w}T\end{array}$$

where $R$ is the universal gas constant, $M_w$ the molar mass, and $T$ the temperature.

The thermodynamic pressure is stored and updated globally:

double Pth;       ///< Current thermodynamic pressure [Pa]
double PthOld;    ///< Previous time step value [Pa]
double Pth0;      ///< Initial thermodynamic pressure [Pa]
double Poutlet;   ///< Outlet hydrodynamic pressure [Pa]

---

Density from ideal gas law

In the TC variant, the lattice density is not a conserved quantity but is recomputed each timestep from the ideal gas law:

void CalculateDensityTC(Temperature& Tx, Composition& Cx,
                         const BoundaryConditions& BC);
double CalculateIdealGasDensity(double p, double Rm, double Temp);

$$\begin{array}{}\text{(3)}& \rho =\frac{p_{\text{th}}{M}_w}{RT}\end{array}$$

This replaces the density moment of the populations in the incompressible case.

---

Density gradient

The gradient of density is needed for the compressible collision operator:

void CalculateDensityGradient(PhaseField& Phase, Velocities& Vel);

This gradient is stored in:

Storage3D<dVector3, 1> GradRho;

Three gradient schemes are supported:

Flag	Scheme
GradRho_Central	Second-order central differences
GradRho_Upwind	First-order upwind
GradRho_VanLeer	Van Leer flux limiter (recommended)

---

Thermally compressible equilibrium distribution

The TC variant uses a modified equilibrium distribution that accounts for the variable density and velocity divergence:

D3Q27 EquilibriumDistributionTC(dVector3 lbvel, double lbph, double lbnut,
                                  double lbDivVel, dVector3 lbFb,
                                  const double weights[3][3][3],
                                  double lbrho, dVector3 lbGradRho,
                                  size_t n);

The divergence of velocity $\mathrm{\nabla }\cdot \mathbf{u}$ enters because density is no longer conserved:

$$\begin{array}{}\text{(4)}& \frac{\partial \rho }{\partial t}+\mathrm{\nabla }\cdot (\rho \mathbf{u})=0\Rightarrow \mathrm{\nabla }\cdot \mathbf{u}=-\frac{1}{\rho }\frac{D\rho }{Dt}\end{array}$$

The divergence is stored in:

Storage3D<double, 1> DivVel;

---

Compressible collision step

The TC collision operator differs from the standard BGK by including:

• Variable kinematic viscosity (Sutherland's law)
• Non-zero velocity divergence
• Density gradient correction

void CollisionTC(Velocities& Vel);

Sutherland viscosity

The kinematic viscosity of a gas depends on temperature. Sutherland's law is used:

double CalculateSutherlandViscosity(double p, double mw, double Temp);

$$\begin{array}{}\text{(5)}& \mu (T)={\mu }_{\text{ref}}{\left(\frac{T}{T_{\text{ref}}}\right)}^{3/2}\frac{T_{\text{ref}}+S}{T+S}\end{array}$$

where $S$ is Sutherland's constant. The local viscosity is stored in:

Storage3D<double, 1> nut;       ///< Kinematic viscosity [lattice units]

---

Hydrodynamic pressure

The local hydrodynamic pressure is maintained separately from the thermodynamic pressure:

Storage3D<double, 1> HydroPressure;

void CalculateHydrodynamicPressureAndMomentum(Velocities& Vel);

This pressure drives the velocity field through the standard LBM mechanism while $p_{\text{th}}$ captures the global compressibility effect.

---

Initialization for TC flow

TC simulations require a special initialization of populations that accounts for the non-constant density field:

void SetInitialPopulationsTC(BoundaryConditions& BC, Velocities& Vel);

This should be called instead of the standard FinalizeInitialization() when Do_ThermalComp = true.

---

Related examples

Example	Description
examples/ThermalCompressibleFlow/CouetteFlow	Shear flow with variable density
examples/ThermalCompressibleFlow/PoiseuilleFlow	Channel flow with thermal effects
examples/ThermalCompressibleFlow/CylinderFlow	Flow around a heated cylinder
examples/ThermalCompressibleFlow/PackedBed2D	2D packed bed with compressibility
examples/ThermalCompressibleFlow/PackedBed3D	3D packed bed with compressibility

---

Numerical considerations

• The Mach number must remain low ($Ma<0.3$) for the low-Mach approximation to hold.
• Choosing GradRho_VanLeer is strongly recommended; central differences can introduce oscillations near steep density gradients.
• The update order matters: CalculateDensityTC must be called before CollisionTC.
• Thermodynamic pressure updates should be applied globally and consistently across MPI ranks to avoid non-physical pressure jumps.