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

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

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

MAC法(Maker And Cell method)を改良した、SMAC法(Simplified MAC法)によるNavier-Stokes方程式の離散化について説明していきます。この手法は流速と圧力を直接未知数として計算する手法です。MAC法ではPoisson方程式の右辺を計算するのに負荷がかかっていましたが、それを簡略化したのがSMACです。なので、基本的にMACとSMACは同じもののようです。ただし、計算時間はSMACのほうが短いです。


以下の内容は『流れ解析のための有限要素法』のpp.160-163を参考にしています。

流れ解析のための有限要素法入門

流れ解析のための有限要素法入門


簡単のため二次元直交座標系で考えます。まず保存型の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)^2}{\partial x} + \frac{\partial (u^nv^n)}{\partial y} &=& -\frac{\partial p^{n+1}}{\partial x} + \frac{1}{Re} \left( \frac{\partial^2 u^n}{\partial x^2} + \frac{\partial^2 u^n}{\partial y^2}\right) \\
\frac{v^{n+1}-v^n}{\Delta t} + \frac{\partial (u^nv^n)}{\partial x} + \frac{\partial (v^n)^2}{\partial y} &=& -\frac{\partial p^{n+1}}{\partial y} + \frac{1}{Re} \left( \frac{\partial^2 v^n}{\partial x^2} + \frac{\partial^2 v^n}{\partial y^2}\right) \\
\frac{\partial u^{n+1}}{\partial x} + \frac{\partial v^{n+1}}{\partial y} &=& 0
\end{eqnarray}

圧力を陰的に扱う(n+1ステップの値を用いる)ところがMACとの違いです。このままではn+1ステップの値があり解けないので、以下のようにnステップの値を用いて中間速度(予測子) u^* v^* を計算します。

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

次に対応する運動方程式をそれぞれ引き算すると

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

を得ます。これが修正子にあたります。ただし

 \displaystyle
\begin{eqnarray} 
p^{n+1} - p^n = \delta p
\end{eqnarray}

としています。これが圧力の補正量です。ここから  u^{n+1} v^{n+1} を連続式に代入します(n+1ステップにおいて連続式を満たすとする)。

 \displaystyle
\begin{eqnarray} 
\frac{\partial^2}{\partial x^2} \left(\delta p \right) + \frac{\partial^2}{\partial y^2} \left(\delta p \right) = \frac{1}{\Delta t} \left( \frac{\partial u^{*}}{\partial x} + \frac{\partial v^{*}}{\partial y} \right)
\end{eqnarray}

これがSMACにおける圧力の補正量に関するポワソン方程式です。圧力そのものではなく、圧力の補正量になっている点がMACとの違いです。

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

1. nステップでの流速  \boldsymbol{(u^n, v^n)} と 圧力  \boldsymbol p^n を用いて中間流速(予測子の式) \boldsymbol{(u^*, v^*)} を得る
2. 中間流速  \boldsymbol{(u^*, v^*)} を用いて圧力の補正量のPoisson方程式を解き、圧力の補正量  \boldsymbol{\delta p} を求める
3. 圧力の補正量  \boldsymbol{\delta p} を用いて、圧力と流速を補正(修正子の式)し、次のタイムステップの値  \boldsymbol{(u^{n+1}, v^{n+1}, p^{n+1})} を得る

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


SMAC法+中心差分による数値計算例です。

SMAC法+風上差分