コーディングを意識したかなり丁寧なSIMPLE法：スタッガード格子版の離散化の解説 - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

以前にSIMPLE法のプログラミングに関する以下の記事を書いたので離散化の記事を書いておきます。

SIMPLE法の記事は世の中に山ほどありますが、ソースコードとセットで書いている人は少ないのでコーディングをしたい方はぜひ読んで見て下さい。

TL;DR

SIMPLE法は離散化はかなり煩雑
SIMPLE法のポイントは運動方程式の線形化処理
式展開の各段階で近似を行うため計算安定ため緩和係数が必須

はじめに
運動方程式の離散化
連続式の離散化
数値解法
緩和係数の導入
GitHub
参考書籍
参考サイト

はじめに

SIMPLE法は非圧縮性流体の計算方法として最も実務的な方法であると考えられますが、陰解法のため、陽解法と比較するとコーディングは難しいです。

私もそれほど頻繁に使用するわけではないので、久しぶりにプログラムを改良するとき等は結構気が重かったりします。本記事では、自分自身の備忘録も兼ねて、プログラミングを意識して離散化についてまとめました。

なお、簡単のためx,yの二次元で記載しています。

運動方程式の離散化

次のx方向の運動方程式のみを記述しています。

$\begin{align} &\phi = u \\ &\frac{\partial \phi}{\partial t}+\frac{\partial u \phi}{\partial x}+\frac{\partial v \phi}{\partial y} + \frac{1}{\rho} \frac{\partial P}{\partial x} - \frac{\partial }{\partial x} \left( \nu \dfrac{\partial \phi}{\partial x}\right) - \frac{\partial }{\partial y} \left( \nu \dfrac{\partial \phi}{\partial y}\right) = 0 \\ % \frac{\partial \phi}{\partial t}+\frac{\partial}{\partial x}\left(u \phi - \nu \dfrac{\partial \phi}{\partial x}\right) % +\frac{\partial }{\partial y}\left( v \phi - \nu \dfrac{\partial \phi}{\partial y}\right) % + \frac{1}{\rho} \frac{\partial P}{\partial x} % &= 0 \end{align}$

スタッガード格子のため、コントロールボリュームは下図のように設けます。

まずは左辺第2,3項の移流項のみを抽出して離散化します。

$\begin{align} \frac{\partial u \phi}{\partial x} + \frac{\partial v \phi}{\partial y} &= 0 \end{align}$

として式展開を行います。

一次風上差分を適用するとx微分の項は次のように展開されます。

$\begin{align} \frac{\partial u \phi}{\partial x} &= \frac{u^{n+1}_{Xp} \phi^{n+1}_{Xp} - u^{n+1}_{Xm} \phi^{n+1}_{Xm}}{\Delta x} \\ u^{n+1}_{Xp} &= \frac{u^{n+1}_{XP} + u^{n+1}_{C}}{2} \\ u^{n+1}_{Xm} &= \frac{u^{n+1}_{C} + u^{n+1}_{XM} }{2} \\ % \phi^{n+1}_{Xp} &= \left\{ \begin{array}{ll} \phi^{n+1}_{C} & (u^{n+1}_{Xp} \geq 0)\\ \phi^{n+1}_{XP} & (u^{n+1}_{Xp} \lt 0) \end{array} \right. \end{align}$

プログラミング上は以下のように記述します。

$\begin{align} u^{n+1}_{Xp} \phi^{n+1}_{Xp} = - max(-u^{n+1}_{Xp},0) \cdot \phi^{n+1}_{XP} + (max(-u^{n+1}_{Xp},0)+u^{n+1}_{Xp})\cdot \phi^{n+1}_{C} \end{align}$

同様に、

$\begin{align} u^{n+1}_{Xm} \phi^{n+1}_{Xm} = max(u^{n+1}_{Xm},0) \cdot \phi^{n+1}_{XM} - (max(u^{n+1}_{Xm},0)-u^{n+1}_{Xm})\cdot \phi^{n+1}_{C} \end{align}$

y微分の項も同様に展開します。

$\begin{align} \frac{\partial v \phi}{\partial y} &= \frac{v^{n+1}_{Yp} \phi^{n+1}_{Yp} - v^{n+1}_{Ym} \phi^{n+1}_{Ym}}{\Delta y} \\ v^{n+1}_{Yp} &= \frac{v^{n+1}_{YP} + v^{n+1}_{C}}{2} \\ v^{n+1}_{Ym} &= \frac{v^{n+1}_{C} + v^{n+1}_{YM} }{2} \\ \end{align}$

$\begin{align} v^{n+1}_{Yp} \phi^{n+1}_{Yp} &= - max(-v^{n+1}_{Yp},0) \cdot \phi^{n+1}_{YP} + (max(-v^{n+1}_{Yp},0)+v^{n+1}_{Yp})\cdot \phi^{n+1}_{C} \\ v^{n+1}_{Ym} \phi^{n+1}_{Ym} &= max(v^{n+1}_{Ym},0) \cdot \phi^{n+1}_{YM} - (max(v^{n+1}_{Ym},0)-v^{n+1}_{Ym})\cdot \phi^{n+1}_{C} \end{align}$

これらを整理すると、

$\begin{align} a_C \phi^{n+1}_{C} &= a_{XP} \phi^{n+1}_{XP} + a_{XM} \phi^{n+1}_{XM} + a_{YP} \phi^{n+1}_{YP} + a_{YM} \phi^{n+1}_{YM} \\ a_{XP}&= max(-u^{n+1}_{Xp},0) \Delta y \\ a_{XM}&= max(u^{n+1}_{Xm},0) \Delta y \\ a_{YP}&= max(-v^{n+1}_{Yp},0) \Delta x \\ a_{YM}&= max(v^{n+1}_{Ym},0) \Delta x \\ a_{C}&= a_{XP} + a_{XM} + a_{YP} + a_{YM} + (u^{n+1}_{Xp} - u^{n+1}_{Xm})\Delta y + (v^{n+1}_{Yp} - v^{n+1}_{Ym})\Delta x \end{align}$

$a_C$ の最後の2項は連続式を用いると0となり、キャンセルされます。

次に移流項以外の項を含めて離散化します。離散式は次のとおりです。

$\begin{align} &\frac{\partial \phi}{\partial t}+\frac{\partial u \phi}{\partial x}+\frac{\partial v \phi}{\partial y} + \frac{1}{\rho} \frac{\partial P}{\partial x} - \frac{\partial }{\partial x} \left( \nu \dfrac{\partial \phi}{\partial x}\right) - \frac{\partial }{\partial y} \left( \nu \dfrac{\partial \phi}{\partial y}\right) = 0 \\ &\frac{\phi^{n+1}_C - \phi^{n}_C}{\Delta t} +\frac{u^{n+1}_{Xp} \phi^{n+1}_{Xp} - u^{n+1}_{Xm} \phi^{n+1}_{Xm}}{\Delta x} +\frac{v^{n+1}_{Yp} \phi^{n+1}_{Yp} - v^{n+1}_{Ym} \phi^{n+1}_{Ym}}{\Delta y} \\ & + \frac{1}{\rho} \frac{P^{n+1}_{Xp} - P^{n+1}_{Xm}}{\Delta x} - \nu \frac{\phi^{n+1}_{XP} -2 \phi^{n+1}_{C} - \phi^{n+1}_{XM}}{\Delta x^2} - \nu \frac{\phi^{n+1}_{YP} -2 \phi^{n+1}_{C} - \phi^{n+1}_{YM}}{\Delta y^2} = 0 \end{align}$

移流項と同様に式展開を行うと次のとおりとなります。

$\begin{align} a_C \phi^{n+1}_{C} &= a_{XP} \phi^{n+1}_{XP} + a_{XM} \phi^{n+1}_{XM} + a_{YP} \phi^{n+1}_{YP} + a_{YM} \phi^{n+1}_{YM}+b - \frac{\Delta x \Delta y}{\rho} \frac{P^{n+1}_{Xp} - P^{n+1}_{Xm}}{\Delta x} \\ a_{XP}&= \left( max(-u^{n+1}_{Xp},0) - \frac{\nu}{\Delta x} \right)\Delta y \\ a_{XM}&= \left( max(u^{n+1}_{Xm},0) - \frac{\nu}{\Delta x} \right)\Delta y \\ a_{YP}&= \left( max(-v^{n+1}_{Yp},0) - \frac{\nu}{\Delta y} \right)\Delta x \\ a_{YM}&= \left( max(v^{n+1}_{Ym},0) - \frac{\nu}{\Delta y} \right)\Delta x \\ a_{C}&= a_{XP} + a_{XM} + a_{YP} + a_{YM} + \frac{\Delta x \Delta y}{\Delta t} + (u^{n+1}_{Xp} - u^{n+1}_{Xm})\Delta y + (v^{n+1}_{Yp} - v^{n+1}_{Ym})\Delta x \\ b &= \frac{\Delta x \Delta y}{\Delta t}\phi^{n}_C \end{align}$

$a_C$ の最後の2項は連続式を用いると0となり、キャンセルされます。

SIMPLE法では定常解を取り扱うことが多いため、 $\Delta t$ を無限大にとって $\Delta t$ を含む項を0とします。

連続式の離散化

同様に連続式も離散化します。

スタッガード格子のため、コントロールボリュームは下図のように設けます。添字の位置が運動方程式と異なる点に注意して下さい。

$\begin{align} &\frac{\partial u }{\partial x}+\frac{\partial v }{\partial y} = 0 \\ &\frac{u^{n+1}_{Xp} - u^{n+1}_{Xm} }{\Delta x} +\frac{v^{n+1}_{Yp} - v^{n+1}_{Ym} }{\Delta y} \\ &(u^{n+1}_{Xp} - u^{n+1}_{Xm}) \Delta y +(v^{n+1}_{Yp} - v^{n+1}_{Ym}) \Delta x = 0 \end{align}$

数値解法

手順1：仮の流速場の計算

前述の離散式の未知数は、 $u^{n+1},v^{n+1},P^{n+1}$ です。これらは陽的に計算できないため、反復法によって求めます。

$u,v,P$ の仮の値を $u^{*},v^{*},P^{*}$ とし、運動方程式に代入すると次のとおりになります。

$\begin{align} \phi^{*} &= u^{*} \\ a_C \phi^{*}_{C} &= a_{XP} \phi^{*}_{XP} + a_{XM} \phi^{*}_{XM} + a_{YP} \phi^{*}_{YP} + a_{YM} \phi^{*}_{YM}+b - \frac{\Delta x \Delta y}{\rho} \frac{P^{*}_{Xp} - P^{*}_{Xm}}{\Delta x} \end{align}$

手順1ではこの式より $u^{*},v^{*}$ について解きます。この式はTDMAによる方法、反復計算による方法等で計算可能です。（前回のプログラムでは反復計算SOR法を使用）

この解法がSIMPLE法のポイントの一つになりますので詳細に説明します。

この計算では、仮の流速の反復計算の一つ前のステップの値 $\phi^{*old}$ から次のステップの値 $\phi^{*new}$ を求めています。各項の係数 $a$ には $u^{*},v^{*}$ が含まれますが、上式の計算では反復計算の一つ前のステップ $u^{*old},v^{*old}$ のまま固定します。この線形化処理を行うことで上式が計算可能となります。

これらを用いて上式を書き直すと以下になります。

$\begin{align} \phi^{*new} &= u^{*new} \\ a_C \phi^{*new}_{C} &= a_{XP} \phi^{*new}_{XP} + a_{XM} \phi^{*new}_{XM} + a_{YP} \phi^{*new}_{YP} + a_{YM} \phi^{*new}_{YM}+b - \frac{\Delta x \Delta y}{\rho} \frac{P^{*}_{Xp} - P^{*}_{Xm}}{\Delta x} \\ a_{XP}&= \left( max(-u^{*old}_{Xp},0) - \frac{\nu}{\Delta x} \right)\Delta y \\ a_{XM}&= \left( max(u^{*old}_{Xm},0) - \frac{\nu}{\Delta x} \right)\Delta y \\ a_{YP}&= \left( max(-v^{*old}_{Yp},0) - \frac{\nu}{\Delta y} \right)\Delta x \\ a_{YM}&= \left( max(v^{*old}_{Ym},0) - \frac{\nu}{\Delta y} \right)\Delta x \\ a_{C}&= a_{XP} + a_{XM} + a_{YP} + a_{YM} + \frac{\Delta x \Delta y}{\Delta t} + (u^{*old}_{Xp} - u^{*old}_{Xm})\Delta y + (v^{*old}_{Yp} - v^{*old}_{Ym})\Delta x \\ b &= \frac{\Delta x \Delta y}{\Delta t}\phi^{n}_C \end{align}$

なお、係数 $a$ 中の $u^{*},v^{*}$ は反復計算で更新しますので、十分に収束していれば問題はないことになります。

手順2 : 圧力補正量の計算

次に、反復計算による補正量 $u^{'},v^{'},P^{'}$ を用いて、 $u^{n+1},v^{n+1},P^{n+1}$ を次のように定義します。

$\begin{align} u^{n+1} &= u^{*}+ u^{'}\\ v^{n+1} &= v^{*}+ v^{'}\\ P^{n+1} &= P^{*}+ P^{'} \end{align}$

これより、元の運動方程式と手順1の仮の値 $u^{*},v^{*},P^{*}$ を代入した運動方程式の差をとると次式が得られます。

$\begin{align} \phi^{'} &= u^{'} \\ a_C \phi^{'}_{C} &= a_{XP} \phi^{'}_{XP} + a_{XM} \phi^{'}_{XM} + a_{YP} \phi^{'}_{YP} + a_{YM} \phi^{'}_{YM} - \frac{\Delta x \Delta y}{\rho} \frac{P^{'}_{Xp} - P^{'}_{Xm}}{\Delta x} \end{align}$

補正量 $u^{'},v^{'}$ は十分に小さいと仮定して、それらを含む項を無視すると次式となります。

$\begin{align} \phi^{'}_{C} &= - \frac{1}{a_C }\frac{\Delta x \Delta y}{\rho} \frac{P^{'}_{Xp} - P^{'}_{Xm}}{\Delta x} \end{align}$

これが流速の補正量の式となります。

次にこの式を連続式に代入します。まずは、わかりやすくするためにコントロールボリュームにあわせて添字を修正します。

$\begin{align} u^{'}_{Xp} &= - \frac{1}{{a_C}|_{Xp}}\frac{\Delta x \Delta y}{\rho} \frac{P^{'}_{XP} - P^{'}_{C}}{\Delta x} \\ u^{'}_{Xm} &= - \frac{1}{{a_C}|_{Xm}}\frac{\Delta x \Delta y}{\rho} \frac{P^{'}_{C} - P^{'}_{XM}}{\Delta x} \\ v^{'}_{Yp} &= - \frac{1}{{a_C}|_{Yp}}\frac{\Delta x \Delta y}{\rho} \frac{P^{'}_{YP} - P^{'}_{C}}{\Delta y} \\ v^{'}_{Ym} &= - \frac{1}{{a_C}|_{Ym}}\frac{\Delta x \Delta y}{\rho} \frac{P^{'}_{C} - P^{'}_{YM}}{\Delta y} \end{align}$

また、連続式は次のとおりに変形します。

$\begin{align} (u^{n+1}_{Xp} - u^{n+1}_{Xm}) \Delta y +(v^{n+1}_{Yp} - v^{n+1}_{Ym}) \Delta x &= 0 \\ (u^{*}_{Xp} - u^{*}_{Xm})\Delta y + (v^{*}_{Yp} - v^{*}_{Ym})\Delta x + (u^{'}_{Xp} - u^{'}_{Xm})\Delta y + (v^{'}_{Yp} - v^{'}_{Ym})\Delta x &= 0 \end{align}$

上式に前述の流速の補正量の式を代入すると次式が得られます。

$\begin{align} A_C P^{'}_{C} &= A_{XP} P^{'}_{XP} + A_{XM} P^{'}_{XM} + A_{YP} P^{'}_{YP} + A_{YM} P^{'}_{YM}+B \\ A_{XP} &= \frac{1}{{a_C}|_{Xp}}\frac{\Delta y^2}{\rho} \\ A_{XM} &= \frac{1}{{a_C}|_{Xm}}\frac{\Delta y^2}{\rho} \\ A_{YP} &= \frac{1}{{a_C}|_{Yp}}\frac{\Delta x^2}{\rho} \\ A_{YM} &= \frac{1}{{a_C}|_{Ym}}\frac{\Delta x^2}{\rho} \\ A_{C}&= A_{XP} + A_{XM} + A_{YP} + A_{YM} \\ B &= - (u^{*}_{Xp} - u^{*}_{Xm})\Delta y - (v^{*}_{Yp} - v^{*}_{Ym})\Delta x \end{align}$

手順2ではこの式より $P^{'}$ を解きます。この式はTDMAによる方法、反復計算による方法等で計算可能です。（前回のプログラムでは反復計算SOR法を使用）

手順3 : 流速、圧力の更新

手順2で求めた $P^{'}$ により、 $P^{*new} = P^{*old} + P^{'}$ として、 $P^{*}$ を更新する。
手順2で求めた $P^{'}$ から、手順2に示す流速の補正量の式より $u^{'},v^{'}$ を求める。
$u^{'},v^{'}$ より、 $u^{*new} =u^{*old} + u^{'}$ 、 $v^{*new} = v^{*old} + v^{'}$ として、 $u^{*},v^{*}$ を更新する。

手順4 : 反復

$u^{'},v^{'},P^{'}$ が十分に小さくなるまで、手順1～3を繰り返す。

緩和係数の導入

SIMPLE法では各所で式の近似を行っているため、計算に不安定が生じやすいです。計算の安定性の観点から緩和係数を導入し、不足緩和を行います。

流速の不足緩和

全体：手順1～3の反復計算において不足緩和を行います。

具体的には反復計算の過程で得られた値 $\phi^{*new}$ を不足緩和係数 $\alpha_U$ を用いて、次のように定義する。

$\begin{align} &\phi^{*old} + \alpha_U(\phi^{*new} - \phi^{*old}) = \alpha_U \phi^{*new} + (1-\alpha_U) \phi^{*old} \rightarrow \phi^{*new} \end{align}$

$\phi^{*old}$ ：反復計算の一つ前のステップでの値

よって、手順1の式は次のように変形できます。

$\begin{align} \phi^{*new} &= u^{*new} \\ \dfrac{a_C}{\alpha_U} \phi^{*new}_{C} &= a_{XP} \phi^{*new}_{XP} + a_{XM} \phi^{*new}_{XM} + a_{YP} \phi^{*new}_{YP} + a_{YM} \phi^{*new}_{YM}+b - \frac{\Delta x \Delta y}{\rho} \frac{P^{*}_{Xp} - P^{*}_{Xm}}{\Delta x} \\ b &= \frac{\Delta x \Delta y}{\Delta t}\phi^{n}_C + (1-\alpha_U) \dfrac{a_C}{\alpha_U} \phi^{*old}_{C} \end{align}$