ゼロ除算

数値計算ではゼロ除算を防止するために，分数の分母にゼロに限りなく近い小さい値(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でもあったような，と頭の片隅に入れておくぐらいでいいのかもしれない．

参考文献

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