用并不存在的流体填满界面 — Ghost Fluid Method

在水与空气相遇的那个网格里，倒进一团本不在那里的假空气。这种鬼影流体(ghost fluid，越过界面外推得到的虚拟流体)在网格上任何位置都不存在。可正是这团鬼影，让冲击波穿过界面时计算保持稳定。要复制的量和要外推的量不是同一件事。原始GFM发响，就响在这个分工弄错的地方。

把密度跳跃拉高，台阶变陡，压力仍沿红色虚线走过去。「show ghost」开关的是界面另一侧的延长线。

打开虚线后，每种流体带着自己的熵伸到对面。挪动界面位置，压力曲线还是一条。断开的是密度。

当两种物质在一个网格里相撞

气体和液体、炸药和金属——两种物性差异巨大的物质在同一网格上相遇。直接用平均物性，界面就会在数值上被抹糊。在密度比1000倍的水-空气界面上，这种抹糊很快蔓延成非物理的振荡。

GFM的思路不同。它不把界面模糊地混合。它把每种流体当作假流体延伸到自己区域之外，然后像只有一种物质那样把单物质求解器跑两遍。界面以水平集(level set，用符号划分区域的距离函数)$\phi$ 的零点来追踪。

这里 $\phi$ 是符号距离函数，其零点 $\phi=0$ 就是界面。

越过界面，什么断裂什么延续

要正确填充鬼影，必须确切知道越过界面什么是连续的。压力 $p$ 与法向速度 $u_n$ 在界面处连续。而熵 $s$ 与切向速度 $u_t$ 则不连续。

方括号 $[\cdot]$ 表示越过界面的跳跃。这四行就是GFM的全部。连续量复制，不连续量外推。上面那张图里只有压力光滑、只有密度呈台阶，原因就在这四行。

在一个网格里把两相平均掉，会出另一种故障。六方程 two-fluid 模型在密度比很大的界面上特征值变成复数，网格一加密，发散来得更早。GFM不做这个平均。一个网格只装一种物质。

把水当作气体 — Stiffened EOS#

水几乎不可压缩，所以理想气体状态方程无法处理它。为了把两种物质归入同一形式的Mie–Grüneisen族，GFM采用刚化气体状态方程(stiffened-gas EOS)。

$\rho$ 是密度，$e$ 是单位质量内能，$\gamma$ 是比热比，$p_\infty$ 是表征物质刚度的常数。理想气体取 $p_\infty = 0$；水约为 $\gamma \approx 4.4$、$p_\infty \approx 6 \times 10^8\,\mathrm{Pa}$。声速由此得到：

正因为 $p_\infty$，相同压力下水的声速远大于空气。把状态方程统一成一行，两种流体就能放进同一个求解器。

填满鬼影网格的两道手法

这是核心步骤。求解流体A的区域时，用A的鬼影填满B占据的网格。手法分两道。

1. 复制：连续量——压力和法向速度——直接取自界面对面真实B网格的值。
2. 外推：不连续量——熵(即密度)和切向速度——从最近的真实A网格越过界面做等熵外推。

等熵外推保持熵 $s = p / \rho^\gamma$ 不变，重新解出密度。

$s_A$ 是真实A流体的熵，$p_\text{copy}$ 是复制过来的压力。外推通常用快速行进(fast marching)方法沿界面法向快速铺开。一次扫掠完成A的鬼影场后，就忘掉B的存在，只解A的单物质Euler方程。对B重复同样的步骤。

只用代码里的初值就能算出鬼影密度，外推把真实密度改成了什么，数字上是清楚的。左侧密度1、压力1、比热比1.4，熵 $s = p / \rho^\gamma$ 等于1。把右侧压力0.1抄过来延伸左侧流体，鬼影密度是 $0.1^{1/1.4}$，约0.193。那不是真实的右侧密度0.125。单物质求解器看到的密度比不是8倍，而是约5.2倍。

反方向用同一条式子。右侧比热比1.67，熵是0.1除以 $0.125^{1.67}$，约3.22。把左侧压力1放到这个熵上，鬼影密度约0.496，接近真实左侧密度1的一半。压力交过去，每种流体的熵留在自己的鬼影里。

声速也直接从状态方程出来。刚度常数为0时，声速是比热比、压力、密度组成的平方根。左侧约1.183，右侧约1.156。400个点、长度1，间距是 $1/399$，约0.00251。CFL取0.4、用左侧声速，时间步长约0.00085。终止时刻0.14除下去，大约165步。声波穿过长度1大约要0.85，0.14还不到它的五分之一。和水-空气密度比1000比，这根激波管的8倍算温和。鬼影密度0.193和0.496已经和真实值分开了。

原始GFM振荡的地方#

当界面一侧非常刚硬(stiff，像水那样 $p_\infty$ 很大)时，原始GFM就会出问题。强冲击波或爆轰波穿过界面时，靠简单复制和外推构造的鬼影状态与真实Riemann解发生偏差。这种偏差在界面附近表现为压力振荡，振荡往往导致计算发散。

在下面的模拟里，让一道压力脉冲穿过刚硬流体(粉色区域)。stiffness ratio越高，越过界面之后的涟漪就越大。

把stiffness ratio推到9，界面之后立刻出现清晰的高频振荡。加密网格也消不掉它，因为根源是边界条件而非网格。

修正GFM — 用界面Riemann问题确定取值#

解法是把鬼影状态解出来，而不是去猜。修正GFM(Modified GFM，MGFM或IGFM)把界面两侧的真实状态作为左右态，求解Riemann问题。

$W_A, W_B$ 是界面两侧的真实状态，$\mathcal{R}$ 是Riemann解法，$p^*, u^*$ 是接触间断的压力与速度。把同一个 $p^*, u^*$ 赋给两侧鬼影，压力-速度连续条件就被精确满足。在上面的模拟里打开「use modified GFM」，相同刚度下振荡随即消失。

Python — 气-气GFM激波管#

用等熵外推GFM求解两种理想气体($\gamma_A = 1.4$、$\gamma_B = 1.67$)相遇的激波管。采用HLL通量，每步做两次鬼影扫掠后用水平集合并。

import numpy as np
 
GAMMA = (1.4, 1.67)   # 左(空气)、右(单原子气体)的比热比
 
def primitive(U, g):
    rho = U[0]; u = U[1] / rho
    e = U[2] / rho - 0.5 * u * u
    p = (g - 1.0) * rho * e          # 刚化气体EOS, p_inf = 0 (理想气体)
    return rho, u, p
 
def conserved(rho, u, p, g):
    E = p / (g - 1.0) + 0.5 * rho * u * u
    return np.array([rho, rho * u, E])
 
def hll(UL, UR, g):
    rL, uL, pL = primitive(UL, g); rR, uR, pR = primitive(UR, g)
    aL = np.sqrt(g * pL / rL); aR = np.sqrt(g * pR / rR)
    sL = min(uL - aL, uR - aR); sR = max(uL + aL, uR + aR)
    FL = np.array([rL * uL, rL * uL * uL + pL, uL * (UL[2] + pL)])
    FR = np.array([rR * uR, rR * uR * uR + pR, uR * (UR[2] + pR)])
    if sL >= 0: return FL
    if sR <= 0: return FR
    return (sR * FL - sL * FR + sL * sR * (UR - UL)) / (sR - sL)
 
def extrapolate_ghost(rho, u, p, mat, side):
    # 用 side 流体填满整个区域: 压力、速度复制, 密度做等熵外推
    rg, ug, pg = rho.copy(), u.copy(), p.copy()
    g = GAMMA[side]
    idx = np.where(mat == side)[0]
    lo, hi = idx[0], idx[-1]
    ref = lo if side == 1 else hi      # 离界面最近的真实网格
    s_ref = p[ref] / rho[ref] ** g     # 需保持的熵
    for i in range(len(rho)):
        if mat[i] != side:
            j = min(max(i, lo), hi)            # 最近的真实网格
            ug[i] = u[j]; pg[i] = p[j]         # 复制
            rg[i] = (pg[i] / s_ref) ** (1.0 / g)   # 等熵外推
    return rg, ug, pg
 
def run_ghost_fluid_tube(N=400, t_end=0.14, cfl=0.4):
    x = np.linspace(0, 1, N); dx = x[1] - x[0]
    x0 = 0.5                                   # 初始界面位置
    mat = np.where(x < x0, 0, 1)
    rho = np.where(x < x0, 1.0, 0.125)
    u   = np.zeros(N)
    p   = np.where(x < x0, 1.0, 0.1)
    t = 0.0
    while t < t_end:
        a = np.sqrt(np.where(mat == 0, GAMMA[0], GAMMA[1]) * p / rho)
        dt = min(cfl * dx / np.max(np.abs(u) + a), t_end - t)
        new = {}
        for side in (0, 1):                    # 两次鬼影扫掠
            rg, ug, pg = extrapolate_ghost(rho, u, p, mat, side)
            g = GAMMA[side]
            U = np.array([conserved(rg[i], ug[i], pg[i], g) for i in range(N)]).T
            F = np.zeros((3, N + 1))
            for i in range(1, N):
                F[:, i] = hll(U[:, i - 1], U[:, i], g)
            F[:, 0] = F[:, 1]; F[:, N] = F[:, N - 1]
            new[side] = U - dt / dx * (F[:, 1:] - F[:, :-1])
        U = np.where(mat[None, :] == 0, new[0], new[1])   # 水平集合并
        for i in range(N):
            rho[i], u[i], p[i] = primitive(U[:, i], GAMMA[mat[i]])
        x0 += dt * u[np.argmin(np.abs(x - x0))]           # 界面移动
        mat = np.where(x < x0, 0, 1)
        t += dt
    return x, rho, u, p
 
if __name__ == "__main__":
    x, rho, u, p = run_ghost_fluid_tube()
    print(f"t=0.14:  max p = {p.max():.3f}, "
          f"接触面附近密度 {rho[180]:.3f} -> {rho[220]:.3f}")

输出显示越过界面的密度跳跃得以保留，而压力平滑相连。两次单物质扫掠就是GFM的脊梁。

写代码时不该踩的陷阱

外推如果顺着网格轴走，斜界面上的熵会混到一起。铺开的方向只能是界面法向。刚度比很大的水-空气，原始的复制加外推几乎总会在压力上留下涟漪。那种情况从一开始就用界面Riemann解填鬼影。水平集不拉回符号距离的话，零点还在，距离值已经乱了，界面速度跟着偏。

把气体-气体上好用的外推原样放到水上时#

$p_\infty$ 为0的两种理想气体($\gamma = 1.4$ 与 $1.67$)只做等熵外推，接触面上的密度台阶还在。水是 $\gamma \approx 4.4$、$p_\infty \approx 6 \times 10^8\,\mathrm{Pa}$。同样的压力下，声速分子里加上这个常数，演示里把 stiffness ratio 推到9时看到的涟漪，只是同一处偏差的缩小版。网格加密也清不掉，它留在边界条件里。

鬼影里的压力和速度继续用抄来的值，还是改成界面Riemann解的 $p^*, u^*$，要在数密度跳跃之前定下来。外推密度掉到负数是另一件限制器的事。保正通量限制器处理的是质量对了、格子密度却为负的情况。刚硬界面看完，接着看那里比较值。