[论文综述] 选择跳跃更小的那一个 — Deng (2018) MUSCL-THINC-BVD 界面重构

比较两个候选重构函数在单元边界上的变差，挑选更小者的BVD原理

东京工业大学 Xiao 课题组在 2018 年做了一项很少有人喜欢公开的测量。他们对可压缩两相流模拟跑了 1000 步，再去测界面厚度——原本只有一个网格的跳跃，已经扩散到了八个网格。即便换用 WENO，结果也差不多。每步看似耗散很小的格式，长时间累计后仍会把界面磨平。本文整理他们对这个问题的回答 —— MUSCL-THINC-BVD 重构。思路很朴素：对每个网格同时构造两个候选，选择在网格边界上跳跃更小的那一个。

一页摘要#

作者 / 期刊: Deng, Inaba, Xie, Shyue, Xiao. Journal of Computational Physics 371 (2018) 945–966.
目标问题: 用五方程模型(five-equation model)模拟可压缩两相流时，物质界面随时间逐渐被磨平。
方案: 在每个网格同时构造 MUSCL(光滑区适用)和 THINC(跳跃区适用)两个候选重构函数，比较它们在网格边界上的跳跃和，取较小者，即 BVD (Boundary Variation Diminishing) 算法。
新意所在: 不需要后处理 anti-diffusion 或人工压缩。同一条 BVD 规则适用于体积分数和其他守恒量，变量间一致性自动得到保证。

两个相互冲突的需求#

可压缩两相流的数值方法同时面临两项要求：光滑区要精度高、耗散低；跳跃区(界面、激波)要单调、且厚度保持薄。让一个函数同时把两件事都做好并不容易。

MUSCL 是分段线性重构，保单调但导数只有一阶精度。每次穿越界面都会刮掉一点，长时间累计就让跳跃变胖。反之，THINC 在每个网格内拟合一个 tanh 跳跃，能把界面压缩到一两格内，但用在光滑区会制造假的台阶。

论文的出发点是一个决定：不要把两者糅合成一个函数，而是在每个网格里挑一个。

MUSCL — 可靠的扩散器#

基线候选是带 minmod 斜率限制器的 MUSCL。

\tilde{q}_i^{\,\text{MUSCL}}(x) = \bar{q}_i + \sigma_i (x - x_i), \quad \sigma_i = \text{minmod}(\bar{q}_i - \bar{q}_{i-1},\; \bar{q}_{i+1} - \bar{q}_i)

其中 $\bar{q}_i$ 为网格平均， $\sigma_i$ 为网格中心斜率，minmod 在两差分同号时返回绝对值较小者，异号时返回零。

候选的边界值：

q^{L,\text{MUSCL}}_{i+1/2} = \bar{q}_i + \tfrac{1}{2}\sigma_i, \quad q^{R,\text{MUSCL}}_{i-1/2} = \bar{q}_i - \tfrac{1}{2}\sigma_i

在任何位置都安全可用，没有振荡，但跳跃厚度会缓慢扩张。

THINC — 模仿跳跃的单调函数#

THINC (Tangent of Hyperbola for INterface Capturing) 在每格内用双曲正切拟合跳跃。

\tilde{q}_i^{\,\text{THINC}}(x) = \bar{q}_{\min} + \tfrac{\bar{q}_{\max}}{2}\left[1 + \theta \tanh\!\left(\beta\!\left(\tfrac{x - x_{i-1/2}}{\Delta x} - \tilde{x}_i\right)\right)\right]

$\bar{q}_{\min}, \bar{q}_{\max}$ 为从邻格平均取出的最小值与振幅， $\theta = \mathrm{sgn}(\bar{q}_{i+1}-\bar{q}_{i-1})$ 表示跳跃方向， $\beta$ 控制跳跃厚度， $\tilde{x}_i$ 是为保持网格平均而求解的跳跃中心。

$\beta$ 在 1.4–2.0 之间稳定，1.6 是标准值。 $\beta$ 越大跳跃越压缩到单格内，看起来很漂亮，但若被错误地用到光滑区，就会产生锯齿状假跳跃。

BVD — 在两个候选之间度量变差#

核心思想出现了。既然有两个候选，就可以问：在与邻格重构相遇的边界上，哪一个产生的跳跃更小。第 $i$ 格的 总边界变差 (Total Boundary Variation, TBV) 定义为

\text{TBV}_i^{P} = |q^{L,\text{MUSCL}}_{i-1/2} - q^{R,P}_{i-1/2}| + |q^{L,P}_{i+1/2} - q^{R,\text{MUSCL}}_{i+1/2}|

其中 $P$ 为候选(MUSCL 或 THINC)。两侧邻格固定为 MUSCL，仅中间格在切换。

选择规则：

\tilde{q}_i^{\,\text{BVD}}(x) = \begin{cases} \tilde{q}_i^{\,\text{THINC}}(x) & \text{if } \delta < C_i < 1 - \delta,\ (\bar{q}_{i+1}-\bar{q}_i)(\bar{q}_i-\bar{q}_{i-1}) > 0,\ \text{TBV}^{\text{THINC}}_i < \text{TBV}^{\text{MUSCL}}_i \\ \tilde{q}_i^{\,\text{MUSCL}}(x) & \text{otherwise} \end{cases}

$C_i = (\bar{q}_i - \bar{q}_{\min})/\bar{q}_{\max}$ 是跳跃位置指标， $\delta \approx 10^{-4}$ 。前两个条件检查"这里是否真有跳跃"，第三个条件才是正式判定。

五方程模型上的一致性#

这篇论文不止是一次重构对比，关键就在这里。五方程模型同时求解体积分数 $\alpha_1$ 、相密度 $\alpha_1\rho_1, \alpha_2\rho_2$ 、动量 $\rho u$ 与总能 $E$ 。同一条 BVD 规则被应用到这五个变量上。

之前的方法只对体积分数做锐化，再对其他变量打补丁来抑制界面处的压力振荡。BVD 强制五个变量在每格遵循同一判定 —— 该格被判为"界面"，五个变量都用 THINC；判为"光滑"，五个变量都用 MUSCL。变量一致性自然落地，不需要单独的 anti-diffusion 步骤。

NumPy 中看 BVD 决策树#

import numpy as np
 
def muscl_minmod_edges(q):
    """ 接收 (q[i-1], q[i], q[i+1])，返回 MUSCL 左右端点值 """
    qm, q0, qp = q
    a, b = q0 - qm, qp - q0
    if a * b <= 0.0:
        slope = 0.0
    else:
        slope = np.sign(a) * min(abs(a), abs(b))
    return q0 - 0.5 * slope, q0 + 0.5 * slope  # qL, qR
 
def thinc_jump_edges(q, beta=1.6, eps=1e-20):
    """ tanh 拟合跳跃: THINC 左右端点值 """
    qm, q0, qp = q
    qmin = min(qm, qp)
    qmax = max(qm, qp) - qmin
    if qmax < 1e-12:
        return q0, q0
    theta = np.sign(qp - qm)
    C = (q0 - qmin + eps) / (qmax + eps)
    if C <= 1e-6 or C >= 1 - 1e-6:
        return q0, q0
    B = np.exp(theta * beta * (2.0 * C - 1.0))
    A = (B / np.cosh(beta) - 1.0) / np.tanh(beta)
    qR = qmin + 0.5 * qmax * (1.0 + theta * A)
    num = 1.0 + theta * A * np.tanh(beta) + theta * np.tanh(beta)
    den = 1.0 + A * np.tanh(beta)
    qL = qmin + 0.5 * qmax * (num / den)
    return qL, qR
 
def bvd_decide(q_window, beta=1.6, delta=1e-4):
    """
    q_window: 长度为 5 的数组 [q[i-2], q[i-1], q[i], q[i+1], q[i+2]]
    返回: (qL_i, qR_i, picked_thinc)
    """
    qm2, qm1, q0, qp1, qp2 = q_window
    monotone = (qp1 - q0) * (q0 - qm1) > 0.0
    qmin = min(qm1, qp1)
    qmax = max(qm1, qp1) - qmin
    C = 0.5 if qmax < 1e-12 else (q0 - qmin) / qmax
    eligible = monotone and (delta < C < 1.0 - delta)
 
    qL_M, qR_M = muscl_minmod_edges([qm1, q0, qp1])
    if not eligible:
        return qL_M, qR_M, False
 
    qL_T, qR_T = thinc_jump_edges([qm1, q0, qp1], beta=beta)
    # 邻格固定为 MUSCL 时的 TBV 比较
    _, qR_left = muscl_minmod_edges([qm2, qm1, q0])
    qL_right, _ = muscl_minmod_edges([q0, qp1, qp2])
    tbv_M = abs(qR_left - qL_M) + abs(qR_M - qL_right)
    tbv_T = abs(qR_left - qL_T) + abs(qR_T - qL_right)
    if tbv_T < tbv_M:
        return qL_T, qR_T, True
    return qL_M, qR_M, False
 
# 测试: 光滑 + 跳跃混合剖面
N = 80
x = (np.arange(N) + 0.5) / N
q = np.where((x >= 0.6) & (x <= 0.8), 1.0, 0.0)
q += 0.5 * np.exp(-((x - 0.25) / 0.04) ** 2)
 
picks = 0
for i in range(2, N - 2):
    _, _, used_thinc = bvd_decide(q[i-2:i+3])
    picks += int(used_thinc)
print(f"THINC 采用率: {picks}/{N-4}")
# 期望值: 仅跳跃两侧的少数几格

在光滑高斯区，minmod 捕获曲率，BVD 倾向 MUSCL；THINC 只在方波两侧 2–3 格胜出。95% 以上的格保持 MUSCL —— BVD 不改变原本的剖面。

看每一格谁胜出#

下面的模拟可亲自调节。把一维标量平流用三种重构并行求解，每步在画布下方的色带上高亮 BVD 选用 THINC 的那些格。

N = 80β = 1.60CFL = 0.40

MUSCL THINC BVDt = 0.00THINC picked: 0% of cells

把 $\beta$ 从 1.4 拉到 2.2，THINC 选中的格里跳跃变得更陡，但光滑高斯两端会出现锯齿。1.6 附近是甜蜜点 —— 这就是论文采用它的原因。还要看色带：青色格只在方波周围少数几个亮起，光滑山从不触发 THINC。这正是 BVD 的自定位特性。

能信到什么程度#

收敛阶低于二阶。MUSCL 限制了光滑区的精度。把候选换成 WENO 可以恢复五阶，但本论文卖的是简洁。
TBV 定义假设邻格使用 MUSCL。相邻两格同时想用 THINC 的多选问题被简化掉了。一遍扫描的判定方式在论文示例中已足够好，后续工作做了推广。
多维情形下 $\beta$ 与界面法向有关。界面法向用 Young 算法估计，在拐角与折点法向估计劣化时，跳跃可能歪斜。
与高阶时间积分的耦合。论文使用 SSPRK3，每个 RK 段都要重做 BVD 判定。如果判定在 RK 段间抖动，界面会颤动 —— 这是后续研究中讨论的话题。

三行总结#

在每格同时构造光滑用 (MUSCL) 与跳跃用 (THINC) 两个候选，取边界跳跃和较小者。
同一条 BVD 规则适用于全部守恒量，变量一致性问题(伪压力振荡)自动消除。
没有 anti-diffusion 也没有人工压缩，界面仍维持单格厚度 —— 简洁就是优势。