離散化スキームの設定2021年1月23日 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
はじめに離散化スキームの設定について。 使用バージョンOpenFOAM v2006、v8 離散化スキームの設定有限体積法による離散化に必要な離散化スキームの設定は system/fvSchemes で行う。
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
ddtSchemes
{
default steadyState;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) bounded Gauss upwind;
div(phi,k) bounded Gauss upwind;
div(phi,epsilon) bounded Gauss upwind;
div(phi,omega) bounded Gauss upwind;
div(phi,R) bounded Gauss upwind;
div(R) Gauss linear;
div((nuEff*dev(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
default none;
laplacian(nuEff,U) Gauss linear corrected;
laplacian((1|A(U)),p) Gauss linear corrected;
laplacian(DkEff,k) Gauss linear corrected;
laplacian(DepsilonEff,epsilon) Gauss linear corrected;
laplacian(DomegaEff,omega) Gauss linear corrected;
laplacian(DREff,R) Gauss linear corrected;
laplacian(1,p) Gauss linear corrected;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default corrected;
}
fluxRequired
{
default no;
p ;
}
設定方法は、各微分演算子などの項目 (ddtSchemes など) に対して、それに関する方程式の項について離散化スキームを指定するものになっている。したがって、ある程度ソルバーが解いている方程式の中身を理解している必要がある。"default" については、個別に指定がなければこの設定を使うというものである。 ddtSchemes時間微分の離散化スキームを指定する。simpleFoam などの定常解析ソルバーの場合はつねに "steadyState" を指定する。pimpleFoam などの非定常解析ソルバーの場合は、つぎのようなものを指定できる。
CrankNicolson はパラメタを取り、つぎのように設定する。
default CrankNicolson 1;
数値は 0〜1 の調整パラメタで、1 ならば純粋な CrankNicolson 法、0 ならば Euler 法になる。 CoEulerは時間刻み幅をクーラン数に応じて局所的に変えるもので、非定常解析ソルバーで定常解析を行うときに使用する。つぎのように設定する。
default CoEuler phi rho 0.5;
rho は密度で、ソルバーに対応したものを指定する。最後の数値はクーラン数である。 SLTS (stabilised local time-step) は、移流方程式が対角優位性を保つように時間刻み幅を局所的に変えるもので、非定常解析ソルバーで定常解析を行うときに使用する。CoEuler 同様、つぎのように設定する。
default SLTS phi rho 0.5;
最後の数値は緩和係数である。 通常の非定常解析の場合、時間微分項のスキームは基本的には Euler でよいが、数値拡散を避けたい場合は backward や CrankNicolson を使う。 divSchemes で "bounded" を使う場合は、ddtSchemes のほうでも "bounded" を指定する。
default bounded Euler;
gradSchemes勾配の離散化スキーム。ここで指定するものはセルの界面の補間スキームであり、つぎのようなものを指定できる。
基本的には "Gauss linear" でよいが、メッシュのひずみが大きい場合などには "leastSquares" などのほうがよい場合がある。 勾配の値に制限をかけることができる。つぎのように設定する。
default cellLimited Gauss linear 1;
最後の数値は制限の強さを指定する 0〜1 の値で、1 が制限、0 が制限なしである。制限にはつぎのものが指定できる。
これらは Barth and Jesperson の方法 (公式は Minmod と呼んでいるが) に基づくもので、上のリストは制限が強い順番に並んでいる。制限が強い方が計算は安定するが、精度は落ちる。なるべく cellLimited など制限が小さいものを使用したほうがよい。 faceLimited と cellLimited の違いは、勾配制限関数を計算するときの注目セルと隣接セルの値の大小を比較する手順の違いである。faceLimited はフェイスごとに隣接セルと値を比較するのに対し、cellLimited は一度に全ての隣接セルと値を比較する。"MD" が付いたものは多次元版である。 OpenFOAM 6 から、cellLimited に Venkatakrishnann および Michalak and Ollivier-Gooch の制限関数を利用できる (参考)。
limitedGrad(U) cellLimited<Venkatakrishnan> Gauss linear 1;
limitedGrad(U) cellLimited<cubic> 1.5 Gauss linear 1;
cubic で設定されている数値は、関数を切り替えるためのパラメタである (原著論文では 1.5 が用いられている)。Barth-Jesperson 制限関数は不連続なため、それが収束性を悪化させるらしい。Venkatakrishnan の制限関数は本来連続であるが、OpenFOAM の実装では値をクリッピングしてあるため不連続である。そのため、連続である Michalak and Ollivier-Gooch の制限関数 (cubic) の使用が推奨されている。 勾配制限は divSchemes で limitedUpwind を使う場合に数値的な振動を抑えるために重要になる。通常は、default ではなく特定の勾配に設定する。
default Gauss linear;
grad(U) cellLimited Gauss linear 1;
ソルバー全体に影響を与えたくない場合は、適当な名前をつけて参照するとよい。
gradSchemes
{
default Gauss linear;
limitedGrad(k) cellLimited Gauss linear 1;
}
divSchemes
{
...
div(phi,k) Gauss linearUpwind limitedGrad(k);
}
以下の divSchemes の項も参照のこと。 divSchemes発散の離散化スキームを指定する。これらは計算上特に重要な設定である。項の中で "div(phi,...)" の形で書かれたものは対流項であり、スキームを慎重に選ぶ必要がある。つぎのようなものを指定できる。
blended は linear と upwind を混ぜたもので、高次精度スキームだと発散するが upwind は使いたくないときなどに利用できる。
div(phi,k) Gauss blended 0.1;
最後の数値は linear の混合率 (0〜1) で、0 だと upwind、1 だと linear になる。 linearUpwind は、つぎのように設定する。
div(phi,k) Gauss linearUpwind grad(k);
最後で勾配 (実際は勾配のスキーム) を指定している。このスキームは多少値の振動が起こるので、それを避けたい場合は勾配制限を使う。スキームで指定する分以外の勾配には制限をかけたくない場合は、つぎのように設定する。
gradSchemes
{
default Gauss linear;
limitedGrad(k) cellLimited Gauss linear 1;
}
divSchemes
{
div(phi,k) Gauss linearUpwind limitedGrad(k);
}
limitedLinear は、つぎのように設定する。
div(phi,k) Gauss limitedLinear 1;
最後の数値は 0〜1 で安定性を調整するもので、1 が安定である。 Gamma は、つぎのように設定する。
div(phi,k) Gauss Gamma 1;
最後の数値は 0〜1 で安定性を調整するもので、1 が安定である。オリジナルの Gamma スキームのパラメタ (β) は 0.1〜0.5 をとるが、OpenFOAM では 0〜1 にスケーリングしてある。 ベクトル用のスキームが別途用意してあるものがあり、それは最後に "V" の文字が付く。たとえば limitedLinear のベクトル版は limitedLinearV などである。 それぞれのスキームの違いは主に数値拡散の大きさであり、数値拡散が小さいほうが計算精度はよいが、一方で安定性が減少して計算できなくなる可能性があり、これといって決め手はない。とにかく安定して計算したい場合は upwind を使う。高次精度のスキームでは計算が発散しやすいが、upwind の計算結果から継続すると計算しやすくなる場合がある。 スカラーの変数で、値が物理的にある範囲を超えようがない場合に、高次精度のスキームだと数値的な振動によりその範囲を超えてしまい問題になることがある。その場合は TVD スキームなど振動を起こさないスキームを使用したほうがよい。たとえば、U に対しては linearUpwindV を使い、k や epsilon などには limitedLinear などを使うとよいだろう。 計算安定化処理を使用するために "bounded" という指定ができる。
div(phi,U) bounded Gauss upwind;
非定常解析ソルバーでこれを使う場合は、ddtScheme でも同様に指定する。これは収束性を改善するもので、定常解析ソルバーでは基本的に有効にしておくべきだが、非定常解析ソルバーの場合は反復計算をするかどうかで使い分ける。 OpenCFD 版 OpenFOAM v1712 から、緩和修正法 (deferred correction) を使用することができる。
div(phi,U) Gauss deferredCorrection linear;
高次精度スキームを 1 次精度風上差分とそれ以外に分解し、後者を陽的に扱うことで計算を安定させる。 laplacianSchemesラプラシアンの離散化スキームを指定する。"laplacian(nuEff,U)" や "laplacian(DkEff,k)" などは拡散項、"laplacian((1|A(U)),p)" などは圧力方程式の設定である。
laplacian(nuEff,U) Gauss linear corrected;
"Gauss linear" は gradSchemes で設定したものと同じセルの界面の補間スキームで、上の例では nuEff の補間に用いられる。通常は linear を設定すればよい。 ラプラシアンの計算にはセルの界面の法線方向の勾配 (surface-normal gradient: snGrad) が用いられる。セル中心間と界面が直交していると想定して計算するのが単純だが、一般的には非直交になるので、セル中心間勾配の界面方向へのマッピングが行われる。それでもメッシュのゆがみが大きいときには正確ではないので、さらに補正が行われる。OpenFOAM ではこれを非直交補正 (non-orthogonal correction) と呼んでいる。非直交補正の項は陽的に扱われ、非直交補正ループにより処理される。ただし、メッシュのゆがみが大きすぎる場合、この補正が計算の安定性を損なわせるため、補正に制限をかけられるようになっている。 snGrad スキームは、上の設定例の "corrected" の部分で設定する。以下のものが指定できる。
ややこしいが、uncorrected は非直交補正を行わないだけで、非直交性は考慮している。orthogonal は非直交性自体を考慮しない。 snGrad スキームの選択は、直交メッシュでない限りは、checkMesh で得られる非直交性 (non-orthogonality) で判断する。非直交性が 5 以下 (誤差が約 5% 以下) であれば uncorrected でよい。5〜60 では corrected にする。60 以上だと補正が大きくなってくるので、たとえば limited 0.5 (補正を法線方向勾配の大きさ以下に制限) や limited 0.333 (補正を法線方向勾配の半分以下に制限) などを使うとよい。非直交性が 80 以上の場合は、計算が難しくなってくるので、いっそ uncorrected にしてしまうか、素直にメッシュを作り直したほうがよい。 snGrad スキームには、上記以外に faceCorrected というものもある。上記のものは界面を挟んで隣接するセル中心の値から界面の法線方向の勾配を計算しようとするが、faceCorrected は値を節点に補間した上で界面の辺中心の値を計算し、それらにより界面の勾配を積分で求め、それを用いて界面の法線方向の勾配を計算する。 interpolationSchemesgradSchemes で設定したものと同じセルの界面の補間スキームである。基本的には linear でよい。 snGradSchemeslaplacianSchemes で設定した snGrad スキームの設定と同じである。laplacianSchemes と合わせておけばよい。 fluxRequiredこれは流束を計算させるかどうかの指定で、ソルバーが要求したときに設定する。 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| PENGUINITIS | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||