はじめに
流体の方程式の離散化には差分法や有限体積法が用いられます。特に有限体積法は多くのソルバへ実装されており、シミュレーションを行う上でその理論的背景を知っておくことは重要です。そこで今回から有限体積法の基礎についてまとめていきます。本記事では差分法と有限体積法の特徴を比較し、実際に１次元非粘性バーカース方程式について有限体積法で離散化してみたいと思います。
有限体積法と有限差分法の違い
有限差分法は連続体を有限の格子点に分割し、テイラー展開により微分項を差分化することで離散的な方程式系を構築する手法です。離散化の基本原理としてテイラー展開のみを用いるため精度が明確、高次精度化が容易という利点があります。
一方、有限体積法は連続体を有限のセル（体積）に分割し、セルごとの物理量の保存則から離散的な方程式系を構築する手法です。物理的な保存則（質量保存、エネルギー保存、運動量保存など）が明確に存在し、それが物理量の移動（輸送）と関連しているような問題、すなわち流体力学や熱伝導の問題を解くのに適しています。
有限体積法は非構造格子への拡張が容易という利点がありますが、高次精度化が煩雑化するという欠点を有しています。
1次元非粘性バーカース方程式を有限体積法で離散化する
有限体積法の具体的な手順について確認します。以下の記事では非粘性バーカース方程式について差分法による離散化を行いました。今回は同様の方程式について有限体積法で離散化を行ってみます。
非粘性バーカース方程式は時間・空間の関数である物理量\(q(x,t)\)について、次のように書かれます。
\frac{\partial q}{\partial t}+q\frac{\partial q}{\partial x}= 0\tag{1}
\end{align}
上式について以下のようにセル中心・セル界面が配置された１次元セル上での離散化を考えます。
有限体積法では離散変数がセル内での空間平均値となるため、対象となる偏微分方程式を積分し、単一セル空間・微小時間に関して平均化します。
\frac{1}{\Delta x \Delta t}\int_{t}^{t+\Delta t}\int_{i-\frac{1}{2}}^{i+\frac{1}{2}}
\left(\frac{\partial q}{\partial t}+q\frac{\partial q}{\partial x}\right)dxdt=0
\tag{2}
\end{align}
上式について積分を実行し離散化していきます。
まず左辺第一項から考えます。空間積分を微分の中に入れて以下のように書き換えます。
&\frac{1}{\Delta x \Delta t}\int_{t}^{t+\Delta t}\int_{i-\frac{1}{2}}^{i+\frac{1}{2}}
\frac{\partial q}{\partial t}dxdt
=
\frac{1}{\Delta t}\int_{t}^{t+\Delta t}\left[\frac{\partial }{\partial t}\left(\frac{1}{\Delta x}\int_{i-\frac{1}{2}}^{i+\frac{1}{2}}qdx\right)\right]dt \tag{3}
\end{align}
有限体積法では離散変数\(\bar{q}_{i}(t)\)を以下のようなセル平均値として定義します。
\bar{q}_{i}(t)=\frac{1}{\Delta x}\int_{i-\frac{1}{2}}^{i+\frac{1}{2}}q(x,t)dx \tag{4}
\end{align}
上式にて空間積分を実行したことで\(q(x,t)\)の位置が決定し、時間のみの関数として\(\bar{q}_{i}(t)\)となることに注意しましょう。上式を式\((3)\)に代入すると
\frac{1}{\Delta t}\int_{t}^{t+\Delta t}\left[\frac{\partial }{\partial t}\left(\frac{1}{\Delta x}\int_{i-\frac{1}{2}}^{i+\frac{1}{2}}qdx\right)\right]dt=\frac{1}{\Delta t}\int_{t}^{t+\Delta t}\frac{\partial \bar{q}_{i}}{\partial t}dt \tag{5}
\end{align}
次に時間積分を実行します。
\frac{1}{\Delta t}\int_{t}^{t+\Delta t}\frac{\partial \bar{q}_{i}}{\partial t}dt=\frac{1}{\Delta t}\left[ \bar{q}_{i}\right]_ {\bar{q}_{i}(t)}^{\bar{q}_{i}(t+\Delta t)}=\frac{\bar{q}_{i}(t+\Delta t)-\bar{q}_{i}(t)}{\Delta t}\tag{6}
\end{align}
\(\bar{q}_{i}(t)=\bar{q}_{i}^{n}\)、\(\bar{q}_{i}(t+\Delta t)=\bar{q}_{i}^{n+1}\)と書くことにすれば、最終的に式\((2)\)の左辺第一項は以下のように書けます。
(式(2)の左辺第一項)=\frac{\bar{q}_{i}^{n+1}-\bar{q}_{i}^{n}}{\Delta t}\tag{7}
\end{align}
色々と変形してきましたが、上式は\(q\)がセル平均\(\bar{q}\)になっているものの、元の方程式\((1)\)の第一項を前進オイラー法で離散化した場合と同等の形になります。
引き続き、以下の式\((2)\)の第二項を考えます。まずは空間積分を実行します。空間積分を実行するため係数\(q\)を微分の中に入れてしまいます。
&\frac{1}{\Delta x \Delta t}\int_{t}^{t+\Delta t}\int_{i-\frac{1}{2}}^{i+\frac{1}{2}}
\left(q\frac{\partial q}{\partial x}\right)dxdt \\
&=\frac{1}{\Delta x \Delta t}\int_{t}^{t+\Delta t}\int_{i-\frac{1}{2}}^{i+\frac{1}{2}}
\frac{\partial }{\partial x}\left(\frac{q^2}{2}\right)dxdt \\
&=\frac{1}{\Delta x \Delta t}\int_{t}^{t+\Delta t}\left[\frac{q^2}{2}\right]_{i-\frac{1}{2}}^{i+\frac{1}{2}}dt \\
&=\frac{1}{\Delta x \Delta t}\int_{t}^{t+\Delta t}\frac{1}{2}\left(q^2_{i+\frac{1}{2}}-q^2_{i-\frac{1}{2}}\right)dt\tag{8} \\
\end{align}
空間積分の結果がセル界面での物理量\(q_{i+\frac{1}{2}},q_{i-\frac{1}{2}}\)で表現されるのがポイントです。このように有限体積法ではセル界面での物理量の流入出を考えることで、後々非構造格子への拡張が可能となります。
ここで界面での物理量の流入・流出という観点で方程式を解釈するため、以下の「数値流束」を導入しておきます。
\tilde{f}_{i+\frac{1}{2}}＝\frac{q^2_{i+\frac{1}{2}}}{2} \\
\tilde{f}_{i-\frac{1}{2}}＝\frac{q^2_{i-\frac{1}{2}}}{2} \\
\tag{9}
\end{align}
こちらの記事でも述べたとおり数値流束は「セル界面を単位時間当たりに通過する物理量\(q\)の総量」を表しています。これを用いると、式\((8)\)は以下のように書くことが出来ます。
\frac{1}{\Delta x \Delta t}\int_{t}^{t+\Delta t}\frac{1}{2}\left(q^2_{i+\frac{1}{2}}-q^2_{i-\frac{1}{2}}\right)dt=\frac{1}{\Delta x \Delta t}\int_{t}^{t+\Delta t}\left(\tilde{f}_{i+\frac{1}{2}}-\tilde{f}_{i-\frac{1}{2}}\right)dt
\tag{10} \\
\end{align}
方程式を積分型に書き換えたことで、差分法による離散化の際に導入した「数値流束」の概念が自然に現れるのが面白いです。言い換えると有限体積法を差分法的に解釈した際に現れるのが保存型表示や数値流束であり、流れのような輸送方程式の離散化は本来は保存則を考慮できる有限体積法の利用が適切なのでしょう。
次に時間積分を実行したいのですが、\(\tilde{f}_{i+\frac{1}{2}},\tilde{f}_{i-\frac{1}{2}}\)が具体的にどんな時間の関数なのかが分かりません。
そこで以下のように\(\tilde{f}_{i+\frac{1}{2}},\tilde{f}_{i-\frac{1}{2}}\)は\(t\sim t+\Delta t\)の間は一定とし、現タイムステップ\(n\)での値\(\tilde{f}_{i+\frac{1}{2}}^n, \tilde{f}_{i-\frac{1}{2}}^n\)で近似することとします。
\tilde{f}_{i+\frac{1}{2}}(t)=\tilde{f}_{i+\frac{1}{2}}^{n}\\
\tilde{f}_{i-\frac{1}{2}}(t)=\tilde{f}_{i-\frac{1}{2}}^{n}\\
\tag{11}
\end{align}
これは\(\tilde{f}_{i+\frac{1}{2}},\tilde{f}_{i-\frac{1}{2}}\)を陽的に離散化（＝現在のタイムステップの情報のみで表される形で離散化）したことに相当します。
もちろん以下のように次のステップの情報を用いて陰的に離散化することも可能です。
\tilde{f}_{i+\frac{1}{2}}(t)=\tilde{f}_{i+\frac{1}{2}}^{n+1}\\
\tilde{f}_{i-\frac{1}{2}}(t)=\tilde{f}_{i-\frac{1}{2}}^{n+1}\\
\tag{12}
\end{align}
ただし、この場合は陰解法になるので最終的に解を求める際には連立一次方程式を解く必要が出てきます。また現ステップの情報と次ステップの情報を両方用いるクランクニコルソン法と呼ばれる手法もあります。
なお上記のいずれの離散化に関しても\(\tilde{f}_{i+\frac{1}{2}},\tilde{f}_{i-\frac{1}{2}}\)は定数として扱われるので、時間に関して１次精度の離散化といえます。このことは以下のように\(\tilde{f}_{i+\frac{1}{2}},\tilde{f}_{i-\frac{1}{2}}\)をテイラー展開してみれば明らかです。
\tilde{f}_{i+\frac{1}{2}}(t+\Delta t)&=\tilde{f}_{i+\frac{1}{2}}(t)+\frac{\partial \tilde{f}_{i+\frac{1}{2}} }{\partial t}\Delta t+\frac{1}{2!}\frac{\partial ^2 \tilde{f}_{i+\frac{1}{2}}}{\partial t^2}(\Delta t)^2 \cdots\\
&=\tilde{f}_{i+\frac{1}{2}}(t)+O(\Delta t)\\
\tag{13}
\end{align}
上式から分かるように、より高次の離散化を行うには\(\tilde{f}_{i+\frac{1}{2}}\)の時間微分項を何らかの形で表現する必要があります。この辺は後々出てくるLax-Wendroff法等が関係するところです。
今回はひとまず１次精度の陽的な離散化式\((11)\)を式\((10)\)へ適用します。\(\tilde{f}_{i+\frac{1}{2}},\tilde{f}_{i-\frac{1}{2}}\)は定数\(\tilde{f}_{i+\frac{1}{2}}^n,\tilde{f}_{i-\frac{1}{2}}^n\)として近似されるので積分の外へ出すことが出来ます。
&\frac{1}{\Delta x \Delta t}\int_{t}^{t+\Delta t}\left(\tilde{f}_{i+\frac{1}{2}}-\tilde{f}_{i-\frac{1}{2}}\right)dt\\
&=\frac{1}{\Delta x \Delta t}\left(\tilde{f}^n_{i+\frac{1}{2}}-\tilde{f}^n_{i-\frac{1}{2}}\right)\int_{t}^{t+\Delta t}dt\\
&=\frac{\tilde{f}^n_{i+\frac{1}{2}}-\tilde{f}^n_{i-\frac{1}{2}}}{\Delta x }\tag{14}\\
\end{align}
では最後に残った\(\tilde{f}_{i+\frac{1}{2}}^n,\tilde{f}_{i-\frac{1}{2}}^n\)は具体的にどのように表現すればよいのでしょうか？自然に考えれば有限体積法の離散変数であるセル平均値\(\bar{q}_{i}\)用いて何らかの近似で表現するのが妥当と思います。
このセル平均値\(\bar{q}_{i}\)から流束\(\tilde{f}_{i+\frac{1}{2}}^n,\tilde{f}_{i-\frac{1}{2}}^n\)を計算する式はリーマンソルバと呼ばれており、「リーマンソルバをどう構築するか（数値流束を如何に表現するか）」は有限体積法の最も重要な問題の一つです。
正確な解を得るためにこれまで様々なリーマンソルバが提案されており、高次精度解法ともかかわりが深い部分です。ただこの辺を書き出すとまた長くなるので、今回は深追いしないでおきたいと思います。
ともあれ、最終的に式\((2)\)の左辺第二項は以下のようになりました。
(式(2)の左辺第二項)=\frac{\tilde{f}^n_{i+\frac{1}{2}}-\tilde{f}^n_{i-\frac{1}{2}}}{\Delta x }\tag{15}\\
\end{align}
以上、式\((7),(15)\)を式\((2)\)に代入すると一次元非粘性バーカース方程式に関する有限体積法による離散化式が得られます。
&\frac{\bar{q}_{i}^{n+1}-\bar{q}_{i}^{n}}{\Delta t}+\frac{\tilde{f}^n_{i+\frac{1}{2}}-\tilde{f}^n_{i-\frac{1}{2}}}{\Delta x }=0\\
&\rightarrow
\bar{q}_{i}^{n+1}=\bar{q}_{i}^{n}+\frac{\Delta t}{\Delta x }\left(\tilde{f}^n_{i+\frac{1}{2}}-\tilde{f}^n_{i-\frac{1}{2}}\right)\tag{16}\\
\end{align}
上式は離散変数がセル平均値\(\bar{q}\)になっているものの、保存型表示での有限差分法による離散化式と形が一致します。（詳しくはこちらの記事を参照。）一般に１次元系の場合、有限差分法と有限体積法の最終的な離散化式は一致することが知られています。
おわりに
今回は有限体積法と差分法を比較し、1次元非粘性バーカーズ方程式について有限体積法で離散化を行いました。差分法でわざわざ導入した数値流束が自然と出てくるのは面白いと思います。高次精度化や非構造格子への対応などまだまだ奥が深そうなので、じっくり勉強していきたいと思います。では今回はこの辺で。
参考文献
[1] 肖 鋒, “数値流体解析の基礎- Visual C++とgnuplotによる圧縮性・非圧縮性流体解析 -“, コロナ社, 2020
[2] 土木学会 応用力学委員会 計算力学小委員会 編, “いまさら聞けない計算力学の常識”, 丸善, 2008年
[3] 谷口伸行, “有限体積法による流れ解析スキームの再評価“, 生産研究, 48号, 1996
コメント