Giving Conduction Its Own τ Left Viscous Heating Stuck at Pr = 1 — The Quasi-Equilibrium g* in Two-Population LB

Pr changed, so why didn't the temperature?#

Picture a Couette channel. The bottom wall is at rest and the top wall slides at speed $U$. Both walls have the same temperature. The middle still heats up, because viscosity turns kinetic energy into heat (viscous heating). The peak temperature rise is exactly proportional to the Prandtl number (Pr, the ratio of momentum diffusion to heat diffusion).

This post is about changing Pr in a lattice Boltzmann (LB) thermal model. Split the distribution into two sets, give each its own relaxation time, and Pr looks easy to tune. Run it, though, and conduction changes while viscous heating stays frozen at its Pr = 1 value. You'll see why, and you'll get the one-line fix: the quasi-equilibrium $g^*$.

Move the Pr slider and the gap between the stress clock and the conduction clock changes, along with which curve the viscous-work flux sticks to.

At Pr = 0.71 the orange curve (stress) and the dashed green curve (viscous-work flux in the g* model) stay on top of each other. The red curve (plain two-τ model) follows the purple conduction clock and drifts away from the stress. Slide Pr to 1 and the four curves collapse into two.

One set of populations only knows Pr = 1#

The BGK collision operator pulls every non-equilibrium moment toward equilibrium at one rate $\omega$. Stress and heat flux relax on the same clock, so viscosity and conductivity are tied to the same $\omega$. That is why the continuous BGK equation gives Pr = 1.

Air has Pr = 0.71 and water about 7. A model locked to Pr = 1 can't even get the ratio of boundary-layer thicknesses right. There are two remedies. One relaxes only the heat-flux moment at a different rate inside the collision operator (MRT, Shakhov-type). The other hands energy to a separate distribution. This post follows the second route.

Two sets: f carries mass and momentum, g carries total energy#

In a two-population model, $f_i$ carries density and momentum, and $g_i$ carries the total energy $E=\rho(c_vT+u^2/2)$ ($i=0,\dots,Q-1$). The kinetic equations are:

$\omega_\mu$ sets viscosity (kinematic viscosity $\nu$), and $\omega_\kappa$ sets conduction (thermal diffusivity $\alpha$). $g_i^{eq}$ is the local equilibrium and $g_i^*$ is the quasi-equilibrium.

The two rates map to the diffusivities as:

$c_s^2=1/3$ is the squared lattice sound speed. With $\tau_\mu=1/\omega_\mu$ and $\tau_\kappa=1/\omega_\kappa$, $\mathrm{Pr}=(\tau_\mu-\tfrac12)/(\tau_\kappa-\tfrac12)$.

The code locks $\tau_\mu=1$. Then $\tau_\mu-\tfrac12=\tfrac12$, so $\tau_\kappa=\tfrac12+(\tfrac12)/\mathrm{Pr}$. The four table entries fall out at once. Pr = 0.5 gives $\tau_\kappa=1.5$. Pr = 0.71 adds $0.5/0.71\approx0.704$, so $\tau_\kappa\approx1.204$. Pr = 1 gives $\tau_\kappa=1$. Pr = 2 gives $\tau_\kappa=0.75$. Since $\omega=1/\tau$, heat at Pr = 2 relaxes at $\omega_\kappa\approx1.333$ while viscosity sits at $\omega_\mu=1$. Viscous work belongs on the slow clock. Plain two-τ puts that work on the fast $\omega_\kappa$. Only the Pr = 1 row has both rates equal, and both columns read 1.000. Pass only that row and the bug stays quiet.

Viscous work hides inside the energy flux#

In this toy, $g^{eq}$ fixes only three moments. Zeroth: $E$. First (energy flux): $(E+p)\mathbf{u}$. Second: $c_s^2E\,\mathbf{I}+E\,\mathbf{u}\mathbf{u}$. The non-equilibrium part of the first moment, $\mathbf{Q}$, is the actual heat flux.

At first order in the Chapman–Enskog expansion (solving for non-equilibrium at small Knudsen number), $\mathbf{Q}\approx-\tau\,\nabla\!\cdot\!\mathbf{R}^{eq}$. Here $\nabla(c_s^2E)$ splits into two pieces:

The first term is Fourier conduction. The second is the raw material for the viscous work $-\mathbf{u}\cdot\boldsymbol{\sigma}$. It has exactly the shape of the non-equilibrium stress in $f$, $\boldsymbol{\Pi}^{neq}=\mathbf{P}-\mathbf{P}^{eq}\approx-\tau_\mu\rho c_s^2(\nabla\mathbf{u}+\nabla\mathbf{u}^T)$. The problem is the $\tau$ in front.

If $g$ is pulled only toward equilibrium at $\omega_\kappa$, both pieces come out carrying $\tau_\kappa$. The viscous-work coefficient scales with $\alpha$ instead of $\nu$. Conduction is right, but viscous heating stays at its Pr = 1 value. That's the red curve in the first viz.

Quasi-equilibrium g*: tie only the flux to the stress clock#

The fix is one line in the definition of $g^*$. It shares every moment with $g^{eq}$ except the energy flux:

$\mathbf{P}$ is the actual momentum flux measured from $f$, and $\mathbf{P}^{eq}$ is its equilibrium value.

Taking the first moment of the kinetic equation gives a one-step update for the two non-equilibrium moments:

$\mathbf{Q}$ itself relaxes at $\omega_\kappa$, so conduction follows $\alpha$. At the same time, if $\mathbf{Q}=\mathbf{u}\cdot\boldsymbol{\Pi}^{neq}$, then $\mathbf{Q}'=\mathbf{u}\cdot\boldsymbol{\Pi}'^{neq}$ holds at the next step. The viscous-work piece rides the stress clock $\omega_\mu$. A Chapman–Enskog calculation says the same thing: the viscous work carrying $\tau_\kappa$ cancels against the $(\omega_\kappa-\omega_\mu)$ term, leaving only $\tau_\mu$.

With $\omega_\kappa=\omega_\mu$ the $g^*$ term vanishes. The model reduces to single relaxation, and Pr = 1. $g^*$ is a correction that only works when Pr $\ne 1$. A thermal lattice with one population kept Pr stuck at 1 even when τ was raised. That is the same wall as the two constants locked in place. Splitting the populations does not leave that wall unless the viscous-work clock moves with them. The question is still which moment relaxes at which rate.

Code: one Couette channel, two kinds of g#

The analytical solution is simple. The steady energy equation is $kT''+\mu(u')^2=0$. With a linear velocity profile, temperature is a parabola:

The peak in the middle is $\theta_{max}=\mathrm{Pr}$. The code below runs D2Q9 with a single periodic cell in $x$ and solves only along $y$. Walls use Guo's non-equilibrium extrapolation. For simplicity, $f$ uses an isothermal equilibrium and energy is coupled one way. In low-Mach ($U=0.1$) Couette flow the pressure is uniform, so this simplification doesn't change the result.

import numpy as np
 
# D2Q9, one periodic cell in x → 1D arrays along y
cx = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
cy = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1])
w = np.array([4/9] + [1/9]*4 + [1/36]*4)
CV, T0 = 1.0, 1.0
 
def feq_lattice(rho, ux):
    cu = cx[:, None]*ux
    return w[:, None]*rho*(1 + 3*cu + 4.5*cu**2 - 1.5*ux**2)
 
def geq_energy(E, rho, ux):
    p = rho*(E/rho - 0.5*ux**2)/CV          # p = rho*T
    qx = (E + p)*ux
    cu = cx[:, None]*ux
    return w[:, None]*(E + 3*cx[:, None]*qx + 4.5*E*cu**2 - 1.5*E*ux**2)
 
def couette_heating(Pr, quasi, N=21, U=0.1, tau_mu=1.0, steps=8000):
    om_mu = 1/tau_mu
    om_ka = 1/(0.5 + (tau_mu - 0.5)/Pr)
    rho = np.ones(N); ux = np.zeros(N)
    E = rho*CV*T0*np.ones(N)
    f = feq_lattice(rho, ux); g = geq_energy(E, rho, ux)
    uw = np.zeros(N); uw[-1] = U
    for _ in range(steps):
        rho = f.sum(0); ux = (cx[:, None]*f).sum(0)/rho; E = g.sum(0)
        fe = feq_lattice(rho, ux); ge = geq_energy(E, rho, ux)
        Pxy = ((cx*cy)[:, None]*(f - fe)).sum(0)
        f = f + om_mu*(fe - f)
        if quasi:   # g*: quasi-equilibrium with u·(P - P_eq) added to the energy flux
            gs = ge + w[:, None]*3*cy[:, None]*(ux*Pxy)
            g = g + om_mu*(ge - g) + (om_ka - om_mu)*(gs - g)
        else:       # plain two-τ: g relaxes only to equilibrium at ω_κ
            g = g + om_ka*(ge - g)
        for k in range(9):
            f[k] = np.roll(f[k], cy[k]); g[k] = np.roll(g[k], cy[k])
        for wall, nb in ((0, 1), (N-1, N-2)):   # non-equilibrium extrapolation walls
            r = f[:, nb].sum(); un = (cx*f[:, nb]).sum()/r; En = g[:, nb].sum()
            fn = feq_lattice(np.array([r]), np.array([un]))[:, 0]
            gn = geq_energy(np.array([En]), np.array([r]), np.array([un]))[:, 0]
            Ew = r*(CV*T0 + 0.5*uw[wall]**2)
            f[:, wall] = feq_lattice(np.array([r]), np.array([uw[wall]]))[:, 0] + f[:, nb] - fn
            g[:, wall] = geq_energy(np.array([Ew]), np.array([r]), np.array([uw[wall]]))[:, 0] + g[:, nb] - gn
    rho = f.sum(0); ux = (cx[:, None]*f).sum(0)/rho
    T = (g.sum(0)/rho - 0.5*ux**2)/CV
    theta = (T - T0)/(U**2/(8*CV))
    return theta.max(), ux
 
print(" Pr   exact   2-tau(g_eq)   g* quasi-eq")
for Pr in (0.5, 0.71, 1.0, 2.0):
    a, _ = couette_heating(Pr, False)
    b, u = couette_heating(Pr, True)
    print(f"{Pr:4.2f}  {Pr:6.3f}   {a:9.3f}    {b:9.3f}")
print(f"u(mid)/U = {u[10]/0.1:.4f}")

 Pr   exact   2-tau(g_eq)   g* quasi-eq
0.50   0.500       1.000        0.498
0.71   0.710       1.000        0.710
1.00   1.000       1.000        1.000
2.00   2.000       1.000        1.997
u(mid)/U = 0.5000

Now run that code in your browser. Both models march forward side by side in the same channel.

Move the Pr slider: the green peak follows the dashed analytical value Pr, while the red one sticks to 1.000. Changing $U$ or $\tau_\mu$ doesn't move the red $\theta_{max}$ off 1 either. Try Pr = 3, which isn't in the output block.

What the 1.000 is telling you#

The plain two-τ column reads 1.000 regardless of Pr. That isn't noise, it's structure. The viscous-work coefficient scales with $\tau_\kappa-\tfrac12$, and so does the conductivity. In their ratio, $\theta_{max}$, $\tau_\kappa$ cancels exactly. For air (Pr = 0.71) that overpredicts the temperature rise between the walls by 41%. For a Pr = 7 fluid like water it underpredicts by a factor of 7.

With $U=0.1$ and $c_v=1$, the denominator is $U^2/(8c_v)=0.00125$. Air's true $\Delta T$ is $0.71\times0.00125\approx8.9\times10^{-4}$. Plain two-τ gives $1.25\times10^{-3}$. On the lattice that is 1.0009 versus 1.0013. Two decimal places print both as 1.00. At Pr = 7 the true rise is $0.00875$ and the model gives $0.00125$. 1.009 versus 1.001. Only $\theta$ splits 1 from 7.

The $g^*$ column reads 0.498, 0.710, 1.000, 1.997. 0.498 is 0.4% under 0.500, and 1.997 is 0.15% under 2. 0.710 matches exact. On 21 nodes the rest is the wall and the grid. The midline velocity $u/U=0.5000$ matches because both models share $f$. The temperature error never shows up in the velocity field.

Go compressible and f gets a quasi-equilibrium too#

The $f$ in this post used an isothermal equilibrium. In compressible two-population models, $f$ also gets a quasi-equilibrium $f^*$. D2Q9's third moments have a defect: $c_{ix}^3=c_{ix}$. When temperature varies or Mach number grows, that defect puts an error term into the stress. $f^*$ carries a pressure tensor that corrects it. $f^*$ repairs the momentum flux; $g^*$ attaches viscous work to the energy flux.

The principle is the same. Equilibrium and quasi-equilibrium share the conserved quantities (density, momentum, total energy). The only difference is the one moment that needs its own relaxation rate. Where to put the pressure term for a non-ideal fluid is the same choice of moment.

The one line to check when your thermal LB code sets τκ separately#

Open the collision step for the energy distribution and look at where $\omega_\kappa$ multiplies. If all of $g$ is pulled toward equilibrium at $\omega_\kappa$, viscous heating is stuck at its Pr = 1 value. One Couette run is enough to check. If $\theta_{max}$ doesn't track Pr, the quasi-equilibrium term is missing.

Validate only at $\omega_\kappa=\omega_\mu$ and the bug stays quiet. A thermal LB check needs one viscous-heating case with Pr $\ne 1$. Temperature diffusing in a fluid at rest never turns this term on. Reynolds' analogy between friction and heat starts from the premise that those two stay in a ratio tied to Pr. If $\theta_{max}$ does not follow Pr, that premise is already broken on the lattice.