不分裂方向求解对流 — Bell-Colella-Glaz 非分裂迎风格式

转完一圈后两种格式的质量误差都在舍入量级，峰值却是0.349和0.806。守恒检查通过，不等于圆盘还在。

1980 年代后期，劳伦斯伯克利实验室的 John Bell 和 Phillip Colella 在爆炸模拟中发现了一道奇怪的疤痕。一个相对网格轴线斜放的火焰面，随着时间推移不断沿对角线方向扩散、扭曲。元凶不是物理，而是离散化。逐方向求解的方向分裂(dimensional splitting)扭曲了对角对流。解法是 Bell-Colella-Glaz(BCG)非分裂二阶迎风格式。方向一拆，旋转圆盘就变成菱形。一项横向导数把这个洞补上。

把 CFL 从 0 推到 0.9，橙色点会贴回单元平均。limiter 设成 none，跳跃旁边的直线会冲过邻值。

CFL = 0.50

Cyan = limited PLM slope inside each cell · Amber dot = predicted left edge state at n+1/2 (advecting velocity u>0). Pick none near the jump to see the slope overshoot the neighbor values.

0.9 时斜率权重几乎消失。改回 minmod，那一下冲出就被削掉。

要处理的方程是最简单的守恒型对流方程。一个无散度(divergence-free)的速度场 $\mathbf{u}$ 携带标量 $s$ 。

\frac{\partial s}{\partial t} + \nabla\cdot(\mathbf{u}\,s) = 0

其中 $s$ 是浓度、温度之类的被动标量， $\mathbf{u}$ 是给定的速度场。用有限体积(FV)更新单元平均值，就得到下面的形式。

s_{ij}^{n+1} = s_{ij}^{n} - \frac{\Delta t}{\Delta x}\left(F^x_{i+\frac12,j} - F^x_{i-\frac12,j}\right) - \frac{\Delta t}{\Delta y}\left(F^y_{i,j+\frac12} - F^y_{i,j-\frac12}\right)

$F^x = u\,s^{n+1/2}$ 是面上的时间居中(time-centered)通量。全部精度最终都归结为：你能多准确地估计面上的值 $s^{n+1/2}$ 。

逐方向分裂会扭曲旋转

方向分裂把 2D 问题拆成先做 x 方向的 1D 更新，再做 y 方向的 1D 更新。实现很简单，一个 1D 求解器就够了。麻烦出在对角流动。

当标量以 $45^\circ$ 流动时，先解 x 再解 y 的顺序，与相反顺序会给出不同的答案。两者都漏掉了夹在两个 1D 步之间的横向耦合(cross term)。在旋转流动中，这个误差每一步都累积，圆盘逐渐扭曲成菱形。用 Strang 分裂交替顺序，也只是消去二阶误差的方向不对称，根本原因依旧存在。

BCG 的药方很简单：x 和 y 同时求解。在一步之内，把所有面值都预测到同一时刻 $t^{n+1/2}$ ，并把横向变化预先纳入这个预测。

把时间换成空间的预测子

核心思路是用 Taylor 展开把面值在时间和空间两个方向都推进半步。看单元 $i$ 右面上的左状态。

s^{n+1/2}_{i+\frac12,\,L} = s^n_{ij} + \frac{\Delta x}{2}\,\partial_x s + \frac{\Delta t}{2}\,\partial_t s

把时间导数原封不动留着没有用。代入对流方程 $\partial_t s = -u\,\partial_x s$ ，把时间导数换成空间导数(Cauchy-Kowalewski 技巧)。

s^{n+1/2}_{i+\frac12,\,L} = s^n_{ij} + \frac{1}{2}\left(1 - \frac{\Delta t}{\Delta x}\,u_{i+\frac12,j}\right)\Delta_x s_{ij}

$\Delta_x s_{ij}$ 是单元 $i$ 的限制斜率(limited slope)， $u_{i+1/2,j}$ 是面上的对流速度。括号里的 $(1 - \mathrm{CFL})$ 承载时间推进。CFL 越接近 1，左状态的贡献越小；为 0 时就退化为单纯的空间外推。

从邻居单元看到的右状态也对称地构造。

s^{n+1/2}_{i+\frac12,\,R} = s^n_{i+1,j} - \frac{1}{2}\left(1 + \frac{\Delta t}{\Delta x}\,u_{i+\frac12,j}\right)\Delta_x s_{i+1,j}

在两个候选中挑一个，靠的是 Godunov 的迎风规则。面速度的符号决定风向。

s^{n+1/2}_{i+\frac12} = \begin{cases} s_{i+\frac12,\,L}, & u_{i+\frac12} > 0 \\[2pt] s_{i+\frac12,\,R}, & u_{i+\frac12} < 0 \end{cases}

上面那张图里，CFL 到 0.9 时橙色点贴向平均，none 会让直线越过邻值。这个差别就是上面预测子里的时间推进系数。

横向导数把角落缝起来

方向分裂漏掉的耦合，BCG 把它直接塞进预测子。在 x 面的左状态上，减去半步的横向通量散度。

s^{n+1/2}_{i+\frac12,\,L} \;\mathrel{-}=\; \frac{\Delta t}{2}\,\frac{\partial (v\,s)}{\partial y}\bigg|_{ij}

$v$ 是 y 方向速度， $\partial_y(v s)$ 是单元 $ij$ 处的横向对流变化。要估计这一项，得先求出 y 面上的 1D 预测状态 $\hat s^y$ 。所以 BCG 分两趟：先在每个方向预测 1D 迎风状态，再把这个预测复用到横向修正，然后用同一条迎风规则选出最终面值。这种“角耦合(corner coupling)”能把对角对流无扭曲地搬运过去。

给斜率套上缰绳 — minmod 与 superbee#

预测子的精度来自斜率 $\Delta_x s$ 。直接用中心差分斜率是二阶精度，但在间断附近会产生振荡(overshoot)。于是用斜率限制器(slope limiter)套上缰绳。最保守的 minmod 取两个单侧差分中较小的一个，符号不同时取 0。

\Delta_x s_{ij} = \operatorname{minmod}\!\left(s_{ij}-s_{i-1,j},\; s_{i+1,j}-s_{ij}\right)

superbee 在保持同样单调性的同时，把斜率立得尽可能陡，把接触面拉得更锐利。在上面的 1D 演示里切换 minmod 和 superbee，可以看到 superbee 一侧的 PLM 直线更陡。

把同样的缰绳套在通量上的一维算法，是给 Lax–Wendroff 加上 TVD 限制器。BCG 把限制器套在单元斜率上，不是套在通量上。

用 numpy 转动的旋转圆盘#

现在来跑 2D 非分裂 BCG。玩具问题是刚体旋转(solid-body rotation)。在绕中心旋转的无散度速度场中浮起一个标量圆盘，在同一网格上比较一阶迎风与 BCG。

import numpy as np
 
N = 64
dx = 1.0 / N
omega = 2 * np.pi                      # 单位时间一圈
 
def face_velocities():
    j = np.arange(N)
    i = np.arange(N)
    u = -omega * ((j + 0.5) * dx - 0.5)        # x 面速度 (只依赖 y)
    u = np.broadcast_to(u[None, :], (N, N)).copy()
    v = omega * ((i + 0.5) * dx - 0.5)         # y 面速度 (只依赖 x)
    v = np.broadcast_to(v[:, None], (N, N)).copy()
    return u, v
 
def minmod_pair(a, b):                          # 保单调的斜率
    return np.where(a * b > 0, np.where(np.abs(a) < np.abs(b), a, b), 0.0)
 
def limited_slopes(s):
    slx = minmod_pair(s - np.roll(s, 1, 0), np.roll(s, -1, 0) - s)
    sly = minmod_pair(s - np.roll(s, 1, 1), np.roll(s, -1, 1) - s)
    return slx, sly
 
def donor_fluxes(s, u, v):                       # 一阶迎风 (施主单元)
    Fx = np.where(u > 0, u * np.roll(s, 1, 0), u * s)
    Fy = np.where(v > 0, v * np.roll(s, 1, 1), v * s)
    return Fx, Fy
 
def bcg_fluxes(s, u, v, dt):                     # 非分裂二阶预测子
    slx, sly = limited_slopes(s)
    # 1) 法向 1D 预测状态
    sLx = np.roll(s, 1, 0) + 0.5 * (1 - dt/dx * u) * np.roll(slx, 1, 0)
    sRx = s - 0.5 * (1 + dt/dx * u) * slx
    hatX = np.where(u > 0, sLx, np.where(u < 0, sRx, 0.5 * (sLx + sRx)))
    sBy = np.roll(s, 1, 1) + 0.5 * (1 - dt/dx * v) * np.roll(sly, 1, 1)
    sTy = s - 0.5 * (1 + dt/dx * v) * sly
    hatY = np.where(v > 0, sBy, np.where(v < 0, sTy, 0.5 * (sBy + sTy)))
    # 2) 横向(角)修正 —— 复用预测状态
    divVy = (np.roll(v * hatY, -1, 1) - v * hatY) / dx
    divUx = (np.roll(u * hatX, -1, 0) - u * hatX) / dx
    sLx -= 0.5 * dt * np.roll(divVy, 1, 0)
    sRx -= 0.5 * dt * divVy
    sBy -= 0.5 * dt * np.roll(divUx, 1, 1)
    sTy -= 0.5 * dt * divUx
    # 3) 用同一迎风规则选最终面值
    sX = np.where(u > 0, sLx, np.where(u < 0, sRx, 0.5 * (sLx + sRx)))
    sY = np.where(v > 0, sBy, np.where(v < 0, sTy, 0.5 * (sBy + sTy)))
    return u * sX, v * sY
 
def advance(s, u, v, dt, scheme):
    Fx, Fy = bcg_fluxes(s, u, v, dt) if scheme == 'bcg' else donor_fluxes(s, u, v)
    return s - dt * ((np.roll(Fx, -1, 0) - Fx) / dx + (np.roll(Fy, -1, 1) - Fy) / dx)
 
def make_disk():
    X, Y = np.meshgrid((np.arange(N)+0.5)*dx, (np.arange(N)+0.5)*dx, indexing='ij')
    return ((X - 0.5)**2 + (Y - 0.78)**2 < 0.13**2).astype(float)
 
u, v = face_velocities()
Umax = omega * 0.5 * np.sqrt(2)
dt = 0.6 * dx / Umax                    # CFL = 0.6
steps = int(round(1.0 / dt))            # 正好转一圈
 
for scheme in ('upwind1', 'bcg'):
    s = make_disk()
    mass0 = s.sum()
    for _ in range(steps):
        s = advance(s, u, v, dt, scheme)
    print(f"{scheme:8s}  peak={s.max():.3f}  "
          f"mass_err={abs(s.sum()-mass0)/mass0:.2e}")

输出类似下面这样。

upwind1   peak=0.349  mass_err=3.1e-16
bcg       peak=0.806  mass_err=2.8e-16

把格式从 1st-order upwind 拨到 BCG unsplit，峰值会从 0.349 走到 0.806，和输出里的数字是同一对。

CFL = 0.60

A scalar disk rotates rigidly about the center. Switch to 1st-order upwind and watch the disk melt into a ring within one revolution; BCG unsplit keeps its edge far longer.

停在 1st-order upwind，圆盘转不完一圈就摊成一个环，边缘发糊。换成 BCG unsplit，边界留得久得多。CFL 提到 0.9，BCG 也会稍微变粗，和一阶的差距还在。

两种格式都把质量守恒到机器精度。守恒型通量差分保证了这一点。分开的是峰值。一阶迎风转一圈就把圆盘压到不足原高的一半，BCG 保住原始高度的 80%。

转两圈之后还剩下什么

用代码里的数把一圈换算出来是这样。N=64，间距就是1/64。角点速度是角速度2π乘0.5再乘根号2，约4.44。CFL取0.6时时间步长约0.00211，填满时间1需要474步。圆盘中心在(0.5, 0.78)，离旋转中心(0.5, 0.5)有0.28。那里的速度是2π乘0.28，约1.76。决定时间步长的不是1.76，而是角点的4.44。

半径0.13，直径0.26，大约16.6个格子。面积除以单元面积，大约217个单元带着标量。只看斜率权重，CFL为0.6时左侧是0.2，右侧是0.8，也就是(1-0.6)和(1+0.6)的一半。CFL放到注明的上限0.8，左侧权重降到0.1。输出里的峰值0.349和0.806，是这些权重加上横向项堆了474次之后的高度。

质量误差3.1e-16和2.8e-16都是舍入。分开的不是守恒，而是丢掉的高度0.651和0.194。一阶丢掉的高度大约是BCG的3.4倍。横向修正每步减去半个时间步长，约0.00106，再乘上横向散度。这一项一圈进去474次。数值耗散不由64个格子决定，而由一个格子里还留多少斜率决定。

移植时漏水的地方

CFL 比分裂格式更紧。 非分裂 2D 因为沿对角线传播信息，稳定极限比一维更小。实践中从 $\mathrm{CFL}\lesssim 0.8$ 起步。
速度场一定放在面上。 把单元中心速度插值到面会破坏无散度性，你以为守恒的质量会漏掉。MAC(交错)布置才是正解。
别漏掉横向项。 只实现法向预测，就只是普通的维度分裂 MUSCL。代码能跑，但旋转扭曲依旧。
limiter 在极值处归零。 在峰顶斜率被削为零，局部退回一阶。这正是峰值慢慢下沉的原因。想再保住一些，就考虑保极值的 limiter。

质量对了、圆盘却死掉的选择

峰值0.349是一阶迎风转完一圈剩下的高度。BCG是0.806。质量误差是3.1e-16和2.8e-16，只看守恒日志，两种都过。看圆盘，只剩一边还像圆盘。

对角旋转的标量要不要用一阶迎风，就由这个差距来定。接触面的高度可以丢掉，供体单元就够。要把高度留下，横向项必须进同一步。网格改成非结构之后，minmod不能原样搬走，极值的处理会变。Barth–Jespersen限制器写的是假极值怎样被邻居的最小值和最大值削掉。那是下一步要定的事。