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

纯单格台阶让一侧差分为零,符号条件通不过,THINC进不了边界变差比较。

东京工业大学 Xiao 课题组在2018年量过一个很少有人愿意公开的厚度。可压缩两相流,1000步。原来一格的跳跃铺成了八格。WENO也没拦住。单步耗散看起来很小,时间一攒,界面还是塌了。

他们拿来的重构是 MUSCL-THINC-BVD。每格放两个候选。边界跳跃更小的那个留下。

把 β 滑块从1.4推到2.2,跳跃格里的 tanh 变窄,高斯肩则以锯齿回应。

N = 80β = 1.60CFL = 0.40

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

看青色条是否只在方波边缘亮起。也看1.6附近这些条是否不再变多。

一页摘要

作者 / 期刊: Deng, Inaba, Xie, Shyue, Xiao. Journal of Computational Physics 371 (2018) 945–966.
问题: 五方程模型里的可压缩两相流,物质界面会随时间铺开。
方案: 每格同时做 MUSCL(光滑段)和 THINC(跳跃),留下边界跳跃更小的一侧。这就是 BVD(Boundary Variation Diminishing)。
不同处: 没有后处理 anti-diffusion,也没有人工压缩。体积分数和其余守恒量走同一条规则,变量之间的一致是跟着来的。

两个相互冲突的需求

可压缩两相流同时要两件事。光滑处要准,耗散要小。跳跃,也就是界面和激波,要单调,还要薄。一个函数很少能两头都做好。

MUSCL 在格内是直线。单调性守得住,斜率只到一阶。每过一次界面就削掉一点。1000步把一格磨成八格,就是这笔账。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 在两个差分同号时取绝对值较小者,异号时给0。

边界值这样读。

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$ 是跳跃方向。 $\beta$ 是厚度。 $\tilde x_i$ 是为保住单元平均而解出的跳跃中心。

$\beta$ 在1.4到2.0之间算稳。标准值是1.6。再大,跳跃被关进一格,看起来干净。进了光滑区就会出锯齿。把 tanh 接到激波一侧的记录在把 TENO 与 THINC 叠进同一模板的重构里。那边不挑选候选,而是用权重去混。

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$ 表示跳跃落在格内的什么位置。 $\delta$ 大约是 $10^{-4}$ 。

选择式的第一道门是差分的乘积。相邻两个单侧差分相乘,必须大于0。有一侧是0,乘积也是0,过不了门。0和1贴在一格上的台阶就是这个形状。跳跃前一格只有右侧差分,左侧差分为零。总边界变差还没量,就已经回到 MUSCL。

第二道门是位置尺度。 $\delta$ 在 $10^{-4}$ 附近。单元平均若落进相邻振幅的这个比例之内,THINC 就没有资格。这是用来切掉平坦峰顶的装置。

五方程模型上的一致性

重构比较没有停在这里。五方程模型一起解体积分数 $\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}")
# 期望值: 仅跳跃两侧的少数几格

代码末尾的注释写着采纳 4/76。紧挨上面的那道门,和这个数字,指的不是同一种分布。代码里的方波由条件语句取 0 或 1。一格跳跃使其中一个差分为零。注释里的四格,要等斜率在两三格上以同一符号爬升的斜坡才会出现。1000步之后的八格,已经是那种斜坡。厚度还停在一格时,这个分支不会去叫 THINC。铺开之后才会叫。

直线上则相反。相邻差分相等时,minmod 斜率就等于那个差分。左面值与左邻的右面值相遇,右侧也一样。MUSCL 一侧的总边界变差是0。THINC 的 tanh 在直线上加了一个弯。面上出现缝,变差变大。不等号留在 MUSCL 一侧。高斯比单元宽时,局部上接近这条直线。光滑段拒绝换候选,原因在这里。

光滑高斯在局部上接近直线,所以替换会被拒绝。方波还贴成一格时,门也不会开。

看每一格谁胜出

前面的平流把三种重构一起往前推,BVD 选中 THINC 的格子用条子涂上。把 $\beta$ 从1.4抬到2.2,只是那些格子里的跳跃变窄。太大时,高斯肩上能看见锯齿。1.6附近条子不再变多,那就是标准值停住的地方。方波略微铺开之后,边缘才出现青色。厚度还是一格时,条子是空的。

能信到什么程度

收敛到不了二阶。候选的一侧是 MUSCL,光滑段的精度就停在 MUSCL。换成 WENO 这类高阶候选,五阶也做得到。这篇论文选了简单的那一侧。

TBV 假定邻居是 MUSCL。两侧同时想走 THINC 的情形,一次扫描里被忽略。论文第4节的图在这个简化下仍然过得去。更一般的形式被留到了后面。

到了多维, $\beta$ 按面法向分开。界面法向用 Young 算法估计。尖角和折点上法向会偏,跳跃就向一侧歪。

时间积分是另一件事。论文用的是 SSPRK3。BVD 的判定必须在每个 Runge–Kutta 阶段重做。阶段之间判定摇摆,界面就会抖。后来的工作又把这个抖动捡了起来。

单格跳跃上选择停住的地方

界面还贴成一格时,这个分支不会去叫 THINC。厚度铺到两三格之后,差分乘积才变成正的。把界面留在一格里,和只有铺开之后才许选择,这两件事彼此错开。

这个错位会不会漏到时间步上限,是下一问。界面保持得很薄,平流的上限还在。毛细项拿掉之后 CFL 0.05 仍然留着的界面平流问的就是那个上限在哪里。若 BVD 真把跳跃关在一格里,推动这一格的时间步能大到什么程度?