超音速流れは上流を知らない — Euler方程式の特性曲線と音波

影響領域の楔が垂直線を越える瞬間がマッハ1で、それ以上の出口では課せる条件の数が特性の符号とともに減ります。

超音速戦闘機は、自分のエンジン音より先に到着します。頭上を通り過ぎたあとで、ようやく轟音が追いかけてきます。理由は Euler 方程式の特性速度にあります。小さな擾乱は三本の信号に分かれ、その軌跡が $(x,t)$ 平面の特性曲線(信号が通り抜けた跡)になります。

M のスライダーを 0.4 に置くと、青い左向き音波が source の上流へ伸びます。1 を超えると、その軌跡は画面から消えます。

background flow u₀ (Mach, c₀=1) = 0.40

The blue tint marks the upstream region the left-going sound has reached. Push u₀ past 1 and the blue tint vanishes: the C− signal now has a positive speed, so every signal runs downstream and the source can no longer talk to what lies upstream.

青い線の傾きの符号が変わる瞬間です。三つの信号がすべて右へだけ走るのか、左側が空になるのかを見ます。

圧縮性方程式は信号の伝播である

1次元非粘性圧縮性流れは、三つの保存則で書かれます。

\partial_t \rho + \partial_x(\rho u) = 0

\partial_t(\rho u) + \partial_x(\rho u^2 + P) = 0

\partial_t(\rho e_{\text{tot}}) + \partial_x\big[(\rho e_{\text{tot}} + P)\,u\big] = 0

$\rho$ は密度、 $u$ は速度、 $P$ は圧力、 $e_{\text{tot}}$ は全比エネルギーです。見た目は流体が「流れる」方程式のようです。しかし、この方程式の本質は移動ではなく信号の伝播です。流れのあらゆる動力学は、決まった速度で走るいくつかの信号に分解されます。その速度が何なのかを見るには、方程式をちょっと線形化すればよいのです。

小さな擾乱が波動方程式をつくる

背景状態 $(\rho_0, u_0, P_0)$ の上に、ごく小さな擾乱を重ねます。 $\rho = \rho_0 + \rho_1$ , $u = u_0 + u_1$ と置いて $P = K\rho^\gamma$ を使い、1次の項だけを残します。背景に沿って動く微分 $D_t \equiv \partial_t + u_0\,\partial_x$ (共動微分)を定義すると、二つの擾乱式を合わせて一つの波動方程式が出てきます。

D_t^2\,\rho_1 - c_0^2\,\partial_x^2\,\rho_1 = 0

ここで $c_0$ は音速です。

c_0 = \sqrt{\gamma\,\frac{P_0}{\rho_0}}

$\gamma$ は比熱比です。この方程式に $\rho_1 \sim e^{i(kx-\omega t)}$ を入れると、分散関係が落ちてきます。

\frac{\omega}{k} = u_0 \pm c_0

$\omega$ は角振動数、 $k$ は波数です。擾乱は背景流れ $u_0$ に乗って、音速 $c_0$ の分だけ左右に広がります。つまり一点の擾乱は、 $u_0 + c_0$ と $u_0 - c_0$ という二つの速度に分かれます。

三本に分かれる特性曲線

二つの音波だけが全てではありません。三番目の信号は流体そのものです。流体に色を塗ったとしましょう。 $t=0$ で左を青、右を赤に染めると、その色の境界 $\phi$ は単純に流れに乗って運ばれます。

\partial_t \phi + u_0\,\partial_x \phi = 0

この信号は速度 $u_0$ で移動します。これがエントロピー(あるいは接触面)特性です。整理すると、1次元流れには三つの特性速度があります。

\lambda_- = u_0 - c_0, \qquad \lambda_0 = u_0, \qquad \lambda_+ = u_0 + c_0

$(x,t)$ 平面で $\dfrac{dx}{dt} = \lambda$ を満たす線を特性曲線と呼びます。 $\lambda_\pm$ は音波、 $\lambda_0$ は物質の移動です。各音波特性に沿ってはRiemann不変量 $u \pm \dfrac{2c}{\gamma-1}$ が一定に保存され、物質特性に沿ってはエントロピーが保存されます。

$M = u_0/c_0$ が 0.4 なら $\lambda_-$ は負です。左向きの音波は上流へ届きます。 $M$ が 1 を超えると $\lambda_- = u_0 - c_0$ が正になり、三つの信号は下流にしか走りません。

三つの速度の固有値分解は、任意方向の法線でのオイラー固有構造へ進むと重根と接線成分が付きます。1次元では、その重複がエントロピー特性ひとつに畳まれています。

超音速では信号の円錐が越えていく

同じ話を時空図で見ると、さらに鮮明になります。sourceから伸び出た三本の特性曲線 $C_-,\,C_0,\,C_+$ のあいだの楔は影響領域(domain of influence)です。sourceが信号を送りうる時空事象の集合です。

$M<1$ なら、この楔は垂直線を挟んで両側に開きます。上流・下流のどちらへも信号が届きます。 $M>1$ なら、楔全体が右へ越えていきます。 $C_-$ 特性までもが右へ傾き、上流のどんな事象も影響領域の外です。これが超音速流れで「上流が下流を知らない」理由です。数値的にも、超音速出口に境界条件を課してはならない根拠がここにあります。

下の図で丸い探針をドラッグして動かしてみましょう。

background flow u₀ (Mach, c₀=1) = 0.40

Drag the probe across the plane. At M = 0.4 the wedge straddles the vertical, so events both upstream and downstream are reachable. Slide u₀ above 1 and the whole wedge leans right — drop the probe anywhere to the left of the source and it stays red.

探針を楔の中に置くと緑(到達可能)、外に置くと赤(到達不可)に変わります。 $u_0$ を1より上に上げたあと、探針をsourceの左のどこに置いても、ずっと赤のまま残ります。上流が信号から完全に切り離されたのです。

Pythonで確かめる特性速度#

特性速度が本当に $u \pm c$ と $u$ なのかを確かめてみましょう。1次元Eulerを原始変数 $W=(\rho,u,P)$ で書くと、準線形形 $\partial_t W + A(W)\,\partial_x W = 0$ になり、係数行列 $A$ の固有値がそのまま特性速度です。

import numpy as np
 
def characteristic_speeds(rho, u, p, gamma=1.4):
    """1D Euler 原始変数ヤコビアン A(W), W=(rho,u,p)。
    A の固有値 = 三つの特性速度。"""
    c = np.sqrt(gamma * p / rho)              # 音速
    A = np.array([[u,      rho,        0.0  ],
                  [0.0,    u,          1/rho],
                  [0.0,    rho*c**2,   u    ]])
    lam = np.sort(np.linalg.eigvals(A).real)
    return c, lam
 
rho, u, p = 1.225, 220.0, 1.0e5               # 空気, u = 220 m/s
c, lam = characteristic_speeds(rho, u, p)
print(f"c0         = {c:8.2f} m/s")
print(f"eig(A)     = {np.round(lam, 2)}")
print(f"u-c,u,u+c  = {np.round([u-c, u, u+c], 2)}")
print(f"Mach       = {u/c:5.3f}  -> {'supersonic' if u > c else 'subsonic'}")

出力は次のとおりです。

c0         =   338.06 m/s
eig(A)     = [-118.06  220.    558.06]
u-c,u,u+c  = [-118.06  220.    558.06]
Mach       =  0.651  -> subsonic

固有値がちょうど $u-c,\,u,\,u+c$ と一致します。今度は実際に擾乱を時間前進させ、信号の前面速度を測ってみます。音響Riemann不変量 $w_\pm = u_1 \pm (c_0/\rho_0)\,\rho_1$ はそれぞれ $u_0 \pm c_0$ で運ばれるスカラーなので、風上差分で別々に移してやればよいのです。

N, L = 400, 1.0
x = np.linspace(0, L, N, endpoint=False)
dx = L / N
rho0, c0 = 1.0, 1.0
u0 = 0.6 * c0                                  # 背景 Mach 0.6
 
bump = np.exp(-((x - 0.5) ** 2) / (2 * 0.01 ** 2))
wp = bump.copy()                              # 右向き不変量 (u0 + c0)
wm = bump.copy()                              # 左向き 不変量 (u0 - c0)
 
def upwind_step(w, a, dt):
    if a > 0:
        dwdx = (w - np.roll(w, 1)) / dx
    else:
        dwdx = (np.roll(w, -1) - w) / dx
    return w - a * dt * dwdx
 
CFL = 0.4
dt = CFL * dx / (abs(u0) + c0)
steps = int(0.25 / dt)
for _ in range(steps):                        # 時間前進ループ
    wp = upwind_step(wp, u0 + c0, dt)
    wm = upwind_step(wm, u0 - c0, dt)
 
T = steps * dt
def front_speed(w):
    return (x[np.argmax(w)] - 0.5) / T
print(f"u0+c0: target {u0 + c0:+.3f}, measured {front_speed(wp):+.3f}")
print(f"u0-c0: target {u0 - c0:+.3f}, measured {front_speed(wm):+.3f}")

u0+c0: target +1.600, measured +1.600
u0-c0: target -0.400, measured -0.400

二つの前面がそれぞれ $u_0+c_0$ と $u_0-c_0$ で動きます。右向き不変量は下流へ、左向き不変量は上流へ。 $u_0$ を $c_0$ より上に上げると、2行目の測定値も正に変わります。

空気の数字で間隔がはっきりします。 $\rho=1.225$ 、 $u=220$ m/s、 $P=10^5$ 、 $\gamma=1.4$ で音速は 338.06 m/s です。出力の固有値は $-118.06$ 、 $220$ 、 $558.06$ で、 $u-c$ 、 $u$ 、 $u+c$ と一桁も違いません。 $558.06-(-118.06)$ は 676.12 m/s で、音速 338.06 の 2 倍と同じです。真ん中の 220 は物質を運ぶ速度なので、音速とは無関係です。マッハ 0.651 は 1 より小さいです。だから $\lambda_-$ は負のままで、上流へ向かう信号が一つ残ります。

二つ目の計算は音速を 1 に置いた場合です。 $u_0=0.6$ なら目標値は $+1.600$ と $-0.400$ です。測定値もその数字です。 $-0.400$ は $0.6-1$ です。符号が負だということが、上流への到達を意味します。 $u_0$ が $c_0$ を超えるとこの行が正に変わり、時空の楔が右へ越えた状態と同じになります。

Riemann 不変量 $u \pm \dfrac{2c}{\gamma-1}$ で $\gamma=1.4$ なら $\gamma-1=0.4$ で、 $2/0.4=5$ です。音波特性に沿って一定な量は $u\pm 5c$ です。コードが運ぶ $w_\pm = u_1 \pm (c_0/\rho_0)\,\rho_1$ は、この量を小さな擾乱へ縮めた形です。 $\rho_0=c_0=1$ なので係数は 1 のまま残ります。

超音速出口で圧力を固定すると波が戻ってくる

影響領域の楔が垂直線を越える瞬間が $M=1$ です。それより上では、上流へ向かう特性がありません。出口が超音速なら三つの特性の符号はすべて外向きです。境界で決められる値はありません。圧力を固定すると、入ってこない特性を数値的に作り、反射波が領域の中へ戻ってきます。

亜音速の出口は違います。 $\lambda_-$ だけが内側を向くので、条件は一つで足ります。出口から戻る波を σ で抑える NSCBC は、その入ってくる特性に反射係数を掛けます。符号を数える前に圧力を課すと、マッハ 0.651 の空気でも上流の信号が出口処理に混ざります。