[論文レビュー] 非線形フラックスを線形に見せかける

安定なaはmax|f'|の少し上に置く。余裕を足すほどChapman–Enskog拡散が衝撃をなまらせる。

厳密なRiemann解法器(セル境界で波の構造を解く仕組み)を書き続けると、状態方程式で止まります。非線形フラックス $f(u)$ のヤコビアンを毎ステップ対角化しなければなりません。実在気体では固有構造が閉じた形で出ません。1995年、Shi JinとZhouping Xinは向きを逆転させました。非線形方程式をそのまま押しません。フラックス自体を新しい未知数へ上げ、線形系を作ってから元の式へ戻します。緩和(relaxation)です。Thomannら(2019)はこの骨格を全速度(all-speed)Eulerまで伸ばしました。迂回がどこで硬く、どこで崩れるか。その代価がこのノートの全体です。

論文情報 — Andrea Thomann, Markus Zenk, Gabriella Puppo, Christian Klingenberg, An all speed second order IMEX relaxation scheme for the Euler equations, arXiv:1907.08398 (2019)。

緩和速度のスライダーaを1未満に下げると、左が1、右が $-0.4$ のBurgers衝撃の波面が揺れます。

relaxation speed a = 1.30relaxation time ε = 1e-3

a = 1.30 ≥ max|u| = 1.00 — sub-characteristic condition satisfied. Smaller ε projects v onto f(u) faster (sharper shock, more relaxation diffusion trade-off).

aを1.3付近に置くと波面が立ちます。epsだけ上げると速度はそのままで、厚さだけが増します。

非線形フラックスがRiemann問題を難しくする#

スカラー保存則 $u_t + f(u)_x = 0$ を考えます。 $u$ は保存量、 $f(u)$ はフラックスです。境界ごとにセル一つ分の衝撃波管を解くには、先に波速 $f'(u)$ が要ります。

問題は実際の流れです。圧縮性Eulerではフラックスが圧力を含み、圧力は密度・内部エネルギーの非線形関数です。実在気体(real gas)の状態方程式が絡むと、固有値・固有ベクトルが解析的に求まりません。毎ステップ数値的にヤコビアンを対角化する費用は重くなります。低Mach領域では音波が $1/M$ で増大し、陽的時間ステップが窒息します。

そこで問いが立ちます。非線形性を毎ステップ正面から相手にせず、先送りできないでしょうか。

Jin–Xinの昇格 — フラックスを未知数へ持ち上げる#

Jin–Xinの答えはこうです。フラックス $f(u)$ を値として計算する代わりに、それを保持する**新変数 $v$ **を導入します。そして $v$ がゆっくり $f(u)$ へ引き寄せられるようソース項を付けます。

\begin{aligned} u_t + v_x &= 0, \\ v_t + a^2 u_x &= \frac{1}{\varepsilon}\bigl(f(u) - v\bigr). \end{aligned}

ここで $u$ は元の保存量、 $v$ はフラックスを代替する緩和変数、 $a$ は定数の緩和速度(relaxation speed)、 $\varepsilon$ は緩和時間です。左辺は完全に線形です。係数行列の固有値は $\pm a$ に固定され、非線形の対角化が消えます。

$\varepsilon \to 0$ では右辺のソースが $v = f(u)$ を強制し、第一式は元の保存則へ戻ります。つまりこの $2\times 2$ 線形系は非線形方程式の粘性近似です。波速が定数 $\pm a$ なので、Riemann問題は一度だけ手で解けます。上の図は、この系で進めたBurgers、 $f(u)=u^2/2$ です。

部分特性条件 — aが小さいと系が爆発する#

緩和速度 $a$ は自由に選べません。安定な近似になるには部分特性条件(sub-characteristic condition、Whitham条件)を守る必要があります。

a \ge \max_u |f'(u)|.

$f'(u)$ は元の方程式の真の波速です。緩和系が持ち歩く凍結速度 $\pm a$ が、その真の速度を常に挟み込むという意味です。Burgersなら $f'(u)=u$ なので $a \ge \max|u|$ で十分です。

下の図で緩和速度 $a$ と解の範囲 $\max|u|$ を変えてみましょう。橙色の曲線 $f'(u)$ が $\pm a$ の帯を出る瞬間が不安定の閾値です。

relaxation speed a = 1.40solution range max|u| = 1.00

The orange curve f′(u) stays inside the ±a band over the whole solution range → sub-characteristic condition holds.

$a$ を解の最大波速より下げると、曲線が帯を突き破ります。このとき緩和近似は物理的拡散ではなく負の拡散を生みます。情報が誤った方向へ流れ、系は爆発します。その理由は次節の展開で正確に見えてきます。

Chapman–Enskogで見る隠れた拡散#

$\varepsilon$ が小さいが0でないとき何が起きるでしょうか。 $v = f(u) + \varepsilon v_1 + \cdots$ と展開(Chapman–Enskog)して緩和系へ代入すると、有効方程式が現れます。

u_t + f(u)_x = \varepsilon\,\partial_x\Bigl[\bigl(a^2 - f'(u)^2\bigr)\,u_x\Bigr] + O(\varepsilon^2).

右辺が緩和がひそかに注入する拡散項です。拡散係数は $a^2 - f'(u)^2$ 。これが非負であるためには、まさに $a \ge |f'(u)|$ が必要です。部分特性条件の正体は**「緩和が生む人工粘性が負にならない条件」**だったのです。

ここでトレードオフが見えます。 $a$ を大きく取れば安定性は余裕ですが、 $a^2-f'^2$ が大きくなり解がなまります。 $a$ を部分特性の限界ぎりぎりに寄せると鋭いが危険です。 $\varepsilon$ は拡散の全体量を調節します。先のシミュレーションで $\varepsilon$ を上げると衝撃波が厚くなるのはこのためです。

余裕の幅が粘性の大きさです。コードのRiemannデータは左が $u=1$ 、右が $u=-0.4$ です。出発時の $\max|u|$ は1なので、部分特性の下限も1です。 $a=1.3$ なら $a^2=1.69$ です。左の $f'=1$ では拡散係数が $1.69-1=0.69$ 。右の $f'=-0.4$ では $1.69-0.16=1.53$ と、さらに大きくなります。同じ $\varepsilon=10^{-4}$ を掛けても、右の方がよく広がります。 $a=0.8$ なら $a^2=0.64$ で、左の係数は $0.64-1=-0.36$ です。符号が反転します。波面の情報が後ろへ戻り、振幅は初期値1に留まれません。

Suliciu — フラックス全部でなく圧力だけを緩和する#

Jin–Xinはフラックスの全成分を緩和します。優雅ですが拡散が過剰です。実務でより好まれるのはSuliciu型緩和です。Euler方程式で真に非線形なのは圧力だけなので、圧力だけを新変数 $\pi$ へ緩和します。

(\rho\pi)_t + \nabla\!\cdot(\rho\pi\,u) + a^2\,\nabla\!\cdot u = \frac{\rho}{\varepsilon}\,(p - \pi).

$\rho$ は密度、 $u$ は速度、 $p$ は実際の圧力、 $\pi$ は緩和圧力です。部分特性条件は $a > \rho\sqrt{\partial_\rho p}$ となり、音速を挟む形になります。緩和系の固有値は $u,\ u\pm a/\rho$ で、すべて線形退化(linearly degenerate)——接触不連続のように扱いやすい波だけが残ります。

Thomannら(2019)の貢献はここから始まります。低Mach領域を狙い、圧力を遅い成分と速い音響成分へ分けます。遅い部分は陽的に、速い音響部分は緩和系の上で陰的(implicit)に解きます。新しい速度変数 $\hat u$ を加えてMach数に依らない拡散を確保し、その結果、低Mach極限で非圧縮Eulerへ収束する漸近保存(asymptotic-preserving)の2次IMEX法を得ます。線形退化構造のおかげで、これらすべてが非線形の対角化なしに回ります。

Python — Burgersを緩和系で解く#

緩和アイデアの骨格だけをnumpyで再現しましょう。Burgers方程式をJin–Xin緩和系で解きます。特性変数 $r,s$ へ分解すると左辺は二つの単純な移流になります。

import numpy as np
 
def burgers_flux(u):
    return 0.5 * u * u                     # f(u) = u^2 / 2
 
def relaxed_step(u, v, a, dx, dt, eps):
    # 特性変数分解: r は +a、s は -a の速度で移動
    r = 0.5 * (u + v / a)
    s = 0.5 * (u - v / a)
    nu = a * dt / dx                       # CFL 数 a*dt/dx
    # 1次風上 (周期境界、np.roll)
    r_new = r - nu * (r - np.roll(r, 1))   # 右進行波
    s_new = s + nu * (np.roll(s, -1) - s)  # 左進行波
    u_new = r_new + s_new
    v_new = a * (r_new - s_new)
    # 緩和ソース: v を f(u) へ引き寄せる (指数積分)
    kappa = 1.0 - np.exp(-dt / eps)
    v_new += (burgers_flux(u_new) - v_new) * kappa
    return u_new, v_new
 
def run_relaxation(u0, a, eps, cfl=0.9, t_end=0.3):
    n = u0.size
    dx = 2.0 / n
    u = u0.copy()
    v = burgers_flux(u)                    # 平衡多様体から出発
    dt = cfl * dx / a                      # 部分特性が CFL を決める
    t = 0.0
    while t < t_end:
        u, v = relaxed_step(u, v, a, dx, dt, eps)
        t += dt
    return u
 
# Riemann 初期条件: 左 1.0、右 -0.4 (衝撃波)
n = 400
x = np.linspace(-1, 1, n, endpoint=False) + 1.0 / n
u0 = np.where(x < 0.0, 1.0, -0.4)
 
u_ok = run_relaxation(u0, a=1.3, eps=1e-4)   # a >= max|u|=1.0  -> 安定
u_bad = run_relaxation(u0, a=0.8, eps=1e-4)  # a <  max|u|      -> 違反
 
print("a=1.3  peak |u| =", round(float(np.max(np.abs(u_ok))), 3))   # ~1.0
print("a=0.8  peak |u| =", round(float(np.max(np.abs(u_bad))), 3))  # 発散

出力は部分特性条件をそのまま証言します。 $a=1.3$ では最大振幅が初期値付近に留まり、滑らかな衝撃波が形成されます。 $a=0.8$ では振幅が数倍に跳ね上がります。非線形ヤコビアンは一度も対角化していません。定数速度 $\pm a$ の線形移流を繰り返しただけです。

ステップ数も $a$ が決めます。 $n=400$ 、区間長2なので $\Delta x=0.005$ です。CFL 0.9なら $\Delta t=0.9\times 0.005/a$ です。 $a=1.3$ では $\Delta t\approx 3.46\times 10^{-3}$ で、 $t=0.3$ までの比は約87です。 $a=0.8$ では $\Delta t\approx 5.63\times 10^{-3}$ と粗くなります。 $a$ を下げると時間刻みはかえって伸びます。 $\varepsilon=10^{-4}$ では $\Delta t/\varepsilon$ が30を超え、 $\kappa=1-\exp(-\Delta t/\varepsilon)$ の指数項が消えます。毎ステップ $v$ は $f(u)$ に張り付き、残るぼやけは $\varepsilon(a^2-f'^2)$ だけです。

批判的に見ると

緩和はただではありません。人工拡散 $\varepsilon(a^2-f'^2)$ が必ず付きまとい、 $a$ を安全に大きく取るほど解がなまります。部分特性条件を守るには $\max|f'(u)|$ の大域的上限が必要で、強い衝撃や真空近傍でこの上限を保守的に取ると拡散が過剰になります。スカラーでは綺麗だった論理が、実在気体Eulerでは $\rho\sqrt{\partial_\rho p}$ の局所推定と、圧力・密度の正値保存のための追加装置を要求します。Thomannらのimplicit-explicit結合は音響部で楕円型方程式を解くため、低Machの利得が線形ソルバの費用と相殺される点も存在します。再現してみると、 $\varepsilon$ と $a$ の同時調整は思ったより敏感です。

OpenFOAMやFluentの密度ベースソルバにこの緩和Riemann解法器が直接入っているわけではありませんが、HLLC・AUSM系の近似Riemann解法器が同じ哲学——波構造を単純化して対角化を避ける——を共有します。Suliciu解法器はSU2のようなオープンソースコードに、接触不連続を厳密に捉える選択肢として実装されています。

下限ぎりぎりのaが次の衝撃で崩れる場所#

このデータの下限は1で、コードは1.3を選んでいます。差0.3が、左状態の拡散係数0.69です。余裕を0まで削って $a=\max|u|$ にすると、波面上の係数は0になります。捕獲は鋭くなります。丸めが $a$ をほんの少し下限の下へ押すと、符号はすぐ反転します。

上限を最初の一歩だけで測ると、次の衝撃で崩れます。 $|u|$ が1を超えた瞬間、安全に見えた $a=1.3$ も部分特性の外です。Suliciuへ移すと不等式は $a>\rho\sqrt{\partial_\rho p}$ に変わります。音速がセルごとに跳ねる実在気体で上限が一拍遅れると、音響を陰的に解いても緩和段階の拡散はすでに負です。圧力だけを触るという選択は、移流と圧力を切り離す流束分離と同じ場所にあります。

[論文レビュー] 非線形フラックスを線形に見せかける — Jin–Xin・Suliciu緩和法