I Plugged in Pressure Directly and a Still Interface Started to Shiver — Two Forms of Nonideal LBM Forcing

A still droplet started to flow on its own#

I put the van der Waals equation of state on top of lattice Boltzmann (LBM, Lattice Boltzmann Method) and solved a two-phase fluid. The initial condition was a single droplet at rest. The density field is smooth, and the velocity is zero everywhere.

After a few steps, small velocities sprouted up near the interface. Nobody pushed it, yet the fluid flows. This artificial velocity is called a parasitic current (a ghost flow that appears at the interface with no physical cause).

Narrowing down the cause led me to one spot in the code. It comes down to what you put into the force term. Do you plug in the pressure $p$ directly, or do you plug in the chemical potential $\mu$? The textbooks say the two are the same. This post writes out, as a ledger, how far that "same" holds true. The answer is one line — they agree in the continuum, and they split on the lattice.

Drop temperature below 1 and the isotherm folds. One pressure then maps to three densities.

The horizontal line is the Maxwell pressure, where the two green lobes have equal area. The gap from the blue dot to the pink dot is the density jump the force term will have to hold.

Where does the force enter#

LBM advects and collides distribution functions to recover the macroscopic equations. For an ideal gas the pressure that comes out on its own is just $c_s^2 \rho$ ($c_s$ is the lattice speed of sound). The true pressure of a nonideal fluid like van der Waals differs from that. Filling that difference is the job of the force term $\mathbf{F}$.

The goal is to revive the following momentum equation.

$p$ is the van der Waals pressure, $\boldsymbol{\tau}$ is the viscous stress, and the last term is the Korteweg stress that builds the interface, with $\kappa$ its strength. Since streaming supplies $c_s^2\rho$, the force only needs to fill the remaining difference.

Here two branches appear. A form that uses the pressure directly, and a form that uses the chemical potential derived from the free energy.

The first is the pressure form (today's subject, "plug in $p$ directly"), the second is the free-energy form. The dial above already folded the isotherm. What remains is why the two forces agree in the continuum.

Gibbs-Duhem — pressure and chemical potential say the same thing twice#

Whether the two forms are the same is settled by a single relation. It is the Gibbs–Duhem relation that ties pressure and chemical potential together at constant temperature.

Substitute this into $\mathbf{F}_\mu$. Since $-\rho\nabla\mu = -\nabla p$, the chemical-potential term turns straight into a pressure gradient. The remaining $c_s^2\nabla\rho$ meshes exactly with the term the pressure form carried as $-\nabla(-c_s^2\rho)$. In the end $\mathbf{F}_p = \mathbf{F}_\mu$.

So the only ground for the claim "you can just plug in $p$ directly" is Gibbs–Duhem. Van der Waals satisfies this relation exactly, because the equation of state and the free energy come from the same thermodynamics. It is the same kind of equivalence as the LBM forcing schemes confirmed by Chapman–Enskog expansion — each wearing a different face, yet all recovering the same macroscopic equation.

Coexistence densities and Gibbs-Duhem, verified in Python#

Words alone are not trustworthy. Let's write the van der Waals equation of state in reduced units, find the coexistence densities per temperature with Newton's method, then directly measure the Gibbs–Duhem defect $\lvert \mathrm{d}p/\mathrm{d}\rho - \rho\,\mathrm{d}\mu/\mathrm{d}\rho \rvert$.

import numpy as np
 
A, B, R = 9.0 / 8.0, 1.0 / 3.0, 1.0   # van der Waals, reduced units (rho_c=1, T_c=1)
 
def p_eos(rho, T):                     # van der Waals pressure
    return rho * R * T / (1.0 - B * rho) - A * rho * rho
 
def mu_eos(rho, T):                    # chemical potential mu = df/drho
    return R * T * (np.log(rho / (1.0 - B * rho)) + B * rho / (1.0 - B * rho)) - 2.0 * A * rho
 
def maxwell(T):                        # equal-area rule: densities where p and mu match in both phases
    x = np.array([0.30, 1.90]); h = 1e-8
    for _ in range(80):
        f = np.array([p_eos(x[0], T) - p_eos(x[1], T), mu_eos(x[0], T) - mu_eos(x[1], T)])
        J = np.empty((2, 2))
        for k in range(2):
            y = x.copy(); y[k] += h
            g = np.array([p_eos(y[0], T) - p_eos(y[1], T), mu_eos(y[0], T) - mu_eos(y[1], T)])
            J[:, k] = (g - f) / h
        x -= np.linalg.solve(J, f)
    return x[0], x[1]
 
# Gibbs-Duhem: dp = rho d(mu). The only ground for plugging in pressure directly.
print("T/Tc  rho_vap  rho_liq |  max| dp/drho - rho*dmu/drho |")
for T in (0.95, 0.90, 0.85):
    rv, rl = maxwell(T)
    r = np.linspace(rv, rl, 400)
    dp = np.gradient(p_eos(r, T), r)
    dmu = np.gradient(mu_eos(r, T), r)
    err = np.abs(dp - r * dmu).max()
    print(f"{T:4.2f}  {rv:7.4f}  {rl:7.4f} |  {err:.2e}")

T/Tc  rho_vap  rho_liq |  max| dp/drho - rho*dmu/drho |
0.95   0.5790   1.4617 |  2.94e-04
0.90   0.4257   1.6573 |  9.44e-04
0.85   0.3197   1.8071 |  1.98e-03

The defect is at the $10^{-3}$ level, and even that is because the derivatives were measured by finite differences. The lower the temperature, the wider the gap between the coexistence densities. Gibbs–Duhem holds exactly in the continuum. Up to here the pressure form and the free-energy form are completely identical.

Read the coexistence columns as subtractions. At $T/T_c=0.95$ the vapor is $0.5790$ and the liquid $1.4617$, a gap of $0.8827$. At $0.90$ the gap is $1.6573-0.4257=1.2316$. At $0.85$ it is $1.8071-0.3197=1.4874$. Dropping temperature by $0.10$ stretches the density jump from $0.8827$ to $1.4874$, about $1.7$ times. Over the same span the defect grows through $2.94\times10^{-4}$, $9.44\times10^{-4}$, and $1.98\times10^{-3}$. The identity $\mathrm{d}p=\rho\,\mathrm{d}\mu$ did not break. A 400-point np.gradient just spends more error on a steeper isotherm. At $T=0.90$ the liquid-to-vapor ratio is $1.6573/0.4257\approx 3.89$. That ratio is the density jump the forcing has to hold.

Stack the two force lines and only one slot disagrees. $-\nabla(p-c_s^2\rho)$ in $\mathbf{F}_p$ and $-\rho\nabla\mu+c_s^2\nabla\rho$ in $\mathbf{F}_\mu$ match once $\mathrm{d}p=\rho\,\mathrm{d}\mu$ holds. The Korteweg term $\kappa\rho\nabla(\nabla^2\rho)$ is written the same way on both sides, so it cannot open a gap. $c_s^2=1/3$ is the lattice-sound correction already booked, not a physical push on the interface. In the continuum ledger the two forms agree, and that sentence ends here. On the lattice the sentence fails by the $1.60\times10^{-2}$ printed later.

On the discrete lattice the two split apart#

The problem is the lattice. There is no guarantee that the continuum $\mathrm{d}p = \rho\,\mathrm{d}\mu$ holds under discrete differentiation. A center-differenced $\nabla p$ and $\rho\,\nabla\mu$ diverge from each other where the density bends sharply, as at an interface.

I solved a still planar interface with the Euler–Lagrange condition, then computed the two force forms directly on top of it.

import numpy as np
 
A, B, R = 9.0 / 8.0, 1.0 / 3.0, 1.0
CS2, KAPPA, T = 1.0 / 3.0, 0.02, 0.90
 
def p_eos(rho):  return rho * R * T / (1.0 - B * rho) - A * rho * rho
def mu_eos(rho): return R * T * (np.log(rho / (1.0 - B * rho)) + B * rho / (1.0 - B * rho)) - 2.0 * A * rho
def dmu(rho):    return R * T * (1.0 / (rho * (1.0 - B * rho)) + B / (1.0 - B * rho) ** 2) - 2.0 * A
def diff1(a, dx): return (np.roll(a, -1) - np.roll(a, 1)) / (2.0 * dx)
def lap(a, dx):   return (np.roll(a, -1) - 2.0 * a + np.roll(a, 1)) / dx ** 2
 
def maxwell():
    x = np.array([0.30, 1.90]); h = 1e-8
    for _ in range(80):
        f = np.array([p_eos(x[0]) - p_eos(x[1]), mu_eos(x[0]) - mu_eos(x[1])])
        J = np.empty((2, 2))
        for k in range(2):
            y = x.copy(); y[k] += h
            g = np.array([p_eos(y[0]) - p_eos(y[1]), mu_eos(y[0]) - mu_eos(y[1])])
            J[:, k] = (g - f) / h
        x -= np.linalg.solve(J, f)
    return x[0], x[1]
 
rv, rl = maxwell()
mu_co = mu_eos(np.array([rv]))[0]
 
# (1) Solve the still planar interface and compare the two forcing forms
NX = 240
xs = np.arange(NX)
rho = 0.5 * (rl + rv) + 0.5 * (rl - rv) * (np.tanh((xs - NX / 4) / 6.0) - np.tanh((xs - 3 * NX / 4) / 6.0) - 1.0)
for _ in range(6000):                                    # relax the Euler-Lagrange residual
    rho -= 0.15 * (mu_eos(rho) - KAPPA * lap(rho, 1.0) - mu_co) / dmu(rho)
 
Gp  = -diff1(p_eos(rho) - CS2 * rho, 1.0) + KAPPA * rho * diff1(lap(rho, 1.0), 1.0) - diff1(CS2 * rho, 1.0)
Gmu = -rho * diff1(mu_eos(rho) - KAPPA * lap(rho, 1.0), 1.0) + CS2 * diff1(rho, 1.0) - diff1(CS2 * rho, 1.0)
gd  = np.abs(diff1(p_eos(rho), 1.0) - rho * diff1(mu_eos(rho), 1.0)).max()
 
print(f"coexistence         rho_vap = {rv:.4f}   rho_liq = {rl:.4f}")
print(f"free-energy form    max|G_mu|        = {np.abs(Gmu).max():.2e}   (well-balanced)")
print(f"pressure form       max|G_p|         = {np.abs(Gp).max():.2e}   (spurious force)")
print(f"gap between forms   max|G_p - G_mu|  = {np.abs(Gp - Gmu).max():.2e}")
print(f"discrete Gibbs-Duhem defect          = {gd:.2e}   <- the gap, exactly")
 
# (2) Refine dx on the same physical interface and the defect converges to 0 fast
print("\ncells/interface |  Gibbs-Duhem defect  order")
prev = None
for n in (10, 20, 40, 80):
    L = 40.0; N = int(L * n / 10)
    z = np.linspace(-L / 2, L / 2, N, endpoint=False); dx = z[1] - z[0]
    r = 0.5 * (rl + rv) - 0.5 * (rl - rv) * np.tanh(z / (0.1 * n))
    d = np.abs(diff1(p_eos(r), dx) - r * diff1(mu_eos(r), dx))[N // 4:3 * N // 4].max()
    order = "" if prev is None else f"{np.log(prev / d) / np.log(2.0):5.2f}"
    print(f"{n:9d}       |  {d:.3e}          {order}")
    prev = d

coexistence         rho_vap = 0.4257   rho_liq = 1.6573
free-energy form    max|G_mu|        = 2.02e-16   (well-balanced)
pressure form       max|G_p|         = 1.60e-02   (spurious force)
gap between forms   max|G_p - G_mu|  = 1.60e-02
discrete Gibbs-Duhem defect          = 1.60e-02   <- the gap, exactly
 
cells/interface |  Gibbs-Duhem defect  order
       10       |  1.140e-02          
       20       |  1.428e-03           3.00
       40       |  5.115e-05           4.80
       80       |  1.611e-06           4.99

Raise resolution and the interface spans more cells. The red peak drops with that cell count.

The red peak is the extra force the pressure form piles on the free-energy form, $G_p - G_\mu$. The green dashed line is zero. More cells, and the peak sits down. The ghost force is a leftover of that difference, not of the physics.

The free-energy form sits at $2.02\times10^{-16}$ on the interface. The pressure form leaves $1.60\times10^{-2}$. The gap matches the discrete Gibbs–Duhem defect to the printed decimal. $1.60\times10^{-2}/2.02\times10^{-16}\approx 7.9\times10^{13}$. Machine roundoff and the ghost force share one line of the table.

The grid table is the same ledger. Ten cells per interface give a defect of $1.140\times10^{-2}$; twenty cells give $1.428\times10^{-3}$, observed order $3.00$. Forty cells give $5.115\times10^{-5}$ (order $4.80$), eighty give $1.611\times10^{-6}$ (order $4.99$). Tightening from 10 to 80 cells, a factor of eight, cuts the defect by $1.140\times10^{-2}/1.611\times10^{-6}\approx 7.1\times10^{3}$. The amount by which centered differences miss $\mathrm{d}p=\rho\,\mathrm{d}\mu$ on the interface slope is the gap between the two forces. That force is what shoves the still interface.

Leaving the pressure form on a thin interface sends the bill#

The free-energy form ($-\rho\nabla\mu$) is balanced at the interface by definition. Even on a coarse lattice the ghost force sits at $2.02\times10^{-16}$. Keeping pressure means the discrete $\nabla p$ has to be the same difference as $\rho\,\nabla\mu$. Otherwise discrete Gibbs–Duhem dies and the table's $1.60\times10^{-2}$ stays.

Resolve the interface with 4 or 5 cells or more and the pressure-form defect falls to order $4.99$. Practical interfaces are about 3 cells thick. There the ghost force tracks the square of the density jump. A force term that matches the continuum identity and still misses the discrete moments leaves the same residue; that ledger is the third-moment limit on the force term. The convenience of "using $p$ directly" is billed in interface thickness.