マルチリージョン

2015年2月18日

はじめに

simpleFoam をマルチリージョン化する。

使用バージョン

OpenFOAM 2.3.1

ファイル

multiRegionSimpleFoam (マルチリージョン版 simpleFoam)

multiRegionBuoyantSimpleFoam (マルチリージョン版 buoyantSimpleFoam)

chtMultiRegionBuoyantSimpleFoam (CHT 版 buoyantSimpleFoam)

chtMultiRegionBuoyantBoussinesqSimpleFoam (CHT 版 buoyantBoussinesqSimpleFoam)

マルチリージョンの概要

chtMultiRegionFoam などのマルチリージョン (複数領域) を扱うソルバーは、基本的には領域ごとに単一領域のソルバー同様の方程式を解いているだけである。

具体的には、単一領域のソルバーの createFields.H で作る変数をリストにしてそれぞれの領域のデータを格納し、そのリストのループを回してそれぞれの領域の方程式を解く。その際、速度のリストをたとえば UFluid として、U = UFluid[i] (i はリストのインデックス) などとして各領域の要素を名前付けしてやれば、単一領域のソルバーの UEqn.H などをほぼ使いまわすことができる。

simpleFoam のマルチリージョン化

Make/options に以下を追加。

EXE_INC = \
    ...
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude

EXE_LIBS = \
    ...
    -lregionModels

インクルードファイル "regionProperties.H" の取り込みを追加する。

#include "regionProperties.H"

そして、"createTime.H" のインクルードの後に regionProperties の定義を追加する。

regionProperties rp(runTime);

これは、constant/regionProperties を読み込む。

"createMesh.H" の代りに "createFluidMeshes.H" をインクルードする。

#include "createFluidMeshes.H"

これは各領域のメッシュを読み込む。

createFluidMeshes.H

const wordList fluidNames(rp["fluid"]);

PtrList<fvMesh> fluidRegions(fluidNames.size());

forAll(fluidNames, i)
{
    Info<< "Create fluid mesh for region " << fluidNames[i]
        << " for time = " << runTime.timeName() << nl << endl;

    fluidRegions.set
    (
        i,
        new fvMesh
        (
            IOobject
            (
                fluidNames[i],
                runTime.timeName(),
                runTime,
                IOobject::MUST_READ
            )
        )
    );
}

"createFields.H" などの代りに "createFluidFields.H" をインクルードする。

#include "createFluidFields.H"

これはおおむね "createFields.H" をリスト化したものである。まず、"createFields.H" の名前を "createFluidFields.H" に変える。そして全体を fluidRegions のループで囲む。

forAll(fluidRegions, i)
{
    ...
}

各変数をリストにする。ループの前でリストを定義する。

PtrList<volScalarField> pFluid(fluidRegions.size());

インスタンスの作成をリストへのセットに書き換える。

        pFluid.set
        (
            i,
            new volScalarField
            (
                IOobject
                (
                    "p",
                    runTime.timeName(),
                    fluidRegions[i],
                    IOobject::MUST_READ,
                    IOobject::AUTO_WRITE
                ),
                fluidRegions[i]
            )
        );

mesh は fluidRegions[i] に置き換える (mesh = fluidRegions[i] としたほうが早かったか)。

"createPhi.H" などをインクルードしている箇所も、もとのコードをコピーして書き換える。

New() を使っているものは、New.ptr() とする (例外もあるようだが)。

        turbulence.set
        (
            i,
            incompressible::RASModel::New
            (
                UFluid[i],
                phiFluid[i],
                laminarTransport[i]
            ).ptr()
        );

さて、メインのコードに戻る。SIMPLE ループを runTime のループに書き換える。

while (runTime.loop())
{
    ...

ループの中身を runTime.write() の前まで fluidRegions のループで囲む。

forAll(fluidRegions, i)
{
    ...
}

各領域のデータ (単一領域ソルバーの "createFields.H" で用意していたもの) をリストから取り出す。

            fvMesh& mesh = fluidRegions[i];
            volScalarField& p = pFluid[i];
            volVectorField& U = UFluid[i];
            surfaceScalarField& phi = phiFluid[i];

            incompressible::RASModel& turb = turbulence[i];

            fv::IOoptionList& fvOptions = fluidFvOptions[i];

            const label pRefCell = pRefCellFluid[i];
            const scalar pRefValue = pRefValueFluid[i];

このあとぐらいで simpleControl を用意する。

            simpleControl simple(mesh);

simple.loop() を用いていないので、前イテレーションの値を保存したりはされないので、自分で追加する必要がある。

            {
                p.storePrevIter();

                #include "UEqn.H"
                #include "pEqn.H"
            }

"turbulence->" を "turb." に書き換える ("UEqn.H" なども)。

	    turb.correct();

インクルードしている "initContinuityErrs.H" と "continuityErrs.H" の cumulativeContErr をリストにする (オリジナルのコードからコピーしてくる。

List<scalar> cumulativeContErr(fluidRegions.size(), 0.0);

収束判定

収束判定の機構がないので、追加する。SIMPLE ループ (simple.loop()) は収束したら勝手に止まるので、判定自体は simpleControl 内で行われていることになる。収束したかどうかは、simpleControl のメンバ関数 criteriaSatisfied() が返すが、これは public ではないので、これを使うために simpleControl をコピーして criteriaSatisfied() を public に移す。そのままコードをコピーしただけだとオリジナルとダブっていると見えるようなので、"simpleControl.C" の defineTypeNameAndDebug() の部分をコメントにする。

//namespace Foam
//{
//    defineTypeNameAndDebug(simpleControl, 0);
//}

さて、収束は全領域で判定する必要があるので、少し工夫がいる。まず、ステップループの頭でブール変数 convergence を定義する。

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        bool convergence = true;

領域ループの最後で次のようにする。

        forAll(fluidRegions, i)
        {
            ...

            bool regionConvergence = simple.criteriaSatisfied();

            if (regionConvergence)
            {
                Info<< nl << "Region " << fluidRegions[i].name()
                    << " convergence criteria" << nl << endl;
            }

            convergence = convergence && regionConvergence;
        }

これで全領域が収束しないかぎり convergence は true にならない。ステップループの最後に収束判定を入れる。

        if (convergence)
        {
            Info<< nl << "SIMPLE solution converged in "
                << runTime.timeName() << " iterations" << nl << endl;

            runTime.writeAndEnd();
        }

マルチリージョン化した simpleFoam のテスト

上の図は、icoFoam のチュートリアル cavity のようなケースを 2 つの領域を作って解いてみた。それぞれの領域の上面の速度は変えてある。それぞれ別々に解けているのが確認できる。

buoyantSimpleFoam のマルチリージョン化

ついでに buoyantSimpleFoam もマルチリージョン化してみた。上の図は、温度分布を与えたもので、それぞれの領域の初期温度と境界の温度を変えている。

上の図のように、領域界面の境界条件を設定すれば、熱伝達も表現できる。

固体も含むマルチリージョン化

固体と流体の熱連成のために、multiRegionBuoyantSimpleFoam について固体も含んだマルチリージョン化を行う。固体領域の計算は chtMultiRegionFoam/chtMultiRegionSimpleFoam のコードを流用する。

fluid ディレクトリを作り、その中に流体用のコードを移動する。

solid ディレクトリを作り、chtMultiRegionFoam/solid、chtMultiRegionSimpleFoam/solid から以下のファイルをコピーする。

  • createSolidMeshes.H
  • createSolidFields.H
  • setRegionSolidFields.H
  • readSolidMultiRegionSIMPLEControls.H (chtMultiRegionSimpleFoam)
  • solveSolid.H (chtMultiRegionSimpleFoam)

Make/options で fluid、solid ディレクトリをパスに追加する。

EXE_INC = \
    -I./fluid \
    -I./solid \
    ...

以下をはじめのインクルードに追加。

#include "coordinateSystem.H"
#include "solidThermo.H"

Make/options に以下を追加。

EXE_INC = \
    ...
    -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude

EXE_LIBS = \
    ...
    -lsolidThermo

流体領域のメッシュ、フィールドの準備の後に、固体領域のメッシュ、フィールドの準備をするコードを挿入する。

    #include "createSolidMeshes.H"
    #include "createSolidFields.H"

流体領域の計算の後に、固体領域の計算を行うコードを追加する。

        forAll(solidRegions, i)
        {
            Info<< "\nSolving for solid region "
                << solidRegions[i].name() << endl;
	    
            #include "setRegionSolidFields.H"
            #include "readSolidMultiRegionSIMPLEControls.H"
            #include "solveSolid.H"
        }

右の領域を固体に変えて計算してみた。

buoyantBoussinesqSimpleFoam のマルチリージョン化

buoyantBoussinesqSimpleFoam を固体も含めてマルチリージョン化する。問題は、buoyantBoussinesqSimpleFoam は非圧縮性流体を想定しているが、領域間の熱伝達の境界条件が圧縮性の熱物性を要求することである。そこで、熱伝達境界条件を用いるために、温度のフィールド T のためだけに rhoThermo を用いることにする。

volScalarField& T = thermo.T();

この方法で起りそうな問題は、非圧縮性の物性と圧縮性の物性との整合が取れなくなることだろう。双方の熱拡散率が合うように気を付けておけばよいのだろう。