【Navier-Stokes方程式】SIMPLE法によるNavier-Stokes方程式の離散化

SIMPLE法（Semi-Implicit Method for Pressure. Linked Equations）によるNavier-Stokes方程式の離散化について説明していきます。やっとこさSIMPLE法が理解出来ました。いきなり空間離散化して項を消したりするのでわかりづらいのです。実はこれHSMAC法とほぼ同じなんです！『数値流体工学』のpp.158-162を読んでいてはっとしました。

数値流体工学

作者: 荒川忠一
出版社/メーカー: 東京大学出版会
発売日: 1994/01/01
メディア: 単行本
この商品を含むブログを見る

この手法は流速と圧力を直接未知数として計算する手法です。SIMPLE法では圧力は完全に陰的に、流速は半陰的に取り扱われます。また、昔はよく定常状態を計算するために用いられていました。

以下の内容は『流体解析の基礎』のpp.215-220を参考にしています。

流体解析の基礎

作者: 河村哲也
出版社/メーカー: 朝倉書店
発売日: 2014/03/26
メディア: 単行本
この商品を含むブログ (1件) を見る

簡単のため二次元直交座標系で考えます。まず保存型のNavier-Stokes方程式は

$\displaystyle \begin{eqnarray} \frac{\partial u}{\partial t} + \frac{\partial (u^2)}{\partial x} + \frac{\partial (uv)}{\partial y} &=& -\frac{\partial p}{\partial x} + \frac{1}{Re} \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right) \\ \frac{\partial v}{\partial t} + \frac{\partial (uv)}{\partial x} + \frac{\partial (v)^2}{\partial y} &=& -\frac{\partial p}{\partial y} + \frac{1}{Re} \left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}\right) \end{eqnarray}$

です。連続式は

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

です。Navier-Stokes方程式と連続式を以下のように時間に対して離散化します。

$\displaystyle \begin{eqnarray} \frac{u^{n+1}-u^n}{\Delta t} &+& \frac{\partial (u^n u^{n+1})}{\partial x} + \frac{\partial (v^n u^{n+1})}{\partial y} \\ &=& -\frac{\partial p^{n+1}}{\partial x} + \frac{1}{Re} \left( \frac{\partial^2 u^{n+1}}{\partial x^2} + \frac{\partial^2 u^{n+1}}{\partial y^2}\right) \\ \frac{v^{n+1}-v^n}{\Delta t} &+& \frac{\partial (u^n v^{n+1})}{\partial x} + \frac{\partial (v^n v^{n+1})}{\partial y} \\&=& -\frac{\partial p^{n+1}}{\partial y} + \frac{1}{Re} \left( \frac{\partial^2 v^{n+1}}{\partial x^2} + \frac{\partial^2 v^{n+1}}{\partial y^2}\right) \\ &&\frac{\partial u^{n+1}}{\partial x} + \frac{\partial v^{n+1}}{\partial y} = 0 \end{eqnarray}$

圧力は完全に陰的に（n+1ステップの値を用いる）、流速は半陰的に（nステップとn+1ステップの値を用いる）取り扱われます。半陰的に取り扱うことにより非線型項を線型化し、非線型代数方程式を解くことを回避しています。

さて、上式の運動量方程式をある点Pで離散化（中心差分でも一次の風上差分でもok）すると、以下のように記述できます。つまりP（中心）の周囲の点W（西）、E（東）、N（北）、S（南）の5点の情報を使って離散化できます。

$\displaystyle \begin{eqnarray} a_P u_P &=& a_W u_W + a_E u_E + a_N u_N + a_S u_S + \Delta y (p_E - p_W) + b\\ c_P v_P &=& c_W v_W + c_E v_E + c_N v_N + c_S v_S + \Delta x (p_N - p_S) + d \end{eqnarray}$

ここで、 $u、v$ は未知数（n+1ステップでの値）で、 $a、b、c、d$ は係数です。具体的な係数は別の記事でまとめます。今回はそれよりも「計算の流れ」に注目しましょう。

次に、推測した圧力を $p^*$ が与えられたとして（例えば $p^n$ ）、上式の離散化式を用いて計算した流速を $u^*,v^*$ とします。真の値 $p^{n+1}、u^{n+1}、v^{n+1}$ とのずれ（補正量）をそれぞれ $p'、u'、v'$ とすると

$\displaystyle \begin{eqnarray} p^{n+1} &=& p^* + p'\\ u^{n+1} &=& u^* + u'\\ v^{n+1} &=& v^* + v' \end{eqnarray}$

と書けます。これが圧力の補正式になります。真の値はもとのNavier-Stokes方程式を、 $u^*,v^*$ は上式の離散化式を満たすことに注意して下さい。すなわち

$\displaystyle \begin{eqnarray} a_P u_P^{n+1} &=& a_W u_W^{n+1} + a_E u_E^{n+1} + a_N u_N^{n+1} + a_S u_S^{n+1} + \Delta y (p_E^{n+1} - p_W^{n+1}) + b\\ c_P v_P^{n+1} &=& c_W v_W^{n+1} + c_E v_E^{n+1} + c_N v_N^{n+1} + c_S v_S^{n+1} + \Delta x (p_N^{n+1} - p_S^{n+1}) + d \end{eqnarray}$

$\displaystyle \begin{eqnarray} a_P u_P^* &=& a_W u_W^* + a_E u_E^* + a_N u_N^* + a_S u_S^* + \Delta y (p_E^* - p_W^*) + b\\ c_P v_P^* &=& c_W v_W^* + c_E v_E^* + c_N v_N^* + c_S v_S^* + \Delta x (p_N^* - p_S^*) + d \end{eqnarray}$

です。この2組の式をそれぞれ引き算してやると

$\displaystyle \begin{eqnarray} a_P u_P' &=& a_W u_W' + a_E u_E' + a_N u_N' + a_S u_S' + \Delta y (p_E' - p_W') \\ c_P v_P' &=& c_W v_W' + c_E v_E' + c_N v_N' + c_S v_S' + \Delta x (p_N' - p_S') \end{eqnarray}$

のように補正量の式が得られます。次の手順がSIMPLE法の特徴です。上式の隣接点に関する項（W、E、N、Sの項）を省略してしまいます！すなわち

$\displaystyle \begin{eqnarray} u_P' &=& \frac{\Delta y}{a_P} (p_E' - p_W') \\ v_P' &=& \frac{\Delta x}{c_P} (p_N' - p_S') \end{eqnarray}$

としてしまいます。何故こんなことが可能なのでしょうか？補正量は値が収束すれば0になるので項を落としても大丈夫なのです！上式に仮の値を足すと真の値が求まります。すなわち

$\displaystyle \begin{eqnarray} \bar{u}_P &=& u^* + \frac{\Delta y}{a_P} (p_E' - p_W') \\ \bar{v}_P &=& v^* + \frac{\Delta x}{c_P} (p_N' - p_S') \end{eqnarray}$

です。これが流速の補正式になります。本来ならば左辺は $u^{n+1}、v^{n+1}$ となりますが、隣接項を落としているので $\bar{u}、\bar{v}$ としておきます。このようにして求めた値を連続式に代入します。すると

$\displaystyle \begin{eqnarray} g_P p_P' &=& g_W p_W' + g_E p_E' + g_N p_N' + g_S p_S' + h \end{eqnarray}$

と圧力に関する補正量の式となります。 $g$ は係数で $u^*、v^*$ より計算可能です（ぐちゃぐちゃするので別の記事で具体的な式を紹介します）。この式を解くと圧力の補正値が求まります。

計算手順は以下のように1～3の繰り返しとなります。

1. nステップでの圧力 $\boldsymbol p^n \boldsymbol{(p^*)}$ を用いて仮の流速 $\boldsymbol{(u^*, v^*)}$ を得る
2. 仮の流速 $\boldsymbol{(u^*, v^*)}$ を用いて圧力の補正量の式を解き、圧力の補正量 $\boldsymbol{p'}$ を求める
3. 圧力の補正量 $\boldsymbol{p'}$ を用いて、流速を補正（流速の補正式）し、圧力を補正（圧力の補正式）し、 $\boldsymbol{(\bar{u}, \bar{v}, \bar{p})}$ を求める
4. 補正量が十分に小さければ $\boldsymbol{(u^{n+1} = \bar{u}, v^{n+1} = \bar{v}, p^{n+1} = \bar{p})}$ として次のタイムステップの値を得る、補正量が十分に小さくなければ $\boldsymbol{p^* = \bar{p}, u^* = \bar{u}, v^* = \bar{v})}$ として1に戻り反復を繰り返す