特征线相撞之处,激波诞生 — 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$ 为零,所以 $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$ 。这和后来大约 $0.7$ 的 shock height 对得上。走到 $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 作为答案。

把 $0.90$ 和 $0.20$ 代入 $s = (q_L + q_R)/2$ ,得到 $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 方法捕捉弱解激波。

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 问题的精确解,从而保证弱解。

亲眼看直线相撞

前面拨过的图里,青色曲线和浅色初始虚线在激波两侧切下的面积相等。这就是等面积律。 $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 这类简单的线性化求解器在跨音速 rarefaction(扇形跨过零)区域会出现"entropy glitch"。Godunov、HLLE、精确 Riemann 求解器自动满足熵条件;若使用 Roe,需要显式的 Harten 形式的 entropy fix。

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

堵车要不要当成激波

车挤进一条车道,速度掉下来。前车慢,后车快。保险杠碰上的那一点,就是 (x,t) 平面上的交点。这个头是 Burgers 激波,解开之前按 Rankine–Hugoniot 速度往后挤。

把堵塞当成连续曲线抱着,多值就还在。切成一条线,还是放开成膨胀波,要看左右谁大。台阶附近若开始振荡,先去看通量限制器。