ITとCFD入門サイト 
格子法を用いた流体解析（CFD）の解説を行っています。第５回は、陽解法と陰解法について解説します。
前回は圧力分離型ソルバー（MAC法）について解説しました。まだ見てないかたはこちらからどうぞ。
ざっくり説明すると、CFDではポアソン方程式で圧力項だけ分離させて計算することが一般的であり、その圧力項をナビエ・ストークス方程式に入れることで速度を更新する、という流れをMAC法を例に解説しました。
陽解法（Explicit）
まずは前回の内容を振り返りましょう。前回は、圧力分離型ソルバーであるMAC法について解説しました。
MAC法では、まずポアソン方程式を用いて圧力を計算し、その後にナビエストークス方程式を使って速度を更新するという流れでした。ポアソン方程式は下記のとおりです。
$$\nabla^2 p = \frac{\rho}{\Delta t}\nabla \cdot v$$
ナビエストークス方程式は下記を使用します。
$$\frac {\partial v}{\partial t}+(v^n \cdot \nabla )v^n=\frac {1}{\rho} \nabla p^* + { \nu \nabla^2 v^n}$$
ポアソン方程式で得た圧力を仮の圧力（*）とし、速度は現在の値（n）を使用しています。左辺第一項の非定常項を離散化すると下記のように次の時間ステップが現れます。
$$\frac {\partial v}{\partial t} = \frac{u^{n+1}-u^{n}}{\Delta t}$$
このように、次の時間ステップの速度の値（$v^{n+1}$）と現在の時間ステップの速度の値（$v^n$）の差分が出てきます。ここでは簡単化のために非定常項しか差分化していませんが、全ての項を差分化した後に左辺にn+1の項が集まるように式を整理すれば、次の時間ステップにおける速度の計算ができます。
さて、前回行った計算を振り返りました。ここで今回重要になってくるのが、ナビエストークス方程式で使用した速度の時間ステップです。
$$\frac {\partial v}{\partial t}+(v^n \cdot \nabla )v^n=\frac {1}{\rho} \nabla p^* + { \nu \nabla^2 v^n}$$
ここで、非定常項以外に注目してください。全ての速度変数が現在の時間ステップにおける速度（$v^n$）を使用していることがわかると思います。
このように現在の時間ステップの値（n）を用いて未来の時間ステップの値（n+1）を計算する手法を陽解法といいます。
当たり前のように見えるかもしれませんが、陽解法は安定性が低くなります。実はこれから紹介する陰解法と呼ばれる別の手法のほうが安定性が良いです。
陰解法（Implicit）
陰解法は未来の時間ステップの値（n+1）を用いて未来の時間ステップの値（n+1）を計算する手法です。
ナビエストークス方程式は以下のようになります。
$$\frac {\partial v}{\partial t}+(v^{n+1} \cdot \nabla )v^{n+1}=\frac {1}{\rho} \nabla p^* + { \nu \nabla^2 v^{n+1}}$$
右辺の速度が全て未来の値（n+1）になっていることがわかります。非定常項を離散化して、左辺に未知の変数、右辺に既知の変数で分けてみましょう。
$$u^{n+1}+\Delta t (v^{n+1} \cdot \nabla )v^{n+1} \Delta t { \nu \nabla^2 v^{n+1}}=u^{n+1} + \Delta t \frac {1}{\rho} \nabla p^*$$
左辺は非定常項の未知の速度（n+1）・移流項・粘性項で構成されています。右辺は非定常項の既知の速度（n）・圧力項で構成されています。
この式は未知の変数が多いため単純な計算では解けません。よって、行列による収束計算を行います。
ここでは詳細は説明しませんが、収束計算のイメージとしては、少しずつ変数を修正していって等号が成り立つ条件を探すとお考え下さい。
予想されている通り、陰解法では収束計算を扱うため、同じ条件では陽解法よりも計算時間が長くなります。ただし、陰解法にも利点があります。未知の値を使用して計算しているので安定性が良いというのもありますが、一番の強みはクーラン数の制限が実質無制限になることです。
クーラン数に関するCFL条件とは下記の移流に関する式になります。
$$v \frac{\Delta t}{\Delta x} < 1$$
詳細は別記事をご覧ください。
CFDにはクーラン数と拡散数という制限があり、陽解法では必ずこれらを満たす必要があります。一方で、陰解法ではこれらを無視することができます。
仕組みに関しても少し説明しておきます。クーラン数と拡散数を満たさないといけない理由は、$\Delta t$が低くて$\Delta x$が大きいと情報の伝搬が伝わりきらないからです。よって、これらの情報を満たさないと情報が失われて計算が破綻します。
一方で、陰解法では収束計算の際に全てのセルに情報が伝搬するため、理論上は制限がないことになるのです。この特性は非常に便利なので、CAEソフトウェア等の実用の場面では多く使用されています。
ただ、陰解法でなんでも上手くいくかというとそうではありません。陰解法でも時間ステップは小さくとったほうが良いです。逆に間違った答えでも平気で出力するので、正しく結果を解釈する能力が必要になります。
おわりに
今回は陽解法と陰解法について解説しました。現在の時間ステップの値（n）を用いて未来の時間ステップの値（n+1）を計算する手法を陽解法といい、未来の時間ステップの値（n+1）を用いて未来の時間ステップの値（n+1）を計算する手法を陰解法といいます。
陰解法では収束計算を使用するので、クーラン数に関するCFL条件や拡散数の条件が当てはまらなくなります。時間ステップを多少大きめにとっても計算ができるので、広く使用されています。
次回は移流項の離散化について解説します。格子法において、移流項の離散化が一番の山場です。あまりにも難しい話はしないつもりですが、覚悟して挑んでください。
youtubeでも解説しています。そちらもどうぞ。
