数学とか語学とか楽しいよね

フランス語、ドイツ語、ロシア語、アラビア語、オランダ語、英語、スペイン語、ラテン語とか数学とか数値計算(有限要素法、有限体積法、差分法、格子ボルツマン法、数理最適化、C++コード付き)とか勉強したことをまとめます。右のカテゴリーから興味のある記事を探してください。最近はクラシックの名演も紹介しています。noteにも書いています。https://note.mu/kumag_min

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

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

数値流体工学

数値流体工学


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


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

流体解析の基礎

流体解析の基礎


簡単のため二次元直交座標系で考えます。まず保存型の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に戻り反復を繰り返す


実は圧力は不足緩和しないといけませんし、流速の更新式にも緩和係数を入れるようです。これも別記事で説明します。

反復一回につき、仮の流速の計算と圧力の補正量の計算の際に、合計二回連立方程式を解く必要があります。

空間方向の離散化には、保存型をそのまま用いれば有限体積法を、非保存型を用いれば有限要素法や差分法を用いることができます。普通は有限体積法が用いられています。