サイトアイコン アマチュア無線局JS2IIU

Pythonで攻略する流体力学の数値計算 第14回：2次元流れの解法(2) ― MAC法の基本アルゴリズム

JS2IIU

こんにちはJS2IIUです。Pythonで攻略する流体力学シミュレーション連載の第14回です。

前回は、2次元限定で圧力を消去できる「流れ関数-渦度法」を学びました。非常にスマートな手法でしたが、3次元への拡張が難しく、境界条件の扱いも一癖あるという課題がありました。

今回からはいよいよ、現代の数値流体力学（CFD）において最も標準的であり、商用・オープンソースを問わず多くの解析ソフトの基盤となっている手法、MAC（Marker-and-Cell）法の解説に入ります。

MAC法は1次元や2次元だけでなく、そのまま3次元へと拡張が可能です。その最大の特徴は、圧力と速度を正面から扱い、それらを「チェス盤」のようにずらして配置するスタガード格子という巧妙なアイデアにあります。今回もよろしくお願いします。

1. なぜ「ずらす」のか？ ― スタガード格子の魔法

1.1 コロケート格子とチェッカーボード現象

これまで、速度 \(u, v\) や圧力 \(p\) はすべて同じ格子点上に定義してきました。これをコロケート格子 (Collocated Grid) と呼びます。一見自然な配置に思えますが、実は非圧縮性流体計算において深刻な問題を引き起こします。

非圧縮性ナビエ・ストークス方程式の運動量保存式は次のように書けます：

$$
\frac{\partial u}{\partial t} = -\frac{1}{\rho}\frac{\partial p}{\partial x} + \text{(その他の項)}
$$

この圧力勾配項を、格子点 \(i\) において中心差分で離散化すると：

$$
\left.\frac{\partial p}{\partial x}\right|_{i}
\approx
\frac{p_{i+1} – p_{i-1}}{2\Delta x}
$$

ここで重要なのは、この式に \(p_i\) 自身が現れていないことです。つまり、\(i-1\) と \(i+1\) の圧力値だけから勾配を計算しており、注目している地点 \(i\) の圧力値は何の影響も与えません。

この性質が何を意味するのでしょうか？例えば、圧力場が次のような交互パターンを持っていたとします：

$$
p_0 = 0, \quad p_1 = 1, \quad p_2 = 0, \quad p_3 = 1, \quad p_4 = 0, \quad \cdots
$$

このとき、各点での圧力勾配を計算してみましょう：

$$
\left.\frac{\partial p}{\partial x}\right|_1 = \frac{0 – 0}{2\Delta x} = 0
$$

$$
\left.\frac{\partial p}{\partial x}\right|_2 = \frac{1 – 1}{2\Delta x} = 0
$$

$$
\left.\frac{\partial p}{\partial x}\right|_3 = \frac{0 – 0}{2\Delta x} = 0
$$

驚くべきことに、すべての地点で圧力勾配がゼロと計算されてしまいます！実際には圧力が激しく振動しているにもかかわらず、中心差分はそれを「滑らかな場」と誤認識するのです。これがチェッカーボード現象 (Checkerboard Oscillation) です。

さらに、連続の式（質量保存則）も同様の問題を抱えています：

$$
\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0
$$

コロケート格子で中心差分を用いると：

$$
\frac{\partial u}{\partial x}\bigg|_{i}
\approx
\frac{u_{i+1} – u_{i-1}}{2\Delta x}
$$

この場合も \(u_i\) が式に現れず、速度場に2格子幅の振動モードが発生しても検出されません。

1.2 スタガード格子による根本的解決

この問題を解決するために、1965年にF. H. HarlowとJ. E. Welchによって提案されたのが、スタガード格子 (Staggered Grid)、またの名をMAC格子 (Marker-And-Cell Grid) です。

物理量の配置原理

スタガード格子では、各物理量を次のように配置します：

この配置には深い物理的意味があります。

圧力勾配と速度の自然な結びつき

\(x\) 方向の運動量方程式を考えます。セルの辺に位置する速度 \(u_{i+1/2,j}\) に働く圧力勾配は：

$$
\left.\frac{\partial p}{\partial x}\right|_{i+1/2,j}
\approx
\frac{p_{i+1,j} – p_{i,j}}{\Delta x}
$$

この差分式は、速度 \(u_{i+1/2,j}\) がまさに存在する位置（2つの圧力セルの境界）における勾配を直接的に表現しています。隣接する2つの圧力値だけを使い、その間の位置での勾配を計算する—これは物理的に極めて自然です。

同様に、連続の式における発散も：

$$
\nabla \cdot \mathbf{u}\bigg|_{i,j}
\approx
\frac{u_{i+1/2,j} – u_{i-1/2,j}}{\Delta x} + \frac{v_{i,j+1/2} – v_{i,j-1/2}}{\Delta y}
$$

セル \((i, j)\) への流入・流出を、そのセルを囲む4つの辺上の速度で直接評価できます。この式は、まさに「セルへの質量流入量－流出量」という物理的意味を忠実に表現しています。

チェッカーボード現象の根絶

スタガード格子では、圧力勾配の計算に必ず隣接する2つの圧力値のみを使用します。したがって、もし圧力にチェッカーボード状の振動があれば：

$$
p_i = 0, \quad p_{i+1} = 1 \quad \Rightarrow \quad \frac{p_{i+1} – p_i}{\Delta x} = \frac{1}{\Delta x}
$$

この大きな勾配が速度に強く作用し、振動を即座に減衰させます。つまり、チェッカーボード模様は物理的に許されない状態として自動的に排除されるのです。

1.3 数学的な安定性の視点

より厳密に言えば、スタガード格子は離散化された方程式系の固有値構造を改善します。コロケート格子では、圧力と速度のカップリング行列に虚固有値（振動モード）が存在しますが、スタガード格子ではこれが消失し、すべての固有値が実数となります。これは数値安定性において決定的な違いです。

また、スタガード格子で離散化された連続の式と運動量方程式は、連続体の場合の以下の重要な性質を保存します：

  1. 運動量保存則の厳密な満足：各セル境界での運動量フラックスが厳密に釣り合う
  2. 運動エネルギーの2次精度保存：数値拡散が最小限に抑えられる
  3. 圧力ポアソン方程式の対称性：解の一意性と計算効率が保証される

このように、スタガード格子は単なる「便利なトリック」ではなく、物理法則を離散世界で最も忠実に再現するための必然的な選択なのです。

2. MAC法のアルゴリズムの流れ

MAC法の目的は、「質量保存（連続の式）を満足するように、各ステップで圧力を決定し、速度場を更新する」ことです。この手法は、数学的にはHelmholtz-Hodge分解に基づいており、任意のベクトル場を「発散ゼロの場」と「勾配場」に分解するという定理を応用しています。

2.1 基本的な考え方：射影法（Projection Method）

非圧縮性ナビエ・ストークス方程式は次の2つの式で構成されます：

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

$$
\nabla \cdot \mathbf{u} = 0
$$

ここで重要なのは、圧力 \(p\) は独立した時間発展方程式を持たず、連続の式を満たすための従属変数であることです。つまり、圧力は「速度場が発散を持たないように調整する力」として機能します。

MAC法では、この「調整」を次の3ステップに分解して実行します。

2.2 ステップ1：仮の速度の計算

まず、圧力勾配項を一時的に無視し、移流項と粘性項だけから仮の速度 \(\mathbf{u}^*\) を計算します：

$$
\mathbf{u}^* = \mathbf{u}^n + \Delta t \left[ -(\mathbf{u}^n \cdot \nabla)\mathbf{u}^n + \nu \nabla^2 \mathbf{u}^n \right]
$$

この式は、次のような時間積分スキームに対応しています：

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

なぜ圧力項を省略するのか？

完全なナビエ・ストークス方程式を一度に解こうとすると、圧力と速度が互いに依存し合う連立方程式になります。これを直接解くのは計算コストが非常に高く、また数値的にも不安定です。

そこで、まず「圧力がない状態で流体がどう動きたがっているか」を計算し、その後「発散が生じてしまった分を圧力で修正する」という2段階のアプローチをとります。これが計算効率と安定性の両方を実現する鍵です。

仮の速度の物理的意味

\(\mathbf{u}^\) は、慣性力（移流項）と粘性力だけに従って流体が動いた場合の速度場です。この時点では一般に \(\nabla \cdot \mathbf{u}^ \neq 0\) であり、流体が「圧縮・膨張」してしまっています。非圧縮性条件を満たすためには、この発散を打ち消す必要があります。

2.3 ステップ2：圧力のポアソン方程式の解決

次に、仮の速度を正しい（発散ゼロの）速度に修正する圧力場を求めます。完全な運動量方程式から：

$$
\frac{\mathbf{u}^{n+1} – \mathbf{u}^n}{\Delta t} = -(\mathbf{u}^n \cdot \nabla)\mathbf{u}^n + \nu \nabla^2 \mathbf{u}^n – \frac{1}{\rho}\nabla p^{n+1}
$$

ステップ1の式と比較すると：

$$
\frac{\mathbf{u}^{n+1} – \mathbf{u}^*}{\Delta t} = -\frac{1}{\rho}\nabla p^{n+1}
$$

整理すると：

$$
\mathbf{u}^{n+1} = \mathbf{u}^* – \frac{\Delta t}{\rho} \nabla p^{n+1}
$$

ここで、新しい速度 \(\mathbf{u}^{n+1}\) は必ず連続の式を満たさなければならないので：

$$
\nabla \cdot \mathbf{u}^{n+1} = 0
$$

この条件を上式に代入します：

$$
\nabla \cdot \mathbf{u}^{n+1} = \nabla \cdot \mathbf{u}^* – \frac{\Delta t}{\rho} \nabla \cdot (\nabla p^{n+1}) = 0
$$

ラプラシアン \(\nabla^2 = \nabla \cdot \nabla\) を用いて整理すると、圧力のポアソン方程式が得られます：

$$
\nabla^2 p^{n+1} = \frac{\rho}{\Delta t} \nabla \cdot \mathbf{u}^*
$$

方程式の物理的解釈

右辺の \(\nabla \cdot \mathbf{u}^*\) は、仮の速度場の発散、すなわち「どこで流体が溜まっているか（正の発散）、または不足しているか（負の発散）」を表します。

左辺の \(\nabla^2 p^{n+1}\) は、圧力の2階微分であり、圧力がどのように空間的に変化して流体を「押し戻す」かを決定します。

この圧力分布が、ちょうど発散をゼロにするように自動的に決まるのです。

境界条件の重要性

ポアソン方程式を解くには境界条件が必要です。一般的に壁面では：

$$
\frac{\partial p}{\partial n} = 0 \quad \text{(Neumann条件)}
$$

これは「壁を貫通する流れがない」という条件から導かれます。壁に垂直な速度成分 \(u_n\) はゼロなので：

$$
\frac{\partial u_n}{\partial t} = -\frac{1}{\rho}\frac{\partial p}{\partial n} = 0 \quad \Rightarrow \quad \frac{\partial p}{\partial n} = 0
$$

2.4 ステップ3：速度の修正（射影）

求まった圧力場を用いて、仮の速度を修正します：

$$
\mathbf{u}^{n+1} = \mathbf{u}^* – \frac{\Delta t}{\rho} \nabla p^{n+1}
$$

なぜ必ず連続の式を満たすのか？

この修正式の両辺の発散をとってみましょう：

$$
\nabla \cdot \mathbf{u}^{n+1} = \nabla \cdot \mathbf{u}^* – \frac{\Delta t}{\rho} \nabla^2 p^{n+1}
$$

ここで、ステップ2で解いたポアソン方程式 \(\nabla^2 p^{n+1} = \frac{\rho}{\Delta t} \nabla \cdot \mathbf{u}^*\) を代入すると：

$$
\nabla \cdot \mathbf{u}^{n+1} = \nabla \cdot \mathbf{u}^* – \frac{\Delta t}{\rho} \cdot \frac{\rho}{\Delta t} \nabla \cdot \mathbf{u}^* = 0
$$

完璧に打ち消し合います！これは数学的に保証された結果であり、丸め誤差の範囲内で常に成立します。

Helmholtz-Hodge分解との関係

数学的に言えば、このステップは任意のベクトル場 \(\mathbf{u}^*\) を次のように分解しています：

$$
\mathbf{u}^* = \mathbf{u}^{n+1} + \nabla \phi
$$

ここで \(\mathbf{u}^{n+1}\) は発散ゼロ（ソレノイダル）成分、\(\nabla \phi\) は勾配（ポテンシャル）成分です。\(\phi = \frac{\Delta t}{\rho}p^{n+1}\) と対応しており、ステップ3は「仮の速度場から勾配成分を引き算して、発散ゼロ成分だけを取り出す」射影（Projection）操作なのです。

2.5 アルゴリズムの数学的性質

このMAC法（射影法）のアルゴリズムは、以下の優れた性質を持ちます：

  1. 質量保存の厳密性：ポアソン方程式の解が存在する限り、\(\nabla \cdot \mathbf{u}^{n+1} = 0\) が機械精度で成立
  2. 無条件安定性：圧力項の扱いが陰的であるため、時間刻み幅に関する安定性制約がない（CFL条件は移流項・粘性項からのみ）
  3. 効率性：圧力と速度の連立方程式を解く代わりに、ポアソン方程式という単一のスカラー方程式を解くだけで済む
  4. 3次元への拡張性：アルゴリズムの構造が次元に依存しないため、そのまま3次元に適用可能

この数学的エレガンスと計算効率の両立こそが、MAC法が半世紀以上にわたって流体計算の主流であり続ける理由です。

3. Pythonによる実装：MAC法で解くキャビティ流れ

それでは、前回と同じく「キャビティ流れ」を題材に、MAC法をPython（PyTorch）で実装してみましょう。スタガード格子のインデックス管理に注意しながら、アルゴリズムの骨格を構築します。

Python
import torch
import matplotlib.pyplot as plt

# --- 設定パラメータ ---
nx, ny = 41, 41
lx, ly = 1.0, 1.0
dx, dy = lx / (nx - 1), ly / (ny - 1)
dt = 0.001
rho = 1.0
nu = 0.1
nit = 50 # 圧力ポアソン方程式の反復回数

# --- 変数の初期化 (スタガード格子) ---
# p: (ny, nx) セル中心
# u: (ny, nx+1) x方向の辺の中点（左右の辺）
# v: (ny+1, nx) y方向の辺の中点（上下の辺）
p = torch.zeros((ny, nx))
u = torch.zeros((ny, nx + 1))
v = torch.zeros((ny + 1, nx))

def avg_x(field):
    return 0.5 * (field[:, 1:] + field[:, :-1])

def avg_y(field):
    return 0.5 * (field[1:, :] + field[:-1, :])

def build_rhs(u, v, dx, dy, dt, rho):
    """圧力ポアソン方程式の右辺(RHS)を計算する"""
    rhs = torch.zeros((ny, nx))
    # 連続の式からのズレ(発散)を計算
    # 1/dt * div(u*)
    rhs[1:-1, 1:-1] = (rho / dt) * (
        (u[1:-1, 2:-1] - u[1:-1, 1:-2]) / dx +
        (v[2:-1, 1:-1] - v[1:-2, 1:-1]) / dy
    )
    return rhs

def pressure_poisson(p, rhs, dx, dy):
    """SOR法を用いて圧力のポアソン方程式を解く"""
    for _ in range(nit):
        pn = p.clone()
        p[1:-1, 1:-1] = (
            ( (pn[1:-1, 2:] + p[1:-1, :-2]) * dy**2 +
              (pn[2:, 1:-1] + p[:-2, 1:-1]) * dx**2 -
              rhs[1:-1, 1:-1] * dx**2 * dy**2 ) /
            (2 * (dx**2 + dy**2))
        )
        
        # 圧力の境界条件 (Neumann条件: dp/dn = 0)
        p[:, -1] = p[:, -2] # 右
        p[0, :] = p[1, :]   # 下
        p[:, 0] = p[:, 1]   # 左
        p[-1, :] = 0.0      # 上(基準点として0に固定)
    return p

def mac_step(u, v, p, dx, dy, dt, rho, nu):
    """MAC法の1ステップを実行"""
    un = u.clone()
    vn = v.clone()
    
    # 1. 仮の速度の計算 (簡易化のため中央差分を使用)
    # u-momentum
    v_on_u = avg_y(avg_x(vn))[1:-1, :]
    u_adv_x = un[1:-1, 1:-1] * (un[1:-1, 2:] - un[1:-1, :-2]) / (2 * dx)
    u_adv_y = v_on_u * (un[2:, 1:-1] - un[:-2, 1:-1]) / (2 * dy)
    u_diff = (
        (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, :-2]) / dx**2
        + (un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[:-2, 1:-1]) / dy**2
    )
    u[1:-1, 1:-1] = un[1:-1, 1:-1] + dt * (-u_adv_x - u_adv_y + nu * u_diff)
    
    # v-momentum
    u_on_v = avg_x(avg_y(un))[:, 1:-1]
    v_adv_x = u_on_v * (vn[1:-1, 2:] - vn[1:-1, :-2]) / (2 * dx)
    v_adv_y = vn[1:-1, 1:-1] * (vn[2:, 1:-1] - vn[:-2, 1:-1]) / (2 * dy)
    v_diff = (
        (vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, :-2]) / dx**2
        + (vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[:-2, 1:-1]) / dy**2
    )
    v[1:-1, 1:-1] = vn[1:-1, 1:-1] + dt * (-v_adv_x - v_adv_y + nu * v_diff)

    # 境界条件 (Lid-driven cavity)
    u[-1, :] = 1.0  # Top lid velocity
    u[0, :] = 0.0; u[:, 0] = 0.0; u[:, -1] = 0.0
    v[-1, :] = 0.0; v[0, :] = 0.0; v[:, 0] = 0.0; v[:, -1] = 0.0

    # 2. 圧力のポアソン方程式
    rhs = build_rhs(u, v, dx, dy, dt, rho)
    p = pressure_poisson(p, rhs, dx, dy)
    
    # 3. 速度の修正
    u[1:-1, 1:-1] -= dt / rho * (p[1:-1, 1:] - p[1:-1, :-1]) / dx
    v[1:-1, 1:-1] -= dt / rho * (p[1:, 1:-1] - p[:-1, 1:-1]) / dy
    
    return u, v, p

# --- シミュレーション実行 ---
for n in range(500):
    u, v, p = mac_step(u, v, p, dx, dy, dt, rho, nu)

# --- 可視化 (セル中心での速度に補間) ---
u_center = 0.5 * (u[:, 1:] + u[:, :-1])
v_center = 0.5 * (v[1:, :] + v[:-1, :])

plt.figure(figsize=(8, 6))
plt.streamplot(torch.linspace(0, 1, nx).numpy(), torch.linspace(0, 1, ny).numpy(), 
               u_center.numpy(), v_center.numpy(), color=p.numpy(), cmap='jet')
plt.colorbar(label='Pressure')
plt.title("Cavity Flow Simulation using MAC Method")
plt.savefig('cavity_mac.png')
plt.show()

4. アルゴリズムの詳細解説

スタガード格子におけるインデックスの扱い

コードの中で、u[1:-1, 1:-2] や v[1:-2, 1:-1] といった複雑なスライシングが登場します。これはスタガード格子の宿命です。

圧力ポアソン方程式の役割

この方程式を解くステップがMAC法の「心臓部」です。仮の速度を代入して計算された rhs は、その地点でどれだけ流体が「溜まろうとしているか（または足りないか）」を表します。ポアソン方程式を解いて得られる圧力場は、その溜まりを解消するために流体を押し戻す力を提供します。結果として、流体はどこにも溜まることなくスムーズに流れる（連続の式を満たす）ようになります。

境界条件の重要性

キャビティ流れでは、上端の境界条件 u[-1, :] = 1.0 がすべての動きの源泉です。この速度が粘性を介して下層に伝わり、箱全体に大きな循環流を作ります。

また、セクション2.3で理論的に導出したように、圧力の境界条件として周囲の壁で勾配ゼロ（Neumann条件 \(\frac{\partial p}{\partial n} = 0\)）を課すことが、領域内の質量バランスを保つために不可欠です。これにより、壁を貫通する流れが物理的に禁止され、領域内の質量が厳密に保存されます。

5. MAC法の展望と次回の予告

MAC法は、圧力と速度のカップリングという非圧縮性流体最大の難問に対して、非常に堅牢で物理的にクリアな回答を示しました。しかし、実装を見て感じた通り、ポアソン方程式の反復計算（pressure_poisson）には多くの時間がかかります。

また、今回は「中央差分」を使用しましたが、第11回で学んだ通り、レイノルズ数が高くなると移流項の計算で振動が発生します。実用的なコードでは、ここに上流差分や高精度スキームを組み込む必要があります。

次回、第15回では、このMAC法をさらに洗練させたSMAC法や、計算効率を極限まで高めたプロジェクション法（フラクショナル・ステップ法）を解説します。現代の最速アルゴリズムが、どのような数学的トリックを使っているのかを解き明かしていきましょう。

最後まで読んでいただきありがとうございます。

モバイルバージョンを終了