圧力を方程式から消した男 — Chorinの射影法とフラクショナルステップ

流入と流出がほんの少しずれると、NeumannのPoissonは解を持たない。射影を切ると、その発散がステップごとに残る。

1967年、Alexandre Chorinは圧力を方程式から一時的に消し去りました。非圧縮Navier–Stokes方程式を解くとき、もっとも扱いにくい項が圧力だからです。圧力には時間微分がありません。密度が一定なので、状態方程式から圧力を求めることもできません。Chorinの答えは大胆でした。まず圧力を無視して速度を進め、その結果を強制的に発散ゼロ(divergence-free)の状態へ戻すのです。この記事では、そのフラクショナルステップ — 一般に射影法(projection method)と呼ばれる手法 — がどう働くのかを、Helmholtz–Hodge分解という一行の定理から出発し、実際に動く2Dソルバーまで作って確かめます。読み終える頃には、非圧縮コードの実行時間の大半がなぜ圧力Poisson方程式に費やされるのかが分かります。

左パネルのソース強度を0から14へ上げると発散の色が生き、右の射影は渦だけを残します。

source strength6.0max|∇·u*| = 0.000 → 0.0000

Red = positive divergence (source), blue = negative (sink), dark = zero. Raise the source and the left panel lights up; the right panel stays dark — projection strips the compressible part and keeps only the swirl.

ソースが0なら両側はほぼ同じです。14まで上げると左だけが赤と青に染まり、右は暗いまま残ります。射影が圧縮成分を丸ごと外した結果です。

圧力には時間微分がない

非圧縮Navier–Stokes方程式を書き下すと、問題がそのまま見えてきます。

\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u}, \qquad \nabla\cdot\mathbf{u} = 0

ここで $\mathbf{u}$ は速度、 $p$ は圧力、 $\nu$ は動粘性係数(運動量の拡散率)です。運動量方程式は $\mathbf{u}$ の時間発展を教えてくれます。ところが第二式 $\nabla\cdot\mathbf{u}=0$ は発展方程式ではありません。各瞬間に守るべき拘束条件です。

圧縮性流れでは連続の式が密度を発展させ、密度が状態方程式を通じて圧力を決めます。非圧縮の極限ではその連鎖が切れます。圧力は時間とともに「発展」しません。圧力は、速度が発散しないように毎瞬間みずからを合わせるラグランジュ乗数(拘束を課す未知数)にすぎません。ですから圧力を時間前進させようとする試みそのものが徒労なのです。

すべての速度場は二つに分かれる — Helmholtz–Hodge#

突破口はベクトル解析の古い定理にあります。適当な領域では、任意のベクトル場 $\mathbf{w}$ は発散ゼロの部分と勾配の部分の和に一意に分解されます。

\mathbf{w} = \mathbf{u} + \nabla\phi, \qquad \nabla\cdot\mathbf{u} = 0

$\mathbf{u}$ は渦(ソレノイダル)成分、 $\nabla\phi$ は圧縮(勾配)成分です。この分解から $\mathbf{u}$ だけを取り出す演算を射影演算子 $\mathbb{P}$ と呼びます。手順は単純です。上式の発散をとると $\nabla\cdot\mathbf{u}=0$ なので、

\nabla^2\phi = \nabla\cdot\mathbf{w}

このPoisson方程式で $\phi$ を解き、元の場から $\nabla\phi$ を引けば、発散ゼロの部分だけが残ります。

\mathbf{u} = \mathbb{P}\mathbf{w} = \mathbf{w} - \nabla\phi

予測し、射影する — フラクショナルステップ

Chorinのアルゴリズムは、この分解をそのまま時間積分に乗せます。1ステップを二つに分けます。

予測子(predictor) — 圧力を落としたまま速度を進め、中間速度 $\mathbf{u}^*$ を得ます。

\frac{\mathbf{u}^* - \mathbf{u}^n}{\Delta t} = -(\mathbf{u}^n\cdot\nabla)\mathbf{u}^n + \nu\nabla^2\mathbf{u}^n

$\mathbf{u}^*$ は発散ゼロではありません。圧力を無視したので当然です。

修正子(corrector) — $\mathbf{u}^*$ を発散ゼロへ戻します。圧力Poissonを解き、

\nabla^2\phi = \frac{1}{\Delta t}\nabla\cdot\mathbf{u}^*

勾配を引いて次ステップの速度をつくります。

\mathbf{u}^{n+1} = \mathbf{u}^* - \Delta t\,\nabla\phi

ここで $\phi$ が事実上の圧力の役割を果たします( $p^{n+1}\approx\rho\phi$ )。圧力は発展させる対象ではなく、毎ステップ拘束条件を満たすために解き直す値です。この二行が射影法のすべてです。

圧力Poissonと境界条件の落とし穴#

ここから実務の苦しみが始まります。予測子は安価です。しかし修正子は毎ステップPoisson方程式を一つ解きます。大きな格子ではこの楕円型ソルブが全体コストの70〜90%を占めます。GMRESやマルチグリッドが非圧縮コードの心臓である理由です。

境界条件も落とし穴だらけです。 $\phi$ の境界条件は壁でNeumann( $\partial\phi/\partial n = 0$ )です。Neumann問題は定数分の不定性があり、解が存在するには適合条件 — 領域全体で $\int \nabla\cdot\mathbf{u}^*\,dV = 0$ — を満たす必要があります。流入・流出の流量が合わなければ、Poisson方程式はそもそも解けません。NaNが出たら、まず境界の流量を疑ってください。

collocated(非スタッガード)格子では、圧力と速度を同じ点に置くとチェッカーボード振動が生じます。スタッガード格子を使うか、Rhie–Chow補間で抑える必要があります。その面速度が緩和係数を覚えていると、収束解そのものが変わります。ちなみにSIMPLE系は同じ射影の考えを毎ステップ何度も繰り返して定常状態へ追い込むのに対し、ここでのフラクショナルステップはステップあたり一度の射影で非定常流れを時間精度よく追いかけます。

Python:二重せん断層をFFTで巻き上げる#

周期境界では、圧力PoissonをFFTで一度に解きます。Fourier空間ではラプラシアンが $-|\mathbf{k}|^2$ の掛け算に変わるからです。二重せん断層(二つの噛み合ったジェット)を初期条件に与えると、Kelvin–Helmholtz渦へと巻き上がります。射影法の古典的な検証問題です。

import numpy as np
 
N = 128                # 格子(一辺)
h = 1.0 / N            # 格子間隔
nu = 5e-3              # 動粘性 -> Re = U L / nu ~ 200
dt = 2e-3              # 拡散・移流の安定範囲内
steps = 3000
 
x = (np.arange(N) + 0.5) * h
X, Y = np.meshgrid(x, x, indexing='ij')
 
# 二重せん断層 + 小さな摂動
rho, delta = 1.0 / 30.0, 0.05
u = np.where(Y <= 0.5, np.tanh((Y - 0.25) / rho), np.tanh((0.75 - Y) / rho))
v = delta * np.sin(2 * np.pi * X)
 
# FFT Poisson用の波数
k = 2 * np.pi * np.fft.fftfreq(N, d=h)
KX, KY = np.meshgrid(k, k, indexing='ij')
K2 = KX**2 + KY**2
K2[0, 0] = 1.0         # 平均モードのゼロ割りを回避
 
def divergence(a, b):
    dadx = (np.roll(a, -1, 0) - np.roll(a, 1, 0)) / (2 * h)
    dbdy = (np.roll(b, -1, 1) - np.roll(b, 1, 1)) / (2 * h)
    return dadx + dbdy
 
def projection_correct(a, b):
    # Laplacian(phi) = div  ->  a <- a - grad(phi)
    phi_hat = np.fft.fft2(divergence(a, b)) / (-K2)
    phi_hat[0, 0] = 0.0
    phi = np.real(np.fft.ifft2(phi_hat))
    dpx = (np.roll(phi, -1, 0) - np.roll(phi, 1, 0)) / (2 * h)
    dpy = (np.roll(phi, -1, 1) - np.roll(phi, 1, 1)) / (2 * h)
    return a - dpx, b - dpy
 
def predictor(a, b):
    # 移流(中心差分) + 拡散、明示的Euler 1ステップ
    ax = (np.roll(a, -1, 0) - np.roll(a, 1, 0)) / (2 * h)
    ay = (np.roll(a, -1, 1) - np.roll(a, 1, 1)) / (2 * h)
    bx = (np.roll(b, -1, 0) - np.roll(b, 1, 0)) / (2 * h)
    by = (np.roll(b, -1, 1) - np.roll(b, 1, 1)) / (2 * h)
    lap = lambda f: (np.roll(f, -1, 0) + np.roll(f, 1, 0)
                     + np.roll(f, -1, 1) + np.roll(f, 1, 1) - 4 * f) / h**2
    astar = a + dt * (-(a * ax + b * ay) + nu * lap(a))
    bstar = b + dt * (-(a * bx + b * by) + nu * lap(b))
    return astar, bstar
 
for n in range(steps):
    us, vs = predictor(u, v)                 # 予測:中間速度 u*
    d_before = np.abs(divergence(us, vs)).max()
    u, v = projection_correct(us, vs)        # 射影:発散を除去
    if n % 500 == 0:
        d_after = np.abs(divergence(u, v)).max()
        print(f"step {n:4d}  |div u*|={d_before:.2e} -> |div u|={d_after:.2e}")

出力を見ると、各ステップで $|\nabla\cdot\mathbf{u}^*|$ は $10^{-2}$ の規模ですが、射影後の $|\nabla\cdot\mathbf{u}|$ は $10^{-13}$ まで落ちます。FFT射影が発散を機械精度で消したのです。

一辺 128 なら間隔は 1/128 です。動粘性 5e-3、速度の規模 1、領域の長さ 1 なら、Re は 1 を 5e-3 で割った 200 です。コメントの Re 約 200 と一致します。時間刻み 2e-3 を 3000 回積むと、積分時間は 6 です。せん断層の厚さ 1/30 は約 0.033、間隔 1/128 は約 0.0078 なので、四セルあまりです。擾乱振幅 0.05 はジェット速度 1 の 5% なので、巻き上がりはこの小さな正弦波から始まります。予測直後の発散 $10^{-2}$ と射影後の $10^{-13}$ の比は十一桁です。FFT が平均モードを 0 に固定する一行は、Neumann 問題の定数の自由度を閉じます。流入と流出が合わず適合条件が破れると、その分母だけ直しても解は生まれません。費用の 70% から 90% という話も、同じ行列を指します。一歩のうち予測子が持つ取り分は、残りの 10% から 30% だけです。

時間刻み 2e-3 での補正は、中間速度からその刻みだけ $\phi$ の勾配を引く作業です。発散 $10^{-2}$ を消すには、 $\phi$ のラプラシアンは発散を時間刻みで割った値、およそ 5 になります。圧力勾配が一歩で速度をそれだけ戻します。射影を切ると、この補正は速度に入らず発散として残ります。3000 ステップは、その残りかすが格子を満たすのに十分です。

射影を切ると何が起きるか

修正子がないと、速度場は毎ステップ少しずつ発散を溜め込みます。質量が保存されず、渦はぼやけ、やがて格子を埋める雑音になります。下のシミュレーションで確かめてみましょう。

dt0.60max|∇·u| = 0.000

Double shear layer rolling up. Red/blue = vorticity sign. Turn projection OFF and max|∇·u| climbs while the vortices dissolve into noise.

ProjectionをONにすると、せん断層はきれいな渦のペアへ巻き上がります。OFFに切り替えた瞬間、右上の max|∇·u| が跳ね上がり、渦の模様が崩れます。dtを上げると移流がより激しくなり、崩壊が速く訪れます。この一度のトグルが、「なぜ射影が必要か」への最短の答えです。

流量がひとしずくずれるとPoissonは答を持たない#

非圧縮では、圧力は発展させる変数ではなく、毎ステップ $\nabla\cdot\mathbf{u}=0$ を課すラグランジュ乗数です。射影は圧力を落として $\mathbf{u}^*$ を作り、Poissonで $\phi$ を解いて勾配を引きます。その順序自体は短いです。壊れる場所は境界です。流入と流出が合わないと適合条件がなく、行列は定数の自由度だけを残して解を出しません。

一歩の射影で時間精度を守るなら、対流をその一歩の中で解き直しません。粘性だけ次数を上げ、対流は外へ押し出す選択は、Gearが対流を外挿して誤差を大きくした場所で分かれます。予測子の対流をより正確に見せようとすると、射影が合わせた発散ゼロが揺れます。