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

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

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

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

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

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

力从哪里进来

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$ ")，后者是自由能形。要看两者是否真的相同，先要了解范德瓦尔斯状态方程的样貌。下面就直接把温度降下来看看。

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.

把 temperature 降到1.0以下，等温线便折成S形。一个压力对应三个密度。此时物理上的两相压力由两个绿色瓣面积相等的位置(麦克斯韦等面积法则)决定。蓝点与粉点的间距，就是forcing项必须支撑的密度差。

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

两种形式是否相同，由一个关系式一锤定音。这就是在等温下连接压力与化学势的吉布斯-杜安(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}$ 量级，而且那还是用有限差分测导数造成的。温度越低，共存密度的间距越大。吉布斯-杜安在连续下精确成立。到这一步为止，压力形与自由能形完全相同。

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

问题出在格子。连续的 $\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

三行是核心。自由能形在界面上力为 $10^{-16}$ ，实际上就是0。压力形则留下 $1.6\times10^{-2}$ 的力。而且两种形式的差与离散吉布斯-杜安亏损小数点后都精确一致。压力形漏出的幽灵般的力，其真身正是这个亏损。这个力推动静止的界面，制造出寄生流。

下面来改变界面被铺开的宽度(格子分辨率)看看。

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.

调高 resolution，让界面跨越更多的单元格，红色压力形曲线的峰便会塌向0。绿色自由能形从头到尾都贴着0。幽灵般的力并非物理，而是离散化的副产物。

那么该代入什么

归纳起来选择有两个。第一，使用自由能形( $-\rho\nabla\mu$ )。这种形式在定义上于界面处保持平衡，即便在粗格子上也没有幽灵般的力。第二，若真想直接使用压力，就不要随意离散化 $\nabla p$ ，而要写成与 $\rho\,\nabla\mu$ 一致。这样离散下吉布斯-杜安也成立，平衡得以存活。

若能把界面充分解析到4~5个单元格以上，压力形的误差就会像上表那样迅速消失。但实务中的界面通常薄到3个单元格上下。在那个区域，压力形会产生正比于密度差平方的幽灵般的力。同样的症状，我曾在寄生流与well-balanced界面张力中从表面张力一侧看过。根源只有一个 — 连续下平衡的项，在离散下是否也搬成了平衡。

"直接使用 $p$ "这份便利并非免费。其代价会以界面厚度这种货币来结算。