収束基準 の選択

収束基準 の選択

完全に非圧縮性である物質は存在しませんが、多くの場合、非圧縮性であるとの仮定は良好な近似となります。この仮定を数値解法スキームで使用する場合は、非圧縮性挙動をもたらす物理的メカニズムを与える何らかの方法を考える必要があります。

(1) \frac{d\rho }{dt}+\rho \nabla \cdot \vec{U}=0

非圧縮性は、速度の発散Ñ ·Uで計算します。これは、単位体積あたりの流体体積の時間変化率と等価です。発散は、質量保存の式内の項として出現します。ここで、dr/dtは、流体とともに移動する流体要素の密度変化です。流体が非圧縮性のとき、その密度は常に一定で、式1から速度の発散はゼロでなければならないことが示されます。

\nabla \cdot \overrightarrow{U}=0

非圧縮性の挙動を与える物理的メカニズムは、圧力波の急速な伝搬であり、流体の物質速度よりもはるかに速く流体内を移動しなければなりません。ほとんどの場合、圧力波の数値的な伝搬は、圧力と速度を連成させた何らかの反復計算スキームで求めることができます。このスキームでは、速度の発散の大きさが、 収束基準 と呼ばれる、設定した絶対数値eを下回るまで計算を繰り返すことが目標となります。

収束基準の選択

「どのように 収束基準を選択すべきか」は、よくある質問です。その答えはシンプルながら捉えにくいもので、解を求めるために選択した数値法によって異なります。

式1を離散化した式で、1つのタイムステップdtで生じる密度変化が示されます。

(2) \frac{\delta \rho }{\rho }=\delta t\nabla \cdot \overrightarrow{U}

この結果については、確認すべきことがいくつかあります。まず、右辺は基本的に以下の桁数となります。

\delta t\nabla \cdot \overrightarrow{U}\approx \frac{\delta tU}{\delta x}

ここで、dxは数値解法で使用する離散化された要素の長さまたはサイズです。この量は、しばしば流れの「クーラン数」と呼ばれます。正確な結果を得るには、クーラン数を1未満にすべきであることが知られています。また、陽的計算法を使用する場合も、数値安定性を得るために1未満とすべきです。

関係式2から、クーラン数が1を超えると、密度変化の値が密度自体よりも大きくなる可能性があることがわかります。この点からも、クーラン数が精度や安定性といった問題に緊密に関係している理由がわかります。さらに、空間的な次元数が1から2または3へと上がるにつれて、安定条件がより厳しくなる理由も説明されます。安定性を得るには、全方向でのフラックスから生じる変化の累積によって、密度(または他の量)を現在の値よりも大きく変化させてはなりません。したがって、より多次元ではさらに多くのフラックスが生じるため、要素が必要以上に空になることを避けるために、より小さなタイムステップが必要となります。

関係式2に戻ると、与えられた収束基準eから予期される密度変化量がわかります。たとえば、e = 0.001とすると、1つの時間単位で予期される密度変化は1%の1/10であり、適用される問題の多くで「非圧縮性」とみなすのに十分な小ささです。

e = 0.001は、非圧縮性流れ計算の多くのケースに適しているといわれる値です。自然対流の場合や、極めて小さな密度変動が重要となる問題では、支配方程式をどのように定式化するかに応じて、より小さな値をeに指定しなければならないこともあります。

eの値が大きすぎる場合かつ圧力の解法でいくらかの過緩和を組み込む場合、流体が圧縮性の挙動を示すことがあります。この場合、数値的な圧力波の計算領域にわたる移動が見られることもあります。eの値を小さくするか、または過緩和を減らすと、こうした人工的な波の振幅が小さくなります。

答えるべき質問の1つとして、1つのタイムステップで生成される密度の誤差が、多数のタイムステップにわたって大きな影響を与えるほど累積される可能性があるかどうかという問題があります。これについては、適切な数値法を採用すれば、そうした可能性はありません。

自己修正型数値法の選択

数値反復計算法を用いる場合、一般的には、精度の高い結果を得るために、極めて精密な収束基準が必要であると考えられています。計算式が「自己修正」型で定式化されている場合、その考え方は必ずしも正しくはありません。基本的な考え方は、非圧縮性の流体に対する式を用いて説明することができます。ここでは、非圧縮性条件として速度の発散Dをゼロとします。

(3) D=\nabla \cdot \overrightarrow{U}=0

流体の圧力pは、運動量方程式(ナビエストークス方程式など)の発散から導かれる式を満たす必要があります。すべての項を保持して、以下の圧力式が成り立ちます。

(4) \frac{\partial D}{\partial t}={{\nabla }^{2}}\left( {p}/{\rho }\; \right)+\upsilon {{\nabla }^{2}}D-Q

ここで、

Q=\nabla \cdot \left[ \nabla \cdot \left( \overline{UU} \right) \right]

ここで、rは流体密度、uは動粘度です。これは、圧力のポアソン方程式で、収束基準が必ず設定される反復計算プロセスによって、各タイムステップが最も効率的に計算されます。

多数のタイムステップにわたって累積される反復計算の収束誤差を最小限に抑えるために、各タイムステップで高レベルの収束を設定するか、あるいはMAC (Marker-and-Cell)法(F.H. Harlow、J.E. Welch、Physics of Fluids Volume 8、p.2182 (1965))とともに初めて使用された自己修正型処理を使用します。式4を解くために考案された多くの数値法では、式3から、左辺をゼロとしたいと思われるかもしれません。しかし、自己修正型処理では、左辺は以下で置き換えられます。

\frac{\partial D}{\partial t}\approx \frac{{{D}^{n+1}}-{{D}^{n}}}{\delta t}\to \frac{-{{D}^{n}}}{\delta t}

ここで、nはn番目のタイムステップを表します。

Dn項を保持する修正処理によって、前サイクルから残った速度発散における誤差が除去される傾向があります。たとえば、n番目のタイムステップの最後で、いずれかの要素内にわずかな膨張に相当する少量の正の残差がDに存在していた場合、ステップn+1では、それを修正するための同量の圧縮が生成されるソース項(-Dn/dt)が式に追加されます。

この処理によって、圧力の方程式に比較的粗い収束基準を設定した場合でも、非圧縮性の累積誤差は少なくなります。自己修正型処理に関する詳細については、次の参考文献を参照してください。C.W. Hirt、F.H. Harlow、Journal of Computational Physics Volume 2、p.114 (1967)。

自己修正型機能の最大の特長としては、必ずしも評価されるわけではありませんが、初期条件をそれほど厳密に扱わない点があります。一部の非圧縮性流れソルバでは、非圧縮性条件を満たす初期条件が要求されます。そのため、これらの手法では、使用可能な初期条件を得るために、別の問題に対する解が必要となることがあります。MAC法や自己修正型機能を使用する他の手法では、初期条件における非ゼロの速度発散は、計算の最初のサイクルで自動的に除去されます。ユーザにとって、これは非常に便利な機能であることはいうまでもありません。

FLOW-3Dプログラムは、自己修正型処理だけでなく、解の計算中に発生している状態に応じてeが調整される収束基準の自動設定機能を備えています。後者もまた、さまざまな流れの段階を有する問題に対して極めて有用で、ユーザにとって扱いやすい機能です。

^ back to top