方程少了，反而不炸了 — 四方程模型与界面平衡条件

均匀压力留得住，是因为直接迎风Γ；把守恒的ρΓ再除以ρ，压力就会弯。

删掉一个方程，代码反而不炸了。可压缩多相流(不同流体掺在一起的可压缩流动)通常反过来教：每一相(phase，物质的不同形态)各自带压力、速度和温度的七方程最安全，方程越少越悬。Collis(2025)把这个次序翻过来。最短的四方程反而撑得住。钥匙是界面平衡。

打开Naive开关，把right gas的γ推到1.66或1.09，界面上的红色压力线升高，绿线保持平直。

right gas γ = 1.66CFL = 0.60

Uniform velocity, uniform pressure, only γ jumps across the shaded interface. The naive scheme grows a pressure spike at the interface; the consistent one keeps P flat to machine precision.

红峰随着γ的跳跃一起变高。步数推到900时绿线是否还贴在地板上，就是这个条件活着还是死了。

平衡的阶梯 — 七、六、五、四方程

扩散界面方法(diffuse interface，把交界面摊在几个单元里来解)需要一条规则，来规定抹开的单元里两相各处在什么状态。这条规则就是“把什么当成平衡”。

七方程(Baer–Nunziato)：每一相的速度、压力、温度都不相同。最一般，也最贵。
六方程：速度共享。压力分开，再用松弛(relaxation)对齐。
五方程：速度和压力共享。只把温度分开。
四方程：速度、压力、温度全部共享。最短。

物理界面的厚度是纳米。数值界面只是把那层薄片吹成几个格子。在这么窄的一条带里，热平衡和力学平衡实际上马上就完成。界面附近的物理，四方程最接近。

T^{(1)} = T^{(2)}, \quad P^{(1)} = P^{(2)}, \quad \mathbf{u}^{(1)} = \mathbf{u}^{(2)}

上标 $(1),(2)$ 是两相。四方程在单元内部把这三个等式全部强制住。

方程越少为何越稳健

五方程及以上通常再加一条输运方程，给体积分数(volume fraction，一个单元里某一相占的体积比例)，或者给贴在比热比上的 $1/(\gamma-1)$ 。在PDE层面这些方程是重复的。离散之后它们会彼此错开。同一个量，两条路径，两个值。

四方程把多余的拿掉。只守恒质量、动量和能量。压力用混合物状态方程闭合。方程少，错开的地方也少。花费也下降。

只剩一个问题。那些多余方程原来提供的余量没了，守恒型离散就会在界面上把压力拉开。

界面平衡条件被打破的那一刻

界面平衡条件(IEC，interface equilibrium condition)很直白。压力、温度、速度一开始是均匀的，界面走过去，这份均匀就不该破。

违背它的格式会在界面附近长出原本没有的压力振荡。振荡随时间变大。密度比一大，负密度或负内能冒出来，计算就停。

为什么守恒型会违背它。理想气体混合物里令 $\Gamma \equiv 1/(\gamma-1)$ ，总能量写成

E = P\,\Gamma + \tfrac{1}{2}\,\rho\,u^2

$E$ 是总能量， $P$ 是压力， $\rho$ 是密度， $u$ 是速度。 $\Gamma$ 随流体而变，在界面上跳跃。

在均匀的 $u,P$ 下用一阶迎风推一格。 $\nu = u\,\Delta t/\Delta x$ 是CFL数。守恒型的能量更新把 $+P$ 项消掉，剩下

E_i^{n+1} = E_i - \nu\,(E_i - E_{i-1})

代进总能量的定义，右边就裂开。

E_i^{n+1} = P_0\big[\,\Gamma_i - \nu(\Gamma_i - \Gamma_{i-1})\,\big] + \tfrac{1}{2}\,u_0^2\,\rho_i^{n+1}

新压力要保持为 $P_0$ ， $\Gamma$ 就必须按那个方括号来更新。这是非守恒迎风。守恒型把 $\rho\Gamma$ 整个搬走，再除以 $\rho$ 。两个结果在界面上分开。压力明明均匀，非守恒路径一对不齐压力就晃，这个位置就是Abgrall留下的判据。这个差就是压力振荡。

刚才推动γ时红线变高，就是这个差。换成IEC-consistent，绿线在900步时仍然是平的。

Python — 振荡的真身与药方#

同一套计算。均匀的 $u,P$ ，只有 $\gamma$ 在界面上跳，用一阶迎风推进。两个格式只差这一处。

import numpy as np
 
def advect_interface(scheme, gamma_R=1.66, nu=0.6, steps=300, N=200):
    x = np.linspace(0.0, 1.0, N, endpoint=False)
    ramp = np.clip((x - 0.28) / 0.06, 0.0, 1.0)   # 抹开在几个单元上的界面
    gamma = 1.4 + (gamma_R - 1.4) * ramp          # 空气(1.4) -> 所选气体
    rho = 1.0 + (0.2 - 1.0) * ramp                # 密度也跳变
    u0, P0 = 1.0, 1.0
    Gam = 1.0 / (gamma - 1.0)                      # 混合物刚度 1/(gamma-1)
    mom = rho * u0
    E = P0 * Gam + 0.5 * rho * u0**2               # E = P*Gamma + 0.5*rho*u^2
 
    for _ in range(steps):
        rhoG = rho * Gam                           # 更新前的 rho*Gamma
        r_up, m_up, e_up = np.roll(rho, 1), np.roll(mom, 1), np.roll(E, 1)
        rho = rho - nu * (rho - r_up)              # 质量:守恒型迎风
        mom = mom - nu * (mom - m_up)              # 动量:守恒型
        E   = E   - nu * (E   - e_up)              # 能量:守恒型
        if scheme == "naive":
            rhoG = rhoG - nu * (rhoG - np.roll(rhoG, 1))  # 守恒输运 rho*Gamma
            Gam = rhoG / rho                              # 再除回 rho
        else:  # consistent
            Gam = Gam - nu * (Gam - np.roll(Gam, 1))      # 直接对 Gamma 做非守恒迎风
 
    P = (E - 0.5 * mom**2 / rho) / Gam             # 还原压力
    return float(np.max(np.abs(P - P0)))
 
for s in ("naive", "consistent"):
    print(f"{s:11s}  max|P-P0| = {advect_interface(s):.3e}")

输出分开的不是两位，是十三位。

naive        max|P-P0| = 3.83e-02
consistent   max|P-P0| = 4.44e-15

naive在界面上把压力偏了大约4%。consistent平到浮点舍入。造成差别的是一行。把 $\Gamma$ 用 rho 除回来，还是直接搬运 $\Gamma$ 本身。

把 $\Gamma$ 摊成数字，跳跃的高度就出来了。 $\Gamma = 1/(\gamma-1)$ ，空气 $\gamma=1.4$ 给出 $\Gamma=2.5$ 。氦气 $\gamma=1.66$ 给出 $\Gamma=1/0.66 \approx 1.515$ 。SF6 $\gamma=1.09$ 给出 $\Gamma=1/0.09 \approx 11.11$ 。空气到氦气的差约0.985。空气到SF6约8.61。SF6的跳跃大约是氦气的8.74倍。γ推到1.09时红线变高，就是这个倍数。

用代码里的 $\nu=0.6$ ，只在一个尖锐格子上推一步。右格是 $\rho=0.2$ 、 $\Gamma\approx 1.515$ ，迎风邻居是空气 $\rho=1$ 、 $\Gamma=2.5$ 。把 $\rho\Gamma$ 按守恒形式搬走再除以 $\rho$ ，得到 $\Gamma_{\mathrm{naive}}\approx 2.384$ 。直接推 $\Gamma$ ，得到 $\Gamma_{\mathrm{cons}}\approx 2.106$ 。把上面的能量式倒回压力， $P/P_0 = \Gamma_{\mathrm{cons}}/\Gamma_{\mathrm{naive}} \approx 0.883$ ，所以一步的 $|P-P_0|\approx 0.117$ 。同一个尖锐格子换成SF6， $|P-P_0|\approx 0.692$ 。输出里的 $3.83\times 10^{-2}$ 不是这个尖锐格子。斜坡宽度 $0.06$ 在 $\Delta x=1/200$ 上是12个单元，300步里界面走了 $\nu \cdot 300 \cdot \Delta x = 0.9$ 。把抹开的界面推这么久，最大误差会小于一步跳跃，停在 $0.0383$ 。尖锐一步的 $0.117$ 与这个 $0.0383$ 之比约为 $3.05$ 。摊到12个单元之后，误差仍是几个百分点。consistent的 $4.44\times 10^{-15}$ 是两条路径已经一致之后剩下的舍入。

论文真正做了什么

上面的玩具只推了一个理想气体的 $\gamma$ 。Collis(2025)的场地更宽。

用 NASG状态方程(Noble–Abel Stiffened Gas，能写到液体的状态式)同时闭合水、空气、氦气和SF6。混合压力按Amagat定律收成一个二次方程。
用高阶 ENO重构(WENO/TENO)替换一阶迎风。重构和四方程的热–力学平衡假设咬合，IEC因此达到机器精度。
为了强激波打在大密度比界面上的情形，把 保正限制器(positivity-preserving limiter)扩到四方程。密度和声速平方不被漏成负数。
这些都没有再加体积分数或比热比方程，守恒性和标量滤波也都留着。

限度也在。液相限单组分，气相限理想气体组分。组分一多，压力和温度的平衡就要用迭代求解器。目前仍是预印本(同行评审之前)。

用ρ把Γ除回来，平坦的压力一步就弯

一个尖锐格子上，氦气已经是 $|P-P_0|\approx 0.117$ 。SF6是 $0.692$ 。300步输出的 $0.0383$ 是这个跳跃被12个单元的斜坡分摊之后的结果。只要留下“守恒地搬走 $\rho\Gamma$ 再除以 $\rho$ ”这一行，本来均匀的 $P_0$ 就在界面上立刻变弯。少写几个方程，换来把这一行改成非守恒，这就是四方程的全部选择。

五方程共享速度和压力，仍留下一个多余方程。如果那一项正是把界面声速压垮的东西，那么为了保住IEC而删掉它，是把病也删了，还是只藏进混合物状态方程里，还没定。水和空气各掺一半时声速垮掉的五方程是否还坐在四方程的闭合里面，是下一处分叉。