緩和係数を一つ変えたら収束解が動いた — Rhie–Chow面速度がアンダーリラクゼーションを覚えている場所

標準のRhie–Chowは収束後も圧力減衰項の前にα_uを残します。前回の面速度をセル平均ではなく面に保存した値で使えば、収束解はαに依存しなくなります。

緩和係数だけ変えたのに答えが変わった

定常SIMPLE計算が残差 $10^{-8}$ まで下がりました。運動量の緩和係数 $\alpha_u$ を0.7から0.3に変えて再計算します。反復回数は3倍ほどに増えましたが、やはり収束しました。ところが二つの速度場が一致しません。格子も境界条件も離散化も同じです。

緩和係数は反復の経路だけを変えるべきものです。収束点がそのつまみで動くなら、結果は信用できません。本記事では、原因がRhie–Chow面速度の式の一行にあることを示します。1次元の静水圧カラムで差を数値として測り、Majumdar(1988)の補正がなぜそれを0にするのかを確認します。

前提知識としてはSIMPLEとRhie–Chowの基本構造で十分です。コロケート格子(速度と圧力を同じセル中心に置く配置)を前提とします。

自分で動かす二層カラム

α_uスライダーとstandard/Majumdarの切り替えが、同じカラムでセル速度の棒を黄色の目盛りからどれだけずらすかを変えます。

alpha_u = 0.30 density ratio = 10 iters/frame = 2

iterating | iter = 0 | residual = 1.0e+0 | peak |u_P| = 0.0000 | cells |u|>0.01 = 0 | max |p-p(a=1)| = 0.00e+0

standardで $\alpha_u$ を0.1まで下げると、棒の符号がセルごとに反転して広がります。Majumdarに切り替えると、どの $\alpha_u$ でも棒は目盛りの上に戻ります。

SIMPLEの1反復でαが入る場所#

セルPの離散運動量方程式は次の形です。

a_P u_P = \sum_{nb} a_{nb} u_{nb} + b_P - V_P (\nabla p)_P

$a_P$ は対角係数、 $a_{nb}$ は隣接係数、 $b_P$ は圧力以外のソース、 $V_P$ はセル体積です。

SIMPLEはこの式をそのまま解きません。対角を $a_P/\alpha_u$ に大きくし、その分を前回の反復値 $u_P^{old}$ で右辺に足します。

\frac{a_P}{\alpha_u} u_P = \sum_{nb} a_{nb} u_{nb} + b_P - V_P (\nabla p)_P + \frac{1-\alpha_u}{\alpha_u} a_P u_P^{old}

収束すれば $u_P = u_P^{old}$ なので、追加した二つの項は打ち消し合います。セル中心の運動量方程式だけを見れば、緩和は収束点に痕跡を残しません。問題は面で起こります。

面速度の式に $u^{old}$ が混ざる経路#

緩和した式をセル速度について整理すると次のようになります。

u_P = \hat{u}_P - d_P (\nabla p)_P,\qquad d_P = \frac{\alpha_u V_P}{a_P}

\hat{u}_P = \frac{\alpha_u}{a_P}\Big(\sum_{nb} a_{nb} u_{nb} + b_P\Big) + (1-\alpha_u)\,u_P^{old}

$\hat{u}_P$ は圧力勾配を除いた残りの速度、 $d_P$ は圧力勾配に対する速度の感度です。

Rhie–Chow補間は、面fで $\hat{u}$ と $d$ を補間し、圧力勾配だけを面の両側のセル値から計算し直します。

u_f = \overline{\hat{u}}_f - d_f (\nabla p)_f

$\overline{\,\cdot\,}_f$ は2セル間の線形補間、 $(\nabla p)_f = (p_N - p_P)/\Delta x$ はコンパクトな面勾配です。

問題は $\overline{\hat{u}}_f$ の中に $(1-\alpha_u)\,\overline{u^{old}}_f$ が入っていることです。面速度の「前回値」をセル速度の平均で代用しているのです。セル速度の平均は面速度と同じではありません。その差こそがRhie–Chow減衰項です。

収束するとαが減衰項の前に残る

収束状態で $u^{old} = u$ を代入して整理すると、標準の式は次のようになります。

u_f = \overline{u}_f + \alpha_u \,\overline{\Big(\tfrac{V}{a}\Big)}_f \Big[\,\overline{(\nabla p)}_f - (\nabla p)_f\,\Big]

角括弧はセル勾配の平均とコンパクトな面勾配の差です。滑らかな圧力では $O(\Delta x^2)$ で、チェッカーボード(格子間隔ごとに符号が反転する圧力モード)で最大になります。この項がチェッカーボードを連続の式に見えるようにする仕組みです。

この仕組みの強さが $\alpha_u$ 倍に弱まります。 $\alpha_u = 0.1$ なら減衰は10%しか残りません。収束解は $\alpha_u$ に依存し、 $\alpha_u$ を下げるほどsimpleFoamで圧力残差が落ちない格子模様がそっと戻ってきます。

Majumdarの処方は簡単です。面速度の前回値を面に別途保存し、そのまま使います。

u_f = \overline{\hat{u}^{*}}_f - d_f (\nabla p)_f + (1-\alpha_u)\, u_f^{old}

$\hat{u}^{*}$ は $\hat{u}$ から $(1-\alpha_u)u^{old}$ の項を除いた値です。収束すれば $u_f = u_f^{old}$ なので $(1-\alpha_u)$ が両辺で整理され、減衰項の係数は1に戻ります。

コード:同じカラム、5つのα#

テスト問題はあえて答えが自明なものを選びました。両端が壁の1次元カラムに、下半分は $\rho = 10$ 、上半分は $\rho = 1$ を満たし、重力をかけます。正解は静止状態 $u = 0$ です。1次元では連続の式がすべての面で $u_f = 0$ を強制するため、圧力補正は1回の積分で終わります。残るセル速度 $u_P$ がそのまま補間誤差です。

# 1次元静水圧カラム: 下半分 rho=10、上半分 rho=1、両端は壁。正解は u = 0。
def density_column(n, rho_lo=1.0, rho_hi=10.0):
    return [rho_hi if (i + 0.5) / n < 0.5 else rho_lo for i in range(n)]
 
def cell_grad(p, h):
    # 面圧力 = 隣接平均、壁面は線形外挿
    n = len(p)
    pf = [1.5 * p[0] - 0.5 * p[1]]
    pf += [0.5 * (p[i] + p[i + 1]) for i in range(n - 1)]
    pf += [1.5 * p[-1] - 0.5 * p[-2]]
    return [(pf[i + 1] - pf[i]) / h for i in range(n)]
 
def simple_hydrostatic(n, alpha, majumdar, a=1.0, g=1.0, tol=1e-13):
    h = 1.0 / n
    f = [-r * g for r in density_column(n)]
    p, u, uf = [0.0] * n, [0.0] * n, [0.0] * (n + 1)
    d = alpha / a                      # d = alpha*V/a_P (V = 1 単位)
    for it in range(1, 5000):
        gp = cell_grad(p, h)
        # 緩和した運動量予測子
        us = [alpha * (f[i] - gp[i]) / a + (1 - alpha) * u[i] for i in range(n)]
        # Rhie-Chow: 圧力項を除いた u_hat を面へ補間
        uh = [us[i] + d * gp[i] - majumdar * (1 - alpha) * u[i] for i in range(n)]
        ufs = [0.0] * (n + 1)
        for i in range(1, n):
            ufs[i] = 0.5 * (uh[i - 1] + uh[i]) - d * (p[i] - p[i - 1]) / h
            ufs[i] += majumdar * (1 - alpha) * uf[i]
        # 圧力補正: 1次元では連続の式が面速度0を直接要求する
        pc = [0.0] * n
        for i in range(1, n):
            pc[i] = pc[i - 1] + h * ufs[i] / d
        p = [p[i] + pc[i] for i in range(n)]
        gc = cell_grad(pc, h)
        un = [us[i] - d * gc[i] for i in range(n)]
        uf = [0.0] * (n + 1)          # 補正後の面速度は厳密に0
        res = max(abs(un[i] - u[i]) for i in range(n))
        u = un
        if res < tol and it > 5:
            return u, p, it
    return u, p, it
 
n = 80
_, p_ref, _ = simple_hydrostatic(n, 1.0, False)
print("alpha | iters | peak|u|  cells>0.01 | max|p-p(1)|   <- standard / majumdar")
for alpha in (1.0, 0.7, 0.5, 0.3, 0.1):
    row = []
    for mj in (0, 1):
        u, p, it = simple_hydrostatic(n, alpha, mj)
        wide = sum(abs(v) > 0.01 for v in u)
        dp = max(abs(p[i] - p_ref[i]) for i in range(n))
        row.append(f"{it:4d} {max(map(abs, u)):.4f} {wide:3d} {dp:.2e}")
    print(f" {alpha:.1f}  | " + "  /  ".join(row))
 
u, _, _ = simple_hydrostatic(n, 0.1, 0)
print("standard, alpha=0.1, cells 34..45:", " ".join(f"{v:+.3f}" for v in u[34:46]))

alpha | iters | peak|u|  cells>0.01 | max|p-p(1)|   <- standard / majumdar
 1.0  |    6 2.2500   2 0.00e+00  /     6 2.2500   2 0.00e+00
 0.7  |   26 2.0499   6 4.59e-03  /    27 2.2500   2 8.88e-16
 0.5  |   42 1.8640   6 8.24e-03  /    45 2.2500   2 0.00e+00
 0.3  |   79 1.5925  10 1.27e-02  /    84 2.2500   2 8.88e-16
 0.1  |  247 1.0811  16 1.92e-02  /   271 2.2500   2 0.00e+00
standard, alpha=0.1, cells 34..45: +0.041 -0.079 +0.152 -0.292 +0.562 -1.081 +1.081 -0.562 +0.292 -0.152 +0.079 -0.041

1.0811と16セルが語ること#

$\alpha_u = 1$ では二つの式は一致します。界面の両側の2セルに $\mp 2.25$ が残ります。これは緩和とは無関係な別の誤差です。体積力 $\rho g$ の段差をRhie–Chowの面が見ていないために生じます。密度比 $r$ での大きさは $(r-1)g/4$ で、ウィジェットの密度比スライダーで確認できます。Denner(2013)のようなbalanced-force処理は、この項を面で圧力と同じ方法で扱って消します。

標準の式では、 $\alpha_u$ を下げるほど最大速度が2.25から1.0811に減ります。数字だけ見ると良くなったようです。しかし $|u| > 0.01$ のセルは2個から16個に増えました。最終行を見ると符号がセルごとに反転し、振幅は1セルあたり約1.9倍ずつ小さくなります。減衰が弱まった分、誤差が局所にとどまれず、振動モードとして広がったのです。

圧力も動きました。 $\alpha_u = 0.1$ では $\alpha_u = 1$ の解との差が $1.92\times10^{-2}$ です。Majumdar側は5ケースすべてで $10^{-15}$ 以下で、機械精度で同じ解です。反復回数はMajumdarのほうがわずかに多いですが、収束経路の違いにすぎません。

反復回数を式ごとに並べると、緩和が経路だけを伸ばすことと、解そのものを動かすことが分かれます。 $\alpha_u = 1$ ではどちらも6回です。0.7なら26と27、0.5なら42と45、0.3なら79と84、0.1なら247と271です。0.7から0.3へ下げると標準側は $79/26 \approx 3.0$ 倍です。冒頭の「3倍」はこの欄です。Majumdarは一拍遅いだけで、比率はほぼ同じです。

解が割れる列は peak $|u|$ です。標準は 2.2500、2.0499、1.8640、1.5925、1.0811 と下がります。 $\alpha_u = 0.1$ では減衰項の前の係数が 0.1 なので、チェッカーボードを押さえる力は 10% しか残りません。しかし最大速度は 2.25 の 10% である 0.225 ではなく 1.0811 です。誤差が減ったのではありません。界面の $\mp 2.25$ が隣へにじんだのです。 $|u| > 0.01$ のセルが 2、6、6、10、16 と増えたのがそのにじみです。

34番から符号が反転する列の振幅比はほぼ一定です。 $1.081/0.562 \approx 1.92$ 、 $0.562/0.292 \approx 1.92$ 、 $0.292/0.152 \approx 1.92$ です。1セルあたり約 1.9 倍。圧力差も $4.59\times10^{-3}$ 、 $8.24\times10^{-3}$ 、 $1.27\times10^{-2}$ 、 $1.92\times10^{-2}$ と、 $\alpha_u$ が小さいほど大きくなります。Majumdar側は五つの $\alpha$ すべてで peak $|u| = 2.2500$ 、広いセルは 2 個、圧力差は 0 か $8.88\times10^{-16}$ です。2.25 は $(r-1)g/4$ に $r = 10$ 、 $g = 1$ を入れた値です。 $\alpha_u$ は反復回数にしか残りません。

Majumdar補正が保存するもの#

実装コストは面ごとに実数1つです。圧力補正が終わった後、補正された面フラックスを $u_f^{old}$ として保存すればよいのです。ほとんどのコードはすでに質量フラックスを面に保持しているので、実質的にタダです。

注意点が一つあります。保存するのは補間直後の $u_f^{*}$ ではなく、圧力補正まで終えた面速度です。補正前の値を保存すると、連続の式を満たさないフラックスが次の反復の記憶になります。

時間項も同じ落とし穴を掘る

非定常計算では時間項が同じ役割を果たします。陰的時間離散化の $\rho V u_P^n/\Delta t$ が $\hat{u}_P$ に混ざってセル平均で補間されると、到達した定常解が $\Delta t$ に依存します。Choi(1999)らがこの依存性を指摘しました。OpenFOAMの ddtCorr は、前の時間ステップの面フラックスを直接使ってこの差を補正する項です。

原理は緩和と同じです。「前回値」が必要なあらゆる場所で、面には面の値を使います。複数の時間レベルを使うスキームでは、保存する面フラックスもその分増えます。密度の段差が残す $\mp 2.25$ のように、静止すべき解が揺れる場所は極座標で well-balanced が角運動量を落とすときと同じ種類の釣り合い誤差です。

緩和係数を収束速度のつまみとしてだけ残すには

検査は簡単です。同じケースを2つの $\alpha_u$ で収束させ、結果を差し引いてみます。差が反復許容誤差のレベルなら、面速度の式はきれいです。差が離散化誤差と同程度なら、面補間が $u^{old}$ をセル平均で混ぜている可能性が高いです。

非定常解析なら同じ検査を $\Delta t$ で繰り返します。定常に収束した解が $\Delta t$ で変わるなら、時間項側の補正が抜けています。どちらの検査も、疑うべきコード1行を指し示すには十分です。