任意方向のオイラー固有構造 — 三つの波、二重の重複、一つの特異点

出張先のホテルでノートパソコンを開きました。三日前から、格子を任意方向に回した瞬間、ソルバーが負の密度を吐き続けていました。$n_x=1, n_y=0$ の整列ケースでは無事です。$n_x=0.01$ になった途端、発散が始まります。ログを逆向きにたどると、一行の符号が反転していました。それは左固有ベクトル行列の最終行でした。三日間、その一行が何を隠しているのかを突き止める作業でした。

本稿はその三日間の整理です。任意の単位法線 $\hat n = (n_x, n_y, n_z)$ 上で 3D オイラー方程式の固有値と左右の固有ベクトルを一度に書き下し、左固有ベクトルがある規約で必ず崩れる理由を示します。最後に同じ式で Roe 型流束を組み立てます。

保存形オイラーを任意方向に書く#

3次元非粘性オイラー方程式は、ベクトル形式の保存則です。

ここで $\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" は左右対称の膨張波なので、二つの音響バーが同じ符号で大きく出て、接触はほぼ 0 になります。

右固有ベクトル — 一度に書く#

$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$ の面が一度でも現れた途端、そこでの左固有ベクトルは発散します。これが出張先で三日を奪ったバグでした。

解決は単純です。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, or 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\}$ の四つに減り、左固有ベクトル行列の最終行は二項が綺麗に相殺します。たとえば

のように、形式上分母に $n_x$ があっても分子が $n_x$ で割り切れます。だから 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$ の面で発散する。面ごとに $\max |n_i|$ のインデックスで規約を選ぶ一行が、非整列コードを生かす。