任意方向上的欧拉特征结构 — 三族波、两重重根、一处奇点

出差时在酒店打开了笔记本。三天前把网格旋转到任意方向开始,求解器就反复吐出负密度。$n_x=1, n_y=0$ 的对齐情况一切正常。一旦 $n_x=0.01$,数值立刻发散。倒着翻日志,锁定到一行符号反掉了。那一行,是左特征向量矩阵的最后一行。整整三天,都在追那一行藏着什么。

这篇文章是那三天的整理。在任意单位法向量 $\hat n = (n_x, n_y, n_z)$ 上,把 3D 欧拉方程的特征值与左右特征向量一次写下,看清左特征向量为何在某一种约定下必然崩塌。最后用同一套式子写一个 Roe 型通量。

拖动 $M = v_n/a$,三条射线在哪里分开一目了然。

$|M|<1$ 时声波分在原点两侧,$|M|>1$ 时三条射线都贴到同一个 $x$ 方向。入流出流边界为何随马赫数改指定量个数,答案在这张图里。

守恒型欧拉沿法向写

3D 无粘欧拉方程是一个向量守恒律。

其中 $\mathbf{Q} = (\rho,\rho u,\rho v,\rho w,\rho e_o)^T$ 是守恒变量,$\mathbf{F}\cdot \hat n$ 是每个面上的法向通量。单元面上自然的变量是法向速度 $v_n = u n_x + v n_y + w n_z$。声速 $a^2 = \gamma p/\rho$,单位法向意味着 $n_x^2 + n_y^2 + n_z^2 = 1$。

这种坐标的优势很直接:所有波族都对齐到单一变量 $v_n$。

雅可比矩阵 A(n̂) — 5×5 里藏着五个波#

把 $\mathbf{F}\cdot \hat n$ 对 $\mathbf{Q}$ 求导,得到变换矩阵 $A(\hat n)= \partial(\mathbf F\cdot\hat n)/\partial\mathbf Q$。不需要把全局 $(x,y,z)$ 和局部 $(\xi,\eta,\zeta)$ 分别写两份,一个矩阵显式带着 $\hat n$ 即可。这是 1983 年 Yee–Kutler 之后的标准写法。

特征方程 $\det(A - \lambda I) = 0$ 的根有五个。

流速 $v_n$ 出现了 三次,其中两次是重根。这一重数,是后面让左特征向量动摇的第一条线索。

三个家族,两个重根 — 声波与熵波

五个特征值分成 三个家族。流速 $v_n$ 本身是 熵波,它以同一速度搬运熵和切向动量两种独立信号。所以重数为二。

$v_n \pm a$ 是 声波 两支,携带压力信号和法向速度信号。在 Riemann 扇形中,中央的直线是接触间断,外侧的扇形或激波是声波。

把这个结构记在脑里,调试就快了。如果只有密度跳跃,$\Delta p = \Delta u = 0$,那么只有熵波强度该存活,两个声波强度应严格为 0。下图能在线确认这一点。

从 "no jump" 预设开始,稍微改动 $\rho_R$。中间的灰色柱子开始增长,左右红蓝柱子保持在 0。换到 "Sod"。右声波柱大幅竖起 — 那就是激波。再换到 "123 problem",左右对称的稀疏波,两个声波柱同号大幅出现,接触柱近乎为零 — 那是稀疏-稀疏的指纹。

右特征向量 — 一次写完

解 $A R_i = \lambda_i R_i$,得到一个 5×5 矩阵。一种约定 (R-1) 写为

各符号的意义如下。$e_k = \tfrac12(u^2+v^2+w^2)$ 是单位质量动能,$h_o = h + e_k$ 是单位质量总焓。前三列对应 $\lambda_1, \lambda_2, \lambda_3$;最后两列对应重根 $\lambda_4=\lambda_5=v_n$,是两个切向方向向量,具体取哪两个由约定决定。

第三列与第五列其实是"选了哪一对切向方向"的痕迹。R-2 与 R-3 不过是用另外一对基张同一个子空间。

左特征向量的陷阱 — n_x → 0 时崩掉的一行#

计算 $L = R^{-1}$,最后两行会出现 $1/n_x$ 这样的项。R-1 约定里总有一行要除以 $n_x$。一旦遇到像 $\hat n = (0, 1, 0)$ 这样 $n_x = 0$ 的面,那一行的左特征向量瞬间发散。出差时我追了三天的 bug 就是这个。

修法很短。R-2 用 $n_y$ 作分母,R-3 用 $n_z$。每个面选 $|n_x|, |n_y|, |n_z|$ 中 最大分量 作分母的约定。一行代码:

def pick_convention(n):
    # n: 形状 (3,) 的单位法向量
    k = int(np.argmax(np.abs(n)))  # 0、1 或 2
    return k  # 选 R-1、R-2 或 R-3

三种约定里,至少有一种在该面上是安全的。非结构网格的代码若不在每个面上做这种分支,几乎活不了多久。熵修正为什么放在最后,写在 Roe 求解器的熵修正 里。

落到 2D 后悄悄消失的奇点#

令 $w=n_z=0$,5×5 化为 4×4。五个特征值缩为 $\{v_n - a, v_n, v_n + a, v_n\}$ 四个,$L$ 最后一行的形式上的奇点会两两相消。例如

分母上的 $n_x$ 被分子整除。所以一辈子只写 2D 代码的人,可能终生不会遇到这个陷阱。它只在从 2D 跨到 3D 那一刻才出现。

Python — 一个函数算任意方向的 Roe 型通量#

把三族特征值与 1D Roe 平均的波强度集中到一个函数里,返回任意 $\hat n$ 方向的 Roe 型通量。玩具问题是 倾斜的激波管: 把左右状态投到 $\hat n = (\cos 30°, \sin 30°, 0)$ 的面上。

import numpy as np
 
GAMMA = 1.4
 
def primitive_to_state(rho, vel, p):
    """rho: 标量, vel: (3,), p: 标量 -> 守恒量 Q (5,)"""
    rE = p / (GAMMA - 1) + 0.5 * rho * np.dot(vel, vel)
    return np.array([rho, rho * vel[0], rho * vel[1], rho * vel[2], rE])
 
def flux_normal(Q, n_hat):
    rho = Q[0]
    vel = Q[1:4] / rho
    rE = Q[4]
    p = (GAMMA - 1) * (rE - 0.5 * rho * np.dot(vel, vel))
    vn = np.dot(vel, n_hat)
    H = (rE + p) / rho
    return np.array([
        rho * vn,
        rho * vel[0] * vn + p * n_hat[0],
        rho * vel[1] * vn + p * n_hat[1],
        rho * vel[2] * vn + p * n_hat[2],
        rho * H * vn,
    ])
 
def roe_normal_flux(QL, QR, n_hat):
    rL, rR = QL[0], QR[0]
    vL = QL[1:4] / rL
    vR = QR[1:4] / rR
    pL = (GAMMA - 1) * (QL[4] - 0.5 * rL * np.dot(vL, vL))
    pR = (GAMMA - 1) * (QR[4] - 0.5 * rR * np.dot(vR, vR))
    HL = (QL[4] + pL) / rL
    HR = (QR[4] + pR) / rR
 
    sL, sR = np.sqrt(rL), np.sqrt(rR)
    inv = 1.0 / (sL + sR)
    rho = sL * sR
    vel = (sL * vL + sR * vR) * inv
    H   = (sL * HL + sR * HR) * inv
    vn  = np.dot(vel, n_hat)
    a2  = (GAMMA - 1) * (H - 0.5 * np.dot(vel, vel))
    a   = np.sqrt(max(a2, 1e-12))
 
    eigs = np.array([vn - a, vn, vn + a])
    dvel = vR - vL
    dvn  = np.dot(dvel, n_hat)
    drho = rR - rL
    dp   = pR - pL
 
    alpha_contact = drho - dp / a2
    alpha_minus   = (dp - rho * a * dvn) / (2 * a2)
    alpha_plus    = (dp + rho * a * dvn) / (2 * a2)
 
    # 平均通量 + |lambda| * alpha * R
    Favg = 0.5 * (flux_normal(QL, n_hat) + flux_normal(QR, n_hat))
 
    # 当场构造三个右特征向量
    R_minus = np.array([1, vel[0] - a * n_hat[0], vel[1] - a * n_hat[1],
                        vel[2] - a * n_hat[2], H - a * vn])
    R_cont  = np.array([1, vel[0], vel[1], vel[2], 0.5 * np.dot(vel, vel)])
    R_plus  = np.array([1, vel[0] + a * n_hat[0], vel[1] + a * n_hat[1],
                        vel[2] + a * n_hat[2], H + a * vn])
 
    diss = (np.abs(eigs[0]) * alpha_minus   * R_minus
          + np.abs(eigs[1]) * alpha_contact * R_cont
          + np.abs(eigs[2]) * alpha_plus    * R_plus)
 
    # 切向速度跳跃(由接触波 lambda = vn 携带)
    dvel_tan = dvel - dvn * n_hat
    R_tan = np.zeros(5); R_tan[1:4] = dvel_tan
    R_tan[4] = rho * np.dot(vel, dvel_tan)
    diss = diss + np.abs(vn) * R_tan
 
    return Favg - 0.5 * diss
 
if __name__ == '__main__':
    n_hat = np.array([np.cos(np.pi/6), np.sin(np.pi/6), 0.0])
    QL = primitive_to_state(1.0,   np.array([0.0, 0.0, 0.0]), 1.0)
    QR = primitive_to_state(0.125, np.array([0.0, 0.0, 0.0]), 0.1)
    F = roe_normal_flux(QL, QR, n_hat)
    print('Sod tilted by 30 deg, Roe flux =', F)

这是手写形式,并非最快实现。每个面上的每一项都在可视范围内,这正是调试需要的样子。生产代码要把 pick_convention 按面分支以切换左特征向量约定,并在末尾加上熵修正(Harten 或部分 LLF)。

忘了 np.argmax(np.abs(n_hat)) 这一行,只埋死了 R-1 — 你就会丢三天。在那间酒店里我损失的不是早餐费,而是这三天。

面上该用哪个分母

$n_x=0.01$ 时,R-1 所除的 $1/n_x$ 等于 100。对齐面 $n_x=1$, $n_y=0$ 的分母是 1,所以安静。倾到 0.01,分母变成 100,那一行花了三天。代码里 $\gamma=1.4$,声速是 $a^2=\gamma p/\rho$。变成 0 时最后两行发散。代码里的 30° 面是 $n_x=\cos(\pi/6)$, $n_y=\sin(\pi/6)=1/2$, $n_z=0$。因为 $\pi/6<\pi/4$,$|n_x|$ 大于 $|n_y|$。pick_convention 返回 0,选用 R-1。法向 $(0,1,0)$ 的索引是 1,分母是 $n_y$,约定是 R-2。

降到二维后,和上面的恒等式一样,分子按 $n_x$ 裂开,这一行保持有限。网格变成三维非结构的那一天,每个面才需要 $\max |n_i|$。斜激波管的左侧是 $(\rho,p)=(1,1)$,右侧是 $(0.125,0.1)$,速度为 0。密度比 $1/0.125=8$,压力比 $1/0.1=10$。$\Delta p$ 和 $\Delta\rho$ 都不为零,所以接触强度和两侧声波强度同时还在。这和右声波柱长起来的 Sod 预设是同一结构。边界上声波射线随马赫数进来几条,画在 特征线与声波 里。