1950年,冯·诺依曼故意把激波抹糊了 — 人工粘性与奇偶解耦

洛斯阿拉莫斯,1950 年。John von Neumann 与 Robert Richtmyer 在 Journal of Applied Physics 发了一篇短文,题为 "A Method for the Numerical Calculation of Hydrodynamic Shocks"。激波若不能干净地落在一格里,就先把它抹开。技巧就是这一行。

下面是 Sod 管,$t = 0$ 时抽掉隔膜。左侧 $\rho = 1, p = 1$,右侧 $\rho = 0.125, p = 0.1$。把 ξ 滑块从 0 拨到 5,激波后面的锯齿会熄掉,间断则铺到好几格。

ξ = 0 时激波背面是锯齿，ξ = 2 附近振荡消失、激波剩两三格，ξ = 5 时平台区也发软。冯·诺依曼推荐的 $\xi \approx 2$-$3$,正是这场权衡的甜区。

人工粘性为什么生效、留下多厚、以及同位网格上的奇偶解耦,是这篇笔记的范围。末尾用 NumPy 抓住一个 Mach 2 激波管。

网格上的激波会振荡

理想 Euler 方程不含粘性,真实激波被压缩到分子平均自由程量级。CFD 网格的一格大约比那个宽一百万倍。密度和压力要在一格之内翻倍。

把这种间断送进任何守恒型差分,几乎必然出现一种现象:激波后面冒出的锯齿振荡,通常叫 wiggle。von Neumann 在 1944 年曼哈顿计划时期就遇上了。爆炸模拟的激波后区一旦起伏,TNT 当量计算就整体偏离真值。

原因不复杂。一阶迎风带来强数值粘性,把激波按住但顺手把其他流动也抹平。二阶以上方案精度高,但在间断附近像 Gibbs 现象那样振荡。自然依靠分子粘性在激波内部产生熵,而网格并无此通道。把一格里的激波管当成黎曼问题来解,是 Godunov 在每一格里放一根激波管。人工粘性走另一条路:把间断拉开大约三格,让振荡死去。

1950年的一行式#

冯·诺依曼与 Richtmyer 的回答是模仿自然。既然分子粘性无法直接分辨,那就在网格尺度上加一个假的体积粘性项,并入压力。

$Q_i$ 是单元 $i$ 上的人工体积压力,$\xi$ 是调谐参数,决定激波铺多少格(通常取 2-3),$u_{i+1} - u_{i-1}$ 是速度一阶差分,$\rho_i$ 是单元密度。在动量与能量方程中把每个 $P$ 替换为 $P + Q$,工作就此完成。

$Q$ 与 $(\Delta u)^2$ 成正比。光滑流里几乎为零,只在激波附近变大。开关是 $u_{i+1} \le u_{i-1}$,只在速度局部汇聚的格子上点着。膨胀波和声波不动。乘上 $\rho_i$,重介质得到更大的压力推动。

正文里的 ξ 和代码里的 0.25 不是同一个厚度#

缩写成 $Q = \xi^2 (\Delta u)^2 \rho$。只在 $\Delta u < 0$ 时打开。下面的列表乘的是 0.25 * xi**2 * du**2 * rho。du 跨过一格。改成单格斜率时,分母里有 $2\Delta x$,平方之后带上 $1/4$。代码里的 0.25 就是这个中心差分。上面的式子把这个系数丢掉,只留下 $\xi^2$。

$\xi = 2.5$ 时 $\xi^2 = 6.25$。代码实际乘上的数是 $0.25 \times 6.25 = 1.5625$。去掉 0.25、又用同一个 $\xi$,$Q$ 变为四倍。$Q$ 正比于 $\xi^2$,四倍等于把 $\xi$ 加倍。滑块上 $\xi = 2$ 时安静下来的激波,若按正文公式原样实现,对应的是 $\xi = 4$。$\xi = 5$ 时平台已经发糊。把 $\xi = 2.5$ 不带 0.25 代进去,就会落到那一侧。

200 个格子、时间 0.2、CFL 0.45 都写在列表里。输出那句话只说峰值密度落在 Mach 2 的 Rankine–Hugoniot 值 2.667 的几个百分点之内。拿去掉 0.25 的公式去对这句话,量到的激波厚度已经不是同一个。

用NumPy复刻这一招#

五十行装下核心。donor-cell 对流加人工粘性项,仅此而已。

import numpy as np
 
def vnr_shock_solver(N=200, T=0.2, xi=2.5, gamma=1.4):
    # 初始条件:Sod激波管
    dx = 1.0 / N
    rho = np.where(np.arange(N) < N // 2, 1.0, 0.125)
    mom = np.zeros(N)
    p   = np.where(np.arange(N) < N // 2, 1.0, 0.1)
    E   = p / (gamma - 1)
    t = 0.0
    while t < T:
        u = mom / rho
        c = np.sqrt(gamma * p / np.maximum(rho, 1e-9))
        dt = 0.45 * dx / np.max(np.abs(u) + c)
        # von Neumann-Richtmyer 的 Q
        du = np.zeros(N)
        du[1:-1] = u[2:] - u[:-2]
        Q = np.where(du < 0, 0.25 * xi**2 * du**2 * rho, 0.0)
        # donor-cell 对 rho, mom, E 做对流
        uf = 0.5 * (u[:-1] + u[1:])
        for q_arr, q_new in [(rho, 'rho'), (mom, 'mom'), (E, 'E')]:
            flux = np.where(uf > 0, q_arr[:-1] * uf, q_arr[1:] * uf)
            q_arr[1:-1] -= dt / dx * (flux[1:] - flux[:-1])
        # 压力更新与 P+Q 的体积力
        p = (gamma - 1) * (E - 0.5 * rho * (mom / rho) ** 2)
        Ptot = p + Q
        mom[1:-1] -= dt / (2 * dx) * (Ptot[2:] - Ptot[:-2])
        E[1:-1]   -= dt / (2 * dx) * (Ptot[2:] * u[2:] - Ptot[:-2] * u[:-2])
        t += dt
    return rho, mom, E
 
rho, mom, E = vnr_shock_solver(xi=2.5)
print(f"shock peak rho ≈ {rho.max():.3f}   (Rankine-Hugoniot exact ≈ 2.667 for Mach 2)")

运行后峰值密度落在 Rankine-Hugoniot 解的百分之几以内。把同一段代码改成 ξ = 0,振荡会让最大值与最小值非对称地偏离真实平台。

奇偶解耦:另一个陷阱

抓住激波并不意味着收工。同一个 collocated 网格(所有变量放在单元中心)还藏着冯·诺依曼注意到的另一个毛病。看这个静止状态:速度为 0、比内能恒为 1,但密度按锯齿排列。

压力 $p_i = \rho_i$(等温假设)使单元 $i = 2$ 的压力是邻居的两倍,直觉上应当立刻松弛。然而单元 $i$ 的动量更新为

原因在于 $p_{i+1}$ 与 $p_{i-1}$ 来自同一奇偶子网格,差正好抵消。网格分裂成两个互不相见的独立子网格(偶单元、奇单元)。此锯齿模式既不增长也不衰减。

把上面的模式从 collocated 切到 staggered。staggered 网格把动量放在单元面上。面 $i+1/2$ 处的压力差是 $p_{i+1} - p_i$,奇偶被迫在相邻单元相遇。锯齿立刻落入一阶差分的射程并被阻尼。

失与得

冯·诺依曼 1950 年的技巧同时展示了两件事。网格极限可以靠模仿弥补。每一次模仿又制造一种新的网格病。两条经验在今日 CFD 中仍然活跃。

ZEUS、Athena、PLUTO 等天体物理代码沿用 staggered 网格加人工粘性;OpenFOAM、ANSYS Fluent 则在 collocated 网格上叠一层 Rhie-Chow 插值绕开同一个解耦。同一个问题的两种侧攻。

激波之后的锯齿未必是 bug。它是一种自 1944 年就被知晓的网格病的另一面。

去掉 0.25 的 ξ=2.5 会落到滑块的 ξ=5 一侧#

拿掉 0.25,把 $\xi = 2.5$ 代进正文公式,$Q$ 的大小就和代码在 $\xi = 5$ 时一样。滑块上那个位置,激波已经钝了。$\xi \approx 2$-$3$ 是代码还带着 0.25 时的位置。系数不带上、只搬数字,厚度的说法就变了。

偶数格和奇数格互相看不见,这和人工粘性是两件事。速度放在单元面上,压力才看得到相邻一格。这种布置因激波后锯齿以外的原因出错,记在 交错网格上旋度约束被破坏的情形。