こんにちはJS2IIUです。Pythonで攻略する流体力学シミュレーション連載の第13回です。
前回までは、1次元のバーガース方程式などを題材に、非線形移流や衝撃波の扱いといった「流体計算の難所」を学んできました。いよいよ今回から、計算対象を2次元へと拡張します。
2次元以上の流体計算において、私たちが最初に直面する最大の壁は「圧力 \(P\)」の扱いです。非圧縮性ナビエ・ストークス方程式には、圧力に関する時間発展の式が存在しません。圧力は「流体の体積が変わらないように（連続の式を満たすように）」一瞬で空間全体に伝播し、速度場を拘束する役割を持っています。この速度と圧力のカップリングを解くのは非常に手間がかかります。
そこで今回は、2次元限定ではありますが、「圧力を計算から消去して、計算を劇的にシンプルにする」という魔法のような手法、流れ関数-渦度法（Stream Function-Vorticity Method）を解説します。今回もよろしくお願いします。
1. 「圧力」を消し去るアイデア
2次元の非圧縮性流体の運動を記述する際、通常は速度成分 \(u, v\) と圧力 \(P\) の3つの変数を用います。しかし、非圧縮性ナビエ・ストークス方程式では圧力の時間発展式が存在せず、速度場と圧力場が連続の式によって強く結びついているため、数値的に圧力を扱うことが大きな課題となります。
この問題を解決するために、2次元流れに特有の「流れ関数-渦度法（Stream Function-Vorticity Method）」が考案されました。この手法では、従来の速度・圧力の代わりに、渦度（Vorticity） \(\omega\) と 流れ関数（Stream Function） \(\psi\) という2つのスカラー場を導入します。これにより、圧力を明示的に計算する必要がなくなり、計算手順が大幅に簡素化されます。
渦度 \(\omega\)（オメガ）の定義と役割
渦度は、流体の回転の強さや局所的な回転運動を定量化する物理量です。2次元直交座標系では、次式で定義されます：
$$
\omega = \frac{\partial v}{\partial x} – \frac{\partial u}{\partial y}
$$
この式は、ある点に小さな水車を置いたときの回転の激しさを表現しています。ナビエ・ストークス方程式に回転（\(\text{rot}\)）を作用させることで、圧力項が消去され、渦度の時間発展を直接記述する渦度輸送方程式が得られます。これにより、圧力を介さずに流体の運動を追跡できるようになります。
流れ関数 \(\psi\)（プサイ）の定義と利点
流れ関数は、2次元非圧縮性流れにおいて速度場を一意に記述できるスカラー場です。次の関係式で速度成分と結びつきます：
$$
u = \frac{\partial \psi}{\partial y}, \quad v = -\frac{\partial \psi}{\partial x}
$$
この定義の最大の利点は、流れ関数を用いることで連続の式（\(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0\)）が自動的に満たされる点です。すなわち、流れ関数を導入するだけで非圧縮性条件が保証され、数値計算の安定性と簡便性が大きく向上します。
2つの方程式への集約
このようにして導入された渦度 \(\omega\) と流れ関数 \(\psi\) を用いると、2次元非圧縮性流体の運動は以下の2つの方程式に集約されます：
- 渦度輸送方程式（放物線型/双曲型）：渦度が流れに沿って移流され、粘性拡散によって広がる現象を記述します。
- ポアソン方程式（楕円型）：渦度分布から流れ関数を決定し、そこから速度場を再構成します。
このアプローチにより、圧力の直接計算を回避しつつ、2次元流体の本質的なダイナミクスを効率的にシミュレーションできるのです。
2. アルゴリズムの全貌
流れ関数-渦度法による2次元非圧縮性流体の数値シミュレーションは、各タイムステップごとに次の4つの主要な手順を順番に実行することで進行します。これにより、圧力を明示的に計算することなく、流体の運動を安定かつ効率的に追跡できます。
- 渦度輸送方程式の数値解法（移流拡散ステップ）
まず、現在の速度場 \(u, v\) を用いて、渦度 \(\omega\) の時間発展方程式（渦度輸送方程式）を数値的に解きます。この方程式は、渦度が流れに沿って運ばれる（移流）と同時に、粘性によって広がる（拡散）現象を記述しています。差分法や有限体積法などを用いて、次の時刻の渦度分布を計算します。 - 流れ関数のポアソン方程式を解く
渦度分布が得られたら、次にポアソン方程式 \(\nabla^2 \psi = -\omega\) を解いて新しい流れ関数 \(\psi\) を求めます。ポアソン方程式は楕円型偏微分方程式であり、数値的には反復法（SOR法やガウス・ザイデル法など）や行列計算によって解かれます。これにより、流れの全体的な形状が決定されます。 - 速度場の再構成
新たに得られた流れ関数 \(\psi\) を空間微分することで、速度成分 \(u, v\) を再計算します。具体的には、\(u = \frac{\partial \psi}{\partial y}\)、\(v = -\frac{\partial \psi}{\partial x}\) の関係式を用いて、格子点ごとに速度場を更新します。これにより、非圧縮性条件を厳密に満たした速度場が得られます。 - 境界条件の適用と更新
最後に、物理的な壁や蓋などの境界に対して「滑りなし条件（no-slip condition）」を適用します。速度場の境界条件を渦度の値に変換し、境界格子点の渦度を適切に設定します。これにより、壁面での物理的な流れの拘束が正しく反映され、数値解の安定性と物理的妥当性が確保されます。
この4つの手順を1ステップとして繰り返すことで、2次元流体の時間発展を高精度かつ効率的にシミュレーションすることが可能となります。
3. Pythonによる実装：キャビティ流れの計算
ここからは、流体力学シミュレーションの代表的なベンチマーク問題であるキャビティ流れ（Lid-driven cavity flow）を題材に、流れ関数-渦度法の具体的なPython実装例を詳しく解説します。
キャビティ流れとは、四角い閉じた箱の上面（蓋）だけを一定速度で動かし、内部の流体がどのように回転・循環するかを調べる問題です。非圧縮性流体の基礎的な挙動や数値解法の検証に広く用いられています。
この実装では、
- 渦度輸送方程式の空間離散化に中心差分、時間積分に前進オイラー法を用いています。
- 流れ関数のポアソン方程式の解法には、第9回でも紹介したSOR法（Successive Over-Relaxation：逐次加速緩和法）を採用しています。
- 境界条件の処理や速度場の再構成も含め、流れ関数-渦度法のアルゴリズム全体をコンパクトにまとめています。
以下に、PyTorchを用いた実装例を示します。
import torch
import matplotlib.pyplot as plt
# 計算パラメータの設定
nx, ny = 51, 51 # 格子点数
lx, ly = 1.0, 1.0 # 領域サイズ
dx = lx / (nx - 1)
dy = ly / (ny - 1)
dt = 0.0005 # タイムステップ
re = 100.0 # レイノルズ数
nu = 1.0 / re # 動粘性係数
# 場の初期化 (PyTorchを使用)
psi = torch.zeros((ny, nx))
omega = torch.zeros((ny, nx))
u = torch.zeros((ny, nx))
v = torch.zeros((ny, nx))
def update_boundary_vorticity(omega, psi, dx, dy):
"""
壁での滑りなし境界条件から、境界の渦度を計算する（Woodの式など）
"""
# 下壁 (y=0)
omega[0, :] = -2.0 * psi[1, :] / dy**2
# 上壁 (y=Ly, 蓋が速度1で動く)
omega[-1, :] = -2.0 * psi[-2, :] / dy**2 - 2.0 / dy
# 左壁 (x=0)
omega[:, 0] = -2.0 * psi[:, 1] / dx**2
# 右壁 (x=Lx)
omega[:, -1] = -2.0 * psi[:, -2] / dx**2
return omega
def solve_poisson(psi, omega, dx, dy, max_iter=500, rel_factor=1.5):
"""
SOR法を用いてポアソン方程式 ∇^2 psi = -omega を解く
"""
for _ in range(max_iter):
psi_old = psi.clone()
# 内部点の更新
psi[1:-1, 1:-1] = (1.0 - rel_factor) * psi[1:-1, 1:-1] + \
(rel_factor / (2.0/dx**2 + 2.0/dy**2)) * (
(psi[1:-1, 2:] + psi[1:-1, :-2]) / dx**2 +
(psi[2:, 1:-1] + psi[:-2, 1:-1]) / dy**2 +
omega[1:-1, 1:-1]
)
# 収束判定
if torch.max(torch.abs(psi - psi_old)) < 1e-4:
break
return psi
# --- メインループ ---
nt = 2000 # 総ステップ数
for n in range(nt):
# 1. 境界の渦度を更新
omega = update_boundary_vorticity(omega, psi, dx, dy)
# 2. 渦度輸送方程式を解く (中心差分+前進オイラー法)
omega_n = omega.clone()
# 移流項 (u*d_omega/dx + v*d_omega/dy)
adv = u[1:-1, 1:-1] * (omega_n[1:-1, 2:] - omega_n[1:-1, :-2]) / (2*dx) + \
v[1:-1, 1:-1] * (omega_n[2:, 1:-1] - omega_n[:-2, 1:-1]) / (2*dy)
# 拡散項 (nu * ∇^2 omega)
diff = nu * ((omega_n[1:-1, 2:] - 2*omega_n[1:-1, 1:-1] + omega_n[1:-1, :-2]) / dx**2 +
(omega_n[2:, 1:-1] - 2*omega_n[1:-1, 1:-1] + omega_n[:-2, 1:-1]) / dy**2)
omega[1:-1, 1:-1] = omega_n[1:-1, 1:-1] + dt * (-adv + diff)
# 3. ポアソン方程式を解いて流れ関数を更新
psi = solve_poisson(psi, omega, dx, dy)
# 4. 速度場を更新 (u=d_psi/dy, v=-d_psi/dx)
u[1:-1, 1:-1] = (psi[2:, 1:-1] - psi[:-2, 1:-1]) / (2*dy)
v[1:-1, 1:-1] = -(psi[1:-1, 2:] - psi[1:-1, :-2]) / (2*dx)
# --- 可視化 ---
X, Y = torch.meshgrid(torch.linspace(0, lx, nx), torch.linspace(0, ly, ny), indexing='ij')
plt.figure(figsize=(8, 6))
plt.contourf(X, Y, u.T, levels=20, cmap='RdBu_r')
plt.colorbar(label='u-velocity')
plt.streamplot(X[0,:].numpy(), Y[:,0].numpy(), u.numpy(), v.numpy(), color='black')
plt.title(f"Cavity Flow (Re={re})")
plt.xlabel("x")
plt.ylabel("y")
plt.show()4. コードの深掘り解説
SOR法によるポアソン方程式の解決
solve_poisson 関数内のSOR法に注目してください。ポアソン方程式は、ある点の値が周囲の平均値と関係するという「楕円型」の方程式です。一度の計算では全体に情報が行き渡らないため、何度も繰り返し計算（反復）を行います。
ここで使っている rel_factor（加速係数）は、1.0より大きく2.0未満の値に設定することで、単純なヤコビ法よりも劇的に早く収束させることができます。
渦度の境界条件
流体計算で最も厄介なのが、壁での渦度です。速度は壁でゼロですが、渦度は壁で発生します。今回のコードでは、psi の値から壁の勾配を推定する近似式を使っています。蓋が動く omega[-1, :] の部分だけ、蓋の速度 1.0 を反映させる項が加わっているのがポイントです。
可視化：流線（Streamlines）
結果の表示に streamplot を使用しました。流れ関数 \(\psi\) の等高線を描いても同じような図になりますが、streamplot を使うと流れの向きが矢印で強調されるため、箱の中で水が大きな渦を巻いている様子がはっきりとわかります。
5. 流れ関数-渦度法の限界
この手法は非常に効率的で安定していますが、残念ながら現実のすべての問題に使えるわけではありません。
- 3次元への拡張が困難: 3次元になると流れ関数がスカラではなくベクトルになり、数式が極めて複雑になります。
- 境界条件の制約: 入口や出口があるような「開いた領域」では、境界での \(\psi\) や \(\omega\) の値を決めるのが難しくなります。
しかし、今回のように「閉じた空間の中の流れ」を素早く、かつ美しくシミュレーションしたい場合には、現在でも非常に有力な選択肢です。
6. まとめ
第13回では、2次元流体シミュレーションの第一歩として「流れ関数-渦度法」を学びました。
- 渦度 \(\omega\) と 流れ関数 \(\psi\) を使うことで、計算から圧力を消去できる。
- ポアソン方程式の解法として SOR法 を実装。
- ベンチマークである キャビティ流れ をPythonで再現。
この手法で、私たちは「2次元の広がり」と「回転する流れ」を手に入れました。
次回、第14回では、3次元への拡張も可能で、現代の標準的な手法となっているMAC（Marker-and-Cell）法に挑戦します。「圧力を消去する」のではなく「圧力と真っ向から向き合う」ための、スタガード格子という非常に巧妙なアイデアを解説します。
最後まで読んでいただきありがとうございます。