OpenFOAMのinterFoam,codeStreamで断続的に水を注ぐシミュレーションにトライ

豆腐メンタルな豆腐

2024年8月7日 20:31

OpenFOAMで遊んだ内容を書き留めただけのチラシの裏です。

#OpenFOAM #オープンCAE

概要

私がOpenFOAMを覚えるために取り組んでいた"3 weeks"seriesのDay11の一部内容の紹介です。

今回のケースファイルおよび解説資料はこちらからダウンロード可能です。
Day11- Programming in OpenFOAM - Implementing boundary conditions using condeStream

下記のような箱の中に左側の壁から水を注いだり止めたりする内容になっています。

interFoamを使用(interFoamやVOF法についてはday9にて解説あり)
箱は1辺1.0mで上面は解放
重力はy方向-9.81[m/s^2]
箱には初期状態としてy <= 0.2m の範囲に水が溜まっている
注水口は上図左側の壁の 0.5m <= z <= 0.7m, 0.5m <= y <= 0.7m の範囲
シミュレーション時間は3.0秒で下記のように水を注いだり止めたりする
　0.0 ~ 1.0秒：注水 ⇒ 1.0 ~ 2.0秒：止水 ⇒ 2.0 ~ 3.0秒：注水

（全く同じだと面白くないので少しだけ値を変更しています）

この時間tに依存して変化する境界条件を作りこむのに、codeStream, codeFixedValueが使用されています。

いつも目にしている設定ファイルに直接コードを書き込むだけなので楽と言えば楽ですが、あまりC++経験のない私にはスラスラとAPIGuideを読むことができずしんどかったです。

とはいえ時間経過とともに水が出たり止まったりするのは見ていて楽しく、
特に印象に残った回だったので今回備忘も兼ねてこのメモにまとめてみました。

alpha.waterの設定

Day11- Programming in OpenFOAM - Implementing boundary conditions using condeStream
のスライド内容通りに設定ファイルを作成していきます。
（一部バージョン違いによる記述修正あり）

internalField    #codeStream
{
    code
    #{
        const IOdictionary& d = static_cast<const IOdictionary&>(dict);
        const fvMesh& mesh = refCast<const fvMesh>(d.db());
        scalarField alpha(mesh.nCells(), 0.0);

        forAll(alpha, i)
        {
            const scalar x = mesh.C()[i][0];
            const scalar y = mesh.C()[i][1];
            const scalar z = mesh.C()[i][2];
            if (y <= 0.2)
            {
                alpha[i] = 1.0;
            }
        }
        alpha.writeEntry("", os);
    #};

    codeInclude
    #{
        #include "fvCFD.H"
        #include <cmath>
        #include <iostream>
    #};

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

boundaryField
{
    leftWall
    {
        type            codedFixedValue;
        value           uniform 0;
        name            inletProfile_alpha_left;

        code
        #{
            const fvPatch& boundaryPatch = patch(); 
            const vectorField& Cf = boundaryPatch.Cf(); 
            scalarField& field = *this; 

            field = patchInternalField();

            scalar min = 0.5;  
            scalar max = 0.7;
            scalar t = this->db().time().value();

            forAll(Cf, faceI)
            {
                if (
                   (Cf[faceI].z() > min) &&
                   (Cf[faceI].z() < max) &&
                   (Cf[faceI].y() > min) &&
                   (Cf[faceI].y() < max) 
                   )
                {
                if ((t < 1.0) or (t > 2.0))
                    {
                        field[faceI] = 1.0;
                    }
                else
                    {
                        field[faceI] = 0.0;
                    }
                }
           }
        #};         

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

        codeInclude
        #{
                #include "fvCFD.H"
                #include <cmath>
                #include <iostream>
        #};
    }

    rightWall
    {
        type            zeroGradient;
    }

    frontWall
    {
        type            zeroGradient;
    }

    backWall
    {
        type            zeroGradient;
    }

    bottomWall
    {
        type            zeroGradient;
    }

    top
    {
        type            inletOutlet;
        inletValue      uniform 0;
        value           uniform 0;
    }
}

APIGuideや公式のプログラマズガイドより、fvMesh, fvPatch, vector・・・の各要素へのアクセス関数が確認できるので見てみます。

下記は公式プログラマズガイドからの抜粋です。

上記のinternalField設定では、forAllループ内でセル中心(.C())のy座標値が0.2m以下の座標のセルだけalpha=1.0としていることが確認できます。
またleftwall設定においては、境界面中心(.Cf())のy,z座標が0.5m ~ 0.7mの面のみ指定の時刻でalpha=1.0としていることも確認できます。

Uの設定

alphaと同様にUも参考資料の内容に沿って作成します。

同様にleftwall設定において境界面中心(.Cf())のy,z座標が0.5m ~ 0.7mの面のみ指定の時刻で流速が(1 0 0)となっています。

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    leftWall
    {
        type            codedFixedValue;
        value           uniform (0 0 0);
        name            inletProfile_u_left;

        code
        #{
            const fvPatch& boundaryPatch = patch(); 
            const vectorField& Cf = boundaryPatch.Cf(); 
            vectorField& field = *this; 
        
            scalar min = 0.5;  
            scalar max = 0.7;
            scalar t = this->db().time().value();

            forAll(Cf, faceI){
                if (
                   (Cf[faceI].z() > min) &&
                   (Cf[faceI].z() < max) &&
                   (Cf[faceI].y() > min) &&
                   (Cf[faceI].y() < max) 
                   )
                {
                if ((t < 1.0) or (t > 2.0))
                    {
                        field[faceI] = vector(1,0,0);
                    }
                else
                    {
                        field[faceI] = vector(0,0,0);
                    }
                }
            }
        #};         

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

        codeInclude
        #{
            #include "fvCFD.H"
            #include <cmath>
            #include <iostream>
        #};
    }

    rightWall
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    bottomWall
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    frontWall
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    backWall
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    top
    {
        type            pressureInletOutletVelocity;
        value           uniform (0 0 0);
    }
}

その他の設定

controlDictでendTimeを3秒に変更しておきます。

これ以外は特に変更する必要のあるファイルは無いはずなので、blockMeshでメッシュを作ってinterFoamで計算できるはずです。

application     interFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         3;

deltaT          0.001;

writeControl    adjustableRunTime;

writeInterval   0.05;

purgeWrite      0;

writeFormat     ascii;

writePrecision  8;

writeCompression uncompressed;

timeFormat      general;

timePrecision   6;

runTimeModifiable yes;

adjustTimeStep  yes;

maxCo           1;

maxAlphaCo      1;

maxDeltaT       1;

シミュレーション結果

paraviewでalpha=0.5の等値面の推移を可視化した結果動画です。
ちゃんと意図した通りに断続的に水が注がれている様子が確認できます。

やっぱり動画で結果が見えると楽しいです＾ｐ＾

参考資料

この記事が気に入ったらサポートをしてみませんか？