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

非線形な自己移流が、なめらかな曲線をどう衝撃に砕くのか

1948 年、オランダの物理学者 J. M. Burgers は乱流のモデルとして一行の方程式を投げた。なめらかな曲線が外力なしに有限時間で垂直の崖に砕けるという、非線形双曲型 PDE の最短ドラマです。本稿では、特性線(characteristic、情報が流れる時空曲線)の幾何から衝撃形成時刻を手計算し、弱解(weak solution、微分可能性が壊れる場所で積分形により定義される解)が多価領域を一本で切り落とす仕組みを整理し、Python 60 行で同じ図を再現します。最後に、(x,t) 平面で直線がぶつかる様子を直接いじることができます。

直線が同じ点で出会う#

双曲型方程式は信号伝搬の話です。情報が一定速度で特性線に沿って流れる限りは平和です。線形移流 $\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^*$ を越えた瞬間に垂直の壁が立ち上がります。

弱解と 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 が答えとして残ります。

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 問題の厳密解を計算し、弱解を保証します。

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

下のシミュレーションで振幅と波数を変えてみてください。

Amplitude A: 0.35Wave number k: 1

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

振幅 $A$ を 0.45 まで上げると $t^*$ は 0.35 付近まで縮み、赤い破線が上に上がります。波数 $k$ を 2 にすると 2 か所で同時に砕け、一度砕けたあとは二つの衝撃が合流してより大きな衝撃へ育ちます。上のパネルで、シアンの曲線と薄い破線(初期曲線)の間の面積が等しいこと — equal-area rule — が直接目で確認できます。

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

(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 速度で後ろへ伝播します。

一行まとめ。非線形双曲型 PDE は、なめらかな初期値からも有限時間で不連続を作る。その事実を受け入れ、弱解という回り道で意味を取り戻すことが、衝撃波の数値解析の出発点である。