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

我们一向被教导：想要二阶精度，就得让Runge-Kutta跑好几级。这篇论文只用一级就做到了。诀窍在于一个GRP(Generalized Riemann Problem,广义黎曼问题)求解器——它在界面上不只给出数值,还连同该值的时间导数一起求出。本文先用标量问题确认GRP为何一步即达二阶,再用代码重现这个时间导数如何驯服Kapila多相流模型的刚性源项。

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

单级二阶 — 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$ 是单元宽度。在下面的模拟中亲自调一调吧。

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

取消勾选 GRP 就回到Godunov:界面值(青色点)随时间流逝纹丝不动。重新勾选,界面值便沿特征线(橙色点)流动,送入通量的值是 $\Delta t/2$ 处的那个(橙色圆环)。把斜率 $u_x$ 置零,GRP就精确退化为Godunov。

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

时间导数怎么求?答案就在控制方程本身里。这正是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)

Godunov停留在一阶(0.9~0.97),而GRP收敛到2.0附近。我们没有增加任何一级,只是加入了斜率和时间导数。同样的网格上,误差小了一百倍。

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声速的压缩/膨胀系数。问题在于这个源项可能把 $\alpha_1$ 推出 $[0,1]$ 。显式积分一旦刚性增大就会越过边界而发散。

GRP的时间导数在此发光。因为我们能在新时间层 $t^{n+1}$ 知道界面值,Gauss-Green公式便能在新时间估计 $\nabla\cdot\mathbf{u}$ 的单元平均。于是可以用Crank-Nicolson(半隐式)积分源项,甚至能建立保持体积分数有界的时间步限制。我们用一个简单模型比较一下刚性松弛。

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

把 $\lambda\Delta t$ (刚性×时间步)调到1.6以上,显式Euler(红色)就冲破 $[0,1]$ 带并振荡。同样条件下Crank-Nicolson(青色)始终停在带内。即使把刚性升到3,青色仍单调收敛到平衡。这就是论文所说的"任意 $\Delta t$ 下的体积分数有界保持"。

复现时令人生疑之处#

标量GRP干净地给出了二阶,但真实的Kapila系统没那么温顺。第一,无限制的中心斜率在激波处会振荡。能出二阶是因为测试是光滑的;在间断处必须用minmod一类的限制器,届时阶数会局部跌回一阶。

第二,Wood声速对体积分数是非单调的。这正是界面数值扩散带内马赫数振荡的根源。GRP的时间导数能很好地捕捉压缩与膨胀,但并不能消除这种非单调性本身。论文对这一部分也是用稳健性绕过,而非正面求解。

第三,工程代码(OpenFOAM等)大多建立在Strang分裂加Runge-Kutta结构上,单级GRP很难直接移植。要享受GRP的好处,就得把通量函数本身改写成Lax-Wendroff型。移植成本不小。

接下来读什么#

GRP的本源是Ben-Artzi与Falcovitz的气体动力学GRP。想看清本论文的线性化声学近似从何而来,溯源到那里最合适。多相方面,Saurel的松弛-投影方法(本博客2026-05-09的文章)用另一种方式解决同样的刚性。把两者并置,设计上的分岔就清晰起来:是把源项分离出来,还是把它溶进通量里?

要点只有一个。只要在界面上连同数值一起求出时间导数,时间二阶精度与刚性源项就都能在一次评估中同时处理。