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

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

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

守恒型欧拉沿法向写#

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$ 出现了 三次,其中两次是重根。这一重数,是后面让左特征向量动摇的第一条线索。

拖动 $M = v_n/a$。三条射线呈扇形展开。$|M|<1$ 时,声波分布在原点两侧;一旦 $|M|>1$,所有三条射线偏向同一个 $x$ 方向。入流出流边界为何随马赫数切换给定量个数,答案就在这一帧里。

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

五个特征值分成 三个家族。流速 $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

三种约定里,至少有一种在该面上是安全的。非结构网格的代码若不在每个面上做这种分支,几乎活不了多久。

落到 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 — 你就会丢三天。在那间酒店里我损失的不是早餐费,而是这三天。

关上调试日记#

• 用任意单位法向 $\hat n$ 写出的 5×5 雅可比把全局 $(x,y,z)$ 与局部 $(\xi,\eta,\zeta)$ 用一个式子统一。
• 特征值分三族 — 声波 $v_n \pm a$,熵波 $v_n$(二重)。这一重数在左特征向量里造出子空间。
• R-1 左特征向量约定在 $n_x=0$ 面上发散。每个面按 $\arg\max |n_i|$ 选约定的一行代码,是非结构求解器活下去的关键。