特征线相撞之处,激波诞生 — Burgers 方程与 (x,t) 平面

1948 年,荷兰物理学家 J. M. Burgers 抛出了一个一行的方程作为湍流模型。一条光滑曲线在没有外力的情况下,有限时间内被撕成一道垂直悬崖 — 这是非线性双曲型 PDE 中最短的戏剧。本文用特征线(characteristic,信息在时空中流动的曲线)的几何手算激波形成的时刻,梳理弱解(weak solution,在可微性失效处由积分形式定义的解)如何用一条线截断多值区域,然后用 60 行 Python 重现同样的图像。文末你可以亲手在 (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$ 为零,所以 $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 方法捕捉弱解激波。

import numpy as np
 
def burgers_breaking_time(q0_prime_min):
    # q0_prime_min: q0'(x) 在 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,破裂在两处同时发生;一旦破裂,两道激波合并成更大的一道。上面板里青色曲线与浅色虚线(初始曲线)之间的面积相等 — 等面积律 — 可以肉眼直接确认。

数值激波捕捉失败的两个原因#

(1) CFL 越界. 违反 $\Delta t \le \Delta x / \max|q|$ 时,数值域追上物理域。若振幅看似随时间增长,先怀疑 $dt$ 。像上面代码那样每步用 np.max(np.abs(q)) 重新计算 $dt$ 更稳。激波形成后 $q$ 可能局部突起,这一点尤为重要。

(2) 熵处理错误. 仅满足 Rankine–Hugoniot,非物理的 expansion shock 也会成为解。Roe 这类简单的线性化求解器在跨音速 rarefaction(扇形跨过零)区域会出现"entropy glitch"。Godunov、HLLE、精确 Riemann 求解器自动满足熵条件;若使用 Roe,需要显式的 Harten 形式的 entropy fix。

可压缩 LES 代码出现激波位置振荡或负压时,几乎都是这两种之一。

下次堵车开始时#

车并入一条道,速度下降。前车慢,后车快。后车碰到前车保险杠的瞬间,堵车的头形成。这就是 (x,t) 平面上两条轨迹相遇的点。道路堵塞实际上就是 Burgers 激波,在解开之前以 Rankine–Hugoniot 速度向后传播。

一行总结。非线性双曲型 PDE 即使从光滑初值出发,也会在有限时间内产生间断。接受这一事实,通过弱解这条迂回路重新获得意义,正是激波数值分析的起点。