当压力振荡成棋盘格 — SIMPLE 与 Rhie-Chow 压力-速度耦合

把 α_p 设成 1、每步整格换掉压力，校正循环会过校正然后炸掉。

速度场看起来没问题。压力场一打开，每个网格都在上下跳。棋盘格（checkerboard）。残差在降。只有压力是锯齿。不是代码错误。压力和速度放在同一网格点的那一刻，结果就已经定了。

在 Naive 和 Rhie-Chow 之间切换，种下的锯齿幅度要么停着，要么沉到光滑曲线上。

Naive interpolation happened to settle — nudge the grid to see the sawtooth return. (sweeps: 0)

Cells N: 32

Naive 下幅度不缩。换成 Rhie-Chow，同一块棋盘格就沉下去。Reseed checkerboard 把锯齿重新种上。

锯齿为什么还在、连续方程怎样变成压力方程，在后面。同一个求解器延伸到可压缩，也在后面。

连续方程看不到邻格的压力

在同位网格（速度与压力都存储在网格中心）上，半离散动量方程写作

a_P \mathbf{u}_P = \sum_{nb} a_{nb}\mathbf{u}_{nb} - V_P (\nabla p)_P = \mathbf{H}(\mathbf{u}) - V_P (\nabla p)_P

其中 $a_P$ 是对角系数， $\mathbf{H}(\mathbf{u})$ 汇集邻格贡献， $V_P$ 是网格体积， $(\nabla p)_P$ 是网格中心的压力梯度。

麻烦就在这个中心梯度。用中心差分写成 $(\partial p/\partial x)_P = (p_E - p_W)/2\Delta x$ 。注意其中没有 $p_P$ 。网格看不到自己的压力，只看到隔一格的邻居。于是奇数格与偶数格彼此分离（奇偶解耦）。 $(+,-,+,-,\dots)$ 的锯齿压力藏在离散方程的零空间里，完全不产生残差。

Rhie–Chow：用动量构造面速度#

对策是不要用网格速度的简单平均来求面速度。在面上重新运用动量方程，让该面的邻格压力差直接进入。

u_f = \overline{\left(\frac{\mathbf{H}(\mathbf{u})}{a_P}\right)}_f - \left(\frac{V}{a_P}\right)_f \frac{p_E - p_P}{\Delta x}

上划线是线性插值。后面那一项是在面上重算的压力梯度。现在面速度直接系于 $p_E - p_P$ ，也就是相邻两格的压力差。锯齿模态再也无处可藏。

只要面上的 $p_E-p_P$ 不是零，锯齿就被这一项抓住。从动量里把压力拿掉的另一条路，是 Chorin的投影法。

SIMPLE：把连续方程变成压力方程#

不可压缩流动的根本困难在于没有一个直接求解压力的方程。SIMPLE（Semi-Implicit Method for Pressure Linked Equations）把连续方程 $\nabla\cdot\mathbf{u}=0$ 变成压力方程。把 Rhie-Chow 面速度代入连续方程，就得到关于压力的椭圆型（Poisson）方程。

\nabla\cdot\left[\left(\frac{V}{a_P}\right)_f \nabla p\right] = \nabla\cdot\left(\frac{\mathbf{H}(\mathbf{u})}{a_P}\right)_f

左端是拉普拉斯算子（扩散形式），右端是预测速度的散度。流程为预测-校正。

用猜测的 $p^*$ 求解动量，预测 $\mathbf{u}^*$ 。
求解上述压力方程，更新 $p$ 。
用新的压力梯度校正面速度与网格速度，使连续性成立。
重复 1–3 直到收敛。

PISO 在此基础上每步再多做两三次校正，用于瞬态计算。PIMPLE 为大时间步把外层 SIMPLE 循环与内层 PISO 循环叠合。

没有欠松弛就会发散

有个陷阱。若每步完全替换压力（ $p \leftarrow p^*$ ），耦合会过校正而发散。SIMPLE 只让更新的一部分通过。

p^{\text{new}} = p^{\text{old}} + \alpha_p\,(p^* - p^{\text{old}})

$\alpha_p$ 是压力欠松弛因子（0 到 1）。速度另有一个 $\alpha_u$ 。太小则收敛缓慢，太大则振荡并崩溃。经验法则是 $\alpha_p \approx 0.3$ 、 $\alpha_u \approx 0.7$ ，使两者之和约为 1。

在下面的实验里亲自调节松弛因子。它用带松弛的 Gauss-Seidel 求解压力校正方程。

converging — sweep 0, residual 0.0e+0.

Relaxation factor w (≈ SIMPLE α_p): 1.50

把因子降到 0.3，残差曲线便缓慢地爬着下降。在 1.5 附近下降最快。逼近 2 时，每次扫掠的过校正不断累积，残差反而跳回上去。SIMPLE 的 $\alpha_p$ 也恰好处在同样的平衡之上。

一个求解器通吃所有速度 — 向可压缩延伸

压力基方法真正的魅力在于它不挑马赫数。可压缩时连续方程变为密度方程，用状态方程（EOS）把密度系于压力。引入 $\psi \equiv \partial\rho/\partial p$ （压缩率， $=1/c^2$ ），压力方程便多出一个时间项。

\frac{\partial (\psi\, p)}{\partial t} + \nabla\cdot(\rho\,\mathbf{u}_f) = 0

这个方程同时兼具对流与扩散的性格。当 $M\to 0$ ， $\psi$ 增大使时间项不再主导，压力按椭圆型（瞬时的全局耦合）求解。当 $M$ 较大，则转为双曲型（声波以有限速度传播）。密度基求解器在低马赫下因密度-压力耦合变弱而变得刚性，而压力基方法通过 EOS 把该耦合显式保留。于是一份代码就能跨越亚音速到超音速。

用 Python 复活并抹去棋盘格#

两种模板放到一维周期网格上跑。要看的是棋盘格留下还是消失。

import numpy as np
 
def poisson_sweep(p, f, naive):
    """一维周期网格上 -p'' = f 的一次无松弛 Gauss-Seidel 扫掠。"""
    N = len(p)
    for i in range(N):
        if naive:                       # naive 线性插值: 隔格取值的解耦模板
            p[i] = 0.5 * (p[(i - 2) % N] + p[(i + 2) % N] + f[i])
        else:                           # Rhie-Chow: 联系相邻格的紧凑三点模板
            p[i] = 0.5 * (p[(i - 1) % N] + p[(i + 1) % N] + f[i])
    p -= p.mean()                       # 压力仅相差一个常数 -> 固定均值
    return p
 
def checkerboard_metric(p):
    """(+,-,+,-,...) 分量的大小。为 0 表示无锯齿。"""
    signs = (-1.0) ** np.arange(len(p))
    return abs(np.dot(p, signs)) / len(p)
 
N = 48
x = 2 * np.pi * np.arange(N) / N
f = np.sin(x) + 0.4 * np.sin(2 * x)
f -= f.mean()
 
for naive in (True, False):
    p = 0.8 * (-1.0) ** np.arange(N)    # 以棋盘格初始化
    for _ in range(4000):
        p = poisson_sweep(p, f, naive)
    tag = "naive linear" if naive else "Rhie-Chow  "
    print(f"{tag}  checkerboard = {checkerboard_metric(p):.2e}")
 
# naive linear   checkerboard = 8.00e-01   <- 锯齿残留
# Rhie-Chow      checkerboard = 3.1e-16    <- 消至机器精度

同样的源项、同样的初值、同样的迭代次数。只改了模板。naive 消不掉锯齿，Rhie-Chow 则抹掉。和前面的画面是同一个位置。

初始压力是 $0.8\times(-1)^i$ 。指标若是 $p$ 与同一符号列的内积绝对值除以网格数，一开始就是 0.8。48 个网格、4000 次扫描之后，naive 仍是 $8.00\times 10^{-1}$ 。一位都没削掉。Rhie-Chow 是 $3.1\times 10^{-16}$ 。机器精度。

naive 守住 0.8，是因为模板。锯齿上， $i$ 与 $i\pm 2$ 同号。 $p_i=0.5(p_{i-2}+p_{i+2}+f_i)$ 用平均消不掉这个模态。网格中心的梯度 $(p_E-p_W)/2\Delta x$ 也一样。 $E$ 和 $W$ 都与 $P$ 反号，所以 $p_E-p_W=0$ 。幅度 0.8 的锯齿，在中心处梯度是 0。残差里看不见。面上的 $p_E-p_P$ 反号，跳跃大小是 1.6。这就是 Rhie-Chow 抓住的地方。源项 $f=\sin x+0.4\sin 2x$ 的光滑部分，两种模板都解掉了。剩下的差别只有锯齿。

把 α_p 设成 1，校正先炸掉#

像 $p\leftarrow p^*$ 那样一步把压力整格换掉，耦合就会过校正。式子只放行 $\alpha_p$ 那一份。经验上 $\alpha_p\approx 0.3$ 、 $\alpha_u\approx 0.7$ ，加起来是 1。实验曲线里，0.3 在爬，1.5 附近最快，贴近 2 时每次扫描残差又跳回去。SIMPLE 的 $\alpha_p$ 放在 0 和 1 之间，坐在那条曲线偏慢的一侧。

松弛一改、收敛解也跟着变的地方，是记住欠松弛的面速度。同位网格上若用网格速度的平均来做面速度，相邻压力差就掉了，锯齿留在残差外面。用 $\psi=\partial\rho/\partial p$ 把时间项救回来，是这之后的事。锯齿还活着的压力方程上再加可压缩，只是声速变快，棋盘格还在。