把压力直接代入后静止的界面开始颤动 — 非理想LBM forcing的两种形式

压力形与自由能形的forcing凭吉布斯-杜安关系在连续下相等。但在格子上，压力形会在界面留下幽灵般的力。

静止的液滴自己开始流动了

把范德瓦尔斯状态方程加到格子玻尔兹曼(LBM, Lattice Boltzmann Method)上求解两相流体。初始条件是一个静止的液滴。密度场平滑，速度全部为0。

跑了几步后，界面附近冒出了微小的速度。没人推动，流体却流动起来。这种人为速度被称为寄生流(parasitic current，没有物理原因却在界面产生的幽灵般的流动)。

把原因逐步缩小，最终落到代码的一处。就是往力项里代入什么。是把压力 $p$ 直接代入，还是代入化学势 $\mu$ 。教科书说两者相同。本文把这个"相同"到底成立到哪一步记入账本。答案只有一行 — 连续下相同，格子上分道扬镳。

把 temperature 降到1以下，等温线就会折起来。一个压力会对应三个密度。

temperature0.900 Tc

Drag temperature down from 1.0: the isotherm folds into an S, and the amber tie line drops in so the two green lobes have equal area. The blue and pink dots are the vapor and liquid densities the LBM must hold apart — their gap is what the forcing term exists to sustain.

横线是两个绿色瓣面积相等的麦克斯韦压力。蓝点到粉点的间距，就是力项稍后必须撑住的密度差。

力从哪里进来

LBM让分布函数经过对流与碰撞，复原宏观方程。若是理想气体，自然浮现的压力只有 $c_s^2 \rho$ 一项( $c_s$ 是格子声速)。范德瓦尔斯这类非理想流体的真实压力与之不同。填补那份差额，正是力项 $\mathbf{F}$ 的职责。

目标是复原下面的动量方程。

\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u}\mathbf{u}) = -\nabla p + \nabla \cdot \boldsymbol{\tau} + \kappa \rho \, \nabla (\nabla^2 \rho)

$p$ 是范德瓦尔斯压力， $\boldsymbol{\tau}$ 是粘性应力，最后一项是竖起界面的科特维格(Korteweg) 应力， $\kappa$ 是其强度。由于流式传递给出 $c_s^2\rho$ ，力只需补上剩余的差额。

这里出现两条分支。直接使用压力的形式，以及使用源自自由能的化学势的形式。

\mathbf{F}_p = -\nabla\!\left(p - c_s^2 \rho\right) + \kappa \rho \, \nabla(\nabla^2 \rho)

\mathbf{F}_\mu = -\rho \, \nabla \mu + c_s^2 \nabla \rho + \kappa \rho \, \nabla(\nabla^2 \rho)

前者是压力形(今天的主题"直接用 $p$ ")，后者是自由能形。等温线在何处折起，前面的旋钮已经降过温度。剩下的只是两股力为何在连续下相同。

吉布斯-杜安 — 压力与化学势在讲同一件事两遍#

两种形式是否相同，由一个关系式一锤定音。这就是在等温下连接压力与化学势的吉布斯-杜安(Gibbs–Duhem) 关系。

\mathrm{d}p = \rho \, \mathrm{d}\mu \qquad (T = \text{const})

把这个式子代入 $\mathbf{F}_\mu$ 看看。由于 $-\rho\nabla\mu = -\nabla p$ ，化学势项立刻变为压力梯度。剩下的 $c_s^2\nabla\rho$ 恰好与压力形以 $-\nabla(-c_s^2\rho)$ 形式持有的那一项精确咬合。最终 $\mathbf{F}_p = \mathbf{F}_\mu$ 。

也就是说，" $p$ 可以直接代入"这一主张的唯一依据就是吉布斯-杜安。范德瓦尔斯精确地满足这个关系。因为状态方程与自由能出自同一套热力学。这与用Chapman–Enskog展开确认过的各种LBM forcing格式即便各有不同面孔，最终仍复原同一宏观方程属于同一类等价。

用Python验证的共存密度与吉布斯-杜安#

光靠嘴说不足信。把范德瓦尔斯状态方程以约化单位写出，用牛顿法求各温度下的共存密度，再直接测量吉布斯-杜安亏损 $\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   # 范德瓦尔斯，约化单位 (rho_c=1, T_c=1)
 
def p_eos(rho, T):                     # 范德瓦尔斯压力
    return rho * R * T / (1.0 - B * rho) - A * rho * rho
 
def mu_eos(rho, T):                    # 化学势 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):                        # 等面积法则: 使 p 与 mu 在两相中相等的密度
    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]
 
# 吉布斯-杜安: dp = rho d(mu). 直接使用压力的唯一依据。
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

亏损在 $10^{-3}$ 量级，而且那还是用有限差分测导数造成的。温度越低，共存密度的间距越大。吉布斯-杜安在连续下精确成立。到这一步为止，压力形与自由能形完全相同。

把共存密度逐列相减。 $T/T_c=0.95$ 时气体 $0.5790$ 、液体 $1.4617$ ，差 $0.8827$ 。 $0.90$ 是 $1.6573-0.4257=1.2316$ ， $0.85$ 是 $1.8071-0.3197=1.4874$ 。温度降 $0.10$ ，密度差从 $0.8827$ 拉到 $1.4874$ ，大约 $1.7$ 倍。同一段里亏损依次是 $2.94\times10^{-4}$ 、 $9.44\times10^{-4}$ 、 $1.98\times10^{-3}$ 。不是恒等式 $\mathrm{d}p=\rho\,\mathrm{d}\mu$ 破了。400 点的 np.gradient 在更陡的等温线上多付了差分误差。 $T=0.90$ 的液气比是 $1.6573/0.4257\approx 3.89$ 。这个比就是 forcing 要撑住的密度跳跃。

两行力叠在一起，真正分开的位置只有一处。 $\mathbf{F}_p$ 里的 $-\nabla(p-c_s^2\rho)$ 和 $\mathbf{F}_\mu$ 里的 $-\rho\nabla\mu+c_s^2\nabla\rho$ ，在 $\mathrm{d}p=\rho\,\mathrm{d}\mu$ 下是同一项。科特维格项 $\kappa\rho\nabla(\nabla^2\rho)$ 两边写法相同，造不出差。 $c_s^2=1/3$ 是格子声速已经记上的修正，不是推界面的物理力。连续账本上两种形式相同，这句话到此为止。格子上这句话裂开的大小，是后面表里的 $1.60\times10^{-2}$ 。

在离散格子上两者分道扬镳

问题出在格子。连续的 $\mathrm{d}p = \rho\,\mathrm{d}\mu$ 在离散微分下未必成立。中心差分得到的 $\nabla p$ 与 $\rho\,\nabla\mu$ ，在界面这样密度急剧弯折的地方彼此错开。

用欧拉-拉格朗日条件求解静止的平面界面，再在其上直接计算两种力的形式。

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) 求解静止平面界面，比较两种forcing形式
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):                                    # 欧拉-拉格朗日残差松弛
    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) 同一物理界面只把dx加密，亏损会迅速收敛到0
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

把 resolution 调高，界面跨越的格子变多，红色峰的高度一起掉下去。

resolution (cells / interface)2.5

The red hump is the force the pressure form adds over the free-energy form — the discrete Gibbs–Duhem defect. It lives exactly on the interface and pushes a fluid that should be at rest. Drag resolution up: the hump collapses onto the green zero line. It was never physics, only the price of writing p on too few cells.

红色峰是压力形叠在自由能形上的额外力，也就是 $G_p - G_\mu$ 。绿色虚线是0。格子一多，峰就塌下去。幽灵般的力不是物理，是这个差分留下的。

自由能形在界面上是 $2.02\times10^{-16}$ 。压力形留下 $1.60\times10^{-2}$ 。两者之差与离散吉布斯-杜安亏损小数点后一致。 $1.60\times10^{-2}/2.02\times10^{-16}\approx 7.9\times10^{13}$ 。机器精度和幽灵力的量级写在同一行。

格子表是同一本账。界面每10个格子，亏损 $1.140\times10^{-2}$ ；20个格子是 $1.428\times10^{-3}$ ，观测阶 $3.00$ 。40个格子是 $5.115\times10^{-5}$ （阶 $4.80$ ），80个格子是 $1.611\times10^{-6}$ （阶 $4.99$ ）。从10格收到80格，只加密八倍，亏损就降到 $1.140\times10^{-2}/1.611\times10^{-6}\approx 7.1\times10^{3}$ 分之一。中心差分在界面斜率上没守住 $\mathrm{d}p=\rho\,\mathrm{d}\mu$ 的那一截，就是两种力的差。这股力推着静止界面。

留下压力形，薄界面会来结账

自由能形( $-\rho\nabla\mu$ )在定义上于界面处平衡。即便粗格子，幽灵般的力也贴在 $2.02\times10^{-16}$ 。若仍直接用压力， $\nabla p$ 不能随意差分，必须和 $\rho\,\nabla\mu$ 用同一种差分。离散吉布斯-杜安才活得下去，表里的 $1.60\times10^{-2}$ 才会降下来。

界面解析到4到5个格子以上，压力形亏损会降到阶 $4.99$ 。实务里的界面大约3个格子厚。那个厚度上，幽灵力跟着密度差的平方留下来。力项只满足连续恒等式、离散矩却空着，就会留下同类残差。那本账在力项的三阶矩极限。"直接用 $p$ "的代价，用界面厚度来结。