去掉K之后水被压缩了1.86倍 — 加密网格也不会减小的界面误差

界面单元里ψ+K几乎等于零。K正是用来抹去ψ∇·u的系数，去掉这一项，液体就会按气体的压缩率被压。

可压缩VOF的颜色函数输运方程多了一项： $-(\psi+K)\nabla\cdot u$ 。我把这个 $K$ 写进代码并打印了它的数值。水与空气的界面单元在 $\psi=0.5$ 时得到 $K=-0.499947$ ，于是 $\psi+K=5.3\times10^{-5}$ 。为了让这一项正确而加进去的系数，几乎把这一项本身抹掉了。

本文从Denner的压力基全耦合可压缩多相有限体积法论文出发，用等压封闭直接推导 $K$ 。随后复现把 $K$ 置零或把符号反过来时，一维压缩问题中液体体积分别多减1.86倍和4.23倍的结果。最后看看为什么把网格加密8倍，这个误差只减少7%。

可压缩VOF为什么带着 $-(\psi+K)\nabla\cdot u$ #

不可压缩VOF中颜色函数只是随流而动。因为 $\nabla\cdot u=0$ ，只做输运就够了。一旦考虑可压缩性，流体会膨胀和收缩。颜色函数是体积分数，必须跟上这种膨胀与收缩。

论文中的式子是这样：

\frac{\partial \psi}{\partial t} + \frac{\partial u_i \psi}{\partial x_i} - (\psi + K)\frac{\partial u_i}{\partial x_i} = 0

其中 $\psi$ 是液体体积分数， $u_i$ 是速度， $K$ 是随物性变化的压缩率因子（compressibility factor）。把第三项移过去就得到物质导数形式。

\frac{\mathrm{D}\psi}{\mathrm{D}t} = K\,\nabla\cdot u

若 $K=0$ ，体积分数沿物质路径就不变。看起来合理，其实是错的。界面单元看到的 $\nabla\cdot u$ 是混合物的整体值，而混合物的压缩率由软的一侧主导。

用数字说话。1个大气压下水的体积模量为 $\rho_l a_l^2 = 2.64\times10^9$ Pa，空气为 $\rho_g a_g^2 = 1.40\times10^5$ Pa，比值18,857。混合物压缩率是按体积分数加权的和，所以即使 $\psi=0.5$ ，空气那一项也几乎占满全部。结果混合物声速跌到23.65 m/s，比水慢69倍，比空气慢14倍。

在下面的模拟中直接拖动体积分数试试。

ψ (liquid volume fraction) = 0.500pressure = 1 bar

a_mix = 23.65 m/s  |  ρ_mix = 500.60 kg/m³  |  a_mix / a_gas = 0.069
crossing 2.0 m: liquid 1.23 ms  |  gas 5.85 ms  |  mixture 84.57 ms

把 $\psi$ 滑块从零往上抬一点，声速就从341 m/s垂直跌到100 m/s以下。下方的跑道显示声波脉冲穿过同一段2 m介质所需的时间。液体只占0.1%体积，混合物那条跑道就已明显落后于纯气体。

$K$ 是从哪里来的数#

$K$ 来自一个假设：等压封闭（isobaric closure）——界面单元内两相共享同一压力，各自以自己的体积模量响应。这与论文在ACID中使用的局部单相假设是同一回事。

跟随一个拉格朗日单元，要求相 $k$ 的质量守恒，可以写出

\frac{1}{V_k}\frac{\mathrm{D}V_k}{\mathrm{D}t} = -\frac{1}{\rho_k a_k^2}\frac{\mathrm{D}p}{\mathrm{D}t}

$V_k$ 是该相占据的体积， $\rho_k a_k^2$ 是它的体积模量。两相相加就得到混合物的散度。

\nabla\cdot u = -\left[\frac{\psi}{\rho_l a_l^2} + \frac{1-\psi}{\rho_g a_g^2}\right]\frac{\mathrm{D}p}{\mathrm{D}t}

方括号里就是混合物压缩率，Wood声速由此自动浮现。接着对 $\psi = V_l/V$ 求导并代入上面两式，剩下

\frac{\mathrm{D}\psi}{\mathrm{D}t} = \psi(1-\psi)\left[\frac{1}{\rho_g a_g^2} - \frac{1}{\rho_l a_l^2}\right]\frac{\mathrm{D}p}{\mathrm{D}t}

乘在前面的是两个压缩率之差。两相体积模量相同时，体积分数根本不动。与 $\mathrm{D}\psi/\mathrm{D}t = K\nabla\cdot u$ 对比并消去 $\mathrm{D}p/\mathrm{D}t$ ， $K$ 就掉出来了。

K = \frac{\rho_g a_g^2 - \rho_l a_l^2}{\dfrac{\rho_g a_g^2}{1-\psi} + \dfrac{\rho_l a_l^2}{\psi}}

分母里有 $1/\psi$ 和 $1/(1-\psi)$ ，在纯相单元会发散，所以代码里应该写约分后的形式。

K = \frac{(\rho_g a_g^2 - \rho_l a_l^2)\,\psi(1-\psi)}{\rho_g a_g^2\,\psi + \rho_l a_l^2\,(1-\psi)}

这个形式在 $\psi=0$ 和 $\psi=1$ 处直接为零。分子的符号完全取决于 $\psi$ 指的是哪一相。照抄论文公式时标号一换，符号就可能反掉。后面会测出这个代价。

如果 $\psi+K$ 是 $5\times10^{-5}$ ，两项能一起删掉吗#

取 $\rho_l a_l^2 \gg \rho_g a_g^2$ 的极限。分子趋于 $-\rho_l a_l^2$ ，分母趋于 $\rho_l a_l^2(1-\psi)$ ，于是

K \to -\psi, \qquad \psi + K \to 0

物理上这完全正确。水几乎不可压缩，所以水的体积应该只被输运。把 $(\psi+K)\nabla\cdot u$ 压到零，正是这种行为。

那么干脆两项都不要？纯液体单元是另一回事。 $\psi=1$ 且 $K=0$ ，于是 $(\psi+K)\nabla\cdot u = \nabla\cdot u$ 。液体也应该按自己的压缩率被压缩，而这一项承担的就是这件事。抵消是否彻底也取决于体积模量之比。比值在 $10$ 附近时（超临界流体之间，或气体对气体）， $\psi=0.5$ 处 $\psi+K=0.0909$ ，并不是零。

从1 bar压到10 bar，究竟差多少#

在欧拉网格上直接测。把区间 $[0,1]$ 的左半填水，右半填空气。左端是壁面（ $u=0$ ），右端是允许纯空气流入的开边界。压力在空间上保持均匀，从1 bar线性升到10 bar。采用stiffened-gas状态方程（ $\gamma_l=4.4$ ， $\Pi_l=6\times10^8$ Pa），面速度由上面推导的混合物压缩率积分得到。

水不会到达右边界，所以液体质量守恒。于是精确液体体积直接就是 $\psi_0\,\rho_l(p_0)/\rho_l(p)$ ，参考解白送。

import math
 
P0, PEND = 1.0e5, 1.0e6                    # 1 bar -> 10 bar
GAM_L, PI_L, RHO_L0 = 4.4, 6.0e8, 1000.0   # 水 (stiffened gas)
GAM_G, PI_G, RHO_G0 = 1.4, 0.0,   1.2      # 空气
 
def stiffened_rho(p, gam, pi, rho0):
    return rho0 * ((p + pi) / (P0 + pi)) ** (1.0 / gam)
 
def stiffened_bulk(p, gam, pi):
    return gam * (p + pi)                  # rho*a^2
 
def wood_speed(psi, p):
    rho_m = psi * stiffened_rho(p, GAM_L, PI_L, RHO_L0) \
          + (1 - psi) * stiffened_rho(p, GAM_G, PI_G, RHO_G0)
    comp  = psi / stiffened_bulk(p, GAM_L, PI_L) \
          + (1 - psi) / stiffened_bulk(p, GAM_G, PI_G)
    return math.sqrt(1.0 / (rho_m * comp)), rho_m
 
def compressibility_factor(psi, p, mode):
    bl = stiffened_bulk(p, GAM_L, PI_L)
    bg = stiffened_bulk(p, GAM_G, PI_G)
    k = (bg - bl) * psi * (1 - psi) / (bg * psi + bl * (1 - psi))
    return 0.0 if mode == 'zero' else (-k if mode == 'flip' else k)
 
def squeeze_column(mode, ncell, nstep, psi0=0.5):
    """左壁 u=0，右侧纯气体流入，p(t) 从 1 到 10 bar 线性上升"""
    dx, dt = 1.0 / ncell, 1.0 / nstep
    psi = [1.0 if (i + 0.5) * dx < psi0 else 0.0 for i in range(ncell)]
    for n in range(nstep):
        p    = P0 + (PEND - P0) * (n * dt)
        pdot = PEND - P0
        bl   = stiffened_bulk(p, GAM_L, PI_L)
        bg   = stiffened_bulk(p, GAM_G, PI_G)
        u = [0.0] * (ncell + 1)            # 混合物压缩率 -> du/dx
        for i in range(ncell):
            u[i + 1] = u[i] - dx * pdot * (psi[i] / bl + (1 - psi[i]) / bg)
        nxt = psi[:]
        for i in range(ncell):
            ur, ul = u[i + 1], u[i]
            right  = psi[i + 1] if i + 1 < ncell else 0.0
            fr = ur * (psi[i] if ur > 0 else right)
            fl = (ul * (psi[i - 1] if ul > 0 else psi[i])) if i > 0 else 0.0
            src = (psi[i] + compressibility_factor(psi[i], p, mode)) * (ur - ul) / dx
            nxt[i] = min(1.0, max(0.0, psi[i] - dt * (fr - fl) / dx + dt * src))
        psi = nxt
    exact = psi0 * stiffened_rho(P0, GAM_L, PI_L, RHO_L0) \
                 / stiffened_rho(PEND, GAM_L, PI_L, RHO_L0)
    return sum(psi) * dx, exact
 
print("[1] stiffened-gas properties")
for p in (P0, PEND):
    rl = stiffened_rho(p, GAM_L, PI_L, RHO_L0)
    rg = stiffened_rho(p, GAM_G, PI_G, RHO_G0)
    print(f"  p={p/1e5:4.1f} bar  rho_l={rl:8.3f}  rho_g={rg:6.3f}"
          f"  rho_l*a_l^2={stiffened_bulk(p,GAM_L,PI_L):.3e}"
          f"  rho_g*a_g^2={stiffened_bulk(p,GAM_G,PI_G):.3e}")
 
print("\n[2] K and (psi+K) at 1 bar")
for psi in (0.10, 0.25, 0.50, 0.75, 0.90):
    k = compressibility_factor(psi, P0, 'acid')
    a, _ = wood_speed(psi, P0)
    print(f"  psi={psi:4.2f}  K={k:+9.6f}  psi+K={psi+k:+.3e}  a_wood={a:7.2f} m/s")
 
print("\n[3] liquid-volume change after 1->10 bar compression")
print("  N     exact        ACID         K=0          K flipped    | x exact")
for ncell in (50, 100, 200, 400):
    out = {}
    for mode in ('acid', 'zero', 'flip'):
        v, exact = squeeze_column(mode, ncell, 10 * ncell)
        out[mode] = v - 0.5
    d = exact - 0.5
    print(f"  {ncell:4d}  {d:.4e}  {out['acid']:.4e}  {out['zero']:.4e}  "
          f"{out['flip']:.4e} |  {out['zero']/d:5.3f}  {out['flip']/d:5.3f}")

[1] stiffened-gas properties
  p= 1.0 bar  rho_l=1000.000  rho_g= 1.200  rho_l*a_l^2=2.640e+09  rho_g*a_g^2=1.400e+05
  p=10.0 bar  rho_l=1000.341  rho_g= 6.215  rho_l*a_l^2=2.644e+09  rho_g*a_g^2=1.400e+06
 
[2] K and (psi+K) at 1 bar
  psi=0.10  K=-0.099994  psi+K=+5.891e-06  a_wood=  39.23 m/s
  psi=0.25  K=-0.249982  psi+K=+1.767e-05  a_wood=  27.28 m/s
  psi=0.50  K=-0.499947  psi+K=+5.302e-05  a_wood=  23.65 m/s
  psi=0.75  K=-0.749841  psi+K=+1.590e-04  a_wood=  27.32 m/s
  psi=0.90  K=-0.899523  psi+K=+4.770e-04  a_wood=  39.43 m/s
 
[3] liquid-volume change after 1->10 bar compression
  N     exact        ACID         K=0          K flipped    | x exact
    50  -1.7027e-04  -1.7027e-04  -3.1695e-04  -7.1967e-04 |  1.861  4.227
   100  -1.7027e-04  -1.7027e-04  -3.1559e-04  -7.0606e-04 |  1.853  4.147
   200  -1.7027e-04  -1.7027e-04  -3.1256e-04  -6.7737e-04 |  1.836  3.978
   400  -1.7027e-04  -1.7027e-04  -3.0670e-04  -6.2727e-04 |  1.801  3.684

K = 0、反号、体积模量比滑块会改变液柱缩短的速度。

ψ₀ (initial liquid fraction) = 0.50final p/p₀ = 10log₁₀(ρ_l a_l² / ρ_g a_g²) = 4.3

p/p₀ = 1.00  |  ψ = 0.50000  |  K = 0.00000  |  ψ+K = 5.00e-1
V_liquid: scheme 0.50000 vs exact 0.50000  →  +0.00%

打开 K = 0 再播放，蓝色块比精确解缩得更快。把比值降到 $1.0$ ，三根柱子叠到一起，而 $\psi+K$ 的读数离开 $0$ 。

压力升高十倍的过程中，水实际减少的体积是 $1.70\times10^{-4}$ ，占初始体积的0.034%。ACID封闭把这个值对到小数第五位。 $K=0$ 削掉的却是 $3.13\times10^{-4}$ ，是物理压缩量的1.84倍。

0.034%和0.062%之差听起来很小。但对照的基准不是总体积，而是物理压缩量。格式凭空添加了这个量的84%。在气泡崩溃、空化这类体积变化本身就是答案的问题里，这个比例直接决定结论。

表[2]接近对称。 $\psi=0.10$ 时 $K=-0.099994$ ， $\psi+K=5.891\times10^{-6}$ ，声速 $39.23$ m/s。 $\psi=0.90$ 时 $K=-0.899523$ ， $\psi+K=4.770\times10^{-4}$ ，声速 $39.43$ m/s。两端声速只差 $0.20$ m/s， $\psi+K$ 却拉开约 $81$ 倍。谷底是 $\psi=0.50$ 的 $23.65$ m/s，比 $39.23$ 慢 $15.58$ m/s，接近 $40\%$ 。体积模量比 $2.64\times10^{9}/1.40\times10^{5}$ 等于 $18857$ 。这个比把 $\psi=0.5$ 的 $K$ 贴到 $-0.5$ 旁边，也就是表里的 $-0.499947$ ， $\psi+K$ 为 $+5.302\times10^{-5}$ 。表[3]里 $K=0$ 的倍率从 $N=50$ 的 $1.861$ 降到 $N=400$ 的 $1.801$ ，少了 $0.060$ 。反号的倍率从 $4.227$ 降到 $3.684$ ，少了 $0.543$ 。同样加密 $8$ 倍，符号误差的降幅大约是 $0.060$ 的 $9$ 倍。 $N=400$ 时 $1.801$ 和 $3.684$ 仍然远离 $1$ 。 $N=50$ 上 $3.1695\times10^{-4}$ 与 $1.7027\times10^{-4}$ 之差是 $1.4668\times10^{-4}$ ，四舍五入就是上文的 $1.467\times10^{-4}$ 。 $N=400$ 是 $3.0670\times10^{-4}$ 减去 $1.7027\times10^{-4}$ ，得 $1.3643\times10^{-4}$ 。 $10$ bar时液体密度从 $1000.000$ 升到 $1000.341$ ，增加 $0.341$ ，气体从 $1.200$ 跳到 $6.215$ 。丢掉的 $K$ 向界面索取的，是这份很小的液体压缩的 $1.86$ 倍。

网格加密8倍，误差会消失吗#

表[3]最后两列给出答案。 $N$ 从50加到400，加密8倍， $K=0$ 的比值只从1.861变到1.801。按误差绝对值算，从 $1.467\times10^{-4}$ 到 $1.364\times10^{-4}$ ，减少7%。这谈不上收敛。

原因在于 $K=0$ 不是离散化误差。把 $K$ 置零，所求解的方程本身就变成了 $\mathrm{D}\psi/\mathrm{D}t=0$ 。格式收敛到的是另一个方程的解，没有理由在 $\Delta x\to0$ 时让差距消失。

剩下那7%来自抹开的界面带随网格略微改变形状。误差大小等于在该带内对 $\psi\,\nabla\cdot u$ 的积分，而一阶上风格式的数值扩散在加密时会按单元数把带拉宽。带的物理厚度几乎不变，积分值也就几乎不变。这条带一步能推多远，由VOF界面输运的CFL上限另行规定。

这个性质既是坏消息也是好消息。网格收敛性测试抓不到这个错误。反过来说，跑两套网格而比值维持不变，这份平坦本身就是封闭写错了的信号。

符号反过来会错成什么样

第三列是反号的结果。 $K\to-K$ 时系数变成 $\psi+K \to 2\psi$ ，压缩量跳到物理值的3.68～4.23倍。方向也翻了：明明在压缩，液体体积分数却在下降，画面上看就像气泡在长大。

麻烦在于这一侧同样给出一幅像样的图。界面依旧光滑，体积分数依旧落在 $[0,1]$ 内，也不发散。检查单相极限同样抓不到，因为 $\psi\to0$ 和 $\psi\to1$ 时两种符号下 $K$ 都为零。

所以把论文公式搬进代码时，用注释钉死 $\psi$ 代表哪一相的体积分数，再用上面那个一维压缩算例锁定符号，这样最快。50个单元大约一秒就跑完。

同一个假设在焓里再用一次的地方

论文中的ACID（Acoustically-Conservative Interface Discretisation）把这个等压封闭分别用在密度和焓上。密度是容易的一半。对分密度做线性插值得到

\rho = \rho_g + \psi(\rho_l - \rho_g)

它与等压封闭在数学上完全等价，质量和动量守恒随之成立。

焓不一样。 $h$ 是线性系统中直接求解的一次变量，不能像密度那样换掉插值式。于是方法另外计算一个带有ACID假设的目标面焓。

\hat{h}_f = \frac{1}{\tilde{\rho}_f}\left[\rho_U h_U + \xi_f L_f (\rho_D h_D - \rho_U h_U)\right]

$U$ 和 $D$ 是上游与下游单元， $L_f$ 是插值权重， $\xi_f$ 是限制器。再把它与当前值之差作为延迟修正项加到能量方程右端。

\delta h_f = \hat{h}_f - \tilde{h}_f

$K$ 与 $\delta h_f$ 是同一个假设的两张面孔。在体积分数方程里它代数地闭合，坍缩成一个系数；在焓这个一次变量上它闭合不了，只能延迟成显式源项。延迟的那一半必须在迭代内收敛，所以会反映到非线性迭代次数上。用等压收到同一个压力上，是先用两个压力求解再松弛成一个压力那条路的终点。差别在于这边把那次松弛折进了一个 $K$ 。

Denner能丢掉黎曼求解器的原因#

密度基解法把界面当黎曼问题求解，从而保持接触间断。压力基全耦合解法把 $\rho$ 从求解变量中剔除，把 $u$ 、 $p$ 、 $h$ 绑进同一个线性系统。同位网格上出现的压力-速度解耦靠MWI压住。在极端密度比下，面密度用调和平均加权。用重构保住接触间断的另一条路，是挑选更小振幅的MUSCL-THINC-BVD。

丢掉黎曼求解器，同时也丢掉了自动照料界面物理的那套机制。这个缺口只能用封闭假设手工填上。 $K$ 与 $\delta h_f$ 就是填料。

要隐式跨过23.65 m/s的混合物声速又不把接触间断抹开，封闭就必须准确。声学一旦隐式，CFL的账单就转到输运一侧。把声学和输运拆开再解CFL的分离就是那个位置。1.86倍，就是对这个封闭所做工作量的一次测量。

去掉K之后水被压缩了1.86倍 — 加密网格也不会减小的界面误差

可压缩VOF为什么带着 −(ψ+K)∇⋅u-(\psi+K)\nabla\cdot u−(ψ+K)∇⋅u#

KKK 是从哪里来的数#

如果 ψ+K\psi+Kψ+K 是 5×10−55\times10^{-5}5×10−5，两项能一起删掉吗#