大学院生のための数値流体力学①：Navier-Stokes方程式

Navier-Stokes方程式は，時間項（非定常項），対流項（移流項），圧力項，粘性項（拡散項），外力項の5つに分けられる

\begin{equation}
\frac{\partial u_{i}}{\partial t}+\frac{\partial \left(u_{i}u_{j}\right)}{\partial x_{j}}=-\frac{1}{\rho}\frac{\partial p}{\partial x_{i}}+\nu\frac{\partial}{\partial x_{j}}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)+f_{i}
\tag{2}
\end{equation}

各項の名前と，どのようなはたらきを持っているかはなんとなくでも理解しておくといろいろと便利

時間項

\begin{equation}
\frac{\partial u_{i}}{\partial t}
\end{equation}

時間の経過に伴う流れの変化を表す項

定常計算では0となる

対流項

\begin{equation}
\frac{\partial \left(u_{i}u_{j}\right)}{\partial x_{j}}
\end{equation}

流れが流れ自身の影響で運ばれていく効果（移流）を表す項

速度の2乗が現れる非線形項であるため，離散化の際に慎重に取り扱う必要がある

圧力項

\begin{equation}
-\frac{1}{\rho}\frac{\partial p}{\partial x_{i}}
\end{equation}

圧力差が流れに及ぼす影響を表す項

粘性項

\begin{equation}
\nu\frac{\partial}{\partial x_{j}}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)
\end{equation}

粘性によって流れの分布が平滑になる効果（拡散）を表す項

数値解析の計算を安定化させるはたらきを持つ

外力項

\begin{equation}
f_{i}
\end{equation}

外部からの力を表す項

外部からの力がない場合は0となる

格子解像度について

一切のモデル化を伴わずにNavier-Stokes方程式を解くためには，格子のサイズを流れ場に存在する最小スケールの渦と同程度にする必要がある

この「最小スケールの渦」がとにかく小さいため，現在の最高性能のスーパーコンピューターをもってしてもRe数の低い簡単な流れ場（チャネル乱流など）の計算を行うので精いっぱいである

実際にどのくらいの計算コストが必要なのか，鳥人間コンテストに出てくるような人力飛行機（速度\(U\)，代表長さ\(L\)）で計算してみる

定常飛行をしている機体を考えると，機体による乱れエネルギーの供給（翼や胴体によって空気が乱されること）と空気の粘性による乱れエネルギーの散逸（最小スケールまで分裂した渦が最終的に熱エネルギーに変換されること）は同程度だと考えられる

機体による単位時間当たりの乱れエネルギーの供給\(\varepsilon ~[\mathrm{m^2/s^3}]\)はエネルギー\(U^{2} ~[\mathrm{m^2/s^2}]\)を時間\(L/U ~[\mathrm{s}]\)で割ることによって求められる

\begin{equation}
\varepsilon=\mathcal{O}\left[\frac{U^{3}}{L}\right]
\end{equation}

一方，空気の粘性による単位時間当たりの乱れエネルギーの散逸\(\varepsilon ~[\mathrm{m^2/s^3}]\)は，動粘性係数\(\nu ~[\mathrm{m^2/s}]\)と渦の最小スケール\(\eta ~[\mathrm{m}]\)を用いて次元解析的に求められる

\begin{equation}
\varepsilon=\mathcal{O}\left[\frac{\nu^{3}}{\eta^{4}}\right]
\end{equation}

供給と散逸が等しいとすると，代表長さ\(L\)と渦の最小スケール\(\eta\)の関係が導かれる

\begin{equation}
\frac{L}{\eta}=\mathcal{O}[Re^{3/4}]
\end{equation}

高Re数の流れ場では，乱れが供給される代表長さに比べて乱れが散逸するスケールは非常に小さく，その比はおおむねRe数の3/4乗に比例することがわかる

一般的な人力飛行機では，Re数は\(10^{5}\)程度，代表長さ（主翼のコード長）は\(1 ~[\mathrm{m}]\)程度なので，人力飛行機まわりの流れに存在する最小スケールの渦（乱れ）は次のように求まる

\begin{equation}
\eta=\mathcal{O}\left[\frac{L}{Re^{3/4}}\right]=\mathcal{O}\left[\frac{1 ~[\mathrm{m}]}{\left(10^{5}\right)^{3/4}}\right]=\mathcal{O}\left[10^{-4} ~[\mathrm{m}]\right]
\simeq0.1 ~[\mathrm{mm}]
\end{equation}

人力飛行機まわりの流れをDNSで計算するためには，格子サイズを最低でも0.1 [mm]以下にする必要があり，全長10 [m]，全幅30 [m]，全高5 [m]の計算領域を確保した場合，格子点数は\(1.5\times 10^{15}\)にも及ぶ

プログラムをかじったことがある人なら，1000兆×1000兆の連立一次方程式を数万回解く必要があるプログラムの大変さが理解できると思う

ちなみにこれまで実施された最大規模の計算は，2004年に行われた分割数\(1.74\times 10^{11}\)の平行平板間乱流（チャネル乱流）である

≫https://doi.org/10.1016/j.ijheatfluidflow.2004.02.010

旅客機のRe数が\(10^{8}\)程度であることを考えると，DNSの実用がまだまだ先であることが実感できる

初期条件について

仮に上記の計算コストを支払えるだけの圧倒的技術力と財力があったとしても，DNSの計算を行うための初期条件の用意が非常に難しい

例えばDNSを使って将来の天気を正確に予測するプログラムを作るとする

そのためには，同時刻における，全世界のありとあらゆる場所の，ありとあらゆるスケールの速度，圧力，温度などの情報を，寸分の狂いなく正確に測定する必要がある

少し考えただけでも困難なことがわかる

カオス理論について

超高性能のスーパーコンピューターと，全知全能の神によって初期条件を設定できたとしても，DNSの計算は難しい

数値計算を行うと，離散化に伴う誤差や，計算を行うことによる丸め誤差が必ず生じる

1回に生じる誤差は極めて小さいものでも，それが積もり積もって非常に大きな結果の違いが生まれることになる

この現象はバタフライ効果やカオス理論として知られ，三重振り子の計算などで確認することができる

≫バタフライ効果 - Wikipedia
≫カオス理論 - Wikipedia

先日のカオス振り子動画を見て極端な場合を試してみたくなったので、3重振り子で高精度物理演算してみた。初期の振り子の角度を10の100乗(グーゴル)分の1ずつだけずらしたものを100個同時に再生。完全に一致してるように見えるけど、1分20秒くらいから「それ」がやってくる。 pic.twitter.com/a7E2OASIDe

— So Takamoto (@tkmtSo) May 25, 2017

フィルター操作に関して，以下の仮定をおく

\begin{equation}
\overline{\frac{\partial f}{\partial x}}=\frac{\partial \overline{f}}{\partial x}, \quad
\overline{\frac{\partial f}{\partial t}}=\frac{\partial \overline{f}}{\partial t}
\end{equation}

本来フィルター操作は空間的にも時間的にも変化しうるものなので，この仮定は厳密には正しくない

もちろんこの過程に起因する誤差も生じる可能性があるが，現状はあまり重要視されていない

連続の式

これを踏まえたうえで，連続の式にフィルター操作を施すと次のようになる

\begin{align}
\overline{\frac{\partial u_{i}}{\partial x_{i}}}=0
\quad \longleftrightarrow \quad
\frac{\partial \overline{u}_{i}}{\partial x_{i}}=0
\end{align}

Navier-Stokes方程式

フィルターをかけたNavie-Stokes方程式は次のようになる

\begin{align}
\overline{\frac{\partial u_{i}}{\partial t}}
+\overline{\frac{\partial \left( u_{i}u_{j} \right)}{\partial x_{j}}}
=
-\frac{1}{\rho}\overline{\frac{\partial p}{\partial x_{i}}}
+\overline{\frac{\partial}{\partial x_{j}}\left\{ \nu\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)\right\} }
\end{align}

\begin{align}
\longleftrightarrow \quad
\frac{\partial \overline{u}_{i}}{\partial t}
+\frac{\partial \left(\overline{u_{i}u_{j}}\right)}{\partial x_{j}}
=
-\frac{1}{\rho}\frac{\partial \overline{p}}{\partial x_{i}}
+\frac{\partial}{\partial x_{j}}\left\{\nu\left(\frac{\partial \overline{u}_{i}}{\partial x_{j}}+\frac{\partial \overline{u}_{j}}{\partial x_{i}}\right)\right\}
\end{align}

ここで，両辺から\(\frac{\partial \left(\overline{u_{i}u_{j}}\right)}{\partial x_{j}}\)を引いて\(\frac{\partial \left(\overline{u}_{i}\overline{u}_{j}\right)}{\partial x_{j}}\)を足すと

\begin{align}
\frac{\partial \overline{u}_{i}}{\partial t}
+\frac{\partial \left(\overline{u}_{i}\overline{u}_{j}\right)}{\partial x_{j}}
=
-\frac{1}{\rho}\frac{\partial \overline{p}}{\partial x_{i}}
+\frac{\partial}{\partial x_{j}}\left\{\nu\left(\frac{\partial \overline{u}_{i}}{\partial x_{j}}+\frac{\partial \overline{u}_{j}}{\partial x_{i}}\right)\right\}- \frac{\partial \left(\overline{u_{i}u_{j}}\right)}{\partial x_{j}}+\frac{\partial \left(\overline{u}_{i}\overline{u}_{j}\right)}{\partial x_{j}}
\end{align}

右辺の\(\frac{\partial}{\partial x_{j}}\)について整理して

\begin{align}
\frac{\partial \overline{u}_{i}}{\partial t}
+\frac{\partial \left(\overline{u}_{i}\overline{u}_{j}\right)}{\partial x_{j}}
=
-\frac{1}{\rho}\frac{\partial \overline{p}}{\partial x_{i}}
+\frac{\partial}{\partial x_{j}}\left\{\nu\left(\frac{\partial \overline{u}_{i}}{\partial x_{j}}+\frac{\partial \overline{u}_{j}}{\partial x_{i}}\right)
-\left(\overline{u_{i}u_{j}}-\overline{u}_{i}\overline{u}_{j}\right)\right\}
\end{align}

よって

\begin{equation}
\tau_{ij}=\overline{u_{i}u_{j}}-\overline{u}_{i}\overline{u}_{j}
\end{equation}

とすると，フィルターをかけたNavier-Stokes方程式が導かれる

\begin{equation}
\frac{\partial \overline{u}_{i}}{\partial t}
+\frac{\partial \left(\overline{u}_{i}\overline{u}_{j}\right)}{\partial x_{j}}
=
-\frac{1}{\rho}\frac{\partial \overline{p}}{\partial x_{i}}
+\frac{\partial}{\partial x_{j}}\left\{\nu\left(\frac{\partial \overline{u}_{i}}{\partial x_{j}}+\frac{\partial \overline{u}_{j}}{\partial x_{i}}\right)-\tau_{ij}\right\}
\end{equation}

DNS，LES，RANSの違いを，旅行の思い出を残すための媒体で例えてみる

わかりやすい？