OpenFOAMのCHTを強連成にしようと試みた話

初めに

　OpenFOAMには固液それぞれ複数領域ある場合の熱流体を計算する標準ソルバーが備わっている。
例：chtMultiRegionFoam、chtMultiRegionSimpleFoam
しかしこれらのエネルギー計算は一つの行列式で解くわけではないため、束性が落ちてしまうという特徴を持つ。（弱錬成）
過去にこれを解決しようと全体のエネルギー式を一つの行列で解こうとしたことがあるので紹介する。

• OpenFuelCellをとても参考にさせていただきました。
• 次回のOpenFOAM-v1712で強連成CHTが実装されるとの噂あり。
• OpenFOAM-extendのHookという関数を使えば複数の行列をマージできるらしい。

作ったもの

　以下のソルバーを作成した。
https://github.com/inabower/OpenFOAM-3.0.1-solvers/tree/master/chtGlobalEnergyPimpleFoam
※OpenFOAM-3.0.1で確認済み

例えば以下のようなメッシュを考える。

• yz平面で（10mm×20mm）×5領域。メッシュは1mm四方。
• 領域は"first"から"fifth"までの五つで全て固体。
• 熱伝導率kが異なる。
• 初期温度300Kで両端はfixedValue:400K。

これをchtMultiRegionFoamで計算すると全体の温度が両端の400Kに全体が近づいていく。
中心温度の時間変化を以下に示す。

弱連成の特徴として計算のタイムステップ⊿tに依って結果が変わる。
実際にこのソルバーを使う際には注意が必要である。

一方で、今回作ったものを赤線で重ねあわせると以下の通り。

タイムステップに依存しない。
※時間があったら机上計算してみます。

中身の説明

図解

　わかりにくいが、本来は各領域ごとに計算を完結させてから次の領域へ・・・といった順番で計算される。
今回のはまず流体領域のUを計算し、エネルギーの計算を一つの行列式で行ってから、流体のpを計算している。
この「全体のエネルギー計算」という部分で以下のようにひと工夫している。

事前に 全体が結合したメッシュを用意しておき、

①[流体]各領域のkappa、rho、Cpと速度計算の結果からK、Ug、rhoPhi、rhoCpPhiを全体領域に転写。
②[個体]各領域のkappa、rho、Cpを全体領域に転写。
③転写された値から全体領域でエネルギー計算を行う。
④[流体]全体領域で得られたeを流体領域に転写
⑤[個体]全体領域で得られたeを個体領域に転写

といった流れになっている。
以下では全体メッシュの作成と転写の方法について説明する。

splitMultiRegion

　このソルバーにはsplitMultiRegionで領域分割することが不可欠となる。（自分で各Regionを作ったcaseは計算できない。）
例えばメッシュにsplitMultiRegionを実行すると各Regionのメッシュが生成される。
そのケースディレクトリにはsplitMultiRegionを使わなかった場合には存在しないファイルがいくつかある。

• 分割前のメッシュ（当たり前だが）
• constant/(リージョン名)/polyMesh/pointRegionAdressing
• constant/(リージョン名)/polyMesh/faceRegionAdressing
• constant/(リージョン名)/polyMesh/boundaryRegionAdressing
• constant/(リージョン名)/polyMesh/cellRegionAdressing

この*RegionAdressingにはある点/面/セルが元のメッシュではどの点/面/セルだったかを示している。
これを利用して分割前のメッシュを"全体のメッシュ"として利用しながら計算を行う。

chtGlobalEnergyPimpleFoam

ソルバーの中身にも少し触れる。
chtMultiRegionFoamをベースにそれと異なる点についてかいつまんで説明する。

全体メッシュの作成

今回は全体メッシュのことをcellRegionと呼称している。
複数リージョンに加えて、更にメッシュを読み込む。

chtGrobalEnergyPimpleFoam/AllCell/createCellMeshes.H

Foam::fvMesh cellRegion
(
    Foam::IOobject
    (
        Foam::fvMesh::defaultRegion,
        runTime.timeName(),
        runTime,
        Foam::IOobject::MUST_READ
    )
);

全体メッシュとのリンク

先ほどの*Adressingを利用する。
例えば流体領域のfaceRegionAdressingは

chtGrobalEnergyPimpleFoam/fluid/createFluidMeshes.H.33

PtrList<labelIOList> fluidFaceRegionAddressing(fluidNames.size());
 
forAll(fluidNames, i)
{
    fluidFaceRegionAddressing.set
    (
        i,
        new labelIOList
        (
            IOobject
            (
                "faceRegionAddressing",
                fluidRegions[i].time().findInstance
                (
                    fluidRegions[i].meshDir(),
                    "faces"
                ),
                polyMesh::meshSubDir,
                fluidRegions[i],
                IOobject::MUST_READ
            )
        )
    );
}

と読み込んでいる。
ただややこしい事にこのfaceRegionAddressingに入っている値は
("対応する面番号" + 1) × ( 1 OR -1 )"その方向が同じかどうか"
となっている。
ここから2つの情報を取り出すために以下の処理を行っている。

chtGrobalEnergyPimpleFoam/fluid/createFluidMeshes.H

PtrList<labelList> fluidFaceMap(fluidNames.size());
PtrList<scalarField> fluidFaceMask(fluidNames.size());
 
forAll(fluidNames, i)
{
    // まずPtrListの中に入るlabelListの要素数をラベルと一緒に指定して
    fluidFaceMap.set
    (
        i,
        new labelList
        (
            fluidFaceRegionAddressing[i].size()
        )
    );
 
    // こっちはPtrListの中に入るscalarFieldの要素数をラベルと一緒に指定して
    fluidFaceMask.set
    (
        i,
        new scalarField
        (
            fluidFaceRegionAddressing[i].size()
        )
    );
    
    // labelListとscalarFieldの値を入力する
    forAll(fluidFaceMap[i], j)
    {
        fluidFaceMap[i][j] = mag(fluidFaceRegionAddressing[i][j]) - 1; // 絶対値-1
        fluidFaceMask[i][j] = sign(fluidFaceRegionAddressing[i][j]); // 符号の取得
    }
}

これを用いて例えば流体領域のphiを全体に転写するときは以下のようになっている。

C++:chtGrobalEnergyPimpleFoam/AllCell/mapFluidToCell.H

surfaceScalarField& phiC = phiCell; // 全体のphi
 
forAll(fluidRegions, i)
{
    const surfaceScalarField& phiF = phiFluid[i]; // 領域iのphi
 
    forAll(phiF, faceI) // 全ての面についてループ
    {
        const label faceJ = fluidFaceMap[i][faceI];       // 対応する面の場所
        const scalarField maskJ = fluidFaceMask[i][faceI];  // 面の符号
 
        phiC[faceJ] = phiF[faceI] * maskJ;
    }
}

cellとboundaryに関しては特にいじる必要がないので

chtGrobalEnergyPimpleFoam/fluid/createFluidMeshes.H

PtrList<labelIOList> fluidCellMap(fluidNames.size());
PtrList<labelIOList> fluidPatchMap(fluidNames.size());
 
forAll(fluidNames, i)
{
    fluidCellMap.set
    (
        i,
        new labelIOList
        (
            IOobject
            (
                "cellRegionAddressing",
                fluidRegions[i].time().findInstance
                (
                    fluidRegions[i].meshDir(),
                    "faces"
                ),
                polyMesh::meshSubDir,
                fluidRegions[i],
                IOobject::MUST_READ
            )
        )
    );
 
    fluidPatchMap.set
    (
        i,
        new labelIOList
        (
            IOobject
            (
                "boundaryRegionAddressing",
                fluidRegions[i].time().findInstance
                (
                    fluidRegions[i].meshDir(),
                    "faces"
                ),
                polyMesh::meshSubDir,
                fluidRegions[i],
                IOobject::MUST_READ
            )
        )
    );
}

と読み込んで、

C++:chtGrobalEnergyPimpleFoam/AllCell/mapFluidToCell.H

volScalarField& TC = TCell;
 
forAll(fluidRegions, i)
{
    const volScalarField& TF = thermoFluid[i].T();
 
    forAll(TF, cellI)
    {
        const label cellJ = fluidCellMap[i][cellI];
        TC[cellJ] = TF[cellI];
    }
}

と使う。

問題点

並列計算できない。

Adressingファイルはconstantに入っておりdecomposeParで分割されないため対応がちぐはぐになってしまう。

速度はむしろ遅い。

メッシュ間を値が大量にコピーされるためかと思われる。

乱流計算を組み込んでない。

転写する項目を増やせば可能。

最後に

　とても中途半端で終わりますが皆さんの開発の中で何か参考になる部分があれば幸いです。