特性線がぶつかると衝撃波が生まれる — Burgers 方程式と (x,t) 平面

衝撃速度は両状態の算術平均です。左が右より小さい不連続は、衝撃ではなく膨張波です。

1948 年、Burgers は乱流のモデルとして方程式を一行投げました。外力なしに、なめらかな曲線が有限時間で崖になります。その崖の誕生日は、特性線(characteristic、情報が流れる時空曲線)が交わる時刻です。弱解(weak solution、微分が壊れたあとも積分形で残る解)は、その重なりを一本で切ります。

振幅 $A$ と波数 $k$ を上げると、赤い破線 — 砕ける時刻 — が上へ動きます。

Amplitude A: 0.35Wave number k: 1

q0(x) = 0.55 + A sin(2 pi k x)

その破線は、山から出た直線が前の谷の直線を踏む最初の時刻です。 $A$ を大きくすると交点は早まり、 $k = 2$ では二か所で同時に重なります。

直線が同じ点で出会う

双曲型方程式は信号伝搬の話です。情報が一定速度で特性線に沿って流れる限りは平和です。線形移流 $\partial_t q + u\,\partial_x q = 0$ で $u$ が定数なら、特性線は平行な直線で初期形状はそのまま横に滑ります。この絵には衝撃波の入る余地がありません。

非線形性が入ると一変します。 $u$ が $q$ 自身に依存すると、大きな値は速く、小さな値はゆっくり進みます。速い部分が遅い部分を追い越した瞬間、(x,t) 平面の同じ点に二つの値が同時に届く。関数が多価(multi-valued)になる、すなわち物理的にあり得ない状態になる、その瞬間が衝撃波の誕生です。

自分自身を運ぶ方程式

Burgers 方程式の保存形は次の通りです。

\partial_t q + \partial_x \left(\tfrac{1}{2} q^2\right) = 0

$q$ は保存スカラー、 $f(q) = q^2/2$ がフラックスです。連鎖律で展開すると非保存形が得られます。

\partial_t q + q\,\partial_x q = 0

移流速度は $u(x,t) = q(x,t)$ 。自分自身が自分自身を運ぶ構造です。この一つの非線形性がすべてのドラマの源です。ヤコビアン $\partial f/\partial q = q$ が正の場所では信号は右へ、負の場所では左へ進みます。空間で符号が変わるような分布 — 正弦波がその典型 — では、速い場所が遅い場所を追いかけます。

特性線が描く (x,t) 平面#

特性線は $dx/dt = q$ を満たす曲線です。共動微分 $Dq/Dt$ がそこで 0 なので、 $q$ は出発値で凍結されます。結果として特性線は、(x,t) 平面における 傾き $1/q_0(x_0)$ の直線 です。出発点 $x_0$ における初期値 $q_0(x_0)$ がその直線の速度になります。

初期条件 $q_0(x_0) = 0.55 + A\sin(2\pi k x_0)$ では:

山( $q_0 > 0.55$ )から出る直線は急に右へ傾く(速い)
谷( $q_0 < 0.55$ )から出る直線は緩やかに傾く(遅い)

山の直線はやがて、すぐ前の谷の直線のつま先を踏みます。二本が一点で交わる瞬間、それが最初の衝撃です。

砕ける時刻 t* — 最初の衝突#

隣接する出発点 $x_0$ 、 $x_0 + \Delta x$ の特性線が時刻 $t$ に出会うには次が成立しなければなりません。

x_0 + q_0(x_0)\,t = (x_0 + \Delta x) + q_0(x_0 + \Delta x)\,t

$\Delta x \to 0$ の極限を取ります。

1 + q_0'(x_0)\,t = 0 \;\;\Longrightarrow\;\; t = -\frac{1}{q_0'(x_0)}

負の傾きの場所だけが交わります。最も早い衝突は $q_0'(x_0)$ が最も負になる点で起こります。

t^* = \frac{1}{\max_{x_0}\big(-q_0'(x_0)\big)}

$q_0$ の振幅と波数が大きいほど早く砕けます。上の正弦の場合、 $q_0'(x_0) = 2\pi k A\cos(2\pi k x_0)$ 、最小値は $-2\pi k A$ です。

t^* = \frac{1}{2\pi k A}

$A = 0.35$ 、 $k = 1$ なら $t^* \approx 0.455$ 。それ以前は曲線がなめらかに非対称化していき、 $t^*$ を越えた瞬間に垂直の壁が立ち上がります。

傾きの最小値は $-2\pi k A$ です。 $2\pi \approx 6.28$ なので $2\pi \times 0.35 \approx 2.20$ 、その逆数が $0.455$ です。山は $0.55 + 0.35 = 0.90$ 、谷は $0.20$ です。最初の高さ $0.70$ が、あとで出る shock height 約 $0.7$ と揃います。 $t = 1.2$ まで進んでも跳びは育ちません。角が丸くなるだけです。

$A$ を $0.45$ にすると分母は約 $2.83$ になり、 $t^*$ は $0.35$ 付近まで縮みます。波数も $2$ にすると分母は倍、 $t^* \approx 0.227$ です。振幅が同じでも波長が短いほど早く重なります。

弱解と Rankine–Hugoniot — 多価を切る刃#

$t > t^*$ で素朴な解析は同じ $x$ に三つの $q$ を返します。物理的に不可能です。打開策が弱解です。微分形を捨て、保存則の積分形で定義をゆるめます。

\frac{d}{dt}\int_{x_L}^{x_R} q\, dx \;=\; f\big(q(x_L)\big) - f\big(q(x_R)\big)

これは $q$ が微分不可能 — つまり不連続 — でも意味を持ちます。二つの平坦な状態 $q_L$ 、 $q_R$ が衝撃で結ばれ速度 $s$ で進むとすると次が従います。

s = \frac{f(q_L) - f(q_R)}{q_L - q_R} = \frac{q_L + q_R}{2}

これが Rankine–Hugoniot 条件 です。衝撃速度は両状態の算術平均。この条件が多価領域を一本で切り、両側の切り取り面積が等しくなる(equal-area rule)ように調整します。

加えてエントロピー条件が必要です。物理的な衝撃は $q_L > q_R$ のときだけです。 $q_L < q_R$ なら扇形の rarefaction(膨張波)で解消されます。これを破ると、数学的には弱解だが非物理的な expansion shock が答えとして残ります。

$s = (q_L + q_R)/2$ に $0.90$ と $0.20$ を入れると $s = 0.55$ です。正弦の背景値と同じです。衝撃はその平均速度で滑り、前後の状態だけが残ります。順番を逆にしても速度の式は $0.55$ を返します。変わるのは答えだけで、衝撃が膨張波になります。一つの算術平均は二種類の波を見分けません。

格子幅は $1/256 \approx 0.00391$ です。CFL $0.4$ と初期の $\max|q| = 0.90$ なら $\Delta t \approx 0.00174$ です。 $t = 1.2$ までは約 $690$ ステップです。衝撃のあと $\max|q|$ が局所的に膨らむと、同じ式が歩幅をさらに刻みます。セル境界のその刃は、一セルの衝撃管と同じものです。

Python:特性線追跡から Godunov 衝撃まで#

特性線を直接描き、Godunov 1 次スキームで弱解の衝撃を捕捉します。

import numpy as np
 
def burgers_breaking_time(q0_prime_min):
    # q0_prime_min: x にわたる q0'(x) の最小値(衝撃発生には負である必要)
    if q0_prime_min >= 0:
        return float('inf')
    return -1.0 / q0_prime_min
 
def trace_characteristics(q0_func, x0_arr, t_arr):
    # x[i, j] = x0_arr[i] + q0(x0_arr[i]) * t_arr[j]
    return x0_arr[:, None] + q0_func(x0_arr)[:, None] * t_arr[None, :]
 
def rankine_hugoniot_speed(qL, qR):
    # 衝撃速度 s = (qL + qR) / 2
    return 0.5 * (qL + qR)
 
def godunov_burgers_step(q, dx, dt):
    # 保存形: q^{n+1} = q^n - (dt/dx) * (F_{i+1/2} - F_{i-1/2})
    qL = q                          # 右界面の左状態 = 現在のセル値
    qR = np.roll(q, -1)             # 右界面の右状態 = 次のセル値
    s = 0.5 * (qL + qR)
    f_shock = np.where(s >= 0, 0.5 * qL**2, 0.5 * qR**2)
    f_rare  = np.where(qL >= 0, 0.5 * qL**2,
              np.where(qR <= 0, 0.5 * qR**2, 0.0))
    F_iph = np.where(qL > qR, f_shock, f_rare)   # 右界面のフラックス
    F_imh = np.roll(F_iph, 1)                    # 左界面 = 一つシフト
    return q - (dt / dx) * (F_iph - F_imh)
 
# シミュレーションパラメータ ----------------------------------
N, L = 256, 1.0
A, k = 0.35, 1
x = (np.arange(N) + 0.5) * (L / N)
q0 = lambda xx: 0.55 + A * np.sin(2 * np.pi * k * xx)
q  = q0(x)
 
t_star = burgers_breaking_time(-2 * np.pi * k * A)
print(f"breaking time t* = {t_star:.3f}")
 
t_end, t = 1.2, 0.0
while t < t_end:
    dt = 0.4 * (L / N) / np.max(np.abs(q))    # CFL = 0.4
    q  = godunov_burgers_step(q, L / N, dt)
    t += dt
print(f"shock height (max - min) at t = {t:.2f}: {q.max() - q.min():.3f}")

breaking time t* = 0.455、衝撃高さは約 0.7 と表示されます。Godunov スキームは多価領域を自動的に一本へ切り落とします。その刃が np.where(qL > qR, f_shock, f_rare) の一行です。この一行が各セル界面で Riemann 問題の厳密解を計算し、弱解を保証します。

直線がぶつかる様子を見る

上で動かした図では、シアンの曲線と薄い初期の破線が、衝撃の左右で切り取る面積が等しくなります。equal-area rule です。 $k = 2$ で交点が二か所なら、二つの衝撃はすぐ一つに合流し、跳びが大きくなります。

面積が等しいとは、切ったあとも $\int q\,dx$ が両端のフラックス差と合うということです。山 $0.90$ と谷 $0.20$ の平均 $0.55$ で衝撃が歩くと、切り取った二つの三角形は打ち消し合います。高さ $0.70$ はそのまま、位置だけが $s = 0.55$ で動きます。

数値的な衝撃捕捉が失敗する二つの理由

(1) CFL 違反. $\Delta t \le \Delta x / \max|q|$ を破ると、数値領域が物理領域を追い越します。振幅が時間と共に育つように見えたら、まず $dt$ を疑いましょう。上のコードのように毎ステップ np.max(np.abs(q)) で $dt$ を更新するのが安全です。衝撃形成直後は $q$ が局所的にスパイクするため、特に重要です。

(2) エントロピーの取りこぼし. Rankine–Hugoniot だけでは非物理的な expansion shock も解になり得ます。Roe のような単純な線形化スキームは、ファンが 0 をまたぐ transonic rarefaction で「entropy glitch」を生みます。Godunov、HLLE、exact Riemann solver はエントロピー条件を自動で満たしますが、Roe を使うなら明示的な Harten 形のエントロピー修正が必要です。

圧縮性 LES コードで衝撃位置が振動したり負圧が現れたりする場合、ほぼこのどちらかです。

渋滞を衝撃と見るか

車が一本の車線に集まり、速度が落ちます。前は遅く、後ろは速い。バンパーが触れる点が (x,t) 平面の交点です。その頭は Burgers の衝撃で、解消されるまで Rankine–Hugoniot 速度で後ろへずれます。

渋滞を連続曲線のまま抱えていると、多価が残ります。一本で切るか、膨張でほどくかは、左右の大小が決めます。段差で振動が見えたら、先に疑うのは流束制限です。