離散化スキームの設定

2013年12月15日

はじめに

離散化スキームの設定について。

使用バージョン

OpenFOAM 2.2.2

離散化スキームの設定

有限体積法による離散化に必要な離散化スキームの設定は 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 などの非定常解析ソルバーの場合は、つぎのようなものを指定できる。

ddtSchemes説明
EulerEuler 法 (1 次精度)
CrankNicolsonCrank-Nicolson 法 (2 次精度)
backward後退差分 (2 次精度)
CoEuler時間刻み幅をクーラン数に応じて局所的に変える
SLTS時間刻み幅を局所的に変える

CrankNicolson はパラメタを取り、つぎのように設定する。

    default         CrankNicolson 1;

数値は 0〜1 の調整パラメタで、1 ならば純粋な CrankNicolson 法、0 ならば Euler 法になる。

CoEulerは時間刻み幅をクーラン数に応じて局所的に変えるもので、非定常解析ソルバーで定常解析を行うときに使用する。つぎのように設定する。

    default CoEuler phi rho 1;

rho は密度で、ソルバーに対応したものを指定する。最後の数字はクーラン数である。

SLTS (Stabilised Local Time-Step) は時間刻み幅を局所的に変えるもので、非定常解析ソルバーで定常解析を行うときに使用する。CoCuler 同様、つぎのように設定する。

    default SLTS phi rho 1;

最後の数字は緩和係数である。

基本的には "Euler" でよいが、数値拡散を避けたい場合は "backward" などを使う。

divSchemes で "bounded" を使う場合は、ddtSchemes のほうでも "bounded" を指定する。

    default         bounded Euler;

gradSchemes

勾配の離散化スキーム。ここで指定するものはセルの界面の補間スキームであり、つぎのようなものを指定できる。

gradSchemes説明
Gauss linear線形補間
Gauss midPoint算術平均 (足して 2 で割る)
Gauss pointLinear点の値から補間
leastSquares最小二乗法

基本的には "Gauss linear" でよいが、メッシュのひずみが大きい場合などには "leastSquares" などのほうがよい場合がある。

勾配計算に制限をかけることができる。つぎのように設定する。

    default         cellLimited Gauss linear 1;

最後の数字は制限の強さを指定する 0〜1 の値で、1 が制限、0 が制限なしである。制限にはつぎのものが指定できる。

  • faceLimited
  • faceMDLimited
  • cellLimited
  • cellMDLimited

これらは Barth and Jesperson スキームに基づくもので、上のリストは数値拡散が強い順番に並んでいる。数値拡散が強い方が計算は安定するが、精度は落ちる。なるべく cellLimited など数値拡散が小さいものを使用したほうがよい。

勾配制限は divSchemes で limitedUpwind を使う場合に数値的な振動を抑えるために重要になる。

divSchemes

発散の離散化スキームを指定する。これらは計算上特に重要な設定である。項の中で "div(phi,*)" の形で書かれたものは対流項であり、スキームを慎重に選ぶ必要がある。つぎのようなものを指定できる。

divSchemes説明
linear線形補間 (中心差分、2 次精度)
upwind1 次精度風上差分
linearUpwind線形風上差分 (2 次精度風上差分)
QUICKQUICK スキーム (2 次精度)
Minmodminmod 制限関数 (2次精度 TVD)
SuperBeesuperbee 制限関数 (2次精度 TVD)
vanLeervan Leer 制限関数 (2次精度 TVD)
vanAlbadavan Albada 制限関数 (2次精度 TVD)
UMISTUMIST 制限関数 (2次精度 TVD)
MUSCLMonotonized central-difference 制限関数 (2次精度 TVD)
limitedLinear線形補間に TVD 制限をつけたもの
GammaGamma スキーム (NVD)

linearUpwind は、つぎのように設定する。

    div(phi,k)      Gauss linearUpwind grad(k);

最後で勾配 (実際は勾配のスキーム) を指定している。このスキームは多少値の振動が起こるので、それを避けたい場合は勾配制限を使う。スキームで指定する分以外の勾配には制限をかけたくない場合は、つぎのように設定する。

gradSchemes
{
    default         Gauss linear;
    limitedGrad(k)  cellLimited Gauss linear 1;
}

divSchemes
{
    div(phi,k)      bounded 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 でも同様に指定する。これは収束性を改善するもので、定常解析ソルバーでは基本的に有効にしておくべきだが、非定常解析ソルバーの場合は反復計算をするかどうかで使い分ける。

laplacianSchemes

ラプラシアンの離散化スキームを指定する。"laplacian(nuEff,U)" や "laplacian(DkEff,k)" などは拡散項、"laplacian((1|A(U)),p)" などは圧力方程式の設定である。

    laplacian(nuEff,U) Gauss linear corrected;

"Gauss linear" は gradSchemes で設定したものと同じセルの界面の補間スキームで、上の例では nuEff の補間に用いられる。通常は linear を設定すればよい。ラプラシアンの計算にはセルの界面の法線方向の勾配 (snGrad) が用いられるが、単純な計算だとメッシュのゆがみが大きい時に精度が落ちるため、補正が行われる。これを非直交補正 (non-orthogonal correction) という。最後の "corrected" の部分はこの補正の設定である。ただし、メッシュのゆがみが大きすぎる場合、この補正が計算の安定性を損なわせるため、補正に制限をかけられるようになっている。補正の設定としては以下のものが指定できる。

指定意味
uncorrected/orthogonal補正なし
corrected補正あり
limited 0〜1補正の制限 (0 で uncorrected、1 で corrected と同じ)

補正のかけかたは、checkMesh で得られる非直交性 (non-orthogonality) で判断する。非直交性が 5 以下 (誤差が約 5% 以下) であれば uncorrected でよい。5〜60 では corrected にする。60 以上だと補正が大きくなってくるので、たとえば limited 0.5 (補正を法線方向勾配の大きさ以下に制限) や limited 0.333 (補正を法線方向勾配の半分以下に制限) などを使うとよい。非直交補正が 80 以上の場合は、計算が難しくなってくるので、いっそ uncorrected にしてしまうか、素直にメッシュを作り直したほうがよい。

interpolationSchemes

gradSchemes で設定したものと同じセルの界面の補間スキーム。基本的には linear でよい。

snGradSchemes

laplacianSchemes で設定した非直交補正の設定と同じ。laplacianSchemes と合わせておけばよい。

fluxRequired

これは流束を計算させるかどうかの指定で、ソルバーが要求したときに設定する。