静止的恒星开始流动 — 极坐标 FVM 的 well-balanced 与角动量守恒

要让静止恒星保持静止,压力源须写成面积差 P·ΔS,而不是教科书里的 2P/r·V。

同样的方程,同样的斜率限制器,同样的 Riemann 求解器。只把网格从 cartesian 换成 polar,一颗静止恒星的中心就开始漏气。压力均匀。速度为零。

把 $\Delta r$ 滑条拉大, $2P/r \cdot V$ 会离开面积差 $P(S_{\text{out}} - S_{\text{in}})$ 。

r (cell radius) = 2.00Δr = 0.30P = 1.00

Naive 2P/r · V is only a Taylor approximation of P·(S_out − S_in); the residual grows with Δr, so coarser radial grids produce stronger spurious accelerations.

残差按 $O(\Delta r^{3})$ 变小,却到不了零。网格越粗,这条缝越早露出来。

罪魁是网格本身。曲线坐标给动量方程多加了几何项。其中一行离散写错,静止平衡就垮。把 φ-动量换成角动量的那一招,就挨在这行旁边。

极坐标给动量方程添加了哪些项

在极坐标 $(r, \theta, \varphi)$ 中,径向动量方程为

\frac{\partial \rho u}{\partial t} + \frac{1}{r^{2}} \frac{\partial (r^{2} \rho u^{2})}{\partial r} + \frac{\partial P}{\partial r} + \cdots - \frac{\rho (v^{2} + w^{2})}{r} = 0

其中 $u, v, w$ 分别是 $r, \theta, \varphi$ 方向速度, $P$ 是压力。与 cartesian 形式的两点差异:① 发散项被 $1/r^{2}$ 包裹成保守形式;② 右端出现像伪力一样的 $-\rho(v^{2}+w^{2})/r$ 。

这些伪力并非真正的加速。局部基(local basis)随位置改变方向,同一个动量矢量在坐标系中看起来像在旋转。但数值上它们以"力"的形式加入网格更新,加入的方式不同,结果就不同。

2P/r 的来源 — 面积差造出的假力#

把压力梯度 $\partial P/\partial r$ 塞进保守发散里,离散径向动量方程变为

\frac{\partial (\rho u V)}{\partial t} + \Delta_{r}(F_{rr} S_{r}) + \Delta_{\theta}(F_{r\theta} S_{\theta}) + \Delta_{\varphi}(F_{r\varphi} S_{\varphi}) = \left[ \frac{2P}{r} + \frac{\rho(v^{2}+w^{2})}{r} \right] V

左端的 $F_{rr} = \rho u^{2} + P$ 就是 Riemann 求解器直接返回的通量。右端的 $2P/r \cdot V$ 是今天的主角。

它源自把 $\partial P / \partial r$ 改写成 $\frac{1}{r^{2}} \partial (P r^{2})/\partial r - 2P/r$ :前一项被通量散度自然吸收,后一项 $-2P/r$ 移到右端成了源。逐网格离散后

\frac{2 P_{i}}{r_{i}} \cdot V_{i}

看起来无害。可是把"压力恒为常数、速度全为零"的静止解代入代码后, $\Delta_{r}(F_{rr} S_{r}) = P (S_{i+1/2} - S_{i-1/2})$ 并不为零。因为极坐标网格内外两面的面积本就不同。本应抵消它的 $2P/r \cdot V$ 只是 Taylor 近似,二者并不严格相等。残差留下来,每一步都化为加速。

回到打印出来的数。 $N=40$ ,半径从 1 到 4,步数 200。naive 的 $\max|u|$ 是 $1.834\times 10^{-3}$ 。balanced 是 $0.000\times 10^{0}$ 。初值是 $P=1$ 、 $\rho=1$ 、 $u=0$ 、 $\gamma=1.4$ 。 $\Delta r=(4-1)/40=0.075$ 。时间步长按代码写成 $\Delta t=0.4\,\Delta r/c_0$ , $c_0=\sqrt{\gamma P/\rho}$ 。这个步长走 200 步,假速度堆到 $1.834\times 10^{-3}$ 。平衡形式在同一循环里不动。

网格加到 200,速度残差按 $(\Delta r)^{2}$ 下降,却不是零。前面图里的残差是 $O(\Delta r^{3})$ 。差了一阶。单格源项误差是面积差的 Taylor 余项,相对体积更小。除以 $\rho V$ 变成加速度,再反复乘 $\Delta t$ 变成速度。速度残差看起来低一阶,原因在这里。换成 $P(S_{i+1/2}-S_{i-1/2})$ 之后,阶数比较就结束了。左端通量差和右端是同一个式子。

体积 $V=(r_R^{3}-r_L^{3})/3$ ,面积 $S=r^{2}$ 。 $\sin\theta\,\Delta\theta\,\Delta\varphi$ 按代码注释从两边拿掉。系数相同,抵消还在。只改面积定义, $P\cdot\Delta S$ 就和通量差对不上。用外积去填曲线单元的真实面积, $\partial x/\partial x$ 也会离开 1。这和极坐标里的 $2P/r$ 是同一种洞。离散散度和面积不成对,静止解就站不住。

改一行,平衡回来

办法很简单。把源项换成 和通量差完全相同的表达式:

\frac{2P}{r} V \;\longrightarrow\; P_{i} \left( S_{(r), i+1/2} - S_{(r), i-1/2} \right)

这样在静止解中,通量贡献与源项以机器精度逐位抵消。恒星不再漏。

同样的技巧适用于 $\theta$ -动量方程的 $\cos\theta P / (r \sin\theta)$ 源:

\frac{\cos\theta}{\sin\theta} \frac{P}{r} V \;\longrightarrow\; P_{i,j} \left( S_{(\theta), i, j+1/2} - S_{(\theta), i, j-1/2} \right)

$\varphi$ -动量方程则没有这种几何压力源 — 因为两个方位面积本就相等。

这是 well-balanced 方法最简单的范例。一句话总结:"用同样的离散表达式,把离散发散造出的假加速精确减回去。"

60 行 NumPy — 静止平衡撑得住吗#

在一维径向网格上以 $P = 1$ , $\rho = 1$ , $u = 0$ 起步,看看两种方式各自的动量残差。

import numpy as np
 
def polar_cell_volume(rL: float, rR: float) -> float:
    """球壳体积 (省略 sin θ Δθ Δφ)。"""
    return (rR**3 - rL**3) / 3.0
 
def polar_r_surface(r: float) -> float:
    """球坐标下 r²-加权面积。"""
    return r * r
 
def naive_pressure_source(P: float, r_mid: float, V: float) -> float:
    """教科书式 2P/r · V — 对 P=const 也会在机器精度下出错。"""
    return 2.0 * P / r_mid * V
 
def balanced_pressure_source(P: float, S_in: float, S_out: float) -> float:
    """与通量差分表达式相同,因此能精确抵消。"""
    return P * (S_out - S_in)
 
def radial_static_run(N: int = 40, scheme: str = "balanced", steps: int = 200):
    rIn, rOut = 1.0, 4.0
    dr = (rOut - rIn) / N
    rho = np.ones(N)
    u = np.zeros(N)
    P = np.ones(N)
    gamma = 1.4
    c0 = np.sqrt(gamma * P[0] / rho[0])
    dt = 0.4 * dr / c0
 
    r_face = np.linspace(rIn, rOut, N + 1)
    S = polar_r_surface(r_face)
    V = np.array([polar_cell_volume(r_face[i], r_face[i + 1]) for i in range(N)])
    r_mid = 0.5 * (r_face[:-1] + r_face[1:])
 
    for _ in range(steps):
        # 通过各面的压力通量 (P 恒定 → 精确值)
        flux_diff = P * (S[1:] - S[:-1])
        # 要检验的源项写法
        if scheme == "naive":
            src = naive_pressure_source(P, r_mid, V)
        else:
            src = balanced_pressure_source(P, S[:-1], S[1:])
        # 动量更新: d(ρ u V)/dt = -flux_diff + src
        du = dt / (rho * V) * (src - flux_diff)
        u += du
    return r_mid, u
 
r1, u_naive    = radial_static_run(N=40, scheme="naive",    steps=200)
r2, u_balanced = radial_static_run(N=40, scheme="balanced", steps=200)
print(f"naive    : max |u| = {np.abs(u_naive).max():.3e}")
print(f"balanced : max |u| = {np.abs(u_balanced).max():.3e}")

40 网格、200 步后的输出大致是

naive    : max |u| = 1.834e-03
balanced : max |u| = 0.000e+00

cells N = 40

Initial state: P = 1, ρ = 1, u = 0 across r ∈ [1, 4]. The pressure flux through each spherical face is P · S_r. The discrete source on the right-hand side should cancel it exactly. Toggle the scheme to see who fails.

停在 naive 上加 step,曲线会离开零轴。换成 well-balanced,曲线贴住零。把 cells 调小,只有 naive 的残差变大。40 个网格、200 步的 $\max|u|=1.834\times 10^{-3}$ 就是这次漂移,balanced 打印 $0.000\times 10^{0}$ 。

把 φ-动量替换成角动量#

几何源不止出现在压力项。 $\varphi$ -动量方程里还有 $\rho u w / r + \cos\theta \rho v w / r$ 等。在恒星、吸积盘等强旋转系统中,这些项同时侵蚀精度和守恒性。

办法是换变量。定义 $l = w r \sin\theta$ — 也就是 角动量 本身 — 于是 $\varphi$ -动量方程里的源项 全部消失:

\frac{\partial \rho l}{\partial t} + \frac{1}{r^{2}} \frac{\partial (r^{2} \rho u l)}{\partial r} + \frac{1}{r \sin\theta} \frac{\partial (\sin\theta \rho v l)}{\partial \theta} + \frac{\partial (\rho w^{2} + P)}{\partial \varphi} = 0

离散上,把 Riemann 求解器返回的 $F_{\varphi r}$ , $F_{\varphi\theta}$ , $F_{\varphi\varphi}$ 乘以 $r \sin\theta$ 得到 $\tilde F$ ,再做差分。无须修正、无须源项的处理。即便网格较粗,角动量也由式子本身保证守恒(Kley 1998, A&A 338, L37)。

但在非均匀旋转(如 Keplerian 盘)中,还需移到自转参考系并配合 FARGO 技巧才能放松 CFL。那是另一篇的题目。

两种写法的选用指南

关注点	naive `2P/r · V`	well-balanced `P · ΔS_r`	角动量形式
静止平衡	❌ $O(\Delta r^{2})$ 残差	✅ 比特级精确	✅ (仅 φ 方向)
额外代码改动	无	替换一行源项	新的通量 + 新变量
旋转系统	⚠ 误差累积	△ 仅修压力源	✅ 天然守恒
调试代价	高(无声崩溃)	低	中
会在哪里碰到	"为什么星核在掏空?"	spherical hydro 标配	吸积盘

再次站在恒星面前

笛卡尔网格永远看不到这个陷阱。曲线坐标系往动量方程里加项,如果这些项的离散写法不与 通量差分形式一致,静止气体就不会保持静止。本文要改的就是这一行:

\frac{2P}{r}\,V \;\longrightarrow\; P\,(S_{i+1/2} - S_{i-1/2})

一行,一次提交。恒星再次安静。

留着 2P/r,还是换成 ΔS#

如果只是压力平衡在晃,把 $2P/r\,V$ 换成 $P(S_{i+1/2}-S_{i-1/2})$ 。40 个网格上 $1.834\times 10^{-3}$ 的 $\max|u|$ 会变成 0。 $\theta$ 方向的 $\cos\theta\,P/(r\sin\theta)$ 同样换成 $P\cdot\Delta S_{\theta}$ 。 $\varphi$ 方向两面面积相同,不需要这一行。

旋转很强时,压力源项不够。 $\rho u w/r + \cos\theta\,\rho v w/r$ 还留在 $\varphi$ -动量里。改成 $l = w r \sin\theta$ ,这项源就从方程里消失。这是 Kley 1998 (A&A 338, L37) 的形式。旋转随半径变化,还得再加 FARGO。静止界面自己流起来,和寄生电流是同一种病。离散散度造出的假力,不用同一个离散式减回去,平衡就不会回到原处。