MAC系解法たち

Last updated at 2024-07-01Posted at 2024-04-23

初めに

数値シミュレーションには，幾つかの時間積分法がある．

本稿では，非圧縮性流体の数値シミュレーションにおける時間積分法に焦点を当て，

MAC法¹（Marker and Cell method, Harlow&Welch1965, Welch+1965）
射影法²（Projection method, Chorin1968）
SMAC法（Simplified MAC method, Amsden&Harlow1970）

を用いた計算を行ってみる（HSMAC法³: Highly Simplified MAC method, Hirt+1975も比較対象に入れたかったが，良く理解できなかったのと，これを用いた文献を見ないので，恐らく分野に浸透していないのだろうとして，議論の対象から省く）．レビュー論文として，例えば，McKee+2008がある（Jet buckling（例えば，Cruickshank&Munson1981, Cruickshank1988）やHydraulic jumpにも触れていて個人的に面白い）．上に列挙した手法は，（射影法を含めて）現在ではMAC系解法やMAC系の方法などと呼ばれている⁴．

陽解Euler法で速度を更新するので，時間方向に1次精度である（「殆ど変わらないじゃないか」と思う気持ちをぐっと堪えて，やってみるまで分からないことも多いので）．また，より高次の方法（例えば，Crank-Nicolson法，Adams-Bashforth/Moulton法など）もあるが，今回は議論の対象から省く（これらの手法をどう実装するのか未だ良く分かっていない，というのが本音）．1次精度の方法についての議論はShen1992が，2次精度の方法についてはShen1996が中々面白そうだ，が，まだじっくり読んでいない（宿題）．

Methods

支配方程式

無次元化した，非圧縮性のNavier-Stokes方程式を記す．

\begin{align}
    \nabla \cdot \mathbf{u} 
    &= 0 \\
    \frac{\partial \mathbf{u}}{\partial t} 
    + \left( \mathbf{u} \cdot \nabla \right) \mathbf{u} 
    &= - \nabla p 
    + \frac{1}{\mathrm{Re}} \nabla^2 \mathbf{u} 
    + \mathbf{f}
\end{align}

ここで，$\mathbf{u}$は速度，$p$は圧力，$\mathbf{f}$は外力，$\mathrm{Re}$はReynolds数である．空間次元は$M$次元とする（i.e., $\mathbf{x} \in \Omega \subset \mathbb{R}^M$）．

MAC法

陽的Euler法から，

\begin{align}
    \frac{\mathbf{u}^{(n+1)} - \mathbf{u}^{(n)}}{\Delta t} 
    + \left( \mathbf{u}^{(n)} \cdot \nabla \right) \mathbf{u}^{(n)} 
    &= - \nabla p^{(n+1)} 
    + \frac{1}{\mathrm{Re}} \nabla^2 \mathbf{u}^{(n)}
    + \mathbf{f}^{(n)} \\
    \mathbf{u}^{(n+1)} 
    &= \mathbf{u}^{(n)}
    + \Delta t
    \left(
        - \nabla p^{(n+1)}
        - \left( \mathbf{u}^{(n)} \cdot \nabla \right) \mathbf{u}^{(n)} 
        + \frac{1}{\mathrm{Re}} \nabla^2 \mathbf{u}^{(n)}
        + \mathbf{f}^{(n)}
    \right)
\end{align}

を得る．また，$\nabla \cdot \mathbf{u}^{(n+1)} = 0$を要求し，

\begin{align}
    \nabla^2 p^{(n+1)}
    &= 
    \frac{\nabla \cdot \mathbf{u}^{(n)}}{\Delta t}
    - \nabla \cdot \left( \left( \mathbf{u}^{(n)} \cdot \nabla \right) \mathbf{u}^{(n)} \right)
    + \nabla \cdot \left( \frac{1}{\mathrm{Re}} \nabla^2 \mathbf{u}^{(n)} \right)
\end{align}
    + \nabla \cdot \mathbf{f}^{(n)}

を得る．ここで，速度場の発散や移流の発散，拡散の発散など，本来ゼロとなるべきものや，微分を繰り返して十分ゼロに近いだろうという量がわざわざ残されているのは，数値計算の誤差を見越して，アルゴリズムに誤差の吸収を任せるためらしい．

すなわち，

Solve Poisson equation for $p^{(n+1)}$
Update velocity, $\mathbf{u}^{(n+1)}$

というアルゴリズムである．

射影法

圧力勾配の作用を取り除いて仮の速度を求めておき，その後圧力勾配の作用を考慮して速度を修正する（圧力場$p$を，速度場$\mathbf{u}$を発散無し場に引き戻すための回転無し場を生み出すスカラーポテンシャルと解釈する, cf. Chorin&Marsden1993, 過去の投稿）．

\begin{align}
    \frac{\mathbf{u}^{(*)} - \mathbf{u}^{(n)}}{\Delta t} 
    &= - \left( \mathbf{u}^{(n)} \cdot \nabla \right) \mathbf{u}^{(n)} 
    + \frac{1}{\mathrm{Re}} \nabla^2 \mathbf{u}^{(n)}
    + \mathbf{f}^{(n)} \\
    \frac{\mathbf{u}^{(n+1)} - \mathbf{u}^{(*)}}{\Delta t} 
    &= - \nabla p^{(n+1)} 
\end{align}

これらを足し合わせれば，元の式に戻る．

2式目で発散を取り，$\nabla \cdot \mathbf{u}^{(n+1)} = 0$を要求すれば，$p^{(n+1)}$は

\begin{align}
    \nabla^2 p^{(n+1)} = \frac{\nabla \cdot \mathbf{u}^{(*)}}{\Delta t}
\end{align}

の解である．

すなわち，

Predict velocity, $\mathbf{u}^{(*)}$
Solve Poisson equation for $p^{(n+1)}$
Correct velocity, $\mathbf{u}^{(n+1)} = \mathbf{u}^{(*)} + \Delta t (- \nabla p^{(n+1)})$

というアルゴリズムである（predictor–corrector methodの一つ）．

SMAC法

あるスカラーポテンシャル$q$を導入して，運動量保存則を分離する．

\begin{align}
    \frac{\mathbf{u}^{(*)} - \mathbf{u}^{(n)}}{\Delta t} 
    &= - \nabla p^{(n)} 
    - \left( \mathbf{u}^{(n)} \cdot \nabla \right) \mathbf{u}^{(n)} 
    + \frac{1}{\mathrm{Re}} \nabla^2 \mathbf{u}^{(n)} 
    + \mathbf{f}^{(n)} \\
    \frac{\mathbf{u}^{(n+1)} - \mathbf{u}^{(*)}}{\Delta t} 
    &= - \nabla q
\end{align}

足し合わせて元の式に戻るために，$p^{(n)} + q = p^{(n+1)}$で$q$を定義する．

2式目の発散を取り，$\nabla \cdot \mathbf{u}^{(n+1)} = 0$を要求すれば，$q$は

\begin{align}
    \nabla^2 q = \frac{\nabla \cdot \mathbf{u}^{(*)}}{\Delta t}
\end{align}

の解である．

すなわち，

Predict velocity, $\mathbf{u}^{(*)}$
Solve Poisson equation for $q$
Correct pressure, $p^{(n+1)} = p^{(n)} + q$
Correct velocity, $\mathbf{u}^{(n+1)} = \mathbf{u}^{(*)} + \Delta t (- \nabla q)$

というアルゴリズムである（predictor–corrector methodの一つ）．

ここまで来ると，これらをMAC系と呼んで取り纏める気持ちも分かるし，特に，射影法とSMAC法が良く似ていると感じる⁵．

Results

2次元キャビティ流れを取り上げる．やり尽くされた問題だが，コードも簡単なので，こういう比較には打って付けだろう．Arakawa B型のstaggered grid（Arakawa&Lamb1977）を用い，空間離散化は以下の通りとした（そろそろ，お馴染みになってきたね，cf. 以前の投稿）．

移流項　　: 3次精度風上差分（Kawamura&Kuwahara1984, Kawamura+1986）
拡散項　　: 2次精度中心差分
圧力勾配項: 2次精度中心差分

圧力は，定数を適当に動かす範囲内で一意に定まることが知られているので（cf. Courant&Hilbert1953），$\hat{p} = p - \bar{p}$を可視化する．ただし，$\bar{p}$は，当該タイムステップにおける圧力の場全体での平均値である．この操作は，カラーレンジを移動させる手間を省くためであり，計算には影響を与えない．

MAC法，射影法，SMAC法は，いずれも時間方向に1次精度なので，精度より，安定性に関係しそうな気がする．陽解法の安定性条件をギリギリまで攻め，違いが現れることを期待してみる．

離散化パラメータを$\Delta x = 1 \times 10^{-2}$に固定する．time marchingにだけ注目したいので，渦粘性（例えば，Smagorinsky1963, Lilly1966）などの特別な処理は入れていない．速度を陽的に更新するので，安定性条件（Courant+1967, Neumann&Richtmyer1950）を考える．Courant数$C$の条件，拡散数$D$の条件は，

\begin{align}
    C &= \sum_{j}^{M} \frac{\left| c_j \right| \Delta t_{C}}{\Delta x_{j}} 
    \le 1 \\
    D &= \sum_{j}^{M} \frac{d \Delta t_{D}}{\left( \Delta x_{j} \right)^{2}} 
    \le \frac{1}{2}
\end{align}

である．ここで，$c_j$は$x_{j}$方向の移流速度，$d$は拡散係数（拡散に方向性は考えないことにする），$\Delta t_{C}$，$\Delta t_{D}$はそれぞれ，Courant数，拡散数の安定性条件を満たす時間増分である．粘性が小さいときは前者が，粘性が大きいときは後者が厳しい条件となる．

拡散数を攻める

$M = 2$，$\Delta x_{j} = \Delta x = 1 \times 10^{-2}$を考えると，

\begin{align}
    D 
    = 2 \frac{d \Delta t_{D}}{\left( \Delta x \right)^{2}} 
    \le \frac{1}{2} 
    \quad
    \Leftrightarrow
    \quad
    \Delta t_{D}
    \le \frac{\mathrm{Re}}{4} \times 10^{-4}
\end{align}

$\mathrm{Re} = 100$を選べば，$\Delta t_{D} \le 2.5 \times 10^{-3} \left( \lt 5 \times 10^{-3} = \Delta t_{C}^{\mathrm{(crit)}} \right)$．

	MAC	Proj	SMAC
$\Delta t = 1.0 \times 10^{-3}$ $(C = 0.1, \ D = 0.2)$
$\Delta t = 1.5 \times 10^{-3}$ $(C = 0.15, \ D = 0.3)$
$\Delta t = 2.0 \times 10^{-3}$ $(C = 0.2, \ D = 0.4)$	(diverged)	(diverged)	(diverged)
$\Delta t = 2.5 \times 10^{-3}$ $(C = 0.25, \ D = 0.5)$	(diverged)	(diverged)	(diverged)

時間増分が違うので，gifのタイミングがズレているが，3回繰り返したら映像が止まるようにした（クリックして拡大表示）．

全ての時間進行スキームで，$D = 0.5$に到達するより早く，$D = 0.4$で発散した．$D = 0.2$程度に抑えておけば，全ての方法で安定にシミュレートできた．$D = 0.4$で発散したのも凡そ当然のように感じる（実用上は，こういうギリギリを攻めた時間増分を選ぶ，という作戦は取らないだろうし，もう少し安全側を採用する気もするが，Courant数だけ満足しても安定でない場合がある，と知ることができた，ということを持ち帰っておこう）．

圧力場は目に見える違いが少ないように感じるが，速度場を比較すると，MAC法のそれは，射影法，SMAC法のそれより若干遅い流れになっていそうだ．射影法，SMAC法の間には，目に見える違いは無いように感じる．

速度場の発散に注目すると，MAC法では，ゼロより大きい値に落ち着いているのが興味深い．Poisson方程式のソースがそもそも，正値の方向にバイアスされているから，これが原因だろうか？射影法では，領域中央付近では，凡そきちんとゼロになっている．SMAC法はさらに良い．境界付近でもかなりゼロに近い値が達成できている．Poisson方程式のソースが，非常にゼロに近い，というのも興味深い点である．

ここから，ソース項に由来して発散ゼロがきちんと満たせるか否かが決まってくる，と分かる（式を見返せば当然のことだ）．しかし，空間スキームは全く同じなのに，時間スキームに起因してこうした違いが出てくるのは面白い．

Courant数を攻める

$M = 2$，$\Delta x_{j} = \Delta x = 1 \times 10^{-2}$，$\boldsymbol{c}^{\mathrm{(crit)}} = \left( c_1, c_2 \right)$とすると，

\begin{align}
    C 
    = \Delta t_{C} \frac{\left| c_1 \right| + \left| c_2 \right|}{\Delta x} 
    \le 1
    \quad
    \Leftrightarrow
    \quad
    \Delta t_{C}
    \le \frac{1}{\left| c_1 \right| + \left| c_2 \right|} \times 10^{-2}
\end{align}

top driving lidをイメージして$\left( c_1, c_2 \right) = \left( 1, 0 \right)$，$\mathrm{Re} = 1,000$を選べば，$\Delta t_{C} \le 1 \times 10^{-2} \left( \lt 2.5 \times 10^{-2} = \Delta t_{D}^{\mathrm{(crit)}} \right)$．

	MAC	Proj	SMAC
$\Delta t = 2.0 \times 10^{-3}$ $(C = 0.2, \ D = 0.04)$
$\Delta t = 4.0 \times 10^{-3}$ $(C = 0.4, \ D = 0.08)$
$\Delta t = 6.0 \times 10^{-3}$ $(C = 0.6, \ D = 0.12)$	(diverged)	(diverged)	(diverged)
$\Delta t = 8.0 \times 10^{-3}$ $(C = 0.8, \ D = 0.16)$	(diverged)	(diverged)	(diverged)
$\Delta t = 1.0 \times 10^{-2}$ $(C = 1.0, \ D = 0.2)$	(diverged)	(diverged)	(diverged)

先ほどと同じように，$C = 1.0$に届くより早く，$C = 0.6$で既に爆発している．やはり，拡散数だけ満足しても安定とならない場合があることが分かる．実現象では外力も作用してくるので，かつ上記の条件は線形安定性解析の結果なので，実用上は$\Delta t = \min \left( \Delta t_{C}, \Delta t_{D} \right)$に，さらに安全係数を乗じたものを選ぶのが自然か（cf. 越塚1997）．

結果を見ると，MAC法は全く粗末なものである．これは非常に面白い．線形ソルバーがpoint Jacobi法なので，ここが原因の可能性は否定できないが，上述の通り，全ての時間進行スキームでソルバーも空間離散化手法も揃えているので，速度と圧力のカップリングの重要性が示唆される．射影法とSMAC法は，先程と同じく，あまり違いが見られない．SMAC法の方が質量保存に優れている（かもしれない）ところまで，先程と同じだ．

終わりに

現在では，射影法やSMAC法が広く使われている気がする．HSMAC法は，個人的には解釈が難しいと感じる．また，これを使っている論文を見たことが未だ無い（ので，省いた）．

今回の数値実験の範囲（比較的静かな流れを考えた，かなり限られた範囲）では，安定性については，どれも同じくらいであることが分かった（特定のスキームを選ぶことで，大胆な時間増分を選ぶことができるようになるわけではないのだろう）．一方，速度と圧力のカップリングが極めて重要であることが分かった（MAC法の劣悪さは面白い）．射影法とSMAC法を見比べると，P1/P2 correction schemeと呼ばれているのが納得できる気がしてくる．数字が大きい方が高級なのだとしたら，圧力自体をソルバーに引き渡す射影法（P1）より，圧力の修正量を考えるSMAC法（P2）の方が高級なのかもしれない．ただし，SMAC法では，変数が増えるので，若干だけメモリの消費があるか？（詳しく調べていない）．

現代では，有限差分法にはArakawa C型のstaggered grid（Arakawa&Lamb1977）を用いて非圧縮性を満足する気がするが，SPH（Smoothed-particle hydrodynamics）法（Gingold&Monaghan77, Lucy77）やMPS（Moving Particle Semi-implicit）法（Koshizuka&Oka1996）などに代表される粒子法では，同じ粒子に複数の物理量が埋め込まれるので，本稿の議論は有益になり得るだろうか．

MAC法は，PIC（Particle in Cell）法（Evans&Harlow1957, Harlow1957）のvariantとして登場したものらしい． ↩
（私の解釈では）文脈によっては，部分段階法（Fractional Step method，Témam1969）と呼ばれる/混同することもあるが（私だけ？），部分段階法の方が枠組みがもう少し広い（と思っている）．非圧縮性流体の時間積分法としては，両者は結果的に同じものに帰着することもあるが，本稿では，射影法と呼ぶこととする（Helmholtzの分解の定理と密接な関係にある，部分段階法の視点からは2段階の分解に当たる）． ↩
HSMAC法のことを，Hirt+1975はSOLAと呼んでいる（どうやらSOLution Algorithmの略らしい）．Hirt+1975: "This report describes a highly simplified MAC code, SOLA, that does not use marker particles and does not have built-in setups for internal obstacles or other complicating refinements." ↩
元々，MAC法は「(1)空間離散化においてLagrange記述のMarkerとEuler記述のCellを用いた流体計算法を導入した」「(2)時間離散化において速度と圧力を分離する分離型解法を導入した」ものだった．当時は，(1)を尊重してMarker and Cell法と呼んだ．一方，現代では，空間離散化については，有限要素法，有限体積法，粒子法などが台頭してきたため，忘れられ，(2)の特性を引き継いだものをMAC系の解法と呼ぶようになったようだ．(1)に最小限だけ触れておくと，MAC法では，質量を持たないMarker（Lagrangian記述）の移動を追跡しつつ支配方程式をCell（Euler記述）上で計算し，もし当該CellにMarkerが存在すればそのCellを流体が詰まったCellとし，Markerが無ければ流体を含まないCellと判定する，という形を取ったらしい．ここまで来て，MPMと似ているか？と思ったが，本題ではないし，今回は放っておく． ↩
射影法のことをP1 correction scheme，SMAC法のことをP2 correction schemeと呼ぶこともあるらしいが（参考，13:30~），これらのキーワードをネットで調べても中々ヒットしなかったので，本稿では，射影法，SMAC法と呼ぶこととする．確かに，中身を見ると非常に似ているし，もし射影法をP1と呼ぶならSMAC法をP2と呼ぶ気持ちも理解できるが，名前の浸透具合には未だ疑問が残る（或いは私の勉強不足のせいかもしれないが）． ↩

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up