[论文评论] 单级即达二阶精度 — 面向Kapila多相流的GRP求解器

界面上同时求出时间导数，不必再加 Runge-Kutta 级，时间二阶就会出现。

时间二阶一度被记成 Runge-Kutta 要跑好几级。这篇论文用一级就做出那个阶。工具是 GRP（Generalized Riemann Problem）求解器，界面上同时交回数值和时间导数。

论文: Tuowei Chen, Zhifang Du, Generalized Riemann Problem Method for the Kapila Model of Compressible Multiphase Flows, arXiv:2310.08241 (2023).

关掉 GRP 勾选，斜率滑条会锁住，界面值在时间上冻住。把 CFL 调高，只有写进通量的半步位置往前走。

advection speed a: 1.00left slope uₓ: 0.50CFL: 0.80GRP (off = Godunov)

斜率放到 0，青色点停下，和 Godunov 叠在一起。把对流速度 $a$ 的符号反过来，提供数值的单元就换了。橙色圆圈是 $\Delta t/2$ 那个时刻。

单级二阶 — GRP要解决的问题#

有限体积的一步由界面通量 $\mathbf{F}_{i+1/2}$ 决定。Godunov 在界面上解黎曼问题，得到一个数。那个数对时间是常数。想要时间二阶，就得把同一个空间算子像 Runge-Kutta 那样重复好几次。

GRP 在半步时刻 $t^{n+1/2}$ 上评估通量。

\mathbf{U}(\mathbf{x}_{i+1/2}, t^{n+1/2}) = \mathbf{U}^{n,*}_{i+1/2} + \frac{\Delta t}{2}\left(\frac{\partial \mathbf{U}}{\partial t}\right)^{n,*}_{i+1/2}

$\mathbf{U}^{n,*}$ 是黎曼解（界面值）， $(\partial \mathbf{U}/\partial t)^{n,*}$ 是该界面的瞬时时间导数。两项按 Taylor 展开加起来，界面值往未来走 $\Delta t/2$ 。评估只有一次。

从分段常数到分段线性 — 黎曼与GRP之别#

黎曼问题的界面两侧是分段常数。GRP 问题的两侧是分段线性。从初始间断伸出的波不再是直线，求解器把黎曼解和时间导数一起交回来。

标量对流 $u_t + a u_x = 0$ 里，把单元 $i$ 的斜率记成 $\sigma_i$ 。 $a>0$ 时，界面上的黎曼解是左单元的外推。

U^{*}_{i+1/2} = u_i + \sigma_i \frac{\Delta x}{2}

$u_i$ 是单元平均， $\sigma_i$ 是单元内斜率， $\Delta x$ 是单元宽度。前面的青色点就是这个式子。斜率取 0、退回 Godunov，也从这个式子里出来。

Lax-Wendroff手续 — 把时间导数换成空间变化#

时间导数从控制方程里取出来。

\frac{\partial u}{\partial t} = -a \frac{\partial u}{\partial x} = -a\,\sigma_i

左边的时间导数换成右边的空间导数。 $\sigma_i$ 在重构里已经知道。半步界面值是

U^{n+1/2}_{i+1/2} = U^{*}_{i+1/2} - a\,\sigma_i \frac{\Delta t}{2}

通量是 $F = a\,U^{n+1/2}$ 。Euler 或 Kapila 里推广成 $\partial u/\partial t = -A\,\partial u/\partial x$ ， $A$ 是雅可比。论文用这个手续把声波线性化，近似地解 GRP。

Python — 用GRP标量求解器测量收敛阶#

把光滑初值 $u_0=\sin(2\pi x)$ 对流一个周期，比较 Godunov 和 GRP 的 $L_1$ 误差阶。区域是周期边界上的线性对流。

import numpy as np
 
# 1D 线性对流  u_t + a u_x = 0,  周期域 [0,1]
# 一阶 Godunov vs 单级二阶 GRP(Lax-Wendroff)对比
 
def cell_gradient(u, dx):
    # 为测量光滑解的阶,使用中心(无限制)斜率
    return (np.roll(u, -1) - np.roll(u, 1)) / (2.0 * dx)
 
def godunov_step(u, a, dx, dt):
    # 分段常数: 迎风界面值 (a > 0)
    flux = a * u                       # F_{i+1/2} = a u_i
    return u - dt / dx * (flux - np.roll(flux, 1))
 
def grp_step(u, a, dx, dt):
    g = cell_gradient(u, dx)           # 空间斜率 sigma_i
    u_star = u + g * dx / 2.0          # 界面黎曼解 U*
    u_half = u_star - a * g * dt / 2.0 # Lax-Wendroff 半步 U^{n+1/2}
    flux = a * u_half                  # F_{i+1/2}^{n+1/2}
    return u - dt / dx * (flux - np.roll(flux, 1))
 
def l1_error(num, exact, dx):
    return np.sum(np.abs(num - exact)) * dx
 
def convergence_study(a=1.0, cfl=0.4, t_end=1.0):
    print(f"{'N':>5} {'Godunov':>16} {'GRP':>16}")
    prev = {}
    for N in [40, 80, 160, 320]:
        dx = 1.0 / N
        x = (np.arange(N) + 0.5) * dx
        u0 = np.sin(2 * np.pi * x)     # 光滑初始条件
        dt = cfl * dx / abs(a)
        nsteps = int(round(t_end / dt))
        dt = t_end / nsteps            # 精确落在 t_end
        ug = u0.copy(); ur = u0.copy()
        for _ in range(nsteps):
            ug = godunov_step(ug, a, dx, dt)
            ur = grp_step(ur, a, dx, dt)
        exact = np.sin(2 * np.pi * (x - a * t_end))
        eg = l1_error(ug, exact, dx)
        er = l1_error(ur, exact, dx)
        og = f"{np.log2(prev['g']/eg):.2f}" if 'g' in prev else "  - "
        orr = f"{np.log2(prev['r']/er):.2f}" if 'r' in prev else "  - "
        print(f"{N:>5} {eg:.2e}({og}) {er:.2e}({orr})")
        prev = {'g': eg, 'r': er}
 
convergence_study()

输出如下。

    N          Godunov              GRP
   40 1.63e-01(  - ) 1.31e-03(  - )
   80 8.76e-02(0.90) 2.69e-04(2.28)
  160 4.54e-02(0.95) 6.32e-05(2.09)
  320 2.31e-02(0.97) 1.55e-05(2.03)

把 $\lambda\Delta t$ 滑条推过 1.6，显式 Euler 会不会冲出 $[0,1]$ 区间就会变。再拨初始体积分数和平衡值，看 Crank-Nicolson 是不是仍留在区间里。

stiffness λΔt: 1.60α₀ (initial): 0.85α_eq (target): 0.10

$\lambda\Delta t$ 到 3 时，青色落到平衡，红色还在振荡。表里的阶 2.03 属于光滑对流。这个滑条说的是，半步时间导数交给源项积分之后的事。

N=40 的误差比是 $1.63\times 10^{-1}/1.31\times 10^{-3}$ ，大约 124 倍。同一网格上小一百倍，指的就是这张网格。N=320 上是 $2.31\times 10^{-2}/1.55\times 10^{-5}$ ，大约 1490 倍。网格从 40 加密到 320，因子是 8，Godunov 的误差只降到 1/7.1，GRP 降到 1/85。8 的平方是 64。印出来的阶从 2.28 到 2.09 再到 2.03，往下坐的时候稍微超过了 64。Godunov 一侧的 0.90、0.95、0.97 也朝 1 贴过去。

N=80 的比值是 $8.76\times 10^{-2}/2.69\times 10^{-4}$ ，约 326 倍。N=160 是 $4.54\times 10^{-2}/6.32\times 10^{-5}$ ，约 718 倍。网格每加密一倍，误差比按 124、326、718、1490 走，每一步略超过两倍。

Kapila模型的刚性体积分数源项#

Kapila 五方程模型（BN 七方程的约化）的体积分数方程带着刚性源项。

\frac{\partial \alpha_1}{\partial t} + \mathbf{u}\cdot\nabla\alpha_1 = \alpha_1 K\,\nabla\cdot\mathbf{u}

$\alpha_1$ 是相 1 的体积分数， $K$ 是从 Wood 声速取出的压缩、膨胀系数。这个 $K$ 怎样随体积分数折下去，和混合物里 Wood 声速跌下去的那笔账是同一处。源项可以把 $\alpha_1$ 推出 $[0,1]$ 。显式积分在刚性变大时越过边界，然后发散。

GRP 的时间导数让界面值在新时刻 $t^{n+1}$ 成为已知。用 Gauss-Green 可以在那个时刻估计 $\nabla\cdot\mathbf{u}$ 的单元平均。源项用 Crank-Nicolson（半隐式）积分，还可以写下保住体积分数边界的时间步限制。前面的 $\lambda\Delta t$ 曲线就是这个限制的模型。

复现时令人生疑之处

标量 GRP 给出了二阶，Kapila 系统比那粗糙。无限制的中心斜率在激波上振荡。阶数出现，是因为测试是光滑的。间断上需要 minmod 这类限制器，那时阶在局部掉回一阶。

Wood 声速对体积分数非单调。界面数值抹平带里马赫数发飘，原因在这里。GRP 的时间导数跟得上压缩和膨胀，去不掉这个非单调本身。论文这里也是用稳健性绕开，而不是正面拆掉。

工程代码多半是 Strang 分裂再加 Runge-Kutta，单级 GRP 塞不进去。要把好处留下来，通量函数得改写成 Lax-Wendroff 形式。移植的成本还在。线性化声学近似的出处，是 Ben-Artzi 与 Falcovitz 的气体动力学 GRP。

源项还在通量外面就指望那个阶的失误

表里的 2.03 是斜率光滑时的数字。激波上打开限制器，那个位置局部回到一阶。看见阶是 2 就删掉各级之前，先得有预算把通量函数重写一遍。

源项不溶进通量、留在外面，同样的刚性就交给另一种时间积分。先用两个压力求解、再收成一个压力的松弛是在分离这一侧守住体积分数边界，还是把 $\Delta t$ 再切细。这组输出里的 2.03 回答不了这件事。