Lax–Wendroff 在阶跃处颤抖 — 用 Flux Limiter 构建二阶 TVD 对流格式

φ(r)落在Sweby区域内，台阶就不会长出新极值。光滑波的精度阶仍按限制器分开。

一次测试里,阶跃的前后冒出了奇怪的凸包。我在用 Lax–Wendroff 求解线性对流方程(linear advection,恒定速度场下的纯平移),初始条件是简单的矩形脉冲。按说二阶格式应该比一阶迎风更锐利,结果却是阶跃前出现下冲、后出现上冲。日志里 CFL = 0.5,稳定性条件明显满足。

问题不是算法不稳定,而是算法不单调。

限制器按钮从 lax-wendroff 开始。CFL 滑块放在 0.5，橙色曲线会在台阶前后鼓出凸包；换成 superbee，凸包就退下去。

CFL0.50steps: 0

橙色是数值解，青色虚线是精确解。minmod 两侧都钝，没有凸包。superbee 甚至让台阶看起来比原来更锐。CFL 推到 0.9 附近，没有凸包的限制器也会摊得更开。稳定条件里面，数值扩散照样咬在一起。

振荡的来源 — Godunov 的墙#

一维线性对流方程为

\partial_t q + u\,\partial_x q = 0

若速度 $u$ 恒定,解就是初始形状原样向右平移。数值上我们在单元界面上定义通量 $f_{i-1/2}$ ,然后写

q_i^{n+1} = q_i^n - \frac{\Delta t}{\Delta x}\bigl(f_{i+1/2} - f_{i-1/2}\bigr)

这叫 flux-conserving form,能把质量、动量、能量等守恒量保持到机器精度。Lax–Wendroff 是在此基础上加了二阶 Taylor 修正的线性格式。

但 Godunov 在 1959 年证明的定理很残酷:对线性常系数对流而言,不存在同时保持单调性(monotonicity)且具有二阶以上精度的线性格式。 两者必须舍一。Lax–Wendroff 选了二阶精度,丢掉了单调性 — 所以阶跃处会振荡。

TVD 条件:Harten 不等式#

Harten(1983)给出的出路是"总变差递减"的格式。

\mathrm{TV}(q) = \sum_i |q_{i+1} - q_i|, \qquad \mathrm{TV}(q^{n+1}) \le \mathrm{TV}(q^n)

左边衡量数值解的总波动量。每一步只要这个值不增加,就不会出现新的极值(上冲、下冲)。满足此条件的格式叫 TVD。

要绕过 Godunov 的墙,格式必须非线性化 — 算法要根据解的光滑程度改变形状。光滑区用二阶,陡梯度附近用一阶迎风。完成这个切换的就是 flux limiter。

一行读懂 Flux Limiter#

Lax–Wendroff 的通量是一阶迎风加上修正项:

f_{i-1/2} = u\,q_{i-1} + \tfrac{1}{2}\,|u|\!\left(1 - \left|\tfrac{u\Delta t}{\Delta x}\right|\right) \varphi(r)\,(q_i - q_{i-1})

其中 $u$ 是对流速度, $\Delta x$ 是单元宽度, $\varphi(r)$ 是限制器, $r$ 是前后斜率之比:

r_{i-1/2} = \frac{q_{i-1} - q_{i-2}}{q_i - q_{i-1}}

光滑时 $r \approx 1$ ,极值附近 $r < 0$ 或 $r \gg 1$ 。限制器根据这个 $r$ 决定修正的强度。

几个特例:

$\varphi(r) = 0$ :donor-cell(一阶迎风)
$\varphi(r) = 1$ :Lax–Wendroff
$\varphi(r) = r$ :Beam–Warming

在 Sweby 图中,只要满足 $0 \le \varphi(r) \le \min(2r,\,2)$ 就能保证 TVD。落在这个梯形区域内的四个代表性限制器对比如下。把 $\varphi(r) = 0$ 关掉的地方就是一阶迎风。那个迎风用变钝换掉振荡，换法写在 FTCS 开始发抖时的交换。

限制器	定义	特点
minmod	$\max(0,\min(1,r))$	最谨慎。能保留阶跃,光滑处略显迟钝
superbee	$\max(0,\min(2r,1),\min(r,2))$	最锐利。会把光滑波削成阶梯
van Leer	$(r+	r
MC	$\max(0,\min(\tfrac{1+r}{2},2,2r))$	默认值,最稳妥

用 Python 直接比较#

基于 roll 的向量化不到 40 行。200 单元, $\mathrm{CFL} = 0.5$ ,500 步(周期边界,波走完约 2.5 圈域)。

import numpy as np
 
def limiter_phi(name, r):
    r = np.where(np.isfinite(r), r, 0.0)
    if name == 'donor':    return np.zeros_like(r)
    if name == 'lw':       return np.ones_like(r)
    if name == 'minmod':   return np.maximum(0, np.minimum(1, r))
    if name == 'superbee':
        return np.maximum.reduce([np.zeros_like(r),
                                  np.minimum(2*r, 1),
                                  np.minimum(r, 2)])
    if name == 'vanleer':  return (r + np.abs(r)) / (1 + np.abs(r) + 1e-30)
    if name == 'mc':
        return np.maximum(0, np.minimum.reduce([0.5*(1+r),
                                                2*np.ones_like(r),
                                                2*r]))
 
def advect_limited_fvm(q, cfl, name):
    """一维有限体积对流,周期边界,u = +1。"""
    dq  = q - np.roll(q, 1)                  # 界面 i-1/2: q_i - q_{i-1}
    r   = np.roll(dq, 1) / (dq + 1e-30)      # (q_{i-1}-q_{i-2}) / (q_i-q_{i-1})
    phi = limiter_phi(name, r)
    flx = np.roll(q, 1) + 0.5 * phi * (1 - cfl) * dq
    return q - cfl * (np.roll(flx, -1) - flx)
 
N, cfl, steps = 200, 0.5, 500
x  = np.arange(N) / N
q0 = ((x > 0.2) & (x < 0.4)).astype(float)
for name in ['lw', 'minmod', 'superbee', 'vanleer', 'mc']:
    q = q0.copy()
    for _ in range(steps):
        q = advect_limited_fvm(q, cfl, name)
    l1 = np.abs(q - q0).sum() / N
    print(f'{name:<9s}  L1={l1:.4f}  overshoot={q.max()-1:.3f}  undershoot={q.min():.3f}')

典型结果:

限制器	$L_1$ 误差	上冲	下冲
lax-wendroff	0.056	+0.181	−0.133
minmod	0.081	0.000	0.000
superbee	0.042	+0.002	0.000
van Leer	0.061	0.000	0.000
MC	0.051	0.000	0.000

只有 lax-wendroff 出现接近 20% 的上冲,其余全为 0。误差排名也有意思 — superbee 的 $L_1$ 最小,说明它最能保留阶跃锐度,代价是把光滑的 sine 波削成阶梯。

另做光滑波测试( $q_0 = \sin(2\pi x)$ ,500 步),minmod 的收敛阶约 1.5,van Leer/MC 约 2.0,superbee 约 1.3。"二阶精度"的标签会在每次极值附近限制器激活时局部回落到一阶。

r=1 时四个限制器给出同一个 φ#

把 $r=1$ 代进表里的定义，minmod、superbee、van Leer、MC 都得到 1。 $\varphi=1$ 就是 Lax–Wendroff。光滑的中段上，四个限制器用同一份修正。差别出现在 $r$ 离开 1 的时候。

再代入 $r=2$ 。minmod 停在 1。superbee 贴在 2，也就是 Sweby 区域的上沿。van Leer 是 $4/3$ 。MC 是 1.5。superbee 把台阶的 $L_1$ 压到 0.042、成为表里最小的一个，靠的就是这条上沿。同一次选择把 sine 测试的精度阶降到大约 1.3。van Leer 和 MC 留在大约 2.0，minmod 大约 1.5。 $r$ 小于或等于 0 的极值上，四个都掉到 0，变成 donor-cell。

代码里的修正系数也已经定了。速度为 1、 $\mathrm{CFL}=0.5$ 时，跳变量前面的系数是 0.25。就算 $\varphi=1$ ，也只加上邻居跳跃的四分之一。500 步走出的距离是 $0.5\times 500=250$ 个单元。网格 200，周期边界上转 1.25 圈。上冲 +0.181 是单位波高的 18%。下冲 −0.133 是地板沉到负数的深度。minmod 的 $L_1$ 0.081 大约是 superbee 0.042 的 1.9 倍。抹掉凸包的代价，是台阶多摊开这么多。

工程清单

激波捕获优先用 minmod 或 van Leer。零振荡是首要目标。
LES / 大规模湍流里 superbee 危险。它会把光滑湍流结构削成阶梯,扭曲能谱。折中选 van Leer 或 MC。
边界单元因为模板不足临时用一阶迎风,会让整体精度降到一阶。要么在边界上保住限制器的模板,要么使用 ghost cell。
OpenFOAM 的 limitedLinear 1 是 van Leer 家族的单参数限制器。系数 0 = 一阶迎风,1.0 = 使用完整 TVD 区域。
验证收敛阶一定要用光滑解。在含间断的测试上测 $L_1$ 斜率,是间断本身在拉低阶数,不是限制器。

台阶上用 superbee，sine 会掉到 1.3 阶#

表里 superbee 的 $L_1$ 是 0.042，最小。上冲 +0.002，在台阶上几乎看不见。把同一个限制器放到 sine 上，精度阶掉到大约 1.3。van Leer 和 MC 留在大约 2.0。保住间断的选择，削掉的是光滑波的阶。

标量如果是密度，lax-wendroff 的下冲 −0.133 就是负密度的入口。堵住这个口的限制器写在正值漏出去的地方。被动标量的台阶，minmod 或 van Leer 就够。密度的话，再开下一扇门。