陰的数値法と陽的数値法
数値解法スキームは、陽的または陰的と呼ばれることがあります。従属変数を既知の量で直接的に計算できるとき、その計算は陽的であるといわれます。従属変数が、結合された複数の方程式で定義され、解を求めるのに行列式または反復計算法のどちらかが必要となる場合、その数値法は陰的であるといわれます。
数値流体力学では、支配方程式は非線形であり、一般に、非常に多くの未知の変数を伴います。こうした状況で、陰的に定式化された式は、ほぼ必ず反復計算法を用いて解きます。
反復計算では、開始状態から収束した最終状態へと連続するステップを介して、解の計算を進めます。求める解が非定常問題における1ステップであっても、最終的な定常状態の結果であっても同様です。どちらのケースでも、反復計算の各ステップは、時間の経過過程に似ています。当然ながら、反復計算におけるステップが、実際の時間依存性挙動に対応することは稀です。実際には、この陰的解法の特徴こそが、定常状態の計算に適する理由です。それは、定常状態の計算解を求めるために要求される反復計算回数が、漸近的に定常状態に近づいていく正確な非定常計算で必要となるタイムステップ数より大幅に少ないためです。
その一方で、この「歪んだ非定常」の特徴は、「時間依存性問題に対して陰的解法と陽的解法を使用した場合のそれぞれの結果はどうなるのか」という質問につながります。その答えは、2つの部分から成ります。1つ目は数値安定性に関するもので、2つ目は数値精度に関するものです。
安定性に関する問題
陰的解法は、より複雑なプログラミングを要し、各計算ステップでより多くの計算量が必要となりますが、これを使用する最大の理由は、大きなサイズのタイムステップを設定できる点にあります。ここでは、シンプルな定性モデルを用いて、その仕組みを説明します。Qをある量とし、時間t=ndtにおけるQの値に関して、時間t=(n+1)dtにおける値Qn+1を計算するとします。すなわち、Qn+1=Qn+dtSとし、SをQの変化率とします。
陽的数値法では、Sは、前のタイムステップnにおける既知の量に関して評価されます。これに対して、陰的解法では、Sの一部またはすべての項が、新しいタイムステップn+1における未知の量に関して評価されます。Qの式の右辺と左辺の両方に新しい量が出現するため、新しいn+1値の陰的な定義であるといえます。通常は、これらの新しい量を計算するために、行列解または反復計算解が必要となります。
数値安定性は、タイムステップdtを増加させたときの解の挙動に関連します。任意の大きな値を設定したタイムステップに対しても解が良好な挙動を示すのであれば、その解法は無条件に安定しているといえます。この状況は、陽的解法では見られることはなく、常に条件的に安定です。Qの式をdtで除算して、dtが無限大に近づくようにすると、これを容易に確認できます。この制限内では、式内に残っているn+1項は一切なく、Qn+1に対して存在する解もありません。つまり、解を得るには、タイムステップのサイズに対して何らかの制限を与える必要があることが示されます。
陰的な定式化では、新しいタイムステップn+1における未知の量に対する解は、どのようなタイムステップであっても求めることができます。当然ながら、大きなタイムステップで得られた解は、陰的な定式化を慎重に行った場合を除き、現実的ではない可能性があります。
Qn+1を求める典型的な反復計算式を立てるには、k回目の反復値に関してk+1回目の反復を計算します。ここで、1回目の反復がQnと等しいものとします。Qk+1の式は、多くの場合、Qk+1=Qk+A(Qn-Qk+dtSk)という形式のニュートン近似(または類似の近似法)です。この式で、Aは緩和係数、Skはk回目の反復に関して評価されたSの近似値です。Aを適切に選択すれば、以降の反復で最終的にはQn+1に収束します。
dtが小さい値および大きい値のときに適切な制限が与えられるよう、緩和係数Aは、A=1/(1+Cdt)という形式をとる必要があります。つまり、極めて小さいサイズのタイムステップでは陽的な方程式は回復しますが、極めて大きなサイズのタイムステップでは、式にはdtに依存しない制限された値が与えられます。量Cは、元の式で陰的に近似されたすべての項(すなわち、S内)の特性を定義する正の係数です。たとえば、Qが陰的な粘性項を伴う運動量方程式によって支配される速度成分であるとき、Cは動粘度を格子サイズの2乗で除算した値に比例します。
精度に関する問題
dtが十分に小さい場合、収束には1回の反復計算のみを要し、Qn+1=Qn+dt/(1+Cdt)Snとなります。これは、タイムステップに乗算される不足緩和係数A=1/(1+Cdt)によって、陽的解法よりも陰的な定式化の方が、1つのタイムステップあたり、Qにより小さな変化が加えられることが示されます。
原則として、条件Cdt≤1は、陽的近似に対する安定条件とほぼ同じであることがわかります。さらに、もう1つの原則として、陽的解法の安定性と精度を得るためのタイムステップサイズは、通常は等価です。したがって、Cdt>1のとき、陽的解法は不安定になりますが、陰的解法では反復計算による解の安定性を維持するために、単に不足緩和係数が小さく調整されます。このタイムステップサイズの増大とともに増大する減衰が、非定常挙動の不正確さをもたらします。
陰的解法で最小の不足緩和(すなわち、小さな減衰)が与えられるには、安定した陽的解法の値よりも大幅にサイズが小さいタイムステップを使用しなければなりません。実際には、上記の解析により、陽的解法の安定性制限値Cdt=1において、陰的近似ではまだA=1/2という大きな不足緩和係数を有します。この不足緩和減衰を減少させるには、陽的解法の安定性制限値よりも大幅に小さいサイズのタイムステップが必要となりますが、陰的解法は要求されないため、ほとんど意味はありません。
物理的現象例
圧力波の伝搬を伴う単純な物理現象問題を使用して、陰的解法と陽的解法の違いを示します。開口端から圧力が加えられ、反対側は閉口端であるオルガンパイプを例にとります。圧力波がパイプを通過して、閉口端で反射することがわかっています。十分な時間が与えられると、圧力波はパイプ内を何度も往復して、最終的に圧力分布は開口部に加えられる一定値に落ち着きます。
定常状態の結果のみを求める場合は、圧力波の減衰が多い陰的解法スキームを使用して、できる限り早く定常状態に到達するようにします。この場合、陰的反復計算法に組み込まれた減衰(すなわち、不足緩和)は、強く望まれるものです。
それに対して、非定常の圧力波を調べる場合は、多くの波反射を正確に追跡できるよう、数値的減衰を最小限に抑えることが望まれます。この状況は、陽的解法で解くのが最良です。
陽的解法では、圧力ステップの進行をタイムステップあたり1計算セル未満に制限するようなサイズのタイムステップを必要とします。ただし、大半の差分方程式は隣接セルのみからの量を考慮するため、この制限は精度に関連します。1つのタイムステップで1セル以上を伝搬する圧力波は、圧力に対して定義された影響が一切ない領域に流入します。これは、物理的に非現実的であるばかりか、数値不安定性にもつながります。
それに対して、陰的解法では、反復計算法を通じてすべてのセルが結合されることで、圧力信号が格子内を伝播します。このように離れた距離にあるセル間でのやり取りの代償は、結合方程式を解くのに必要な不足緩和によって生じる、圧力波の減衰または平滑化です。
陰的解法と陽的解法のどちらを採用するかは、最終的には計算での目標に応じます。時間精度を重視する場合は、陽的解法の方が、陰的解法と比べて少ない計算量でより精度の高い解を得られます。こうした理由から、FLOW-3Dでは適用可能な問題にはできる限り陽的解法を使用しますが、必要に応じて陰的解法オプションも利用できます。