Edited at

OpenFOAMでのゼロに近い値を表す変数

More than 1 year has passed since last update.


ゼロ除算

数値計算ではゼロ除算を防止するために,分数の分母にゼロに限りなく近い小さい値(0.000001など)を足すことがよくある.OpenFOAMでコードをカスタマイズするときに,直にこのような数値をコード中に書いてしまいがちであるが,OpenFOAMではSMALLなどのスカラー変数があらかじめ用意されている.


SMALL, VSMALL, GREAT, VGREAT, etc.

CFD-ONLINEで質問している人がすでにいて,以下のスカラー変数が使える.

GREAT = 1.0e+6;

VGREAT = 1.0e+37;
ROOTVGREAT = 1.0e+18;
SMALL = 1.0e-6;
VSMALL = 1.0e-37;
ROOTVSMALL = 1.0e-18;

ただしこれはfloat型の場合で,double型の場合は,

GREAT = 1.0e+15;

VGREAT = 1.0e+300;
ROOTVGREAT = 1.0e+150;
SMALL = 1.0e-15;
VSMALL = 1.0e-300;
ROOTVSMALL = 1.0e-150;

である.これらの値が定義されているのは,$FOAM_SRC/OpenFOAM/primitives/Scalar/scalar/scalar.H$FOAM_SRC/OpenFOAM/primitives/Scalar/floatScalar/floatScalar.Hを参照すればすぐにわかる.


実例: 非ニュートン流体での粘性係数計算

通常の分数除算ではSMALLで十分かと思うが,非ニュートン流体での粘性係数の計算などはべき乗が使われるのでVSMALLを使うようである.例えば,power law流体の粘性係数は物理式で書くと,

$$

\eta = \max (\eta_{\min}, \min (\eta_{\max}, k\dot\gamma^{n-1}))

$$

で計算されるが,ソース($FOAM_SRC/transportModels/incompressible/viscosityModels/powerLaw/powerLaw.C)では以下のようになっている.

Foam::viscosityModels::powerLaw::calcNu() const

{
return max
(
nuMin_,
min
(
nuMax_,
k_*pow
(
max
(
dimensionedScalar("one", dimTime, 1.0)*strainRate(),
dimensionedScalar("VSMALL", dimless, VSMALL)
),
n_.value() - scalar(1.0)
)
)
);
}

中ほどのmax関数の引数のところでVSMALLが使われており,strainRate()がゼロになるときは,VSMALLの方が大きいのでこの値が利用されることがわかる.


まとめ

特別な数値を表す変数の名前についてはコード(ソフトウェア)によって異なるので,OpenFOAMでもあったような,と頭の片隅に入れておくぐらいでいいのかもしれない.


参考文献


  1. 勉強会の活動を通した非ニュートン流体モデルカスタマイズの取組