放弃精确黎曼解换来了什么 —— Roe 近似黎曼求解器与熵修正#

一维 Euler 的黎曼问题有精确解。星区压力的一个非线性方程,Newton 一次就收工。生产级可压缩代码并不在每个界面上用它。答案在手里。还是丢掉。成本。离开理想气体,精确解本身可能就不存在。

拨动左右密度和速度滑条,界面上的 √ρ 权重和三个波速 $\hat{u}-\hat{c}$、$\hat{u}$、$\hat{u}+\hat{c}$ 会重新算出来。

$\hat{u}-\hat{c}<0<\hat{u}+\hat{c}$ 时,波扇跨在界面上。亚声速。把两侧速度狠狠推向同一边,三个波的符号变得一致,状态进入超声速。

精确解为何从工程现场消失

精确黎曼求解器（Godunov）会在每个界面上求解一个非线性方程的根。当网格有数百万个单元时，这个求根循环每步就要跑数百万次。更糟的是，精确解被绑死在理想气体状态方程上。一旦换成真实气体或两相混合物，精确解本身就不复存在。

近似黎曼求解器绕开了这个难题。它在局部把非线性黎曼问题线性化。一组代数公式即可给出通量，无需任何迭代。代价是牺牲一点精确性。

把雅可比矩阵压缩成一个平均

Roe 的想法很简单。用一个常数矩阵 $\hat{A}$ 作用于界面两侧的两个状态 $q_L$、$q_R$，来近似界面上的通量差。

其中 $q$ 是守恒量向量，$f$ 是物理通量。这个条件具有决定性意义。若两个状态由单一波动（激波或接触间断）相连，那么 $\hat{A}$ 会精确传播这个波动，即便在有限幅度下也成立。也就是说，尽管它是近似求解器，但在纯激波或接触间断面前却是精确的。

难点在于如何选取 $\hat{A}$。随便取个平均就会破坏上面的条件。

√ρ 平均从何而来

Roe 的答案是以密度的平方根加权的平均。

其中 $u$ 是速度，$H$ 是总比焓，尖帽符号表示界面平均值。Roe 声速随之为 $\hat{c}=\sqrt{(\gamma-1)(\hat{H}-\hat{u}^2/2)}$。

为什么偏偏是 $\sqrt{\rho}$？把状态 $q$ 和通量 $f$ 写成参数向量 $W=\sqrt{\rho}\,(H,u,1)^T$ 各分量的函数，二者都会变成完美的二次式。既然是二次式，$q$ 和 $f$ 对 $W$ 的导数就是线性的，构造 $\hat{A}$ 的路径积分也就能精确算出。掉出来的正好就是这个 √ρ 加权平均。

权重的分母是 $\sqrt{\rho_L}+\sqrt{\rho_R}$。一侧密度是另一侧的4倍,平方根就是2倍,那一侧的速度和焓会以两倍的力气拉平均。算术平均在有限跳跃上满足不了 $\hat{A}(q_R-q_L)=f(q_R)-f(q_L)$。没有 √ρ,纯激波前面的波速就已经错了。声速是 $\hat{c}=\sqrt{(\gamma-1)(\hat{H}-\hat{u}^2/2)}$。$\gamma=1.4$ 时 $\gamma-1=0.4$。总焓减去动能 $\hat{u}^2/2$ 之后剩下的,就是声速的平方。三个波速变号的地方,和 任意法向上拆开的欧拉特征结构 是同一组三波。

拆成波，再加回来

组装通量分三步。首先把状态跳跃 $\Delta q = q_R - q_L$ 分解为三个特征向量 $\hat{K}_k$ 之和；每个分量的大小就是波强 $\alpha_k$。然后根据每个波特征值 $\hat{\lambda}_k$ 的符号，把它们沿迎风方向输运。

前一项是中心平均，后一项是以特征值大小加权的迎风耗散。这套结构可以直接搬进代码。测试算例选用 Shu–Osher 问题：一道马赫 3 激波冲入正弦密度场，在其身后留下高频结构。这正是 Roe 低数值耗散大显身手的场合。

import numpy as np
 
gamma = 1.4
 
def phys_flux(U):                       # 守恒量 U=(rho, rho*u, E) -> 物理通量
    rho = U[0]; u = U[1] / rho; E = U[2]
    p = (gamma - 1) * (E - 0.5 * rho * u * u)
    return np.array([rho * u, rho * u * u + p, u * (E + p)])
 
def roe_flux(UL, UR, delta):            # Roe 近似黎曼通量 (+ Harten 熵修正)
    rhoL, rhoR = UL[0], UR[0]
    uL, uR = UL[1] / rhoL, UR[1] / rhoR
    pL = (gamma - 1) * (UL[2] - 0.5 * rhoL * uL * uL)
    pR = (gamma - 1) * (UR[2] - 0.5 * rhoR * uR * uR)
    HL = (UL[2] + pL) / rhoL; HR = (UR[2] + pR) / rhoR
    sL, sR = np.sqrt(rhoL), np.sqrt(rhoR)      # sqrt(rho) 权重
    u = (sL * uL + sR * uR) / (sL + sR)        # Roe 平均速度
    H = (sL * HL + sR * HR) / (sL + sR)        # Roe 平均焓
    c = np.sqrt((gamma - 1) * (H - 0.5 * u * u))   # Roe 平均声速
    rr = sL * sR
    drho, dp, du = rhoR - rhoL, pR - pL, uR - uL
    alpha = np.array([(dp - rr * c * du) / (2 * c * c),   # 波强 alpha_k
                      drho - dp / (c * c),
                      (dp + rr * c * du) / (2 * c * c)])
    lam = np.array([u - c, u, u + c])          # 特征值 (波速)
    K = np.array([[1, u - c, H - u * c],       # 右特征向量
                  [1, u,     0.5 * u * u],
                  [1, u + c, H + u * c]])
    al = np.abs(lam)
    small = al < delta                         # 熵修正: |lambda| 下限
    al[small] = (lam[small] ** 2 + delta ** 2) / (2 * delta)
    diss = (al * alpha) @ K
    return 0.5 * (phys_flux(UL) + phys_flux(UR)) - 0.5 * diss
 
def run_shu_osher(N=400, tmax=1.8, cfl=0.4, delta=0.1):
    x = np.linspace(0, 10, N); dx = x[1] - x[0]
    rho = np.where(x < 1, 3.857143, 1 + 0.2 * np.sin(5 * x))   # 激波 + 正弦密度
    u   = np.where(x < 1, 2.629369, 0.0)
    p   = np.where(x < 1, 10.33333, 1.0)
    U = np.array([rho, rho * u, p / (gamma - 1) + 0.5 * rho * u * u])
    t = 0.0
    while t < tmax:
        r = U[0]; v = U[1] / r; pp = (gamma - 1) * (U[2] - 0.5 * r * v * v)
        dt = cfl * dx / np.max(np.abs(v) + np.sqrt(gamma * pp / r))
        dt = min(dt, tmax - t)
        F = np.zeros((3, N + 1))
        for i in range(1, N):
            F[:, i] = roe_flux(U[:, i - 1], U[:, i], delta)
        F[:, 0] = phys_flux(U[:, 0]); F[:, N] = phys_flux(U[:, N - 1])
        U[:, 1:N - 1] -= dt / dx * (F[:, 2:N] - F[:, 1:N - 1])
        t += dt
    return x, U[0]
 
x, rho = run_shu_osher()
print(f"t=1.8  min rho={rho.min():.3f}  max rho={rho.max():.3f}")   # -> min~0.81  max~4.08

大约四十行代码就构成了一个完整的可压缩求解器。没有 Newton 迭代，也没有精确黎曼求解。delta 就是接下来要讲的熵修正参数。

Shu–Osher 初场在 $x<1$ 处是 $\rho=3.857143$、$u=2.629369$、$p=10.33333$,右侧是 $\rho=1+0.2\sin(5x)$、$u=0$、$p=1$。$N=400$,区间 $[0,10]$,CFL 为 0.4,$t=1.8$,$\delta=0.1$。输出里 $\rho$ 最小约 0.81,最大约 4.08。左侧密度 3.857 略升到 4.08 附近,正弦波的波谷落到 0.8 附近。没有 Newton,这个范围仍被保住,就是这次运行的结果。

Roe 违背了熵条件#

Roe 求解器把每个波都当作跳跃来处理。连膨胀波（rarefaction）也被近似为一串小激波。多数情况下这没有问题。但当膨胀波中包含一个声速点（sonic point）时，情况就变了。在那一点上，某个特征值 $\hat{\lambda}_k$ 会变号，穿过 $0$。

当 $|\hat{\lambda}_k|\to 0$ 时，该波对应的迎风耗散消失了。数值格式不再展开成一道光滑的膨胀扇，而是锁死了一道静止的膨胀激波（expansion shock）。它满足 Rankine–Hugoniot 条件，却违背了熵条件 —— 这是一个非物理解。

解药是 Harten 的熵修正。它在 $0$ 附近给特征值大小铺一个下限。

$|\hat{\lambda}_k|=0$ 时,修正式等于 $\delta/2$,因为 $(0+\delta^2)/(2\delta)=\delta/2$。$\delta$ 为 0,地板也是 0,迎风耗散消失。代码里的 delta=0.1 在这块地板上铺了 0.05。$(|\hat{u}|+\hat{c})$ 的 5% 到 10% 是接触间断被抹得太开之前的宽度。

把 $\delta$ 滑条从 0 往上推,$u_L=-0.5$、$u_R=1.0$ 的 Burgers 膨胀 $u_t+(u^2/2)_x=0$ 里,$x=0$ 的蓝色折点会不会动就看得见。声速点 $u=0$ 在区间正中。

当 $\delta=0$ 时，$x=0$ 处会留下一道蓝色折点 —— 即静止的膨胀激波。把 $\delta$ 调大，数值曲线就会滑落到琥珀色的精确膨胀扇之上。取局部 $(|\hat{u}|+\hat{c})$ 的 5%~10% 通常安全。取得太大，接触间断就会被抹糊。

更省：HLL 与 HLLC#

如果觉得 Roe 太重，那就少留几个波。HLL（Harten–Lax–van Leer）只保留左右两个声波，丢掉中间的接触波。

其中 $S_L$、$S_R$ 是左右最外侧波速的估计值。HLL 鲁棒，能很好地保持密度为正。代价是它无法把接触间断抓得锐利 —— 因为没有中间波。

HLLC 中的 C 代表中央接触波（Contact）。它把被丢掉的接触波重新找了回来：三个波，四个常状态区。归根结底，这三者落在同一条"保留几个波"的谱线上。

工程现场的默认选择通常是 HLLC。它既保住了接触间断，又易于强制密度和压力的正性。Roe 分辨率出色，但熵修正是必需的，而且在与网格平行的强激波上还得提防"红斑"（carbuncle）现象。

δ留在0,膨胀激波就停在界面上#

√ρ 权重来自参数化 $W=\sqrt{\rho}\,(H,u,1)^T$,它把 $q$ 和 $f$ 写成二次式。没有这条路径积分,$\hat{A}$ 的跳跃条件在有限幅度上就不成立。声速点上任由 $|\hat{\lambda}_k|\to 0$,Rankine–Hugoniot 仍成立,熵条件却不成立。不铺 Harten 地板的 Roe,会把膨胀激波锁成静止解。

接下来要决定留几个波。HLL 只留两道声波,接触间断被抹开。HLLC 把接触波找回来,定常区变成四块。Roe 在三维里可以到五道波,但熵修正是必需的,与网格平行的强激波还得盯着红斑。低马赫数下人工耗散的尺度会另分一岔。Hope–Collins 把低马赫 Roe 分成三种尺度的记录就是那一岔。