每一步都把糊掉的界面拉回来 — 亲手实现反扩散锐化(IST)

只把界面削薄，格数还会再漂。每步之后固定ε/a，激波旁边的带宽才留得住。

用五方程模型跑了一个激波–气泡问题。200 步之后，氦气泡的边界从 6 个网格膨胀到了 20 个。物理量是对的，但界面糊成了一团粥。这就是扩散界面法(让界面在有限宽度内弥散的一类方法)的宿命。Nguyen 等(2021)在每个时间步之后用一次后处理，把这种弥散拉回来。今天我们从零实现这套反扩散技术，亲眼看清为什么界面厚度会锁定为一个常数。

打开 numerical diffusion，6 个网格的边界会把格子吞进去；打开 sharpening (IST)，这层扩散停在 ε/a 附近。

compression a: 1.00regularization ε: 0.0040

measured band: 0 cells

target ≈ ε/a: 2 cells

steady tanh thickness scales like ε/a

黄色带是目标厚度。调高 compression a 带子变窄，调高 regularization ε 又变厚。看测量格数是否跟着这个比值走就够了。

论文信息

标题: Numerical modeling of multiphase compressible flows with the presence of shock waves using an interface-sharpening five-equation model
作者: Van-Tu Nguyen, Thanh-Hoang Phan, Warn-Gyu Park (Pusan National University)
来源: International Journal of Multiphase Flow, 2021, 103542
DOI: 10.1016/j.ijmultiphaseflow.2020.103542
一句话总结: 把界面锐化(IST)作为后处理附加到五方程可压缩两相模型上，即使紧邻激波也能让界面保持恒定厚度。

界面为何会自己变胖

五方程模型是准守恒系统：守恒律再加一个颜色函数。体积分数 $\alpha_1$ (单元格内流体 1 所占的比例)满足这样一个对流方程。

\frac{\partial \alpha_1}{\partial t} + \mathbf{u}\cdot\nabla \alpha_1 = \alpha_1 K\, \nabla\cdot\mathbf{u}

其中 $\alpha_1$ 是流体 1 的体积分数， $\mathbf{u}$ 是混合速度，右端的 Kapila 项修正两相压缩率的差异。

麻烦出在用激波捕捉格式求解它的时候。一旦在单元格面上平滑插值，本该锋利的 $\alpha_1$ 台阶每一步都会糊开一点。这种数值扩散会累积。步数一多，界面就从几个网格扩展到几十个，气泡形状、激波反射位置这些核心特征整个消失。

WENO 这类高阶格式能减少扩散，但会带来振荡，在多维中还很昂贵。于是作者转了个方向。他们不动格式本身，而是在每一步结束后单独修整 $\alpha$ 场。

反扩散：锐化方程

核心是 Shukla 等(2010)的正则化方程。我们对虚拟时间 $\tau$ (而非物理时间)让 $\alpha$ 演化若干次。

\frac{\partial \alpha}{\partial \tau} = \nabla\cdot\!\left(\varepsilon\, \nabla \alpha\right) - \nabla\cdot\!\left(a\, \alpha(1-\alpha)\, \hat{\mathbf{n}}\right)

$\varepsilon$ 是调控厚度的扩散系数， $a$ 是收紧界面的压缩强度， $\hat{\mathbf{n}} = \nabla\alpha / |\nabla\alpha|$ 是界面法线方向。

两项在拔河。第一项是普通扩散，把 $\alpha$ 摊开。第二项符号相反，是反扩散(anti-diffusion)，沿法线方向收紧界面。多亏 $\alpha(1-\alpha)$ 因子，压缩只在界面( $0<\alpha<1$ )处作用，在纯流体区( $\alpha=0$ 或 $1$ )则关闭。

在一维求解定常态( $\partial_\tau \alpha = 0$ )，两个通量平衡。

\varepsilon\, \frac{d\alpha}{dx} = a\, \alpha(1-\alpha) \quad\Longrightarrow\quad \alpha(x) = \tfrac{1}{2}\!\left[1 + \tanh\!\left(\frac{a}{2\varepsilon}\,x\right)\right]

界面收敛到双曲正切形状，其厚度正比于 $\varepsilon/a$ 。也就是说，厚度由你来定。这正是作者强调的"恒定厚度"性质。当界面始终保持相同的网格数，激波与接触间断就能被稳定捕捉。

一维直接看：弥散与锐化的拔河

上面那场拔河写成脚本就是下面这段。先模仿对流留下的弥散，再在虚拟时间里把锐化方程跑几遍。

import numpy as np
 
N, dx = 201, 1.0 / 200
x = np.linspace(0, 1, N)
 
def band_width(a, lo=0.05, hi=0.95):
    # 视为界面的网格数（厚度度量）
    return int(np.sum((a > lo) & (a < hi)))
 
def smear_once(a, D=0.16):
    # 对流格式每步把界面变胖的效果
    lap = np.zeros_like(a)
    lap[1:-1] = a[2:] - 2 * a[1:-1] + a[:-2]
    return a + D * lap
 
def sharpen_sweep(a, comp_a, eps, iters=200):
    # 通量形式的反扩散正则化 (Shukla 2010, 式 (38)-(39))
    dtau = 0.9 * min(dx * dx / (2 * eps), dx / comp_a)
    for _ in range(iters):
        af = 0.5 * (a[:-1] + a[1:])          # 面心 alpha
        dA = a[1:] - a[:-1]
        s = np.sign(dA)
        J = comp_a * af * (1 - af) * s - eps * dA / dx   # 面通量
        a[1:-1] = np.clip(a[1:-1] - dtau / dx * (J[1:] - J[:-1]), 0, 1)
    return a
 
# 从锋利台阶开始 -> 让它弥散 20 次
a = 0.5 * (1 + np.tanh((x - 0.5) / 0.012))
for _ in range(20):
    a = smear_once(a)
print("弥散后的带宽:", band_width(a), "格")     # -> 变宽
 
a = sharpen_sweep(a, comp_a=1.0, eps=0.004)
print("锐化后的带宽:", band_width(a), "格")     # -> 收缩到目标厚度

弥散后带宽很宽，锐化后缩到 $\varepsilon/a$ 所定的几个网格。格式原样不动。sharpen_sweep 只接收结果数组，并不知道是哪种黎曼求解器填进去的。

默认参数是 comp_a=1.0、eps=0.004。厚度尺度 $\varepsilon/a$ 就是 0.004。定常剖面的斜率系数是 $a/(2\varepsilon)=125$ 。代码起步的 tanh 分母是 0.012，所以目标比初始台阶更陡。 $dx=1/200=0.005$ 。虚拟时间上限取 $dx^2/(2\varepsilon)=0.000025/0.008=0.003125$ 和 $dx/a=0.005$ 里较小的那个。扩散限制胜出。乘上 0.9 之后 $d\tau$ 是 0.0028125。越过这个上限，反扩散会翻号，界面开始振荡。

面通量里的压缩项正比于 $\alpha(1-\alpha)$ 。这个乘积在 $\alpha=0.5$ 时达到 0.25，在纯流体里是 0，压缩随之关掉。带宽函数只数 0.05 到 0.95 之间的格子。只把 $\varepsilon$ 加倍，目标厚度也加倍；只把 $a$ 加倍，厚度减半。开头的 200 步是物理时间。6 个格子变成 20 个，带宽大约是 3.3 倍。smear_once 里的 D=0.16 和锐化循环的 iters=200 是两层不同的东西。论文说的每步 1 到 3 次，前提是每个物理步的扩散很小。Kapila 右端 $\alpha_1 K\nabla\cdot\mathbf{u}$ 在激波把散度打成尖峰时，会直接推动体积分数。这和数值扩散不是一回事，IST 也不擦掉这个源。跟把压力松弛单独放下的 Saurel(2009) 步骤不同，这里只把界面格数拉回去。

把厚度做成常数意味着什么

把界面做薄，和把厚度钉成一个常数，是两回事。作者与早期方法(Tiwari 等 2013)划清界限的地方正在于此。在源项里加锐化函数能减少扩散误差，但厚度会随时间与空间起伏。相反，用后处理固定 $\varepsilon/a$ ，处处都得到相同的厚度。

在可压缩流里还多一层。修整 $\alpha$ 时，混合密度与能量必须随之改变，热力学才能自洽。于是作者不只修正 $\alpha$ ，而是按照混合律把整个守恒变量向量( $\alpha_1\rho_1$ 、 $\alpha_2\rho_2$ 、动量、能量)重新分配。多亏这种"混合一致正则化"，界面处的速度、压力、温度平衡不会被打破。论文的一维纯界面对流测试正说明了这点：当水和空气以相同速度流动，打开锐化后压力与温度仍保持不变。

二维气泡：保持圆形界面

一维通了，就转到二维气泡。论文的激波–气泡图把几格厚的圆一直保持到最后。打开 numerical diffusion，边缘变模糊；打开 sharpening 并调高 compression a，边缘又收成清晰的圆环。

compression a: 1.20

interface band: 0 cells

teal ring = cells with 0.1 < α < 0.9. Fewer = sharper.

青色圆环收紧时，气泡中心不动。既不变大也不变小。减少的只是环上的网格数。压缩沿着 $\hat{\mathbf{n}}$ ，曲率不会被压皱。

复现时撞上的问题

纯流体里 $|\nabla\alpha|\to 0$ ，法线会被零除。论文因为 $a\,\alpha(1-\alpha)$ 在那里是零就略过，但分母不加一个小量就会写出 NaN。 $d\tau$ 一大，界面就振荡。需要 CFL 型上限才稳。论文说的每步 1 到 3 次，前提是每个物理步的扩散很小。粗网格或强激波下，这个次数可能不够。

还有一点。"后处理"这份优雅是有代价的。锐化独立于控制方程去改变 $\alpha$ ，因此那一瞬的对流并不是原始 PDE 的解。换来恒定厚度的同时，界面附近的局部守恒可能会略有偏移。作者用混合一致再分配把它降到最低，但并非恰好为零。OpenFOAM 的 interfoam 所用的压缩项( $\nabla\cdot(\alpha(1-\alpha)\mathbf{u}_r)$ )属于同一类思路，把这套技术移植到可压缩求解器时，可以拿它作对照基准。

不锁住 ε/a，格数会再漂起来#

把锐化放进源项，扩散误差会减小，厚度仍会在时空里晃。步子结束之后固定 $\varepsilon/a$ ，激波旁边才能留下同样的格数。只改 $\alpha$ 、不按混合律把守恒变量重新分回去，那一格的压力和温度平衡会先坏。

$d\tau$ 一旦大于 $dx^2/(2\varepsilon)$ ，反扩散不是在收紧界面，而是在让它振荡。用常数厚度换来的代价是，那一瞬间的对流不再是原始 PDE 的解。局部守恒的那点偏移，和 Abgrall 准则抓住的非守恒项是同一类。