非结构网格上二阶重构制造的虚假极值 — Barth–Jespersen 与 Venkatakrishnan 限制器

求解器本来好好的。二阶一开，密度就掉成了负值。激波前面的面值高出了单元平均。同一套网格上的一阶迎风没有坏。把精度抬上去的那一侧，把计算杀掉了。非结构单元里的一条直线，会造出任何邻居都没有的极值。Barth–Jespersen 把这个极值削掉。Venkatakrishnan 把削过的折点磨圆。

重构制造出本不存在的极值

有限体积法（FVM，以单元平均为未知量的离散化）里，二阶就是单元内部的一条直线。

$u_i$ 是单元 $i$ 的平均。$\nabla u_i$ 是重构梯度，$\mathbf{x}_f$ 是面中心，$\Phi_i\in[0,1]$ 是限制器。

$\Phi_i=1$ 时限制器是关着的。梯度由最小二乘或 Green–Gauss 做成，它平均了邻居。不连续旁边，这个平均太陡。直线伸到面上，就出现邻居从未持有的新最大或最小值。

把梯度 gain 拉高，无限制重构（Φ=1）会冲破灰色带，也就是邻居的最大和最小，并变成红色。

线留在带内时，该单元的 Φ 仍是 1。只有跳跃旁边的线段会先变短。在同样的 gain 下把限制器换成 Barth–Jespersen，线就不再冲出带子。

虚假极值会污染下一步的通量。振荡长到最后，就是负密度和负能量。重构出的面值不得超出邻居已经持有的范围。 这是局部极值原理（LMP）。$\Phi_i$ 就是为守住这个范围而削减梯度的标量。

Barth–Jespersen：用邻居的最大最小关起来#

Barth 和 Jespersen（1989）先把允许范围固定下来。那是单元 $i$ 和邻居一起给出的最大与最小。

$N(i)$ 是单元 $i$ 的面邻居集合。每个面上先看无限制增量 $u_f-u_i$，再求该面的系数 $\phi_f$。

一个单元的限制器，取所有面里最保守的那个值。

落在范围内时 $\phi_f=1$，二阶还在。想越界的面只留下刚好碰到边界的斜率。极值旁边 $\Phi_i\to 0$，格式退回一阶迎风。

只看正增量时，$\phi_f=\min(1,(u_i^{\max}-u_i)/(u_f-u_i))$。分子一旦比分母小，面值就贴在 $u_i^{\max}$ 上停住。上图里红线在灰色边界折断，就是这个比值。

第二个陷阱：不可微

Barth–Jespersen 守住单调性。放进定常求解器，残差会在某一高度不再下降，只在那里抖。原因是 $\min$ 和除法。$\phi_f$ 对解不可微。哪个面给出最小值，会在迭代之间换来换去。限制器一开一关，残差就被弹回去。隐式（implicit）求解器的雅可比矩阵讨厌这个不连续。

光滑处的 Φ 保持 1 就够了。它只该在真正的不连续上启动，开关之间的边界也该是光滑的。下一步就是把尖锐的折点磨圆。

Venkatakrishnan：平滑地削减#

Venkatakrishnan（1993）用有理函数换掉了 $\min(1,y)$ 的折点。源文档提到的 Michalak–Ollivier-Gooch 限制器也在同一族里。每个面的系数写成

$\Delta_-=u_f-u_i$ 是无限制增量。$\Delta_+$ 按增量的符号取 $u_i^{\max}-u_i$ 或 $u_i^{\min}-u_i$。真正用手选的量是 $\epsilon^2$。

$K$ 是调参，$\Delta x$ 是网格尺度。变化比 $\epsilon$ 小时，$\phi_f\to 1$，限制器关闭。$K$ 变大，阈值升高，限制器松开的区域变宽。过大就会放过激波旁边的振荡。把 $K$ 送到 0，它就回到 Barth–Jespersen。

代码把单元尺寸当成 1，用 $\epsilon^2=K^3$ 来算。默认 $K=0.3$ 时，$\epsilon^2=0.027$。在 $\epsilon=0$ 且 $\Delta_+=\Delta_-=\Delta$ 的地方，也就是重构刚好碰到允许界限时，有理函数等于 $(3\Delta^2)/(4\Delta^2)=3/4$。同一位置上 Barth–Jespersen 给出的是 $\phi_f=1$。Venkatakrishnan 在碰到界限之前先丢掉四分之一的斜率。角点因此早一步变钝。

习惯上试的区间是 $K\approx 0.1$ 到 $K\approx 5$。把两端代进代码里的 $\epsilon^2=K^3$，得到 0.001 和 125。阈值之比是 12 万 5 千。$K$ 只挪一点也要做立方，于是被当成光滑而放过的单元会猛增。反过来把 $K$ 放到 0.1 以下，$\epsilon^2$ 实际上就是 0，Barth–Jespersen 的折点会原样回来。

在上图里打开 Venkatakrishnan 并移动 $K$，Φ 会没有台阶地滑过去。

Python — 用标量对流让三种限制器赛跑#

周期边界。方波和高斯峰一起对流（advection）。无限制（Fromm）、Barth–Jespersen、Venkatakrishnan 并排跑，看最后的最小值。低于零，就说明挖出了那么深的虚假波谷。

单元斜率是 $s=(u_{i+1}-u_{i-1})/2$。送到面上的增量是 $0.5s$，所以左右邻居相差 1 时 $s=1/2$，面增量是 0.25。这 0.25 比允许宽度大，Φ 就降到 1 以下。一维中心差分只看左右邻居。非结构单元上，给同样的梯度配权重是另一件事，移动最小二乘按单元重新拟合的那种。限制器接在梯度已经做好之后。

import numpy as np
 
NX, A, CFL = 200, 1.0, 0.4
dx = 1.0 / NX
dt = CFL * dx / A
 
def init_profile():
    x = (np.arange(NX) + 0.5) / NX
    u = np.where((x > 0.1) & (x < 0.3), 1.0, 0.0)   # 方波
    u += np.exp(-((x - 0.65) / 0.06) ** 2)          # 高斯峰
    return u
 
def cell_slope(u):
    return (np.roll(u, -1) - np.roll(u, 1)) / 2.0    # 中心斜率（Fromm）
 
def barth_jespersen_phi(u, s):
    up, um = np.roll(u, -1), np.roll(u, 1)
    umax = np.maximum(u, np.maximum(up, um))
    umin = np.minimum(u, np.minimum(up, um))
    phi = np.ones_like(u)
    for du in (0.5 * s, -0.5 * s):                   # 两个面
        f = np.ones_like(u)
        pos, neg = du > 1e-12, du < -1e-12
        f[pos] = np.minimum(1.0, (umax[pos] - u[pos]) / du[pos])
        f[neg] = np.minimum(1.0, (umin[neg] - u[neg]) / du[neg])
        phi = np.minimum(phi, f)
    return np.clip(phi, 0.0, 1.0)
 
def venkatakrishnan_phi(u, s, K=0.3):
    up, um = np.roll(u, -1), np.roll(u, 1)
    umax = np.maximum(u, np.maximum(up, um))
    umin = np.minimum(u, np.minimum(up, um))
    eps2 = K ** 3                                    # (K*h)^3，单元尺寸 h=1 单位
    phi = np.ones_like(u)
    for du in (0.5 * s, -0.5 * s):
        d = np.where(du > 0, umax - u, umin - u)
        num = (d * d + eps2) * du + 2 * du * du * d
        den = d * d + 2 * du * du + d * du + eps2
        f = np.where(np.abs(du) < 1e-12, 1.0, num / (du * den))
        phi = np.minimum(phi, f)
    return np.clip(phi, 0.0, 1.0)
 
def muscl_rhs(u, limiter):
    s = cell_slope(u)
    phi = limiter(u, s) if limiter else np.ones_like(u)
    uL = u + 0.5 * phi * s              # 面 i+1/2 左状态
    flux = A * uL                       # 迎风通量 (a > 0)
    return -(flux - np.roll(flux, 1)) / dx
 
def advance_muscl(u, limiter):         # SSP-RK2
    k1 = muscl_rhs(u, limiter)
    u1 = u + dt * k1
    k2 = muscl_rhs(u1, limiter)
    return 0.5 * (u + u1 + dt * k2)
 
def run_limiter_race(steps=160):
    fields = {"none": init_profile(), "bj": init_profile(), "venk": init_profile()}
    lims = {"none": None, "bj": barth_jespersen_phi, "venk": venkatakrishnan_phi}
    for _ in range(steps):
        for k in fields:
            fields[k] = advance_muscl(fields[k], lims[k])
    for k, u in fields.items():
        print(f"{k:5s}  min={u.min():+.4f}  max={u.max():.4f}")
 
run_limiter_race()

输出是这样。

none   min=-0.0417  max=1.0231
bj     min=+0.0000  max=1.0000
venk   min=-0.0004  max=1.0009

无限制的最小值 -0.0417，是高度为 1 的方波被挖到零以下的深度。最大值 1.0231 是越过同一高度的量。Barth–Jespersen 把最小和最大钉在 $[0,1]$。Venkatakrishnan 的最小值 -0.0004 大约是那个波谷的 104 分之 1。超过 1 的宽度从 0.0231 缩到 0.0009。

160 步的移动距离是 $A \cdot dt \cdot 160$。$dx=1/200$ 且 CFL 为 0.4 时，$dt=0.002$，距离就是 0.32。200 个单元里的 64 个。方波区间 $(0.1, 0.3)$ 移到 $(0.42, 0.62)$，高斯中心 0.65 到达 0.97。两者速度相同，相对位置保持不变。输出里的最小和最大不是两列波撞出来的。每个拐角各自产生了一个数。方波长度 0.2 是 40 个单元，高斯尺度 0.06 是 12 个单元。决定最小和最大的，就是拐角附近的那几格。

只看方波后沿就够了。无限制（红）在零以下长出涟漪。两种限制器都保持单调。Venkatakrishnan（黄）在拐角处比 Barth–Jespersen（青）略钝。这和 $\epsilon=0$ 时斜率只留下 3/4 的代数结果是同一个方向。代码里的 $K=0.3$ 并没有把折点抹掉，只是把超出量压到 0.0009 这一档。

现场打开限制器时

邻居怎么捆，决定允许范围。Barth–Jespersen 的 $\max/\min$ 取面邻居还是取顶点（vertex）邻居，结果不一样。在非结构网格上用顶点邻居，也就是共享顶点的全部单元，方向偏差会减小，限制器更均衡。

接下来是 $K$。太小，收敛会像 Barth–Jespersen 那样停住。太大，激波上的振荡就会漏过去。$\epsilon^2$ 带着 $\Delta x^3$，网格一加密，同一个 $K$ 也会让限制器更频繁地启动。

在向量系统（Euler 方程）里，对每个守恒变量分开限制，分量之间对不齐，就会生出新的振荡。先投影到特征变量（characteristic variable）再限制更稳妥。代价更高。

残差停住时，先不要怪梯度

二阶求解器死在负密度上，先看面值是不是已经跑出邻居范围。如果是，问题在于造出 Φ 的邻居集合，不在梯度程序。只用面邻居，会漏掉对角方向的极值。重构还会从那条缝里漏出去。

定常计算里如果只剩残差在原地踏步，就把 $K$ 打开。$K=0$ 就是 Barth–Jespersen 的折点，雅可比矩阵会在每次迭代上跳一下。跳跃处 Φ 降到 0，那个单元的阶数就贴在一阶上。光滑区域还能留在二阶附近吗？能的话，是 $\epsilon^2$ 把那片区域放开了。同一次计算里把这个分叉一起量出来会得到什么数，还留在把激波和光滑涡放进同一次阶数测量的记录里。