圧力が市松模様に振動するとき — SIMPLE と Rhie-Chow の圧力-速度連成

α_p を 1 にして圧力を毎回まるごと置き換えると、補正ループは過補正で破綻します。

速度場は問題なさそうに見えました。圧力場を開くと、セルごとに値が上下していました。市松模様（チェッカーボード）です。残差は下がります。圧力だけが鋸歯です。バグではありません。圧力と速度を同じ格子点に置いた時点で、予告されていた結果です。

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$ がありません。セルは自分の圧力を見ず、二つ隣だけを見ます。こうして奇数セルと偶数セルが分離します（odd-even decoupling）。 $(+,-,+,-,\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 はここに修正段階を 2–3 回追加した予測子-修正子で、非定常計算に使います。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 で市松模様を蘇らせ、消す#

主張だけでは不十分です。二つのステンシルを 1 次元周期格子で走らせ、市松模様が残るか消えるかを測ります。

import numpy as np
 
def poisson_sweep(p, f, naive):
    """1D 周期格子での -p'' = f の緩和なし Gauss-Seidel の 1 スイープ。"""
    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: 隣を結ぶ密な 3 点ステンシル
            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$ です。指標が、同じ符号列との内積の絶対値をセル数で割ったものなら、最初は 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 の鋸歯は、中心で勾配がゼロです。残差に出ません。面の $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$ で時間項を生かすのはそのあとです。鋸歯が生きた圧力方程式に圧縮性を載せると、音速だけが速くなり、市松模様はそのままです。