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

超音速戦闘機は、自分のエンジン音より先に到着します。頭上を通り過ぎたあとで、ようやく轟音が追いかけてくるのです。なぜこんなことが起きるのかは、実はEuler方程式の中にすでに書き込まれています。この記事では、小さな擾乱がどうやって三本の信号に分かれるのかを導きます。その信号が $(x,t)$ 平面に描く特性曲線(信号が通り抜けた軌跡)を操作可能な図で確かめ、最後には三つの特性速度をPythonで直接取り出してみます。

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

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}$ が一定に保存され、物質特性に沿ってはエントロピーが保存されます。

下のアニメーションで背景流れ $u_0$ を直接操作してみましょう。

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.

$M = u_0/c_0$ を0.4に置くと、青い信号(左向きの音波)がsourceの左、つまり上流へ伸びていきます。 $M$ を1より上に上げると、その青い軌跡は消えます。 $\lambda_- = u_0 - c_0$ が正に変わり、三つの信号がすべて下流にしか走らなくなるからです。

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

同じ話を時空図で見ると、さらに鮮明になります。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行目の測定値も正に変わります。

特性曲線が残すもの

圧縮性流れの動力学は、三つの特性速度 $u-c,\ u,\ u+c$ で走る信号の伝播です。二つは音波、一つは物質の移動。
影響領域の楔が垂直線を越える瞬間が $M=1$ です。それより上では、上流へ向かう信号がありません。
超音速境界で特性がすべて一方向なら、課せる境界条件の数もその分だけ決まります。特性の符号を数える習慣が、境界条件のミスを防ぎます。