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

チェッカーボード圧力の原因、Rhie-Chow 補間、SIMPLE 圧力方程式、緩和係数

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

この記事では、その鋸歯がなぜ生じるのか、Rhie-Chow 補間がどう取り除くのか、そして連続の式がどのように圧力方程式へ変わるのかを見ていきます。最後には同じソルバを圧縮性まで伸ばします。

連続の式は隣のセルの圧力を見ていない

コロケート格子（速度と圧力を同じセル中心に保存）での半離散運動量方程式は次のように書けます。

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$ 、つまり隣り合うセルの圧力差に直接結びつきます。鋸歯モードはもう隠れる場所がありません。

下のシミュレーションを直接操作してみましょう。二つの方式で圧力を緩和し、市松模様で初期化した圧力がどうなるかを見ます。

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

Cells N: 32

naive 補間では鋸歯の振幅がそのまま残ります。Rhie-Chow に切り替えると、同じ市松模様が滑らかな曲線へ沈みます。Reseed checkerboard を押してもう一度確かめてください。

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 は完全に消します。ビューアで見たものとまったく同じです。

圧力ソルバの前で失敗しないために

コロケート格子では、面速度をセル速度の単純平均で作ってはいけません。Rhie-Chow（またはスタガード格子）で隣接圧力差を直接効かせれば、鋸歯は生じません。
連続の式には圧力を解く方程式がありません。SIMPLE はそれを圧力 Poisson 方程式に置き換えた予測子-修正子ループです。
発散したら、まず $\alpha_p$ を下げてください。0.3 付近から始め、速度との和が約 1 になるよう $\alpha_u$ を合わせます。
低マッハと超音速を一つのコードで扱いたいなら、EOS の $\psi=\partial\rho/\partial p$ で密度を圧力に結び、圧力方程式の時間項を生かします。