超声速流动感知不到上游 —— Euler 方程的特征线与声波

超声速战机比它自己的引擎声更早抵达。它从头顶掠过之后，轰鸣声才姗姗追来。这件事为什么会发生，其实早已写进了 Euler 方程里。本文推导一个微小扰动是如何分裂成三路信号的。我们会用一张可以动手操作的图，看看这些信号在 $(x,t)$ 平面上画出的特征线（信号经过留下的轨迹），最后再用 Python 亲手把三个特征速度算出来。

可压缩方程描述的是信号的传播

一维无粘可压缩流动写成三个守恒方程。

\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$ 后只保留一阶项。定义随背景一起运动的微分算子 $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$ 移动。这就是熵（或接触面）特征。归纳起来，一维流动共有三个特征速度。

\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$ 。把一维 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):
    """一维 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$ 以上，第二行的测量值也会变成正数。

特征线留下的东西

可压缩流动的动力学，就是沿三个特征速度 $u-c,\ u,\ u+c$ 奔跑的信号的传播。其中两个是声波，一个是物质的搬运。
影响域楔形越过竖直线的那一刻，就是 $M=1$ 。在此之上，没有任何朝上游走的信号。
在超声速边界上，若特征全部指向同一方向，能施加的边界条件数目也就随之确定。养成数一数特征符号的习惯，能避免边界条件上的失误。