Hi, I am trying to understand how thermal LBM works, for that purpose I wrote simple code in Python from formulas I found in nice review paper by Sharma et al. My code use double distribution for momentum and temperature calculation.

It looks like that it is working fine for momentum only or temperature only, but not for both, e. g. I quickly got instability when try to simulate natural convection in 2D.

Just a guess, but it looks like that energy(temperature) is not in sync with mass (momentum), e. g. as seen in picture above, temperature is transferred slower than mass, high temperature zone mostly remains in center, but mass is transferred upwards.

In line 89, Boussinesq force term is calculated:

Fi[:, :, i] = cxs9[i]*g[0]*beta*(T-Tref) + cys9[i]*g[1]*beta*(T-Tref) # Boussinesq

In collision step: new density is calculated by F += -(1.0/tau_f) * (F – Feq) + Fi but new energy G += -(1.0/tau_g) * (G – Geq), so should there be no Boussinesq term in this equation? in Sharma et al. paper above such term is only at mass but not energy calculation in collision step

Maybe somebody got any insights? where could problem be?

Code and original post here: https://www.reddit.com/r/CFD/comments/wtd1te/thermal_lbm_convection_instability