非圧縮性流体の計算アルゴリズム-導入とMAC法 フラクショナルステップ法の導入とMAC法

非圧縮性流体

非圧縮性流体力学の方程式はテンソル表記で下のように表せます。

$$
\begin{align}
&質量保存則 :\hspace{ 30pt } \frac{ \partial u_i }{ \partial x_i } = 0 \tag{1} \label{eq:scale-factor-1} \\\\
&運動量保存則:\hspace{ 24pt } \frac{ \partial u_i }{ \partial t} + \frac{ \partial u_i u_j }{ \partial x_j} = -\frac{1}{\rho} \frac{\partial P}{\partial x_i} + \nu \frac{\partial^2 u_i}{\partial x_j^2} \tag{2} \label{eq:scale-factor-2}
\end{align}
$$
 
ベクトル表記だとこうでした。
 
 
$$
\begin{align}
&質量保存則 :\hspace{ 30pt } \nabla \cdot \mathbf{V} = 0  \\\\
&運動量保存則:\hspace{ 24pt } \frac{ \partial \mathbf{V}}{ \partial t} + (\mathbf{V} \cdot \nabla) \mathbf{V} = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \mathbf{V} \qquad \\\\
\end{align}
$$

 

だいぶ話が戻りましたね。

それでは実際に解法へ移っていきます。僕らは前回までで微分方程式を離散化するという方法を学習してきました。ならば、いかに複雑なNavier-Stocks方程式も恐るるに足りません。各項を離散化してやればいいだけです。よし、解ける、、、となったところで一つ壁に当たります。

圧力項どうしよう、、、。これらは4方程式で速度3つ、圧力1つの4つの未知数なので解けるはずなのですが、質量保存則をうまく使って圧力を処理しなければなりません。

そこで計算アルゴリズムを考えていきます。

まず、圧力だけに注目して、それ以外をFとして書いて1次Adoms-Bashforth法で離散化してみます。

$$
\begin{align}
&\frac{ \partial u_i }{ \partial x_i } = 0  \\\\
&\frac{ \partial u_i }{ \partial t} = F -\frac{1}{\rho} \frac{\partial P}{\partial x_i} \\\\
\rightarrow \ &u_i^{n+1} – u_i^n=  \Delta t F^n – \frac{\Delta t}{\rho} \frac{\partial P^{n \ or \ n+1}}{\partial x_i} \\\\
\tag{3} \label{eq:scale-factor-3}
\end{align}
$$

ここでの添字のiは、空間の離散化した添字ではなく、テンソルのiなことに注意してください。少しややこしいですね。まああまり問題になりませんが、この意味でベクトルを使った方がいい時もあります。

ここで悩むのはまず、離散化したのはいいものの、\(P^{n+1}\)とすればいいのか、\(P^n\)とすればいいのか、また\(u\)の方程式の中で\(P\)をどのように扱ったらいいのかというところでしょう。

まず、\(P^n\)だったらどうなるか考えます。方程式がどうなるか書いてみます。

$$
\begin{align}
&u_i^{n+1} – u_i^n=  \Delta t F^n – \frac{\Delta t}{\rho} \frac{\partial P^{n}}{\partial x_i} \\\\
\tag{4} \label{eq:scale-factor-4}
\end{align}
$$

\(P^{n+1}\)が一生求まる気がしません。\(P^{n+1}\)がどこにも存在せず、一生次のステップを求められないまま終わります。よって、\(P^{n+1}\)にするしかありません。この意味で、圧力項はどの解法でも陰解法になります。

項に対して陽解法とか陰解法とか変えていいの?って思うかもしれませんが、普通に大丈夫です。陽解法の方を先に解いてあげて、陰解法の項だけでそれ使って後から解いたりなど色々その辺は柔軟にできます。だって項同士は線形ですから。

なので、とりあえず\(P^{n+1}\)を採用します。

$$
\begin{align}
&\frac{u_i^{n+1} – u_i^n}{\Delta t}=  F^n – \frac{1}{\rho} \frac{\partial P^{n+1}}{\partial x_i} \\\\
\tag{5} \label{eq:scale-factor-5}
\end{align}
$$

そして、2段分離をします。

$$
\begin{align}
&\frac{u_i^{n+1} – \tilde{u_i} }{\Delta t}=  – \frac{1}{\rho} \frac{\partial P^{n+1}}{\partial x_i} \tag{6} \label{eq:scale-factor-6} \\\\
&\frac{\tilde{u_i} – u_i^n}{\Delta t}=  F^n
\tag{7} \label{eq:scale-factor-7}
\end{align}
$$

こんなことしていいのか、、、と思うかもしれませんが、足して見ると元の式に一致すると思います。ここで、下の式ならば\(\tilde{u_i}\)は算出できます。\( \tilde{u_i}\)を中間速度と呼んだりします。

そして、1個目の圧力項の式を\(x_i\)で微分します。

$$
\begin{align}
&\frac{1}{\Delta t} \Bigl( \frac{\partial u_i^{n+1} }{\partial x_i }- \frac{\partial \tilde{u_i}}{\partial x_i} \Bigr) =  – \frac{1}{\rho} \frac{\partial^2 P^{n+1}}{\partial x_i^2} \\\\
\tag{8} \label{eq:scale-factor-8}
\end{align}
$$

ここですごいことに気づきます。連続条件(質量保存則)が隠れてませんか?

$$
\begin{align}
\frac{\partial u_i^{n+1} }{\partial x_i } = 0
\end{align}
$$

これはn+1時点での連続条件のはずです。つまり、この項を0にしたら、n+1時点の連続条件を満たす\(P^{n+1} \)が求まることになります。

そして式を整理すると、

$$
\begin{align}
& \frac{\partial^2 P^{n+1}}{\partial x_i^2} = \frac{\rho}{\Delta t}  \frac{\partial \tilde{u_i}}{\partial x_i}   \\\\
\tag{9} \label{eq:scale-factor-9}
\end{align}
$$

こうなります!!

これで、Pを\(\tilde{u_i}\)から求めることができるような式ができました。これはよく見るとPoisson方程式になってます。電磁気でよく見るあれですね。

なので、圧力はPoisson方程式を解いて求めます。ここでは、時間ステップの逆行列ではなくて、単純な逆行列(連立方程式)を解くことになります。

なので、このアルゴリズムではまず、\eqref{eq:scale-factor-6}式から\(\tilde{u_i}\)を求め、\eqref{eq:scale-factor-9}式から\(P^{n+1} \)を求め、 \eqref{eq:scale-factor-7}式から、\(u_i^{n+1} \)を求めます。

とりあえず、これで求まりそうなアルゴリズムを一つ構築できました。このような方法をフラクショナルステップ法と言います。

なかなかに優れたアルゴリズムなのですが、このような分離の際には分離誤差が生じてしまったり、上記のように空間的に離散化されていない方程式を基にフラクショナル・ステップ法を構成すると中間速度\(\tilde{u_i}\)の境界条件の設定が曖昧となってしまう、といった問題が生じたりします。今はこの問題には深入りしません。また、本来の形とは若干違うことや一般にフラクショナルステップ法というと、そのまま\(F^n\)にするのではなく、移流項を陽解法、粘性項を陰解法でとく方法を指すので後々もっと詳しくやっていきます。

ここからはいろんな解法をみていきます。まず差分系の解法を大別すると以下のようになっています。

計算アルゴリズム移流項粘性項圧力項
MAC法 系統陽解法陽解法(陰解法)
フラクショナル・ステップ法 系統陽解法陰解法(陰解法)
------------------陰解法陽解法(陰解法)
SIMPLE法 系統陰解法陰解法(陰解法)

MAC法

MAC法は最も知られたアルゴリズムです。やり方もシンプルで分かりやすいです。それでは導出していきます。

まず、Navier-Stocks方程式から移流項と粘性項を陽的に時間離散化して、連続式を新たな時刻(\t^{n+1}\)において満足させると、次式になります。

$$
\begin{align}
&\frac{ \partial u_i^{n+1} }{ \partial x_i^{n+1} } = 0 \tag{10} \label{eq:scale-factor-10} \\\\
&\frac{ u_i^{n+1} – u_i^{n}}{ \Delta t} + \frac{ \partial u_i^n u_j^n }{ \partial x_j} = -\frac{1}{\rho} \frac{\partial P^{n+1}}{\partial x_i} + \nu \frac{\partial^2 u_i^n}{\partial x_j^2} \tag{11} \label{eq:scale-factor-11}
\end{align}
$$

やはり、iやjの添字が空間の離散化した添字ではなく、方程式のテンソルの添字なことに注意です。

これの\eqref{eq:scale-factor-11}式の両辺の発散をとります。テンソルの発散をとるにはベクトルになっている添字で微分すればいいです。そうすれば、縮約で内積みたいな形になり、ナブラとの内積\(\nabla \cdot\) みたいになります。

つまり、両辺を\( \frac{ \partial }{ \partial x_i} \)で微分します。

$$
\begin{align}
&\frac{1}{\Delta t} \frac{\partial u_i^{n+1}}{\partial x_i} – \frac{1}{\Delta t} \frac{\partial u_i^n}{\partial x_i} +  \frac{ \partial^2 u_i^n u_j^n }{ \partial x_i \partial x_j} = -\frac{1}{\rho} \frac{\partial^2 P^{n+1}}{\partial x_i^2} + \nu \frac{\partial^3 u_i^n}{\partial x_i \partial x_j^2} \\\\
\tag{12} \label{eq:scale-factor-12}
\end{align}
$$

ここで、連続条件を適用し、少し整理します。

$$
\begin{align}
& \frac{1}{\rho} \frac{\partial^2 P^{n+1}}{\partial x_i^2}  =  \frac{1}{\Delta t} \frac{\partial u_i^n}{\partial x_i} –  \frac{ \partial^2 u_i^n u_j^n }{ \partial x_i \partial x_j} + \nu \frac{\partial^3 u_i^n}{\partial x_i \partial x_j^2} \\\\
& \frac{\partial^2 P^{n+1}}{\partial x_i^2}  =  \frac{\rho}{\Delta t} \frac{\partial u_i^n}{\partial x_i} +  \rho \frac{\partial }{\partial x_i} \bigl( \nu \frac{\partial^2 u_i^n}{ \partial x_j^2} – \frac{\partial u_i^n u_j^n }{\partial x_j}  \bigr) \\\\
\tag{13} \label{eq:scale-factor-13}
\end{align}
$$

これで圧力のPoisson方程式が完成しました。MAC法ではまずこの式を解いて、連続式\eqref{eq:scale-factor-10}を満たす、\(P^{n+1}\)を\eqref{eq:scale-factor-13}式から求めたのちに求めた圧力を使って、\eqref{eq:scale-factor-11} 式を使って次のステップの\(u^{n+1}\)を求めます。

ベクトルでも書いてみます。次のようになります。まず、移流項と拡散項を陽的に離散化します。

$$
\begin{align}
& \nabla \cdot \mathbf{V^{n+1}} = 0 \tag{14} \label{eq:scale-factor-14} \\\\
&\frac{ \mathbf{V^{n+1}} – \mathbf{V^{n}} }{ \Delta t} + (\mathbf{V^n} \cdot \nabla) \mathbf{V^n} = -\frac{1}{\rho} \nabla p^{n+1} + \nu \nabla^2 \mathbf{V^n} \qquad \tag{15} \label{eq:scale-factor-15}
\end{align}
$$
 
そして、両辺の発散をとります。つまり、両辺に\(\nabla \cdot\)を掛けます。そして連続条件を適用して、整理します。

 

$$
\begin{align}
& \nabla \cdot \nabla p^{n+1} = \frac{1}{ \Delta t} \nabla \cdot \mathbf{V^{n}} + \rho \nabla \cdot \bigl( \nu \nabla^2 \mathbf{V^n} – ( \mathbf{V^n} \cdot \nabla) \mathbf{V^n} \bigr)
\tag{16} \label{eq:scale-factor-16}
\end{align}
$$

だいたい慣れましたでしょうか?どの解法でも、連続条件を使って、圧力のPoisson方程式を作ることは変わりません。

そして、もう一つだけMAC法といえば、、、みたいなことがありまして、一般にMAC法はスタガード格子というものを採用します。もちろん、他の格子を採用することも可能ですが、スタガード格子と相性がいいので、セットみたいになってます。

格子というのは空間の離散化の仕方のことで、流体には圧力と速度の2成分あるため、このあたりをどう扱うかで格子にいろんな種類が存在します。レギュラー格子、スタガード格子、コロケート格子といった3種類が代表的なものです。少し長くなってしまうので、格子の話は次回に書きます。

 



[
前の記事へ]  [非圧縮性流体の目次へ]  [次の記事へ]