Notched Disc In Rotating Flow

はじめに

ESI版OpenFOAM v2406を対象にして、タイトルに記した計算を行う。
この計算は、回転する流速場にあるノッチ付き円盤状の水(以下、$\alpha$と記す)が一周してどのくらい初期の形状を保てるか、というものである。
ソルバはinterIsoFoam/interFoamを使用し、$\alpha$の移流のみ計算する。
ESI版はinterIsoFoamのチュートリアルに本問題があり、それを使用する。
Foundation版OpenFOAM 11でも計算を行う。

解析モデル

設定値一覧

名前	設定値
ドメイン長さ $L [m]$	$1$
ディスク半径 $R [m]$	$0.15$
ノッチ幅 $W [m]$	$0.05$
ノッチ高さ $H [m]$	$0.25$
ディスク初期位置 $(x_0,y_0,z_0)$	$(0.5, 0, 0.75)$

回転する流速場

ESI版OpenFOAMでは、function objectで流速場の初期設定ができる。

速度場$U(u,0,v)$の分布は次式の通りである。

u = \frac{\pi}{3.14} (0.5 - z)

v = \frac{\pi}{3.14} (x - 0.5)

functions
{
    setVelocity
    {
        type            setFlow;
        libs            (fieldFunctionObjects);
        writeControl    writeTime;
        mode            rotation;
        scale           1;
        //reverseTime     1;
        omega           1.000507215; // 2*pi/(2*3.14)
        //omega           6.28318530718;
        origin          (0.5 0 0.5);
        refDir          (1 0 0);
        axis            (0 1 0);
    }
}

system/controlDictのオリジナルからreverseTimeをコメントにし、omegaを$\frac{2\pi}{2✕3.14}$に変更する。つまり、ノッチ付き円盤の周期は、$2✕3.14=6.28 [s]$である。これを計算終了時刻とする。

ノッチ付き円盤の動き

評価式

計算開始時と計算終了時の$\alpha$の分布より、次の誤差を計算する。

\epsilon_{ND2} = \frac{1}{L} \sum_{j=1}^{N} A_{j} |\alpha_{J}^{I} - \alpha_{j}^{F}|

変数一覧

変数名	説明
$L$	界面の長さ $[m]$
$A_j$	$セルjの面積 [m^2]$
$\alpha_j^I$	$計算開始時のセルjの\alpha$
$\alpha_j^F$	$計算終了時のセルjの\alpha$

$L$は計算開始時の界面の長さとすると、約$1.442 [m]$である。$A_j$はセルの分割数から算出した値とする。文献では$\alpha_j^I$を1タイムステップ実行後の値としているが、ここでは初期値とする。

ESI版OpenFOAMの計算結果

誤差$\epsilon_{ND2}$の計算結果

	ケース1	ケース2	ケース3
分割数	$100✕100$	$200✕200$	$400✕400$
$interIsoFoam$	$0.001 97$	$0.000 51$	$0.000 24$
$interFoam$	$0.003 42$	$0.001 39$	$0.000 75$

なお、interIsoFoamのreconstructSchemeはisoAlphaである。

Foundation版の回転する速度場

Foundation版では、0/UのinternalFieldsにcodeStreamを適用することで、初期設定ができる。

internalField   #codeStream
{
    codeInclude
    #{
        //#include "fvCFD.H"
        #include "volFields.H"
    #};

    codeOptions
    #{
        -I$(LIB_SRC)/finiteVolume/lnInclude \
        -I$(LIB_SRC)/meshTools/lnInclude
    #};

    codeLibs
    #{
        -lfiniteVolume \
        -lmeshTools
    #};

    code
    #{
        const IOdictionary& d = static_cast<const IOdictionary&>(dict);
        const fvMesh& mesh = refCast<const fvMesh>(d.db());

        vectorField U_0(mesh.nCells(), Foam::Vector<double>(0.0,0.0,0.0));

        forAll(U_0, i)
        {
            const scalar x = mesh.C()[i][0];
            const scalar z = mesh.C()[i][2];
            const scalar u = constant::mathematical::pi * (0.5 - z) / 3.14;
            const scalar v = constant::mathematical::pi * (x - 0.5) / 3.14;
            U_0[i] = Foam::Vector<double>(u,0.0,v);
        }

        // ESI style U_0.writeEntry("", os);
        writeEntry(os, "", U_0);
    #};
};

バージョン11以前では、codeIncludeで"fvCFD.H"をインクルードする。
また、ESI版で使う場合は、writeEntryの書き方を変える。
このコードは実行時にコンパイルされ、速度場は次のようになる。

この速度場は、ESI版と同じ（のはず）である。

Foundation版OpenFOAMの計算結果

Foundation版のinterFoamでは、$\alpha$の発散スキームとして次の3つが用意されている。

誤差$\epsilon_{ND2}$の計算結果

	ケース1	ケース2	ケース3
分割数	$100✕100$	$200✕200$	$400✕400$
$interfaceCompression$	$0.00411$	$0.00177$	$0.00097$
$PLIC$	$0.00449$	$0.00177$	$0.00097$
$MPLIC$	$0.00449$	$0.00210$	$0.00133$

まとめ

本問題に関しては、ESI版の方が誤差が小さく、計算の準備が容易である。