interFoam のスキームの検討

2012年11月9日

はじめに

interFoam において、初期の液面がメッシュと沿っていない場合に、液面付近で変な速度が出る。この速度を極力抑えるためにスキームを吟味してみる。

使用バージョン

OpenFOAM 2.1.1

問題

interFoam のチュートリアルケース damBreak をベースにする。system/setFieldsDict を次のように変更する。

    boxToCell
    {
        //box (0 0 -1) (0.1461 0.292 1);
	box (0 0 -1) (1 0.292 1);
        fieldValues
        (
            volScalarFieldValue alpha1 1
        );
    }

これで setFields を実行すると、つぎのようになる。

interFoam を実行すると、つぎのようになる。

液面には変化なし。液面上に少し速度がでる。

今度はメッシュを少し変更する。constant/polyMesh/blockMeshDict の vertices の "(4 4 0)" と "(4 4 0.1)" をそれぞれ "(4 3 0)" と "(4 3 0.1)" にする。これで blockMesh を実行するとつぎのようになる。

setFields 実行結果。

メッシュ線が液面に沿っていないため、初期がギザギザになる。

interFoam 実行結果。

ギザギザの液面が平らになろうとするため、液面が左に移動する。その際の液面付近の速度が速すぎるような気がするので、スキームの変更でどうにかなるか検討してみる。

スキームの検討

1 秒後の最大速度でスキームを評価する。スキームは以下を検討する。

• ddtSchemes: Euler
• gradSchemes: linear, midPoint, pointLinear, leastSquares, cellLimited, faceLimited (4 x 3 = 12 ケース)
• divSchemes, div(rho*phi, U): limitedLinearV, upwind (2 ケース)
• divSchemes, div(phi,alpha): vanLeer, vanLeer01 (2 ケース)
• laplacianSchemes: corrected, limited 0.5, limited 0.333, uncorrected (4 ケース)
• interpolateSchemes: linear, midPoint, pointLinear (3 ケース)
• snGradSchemes: laplacianSchemes に合わせる

存在するスキームをすべて検討するとたいへんなので、適当に限定。12 x 2 x 2 x 4 x 3 = 576 ケース。fvSolution のほうは、ちょっと試した感じではあまり効果がなさそうなので、ここでは無視する。

スクリプトを作成。

run

#!/bin/sh

DDT="Euler"
GRAD="linear midPoint pointLinear leastSquares cLinear cMidPoint cPointLinear cLeastSquares fLinear fMidPoint fPointLinear fLeastSquares"
DIVPHI="limitedLinearV upwind"
DIVALPHA="vanLeer vanLeer01"
LAPLACIAN="corrected 0.5 0.333 uncorrected"
INTERPOLATE="linear midPoint pointLinear"

MINVEL=100000
rm -r CASE*

I=1
for DDT1 in $DDT ; do
for GRAD1 in $GRAD ; do
for DIVPHI1 in $DIVPHI ; do
for DIVALPHA1 in $DIVALPHA ; do
for LAPLACIAN1 in $LAPLACIAN ; do
for INTERPOLATE1 in $INTERPOLATE ; do

CASE=CASE$I
SCHEME="$DDT1, $GRAD1, $DIVPHI1, $DIVALPHA1, $LAPLACIAN1, $INTERPOLATE1"

echo -n "$CASE, $SCHEME, "

cp -r damBreak $CASE

# DDT
if [ $DDT1 = "Euler" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "ddtSchemes['default']" Euler
fi

# GRAD
# オリジナルデータのフィールドを 4 つにしておく
if [ $GRAD1 = "linear" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][0]" Gauss
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][1]" linear
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "interpolationSchemes['default']" linear
elif [ $GRAD1 = "midLinear" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][0]" Gauss
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][1]" midLinear
elif [ $GRAD1 = "pointLinear" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][0]" Gauss
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][1]" pointLinear
elif [ $GRAD1 = "leastSquares" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][0]" leastSquares
elif [ $GRAD1 = "cLinear" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][0]" cellLimited
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][1]" Gauss
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][2]" linear
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][3]" 1
elif [ $GRAD1 = "cMidLinear" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][0]" cellLimited
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][1]" Gauss
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][2]" midLinear
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][3]" 1
elif [ $GRAD1 = "cPointLinear" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][0]" cellLimited
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][1]" Gauss
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][2]" pointLinear
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][3]" 1
elif [ $GRAD1 = "cLeastSquares" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][0]" cellLimited
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][1]" leastSquares
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][2]" 1
elif [ $GRAD1 = "fLinear" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][0]" faceLimited
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][1]" Gauss
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][2]" linear
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][3]" 1
elif [ $GRAD1 = "fMidLinear" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][0]" faceLimited
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][1]" Gauss
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][2]" midLinear
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][3]" 1
elif [ $GRAD1 = "fPointLinear" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][0]" faceLimited
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][1]" Gauss
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][2]" pointLinear
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][3]" 1
elif [ $GRAD1 = "fLeastSquares" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][0]" faceLimited
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][1]" leastSquares
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "gradSchemes['default'][2]" 1
fi

# DIVPHI1
if [ $DIVPHI1 = "limitedLinearV" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "divSchemes['div(rho*phi,U)'][0]" Gauss
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "divSchemes['div(rho*phi,U)'][1]" limitedLinearV
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "divSchemes['div(rho*phi,U)'][2]" 1
elif [ $DIVPHI1 = "upwind" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "divSchemes['div(rho*phi,U)'][0]" Gauss
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "divSchemes['div(rho*phi,U)'][1]" upwind
fi

# DIVALPHA
if [ $DIVALPHA1 = "vanLeer" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "divSchemes['div(phi,alpha)'][0]" Gauss
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "divSchemes['div(phi,alpha)'][1]" vanLeer
elif [ $DIVALPHA1 = "vanLeer01" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "divSchemes['div(phi,alpha)'][0]" Gauss
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "divSchemes['div(phi,alpha)'][1]" vanLeer01
fi

# LAPLACIAN1
if [ $LAPLACIAN1 = "corrected" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "laplacianSchemes['default'][0]" Gauss
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "laplacianSchemes['default'][1]" linear
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "laplacianSchemes['default'][2]" corrected
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "snGradSchemes['default']" corrected
elif [ $LAPLACIAN1 = "0.5" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "laplacianSchemes['default'][0]" Gauss
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "laplacianSchemes['default'][2]" limited
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "laplacianSchemes['default'][3]" 0.5
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "snGradSchemes['default'][0]" limited
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "snGradSchemes['default'][1]" 0.5
elif [ $LAPLACIAN1 = "0.333" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "laplacianSchemes['default'][0]" Gauss
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "laplacianSchemes['default'][1]" linear
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "laplacianSchemes['default'][2]" limited
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "laplacianSchemes['default'][3]" 0.333
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "snGradSchemes['default'][0]" limited
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "snGradSchemes['default'][1]" 0.333
elif [ $LAPLACIAN1 = "corrected" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "laplacianSchemes['default'][0]" Gauss
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "laplacianSchemes['default'][1]" linear
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "laplacianSchemes['default'][2]" uncorrected
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "snGradSchemes['default']" uncorrected
fi

# INTERPOLATE 
if [ $INTERPOLATE1 = "linear" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "interpolationSchemes['default']" linear
elif [ $INTERPOLATE1 = "midPoint" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "interpolationSchemes['default']" midPoint
elif [ $INTERPOLATE1 = "pointLinear" ] ; then
    pyFoamWriteDictionary.py $CASE/system/fvSchemes "interpolationSchemes['default']" pointLinear
fi

interFoam -case $CASE > $CASE/log

VEL=`sed -n '23,2290p' $CASE/1/U | sed -e 's/(//' | sed -e 's/)//' | awk 'BEGIN{max = 0} {v = sqrt($1^2 + $2^2 + $3^2); if(v > max) max = v} END{print max}'`

echo $VEL

R=`echo "$VEL < $MINVEL" | bc`
if [ $R -eq 1 ] ; then
    MINCASE=$CASE
    MINSCHEME=$SCHEME
    MINVEL=$VEL
fi

I=`expr $I + 1`

done
done
done
done
done
done

echo
echo "Min Case"
echo $MINCASE, $MINSCHEME, $MINVEL

damBreak (上の問題のケース) をコピーし、設定を変えて実行し、1 秒後の最大速度を取り出す、という操作を繰り返している。設定変更には pyFoam を利用。一部、パラメタの数が変わるので、前もって増やしておく。

system/fvSchemes

gradSchemes
{
    default         Gauss linear 0 0;
}

...

laplacianSchemes
{
    default         Gauss linear corrected 0;
}

実行。

$ ./run > log &

実行結果。

CASE1, Euler, linear, limitedLinearV, vanLeer, corrected, linear, 1.68949
CASE2, Euler, linear, limitedLinearV, vanLeer, corrected, midPoint, 1.68174
CASE3, Euler, linear, limitedLinearV, vanLeer, corrected, pointLinear, 2.13367
CASE4, Euler, linear, limitedLinearV, vanLeer, 0.5, linear, 0.393318
CASE5, Euler, linear, limitedLinearV, vanLeer, 0.5, midPoint, 0.521594

...

CASE572, Euler, fLeastSquares, upwind, vanLeer01, 0.333, midPoint, 0.210771
CASE573, Euler, fLeastSquares, upwind, vanLeer01, 0.333, pointLinear, 0.312346
CASE574, Euler, fLeastSquares, upwind, vanLeer01, uncorrected, linear, 0.27239
CASE575, Euler, fLeastSquares, upwind, vanLeer01, uncorrected, midPoint, 0.417624
CASE576, Euler, fLeastSquares, upwind, vanLeer01, uncorrected, pointLinear, 0.276111

Min Case
CASE365, Euler, cLeastSquares, upwind, vanLeer, 0.5, midPoint, 0.133464

1 秒後の最大速度がもっとも小さいものをよしとすると、CASE365 がよいことになる。ただし、upwind である。limitedLinearV の結果もほしい。limitedLinearV の結果だけを取り出す。

limitedLinear

#!/bin/sh

cat log | grep limitedLinearV | awk -F"," 'BEGIN {min=1e5} {if($8 < min){min=$8; line=$0}} END{print line}'

実行。

$ ./limitedLinear
CASE305, Euler, cPointLinear, limitedLinearV, vanLeer01, 0.5, midPoint, 0.143593

しかし、interpolateShcmes が midPoint なのが気になる。その前後を検討してみる。

CASE304, Euler, cPointLinear, limitedLinearV, vanLeer01, 0.5, linear, 0.159685
CASE305, Euler, cPointLinear, limitedLinearV, vanLeer01, 0.5, midPoint, 0.143593
CASE306, Euler, cPointLinear, limitedLinearV, vanLeer01, 0.5, pointLinear, 0.182687

この 3 ケースでもとの damBreak の計算をしてみたところ、1 秒後の挙動が 実験結果 に近そうなのは CASE306 か？ (特に左端の水面の挙動)。

計算結果。

はじめのケースよりも水面の動きが落ち着いて見える。

damBreak の結果。

左: オリジナル、右: CASE306。あまり変わっていない。

最適 fvSchemes (CASE306)

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         cellLimited Gauss pointLinear 1;
}

divSchemes
{
    div(rho*phi,U)  Gauss limitedLinearV 1;
    div(phi,alpha)  Gauss vanLeer01;
    div(phirb,alpha) Gauss interfaceCompression;
}

laplacianSchemes
{
    default         Gauss linear limited 0.5;
}

interpolationSchemes
{
    default         pointLinear;
}

snGradSchemes
{
    default         limited 0.5;
}

おまけ

ちなみに、つぎのような設定にすると、変な結果になる。

laplacianSchemes
{
    default         Gauss linear limited 0.5;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         linear;
}

水面が波打つ。laplacianSchemes と snGradSchemes のスキームを合わせたのはこのため。

追加検討

divSchemes で次のようなものがあるらしい。

div(rho*phi,U)  Gauss linearUpwindV grad(U);

damBreak で確認したかぎりではよさそうなので、検討してみる。

※"grad(U)" の設定はつぎのようにはできない。

pyFoamWriteDictionary.py $CASE/system/fvSchemes "divSchemes['div(rho*phi,U)'][2]" grad(U)

こうする。

pyFoamWriteDictionary.py $CASE/system/fvSchemes "divSchemes['div(rho*phi,U)'][2]" "grad(U)"

これだと、つぎのようになる。

div(rho*phi,U)  Gauss linearUpwindV "grad(U)";

これで問題ないみたい。

実行。

$ ./run > log &

実行結果。

CASE1, Euler, linear, linearUpwindV, vanLeer, corrected, linear, 1.71076
CASE2, Euler, linear, linearUpwindV, vanLeer, corrected, midPoint, 1.71112
CASE3, Euler, linear, linearUpwindV, vanLeer, corrected, pointLinear, 2.30832
CASE4, Euler, linear, linearUpwindV, vanLeer, 0.5, linear, 0.480488
CASE5, Euler, linear, linearUpwindV, vanLeer, 0.5, midPoint, 0.477912

...

CASE284, Euler, fLeastSquares, linearUpwindV, vanLeer01, 0.333, midPoint, 0.20723
CASE285, Euler, fLeastSquares, linearUpwindV, vanLeer01, 0.333, pointLinear, 0.35786
CASE286, Euler, fLeastSquares, linearUpwindV, vanLeer01, uncorrected, linear, 0.574206
CASE287, Euler, fLeastSquares, linearUpwindV, vanLeer01, uncorrected, midPoint, 0.507353
CASE288, Euler, fLeastSquares, linearUpwindV, vanLeer01, uncorrected, pointLinear, 0.324025

Min Case
CASE281, Euler, fLeastSquares, linearUpwindV, vanLeer01, 0.5, midPoint, 0.17293

faceLimited だと damBreak の結果がはげしくなまるので、cellLimited だけ抽出。

cellLimited

#!/bin/sh

cat log | awk -F"," '$3~/^ c/{print $0}' | awk -F"," 'BEGIN {min=1e5} {if($8 < min){min=$8; line=$0}} END{print line}'

実行。

$ ./limitedLinear
CASE149, Euler, cPointLinear, linearUpwindV, vanLeer, 0.5, midPoint, 0.187634

midPoint が気に入らないので、周辺のケースを見る。

CASE148, Euler, cPointLinear, linearUpwindV, vanLeer, 0.5, linear, 0.192668
CASE149, Euler, cPointLinear, linearUpwindV, vanLeer, 0.5, midPoint, 0.187634
CASE150, Euler, cPointLinear, linearUpwindV, vanLeer, 0.5, pointLinear, 0.191547

CASE150 がよさそう。

CASE150 の計算結果。

damBreak の結果。

左: オリジナル、右: CASE150。大きくは変わらない。

ちなみに、vanLeer01 だと以下の通りで、あまり違わない。

CASE160, Euler, cPointLinear, linearUpwindV, vanLeer01, 0.5, linear, 0.192312
CASE161, Euler, cPointLinear, linearUpwindV, vanLeer01, 0.5, midPoint, 0.188031
CASE162, Euler, cPointLinear, linearUpwindV, vanLeer01, 0.5, pointLinear, 0.192372

CASE162 の計算結果。

damBreak の結果。

左: オリジナル、右: CASE162。