【OpenFOAMとOpenMoelica連携（２）】OpenFOAMとOpenMoelicaの連携：buoyantPimpleFoamソルバ

前回の記事ではOpenModelicaのコマンド実行の方法について紹介しました。

本記事は、OpenFOAMとOpenModelicaを連携させた計算方法を紹介する記事です。

OpenFOAMとOpenModelicaの連携にはFMUI規格に沿った、FMU4FOAMというツールがあります。
以下の記事を１つずつ読み進めていけば、FMU4FOAMを使うことができます。
この記事はOpenFOAMのv2412でも動かせるようにしたものです。
※並列計算できるかなどは確認していません。

あわせて読みたい

ただし、このFMU4FOAMは現時点（2025年11月）で調べた限りはOpenFOAMのマルチリージョンソルバには対応していないようでした。

これが少し不満でどうしてもマルチリージョンソルバに対応したOpenFOAMとOpenModelicaの連携がしたく、上記の勉強会でも報告しアドバイスを頂いましたが、解決にはいたらず。
マルチリージョンソルバは以下の図のように領域ごとに「0, constant, system」ファイルを構成するため、FMU4FOAMの境界条件がこのファイル構造を認識するようにプログラムされていないことが原因のようです。

今回はマルチリージョンソルバへの対応のために、OpenFOAMと外部ソルバ（今回はOpenModelica）との連携ができるプログラムを実装する手順を考えたので紹介します。

やりたいこと

OpenFOAMから外部ソルバ（今回はOpenModelida）を実行し、OpenFOAMとOpenModelicaの連携を行います。
本記事では、まずはOpenModelicaのコマンド実行ができるようにします。

計算の流れは以下のようになります。

この仕組みを使ってFMU4FOAMの例題（heatedRoom）と同じことを行いました。
右下のグラフに結果を載せています。

一番左のグラフが大事で、OpenFOAMのprobesで指定した温度センサの値がオレンジ色のターゲット温度に追従するよう、OpenModelicaのPID制御で入口温度境界に反映させる設定にしています。

OpenFOAMと外部ソルバ（OpenMoelica）を連携

以下がFMU4FOAMによる結果です。

FMU4FOAM

じゃっかん結果が違いますが、やりたいことは達成していますよね。
OpenFOAMの流出温度の評価方法とPID制御のパラメータが少し違うことが原因なのかな？

というわけで前段が長くなりましたが、本記事ではFMU4FOAMを使わずOpenFOAMとOpenModelicaを連携させる方法についてまとめておきます。

OpenFOAM v2412
OpenModelica 1.25.0（Linux版）
WSL2（Ubuntu-22.04）

OpenModelicaの計算実行
OpenFOAMとOpenModelicaの連携の構成図
計算毎時に実行する独自プログラム
1. 特定の境界面の温度、特定の場所の温度を取得

OpenModelicaの計算実行

前回の記事ではOpenModelicaをコマンド実行により計算する内容を紹介しました。

それぞれのファイルの関係性を図にしておきます。

これができるとものすごく便利ですよね。
なぜならOpenModelicaのGUIを起動することなく計算を実行できるため、数値を変更したパラスタやOpenFOAMの計算中に値をOpenModelicaに渡して計算を実行するなども可能になります。

OpenFOAMとOpenModelicaの連携の構成図

OpenFOAMとOpenModelicaを連成させる構成図を整理しました。

流れは以下です。

system/controlDictのcodeExecuteにより特定の境界面の温度、特定の場所の温度を取得します。
codeExecuteの途中でupdateMOS_and_Run.pyのPythonファイルを実行します。
このファイルはHVACSystem002.mosを作成、run.shを実行を行うファイルです。
run.shは前回の記事で説明したHVACSystem002.mosを実行するためのスクリプトです。
HVACSystem002.mosには流入境界の値が記述されたboundarySyncFileも出力しています。
boundarySyncFileに書かれた境界面の値を0/U, 0/Tファイルが読み込む

といった流れになります。

計算毎時に実行する独自プログラム

system/controlDictは計算毎時に独自のプログラムを実行するcodedという機能があります。
設定ファイルの全文を載せておきます。

system/controlDict

application     buoyantPimpleFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         300;

deltaT          1e-4;

writeControl    adjustableRunTime;

writeInterval   10;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

adjustTimeStep  yes;

maxCo           0.5;

// libs
// (
//     externalComm
// );

functions
{
    #include "outletTempAverageValue";
    fieldMinMax1
    {
        type fieldMinMax;
        libs      (fieldFunctionObjects );  // Mandatory entries (unmodifiable)// Mandatory entries (runtime modifiable)
        mode magnitude;
        fields   ( U  p T);  // Optional entries (runtime modifiable)
        location yes;
        writePrecision 8;
        writeToFile yes;
        useUserTime yes;
        region region0;
        enabled yes;
        log yes;
        timeStart 0;
        timeEnd 8000;
        executeControl timeStep;
        executeInterval 1;
        writeControl timeStep;//adjustableRunTime;
        writeInterval 1;
    }
    // FMUSimulator
    // {
    //     type            FMUSimulator;
    //     libs            (pyFMUSim);
    //     pyClassName     HVACSystem;
    //     pyFileName      HVACSystem;
    // }
}

application buoyantPimpleFoam;

startFrom latestTime;

startTime 0;

stopAt endTime;

endTime 300;

deltaT 1e-4;

writeControl adjustableRunTime;

writeInterval 10;

purgeWrite 0;

writeFormat ascii;

writePrecision 6;

writeCompression off;

timeFormat general;

timePrecision 6;

runTimeModifiable true;

adjustTimeStep yes;

maxCo 0.5;

// libs

// (

// externalComm

// );

functions

{

#include "outletTempAverageValue";

fieldMinMax1

{

type fieldMinMax;

libs (fieldFunctionObjects ); // Mandatory entries (unmodifiable)// Mandatory entries (runtime modifiable)

mode magnitude;

fields ( U p T); // Optional entries (runtime modifiable)

location yes;

writePrecision 8;

writeToFile yes;

useUserTime yes;

region region0;

enabled yes;

log yes;

timeStart 0;

timeEnd 8000;

executeControl timeStep;

executeInterval 1;

writeControl timeStep;//adjustableRunTime;

writeInterval 1;

}

// FMUSimulator

// {

// type FMUSimulator;

// libs (pyFMUSim);

// pyClassName HVACSystem;

// pyFileName HVACSystem;

// }

}

大事なのは#include "outletTempAverageValue";の部分です。
ここに独自のプログラムを書いておき、includeで読み込んでいます（直接記述すると煩雑になるためファイルを分けています）

特定の境界面の温度、特定の場所の温度を取得

outletTempAverageValueの全文を載せておきます。

outletTempAverageValue1
{
    type    coded;
    name    codeTempFieldValue;
    libs    ("libutilityFunctionObjects.so");

    // 重要: どのリージョンで動かすか
    // region  fluid;

    // 重要: いつ実行するか
    enabled         yes;
    executeControl  timeStep;
    executeInterval 1;

    codeExecute
    #{

        Info << "===== codeTempFieldValue (coded) start =====" << nl;

        // ---- ここから元の処理 ----
        const volVectorField& C = mesh().C();
        const scalarField&    V = mesh().V();

        word outletBoundaryName = "outlet";
        const volVectorField& U = mesh().lookupObject<volVectorField>("U");
        const volScalarField& T = mesh().lookupObject<volScalarField>("T");
    
        // targetNameのIDの取得
        const label outletID = mesh().boundaryMesh().findPatchID(outletBoundaryName);
         // targetNameのIDから速度場の取得
        const fvPatchVectorField& outletUList = U.boundaryField()[outletID];
         // targetNameのIDから温度場の取得
        const fvPatchScalarField& outletTList = T.boundaryField()[outletID];
        // targetNameのIDから面積ベクトルの取得
        const scalarField& outletmagSf = mesh().magSf().boundaryField()[outletID];


        // プローブ点の温度を取得
        scalar Tprobe = GREAT;                 // 初期化（並列用）
        scalar xs = 5;
        scalar ys = 5;
        scalar zs = 2.5;

        Info << " === probe Temperature === " << endl;
        // probeT = T.internalField()[mesh().findCell(C.findNearest(point(xs, ys, zs)))];
        Tprobe = T[ mesh().findCell(point(xs, ys, zs)) ];
        Info << mesh().findCell(point(xs, ys, zs)) << " Tprobe = " << Tprobe << endl;

        Info << " === probe Temperature end === " << endl;

        // ================================== Start WeightAreaAverageTemp==================================
        // Open the file to write the results
        std::ofstream outletWeightAreaTempFile; 
        outletWeightAreaTempFile.open("constant/outletWeightAreaTempFile");

        if (!outletWeightAreaTempFile.is_open())
        {
            FatalErrorIn("codeExecute") << "Cannot open outletWeightAreaTempFile file for writing!" << exit(FatalError);
        }

        scalar Toutlet;

        // MainTankはfloodとinteから
        Toutlet = (gSum(mag(outletUList)*outletmagSf*outletTList)  )
                        / (gSum(mag(outletUList)*outletmagSf) );

        // Write the header part of the file
        outletWeightAreaTempFile << "Tout," << Toutlet <<";"<< std::endl;
        outletWeightAreaTempFile << "Tprobe," << Tprobe <<";"<< std::endl;
        outletWeightAreaTempFile << "Time,"<< mesh().time().timeName()<<";" << std::endl;
        // Close the file
        outletWeightAreaTempFile.close();
        // ================================== End WeightAreaAverageTemp==================================
        
        // ================================== Python同期処理 Start ==================================
        Info << "===== Execute Python updateMOS_and_Run.py =====" << nl;

        // ロックファイル作成
        const fileName lockFile("constant/updateMOS.lock");
        {
            std::ofstream lock(lockFile.c_str());
            if (lock.is_open())
            {
                lock << "LOCKED" << std::endl;
                lock.close();
                Info << "Lock file created: " << lockFile << endl;
            }
            else
            {
                FatalErrorIn("codeExecute") << "Cannot create lock file!" << exit(FatalError);
            }
        }

        // Pythonスクリプトの実行（同期モードでも非同期でもOK、ここでは同期で実行）
        const fileName pyScript("updateMOS_and_Run.py");
        int ret = std::system(("python3 " + pyScript).c_str());
        if (ret != 0)
        {
            Info << "Warning: Python script exited with nonzero status: " << ret << endl;
        }

        // ロックファイルが削除されるまで待機
        while (Foam::isFile(lockFile))
        {
            Info << "Waiting for Python script to finish..." << endl;
            Foam::sleep(0.01);  // 1秒スリープ
        }

        Info << "Python script finished. Lock removed." << endl;
        // ================================== Python同期処理 End ================================
   
        #};
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

outletTempAverageValue1

{

type coded;

name codeTempFieldValue;

libs ("libutilityFunctionObjects.so");

// 重要: どのリージョンで動かすか

// region fluid;

// 重要: いつ実行するか

enabled yes;

executeControl timeStep;

executeInterval 1;

codeExecute

Info << "===== codeTempFieldValue (coded) start =====" << nl;

// ---- ここから元の処理 ----

const volVectorField& C = mesh().C();

const scalarField& V = mesh().V();

word outletBoundaryName = "outlet";

const volVectorField& U = mesh().lookupObject<volVectorField>("U");

const volScalarField& T = mesh().lookupObject<volScalarField>("T");

// targetNameのIDの取得

const label outletID = mesh().boundaryMesh().findPatchID(outletBoundaryName);

// targetNameのIDから速度場の取得

const fvPatchVectorField& outletUList = U.boundaryField()[outletID];

// targetNameのIDから温度場の取得

const fvPatchScalarField& outletTList = T.boundaryField()[outletID];

// targetNameのIDから面積ベクトルの取得

const scalarField& outletmagSf = mesh().magSf().boundaryField()[outletID];

// プローブ点の温度を取得

scalar Tprobe = GREAT; // 初期化（並列用）

scalar xs = 5;

scalar ys = 5;

scalar zs = 2.5;

Info << " === probe Temperature === " << endl;

// probeT = T.internalField()[mesh().findCell(C.findNearest(point(xs, ys, zs)))];

Tprobe = T[ mesh().findCell(point(xs, ys, zs)) ];

Info << mesh().findCell(point(xs, ys, zs)) << " Tprobe = " << Tprobe << endl;

Info << " === probe Temperature end === " << endl;

// ================================== Start WeightAreaAverageTemp==================================

// Open the file to write the results

std::ofstream outletWeightAreaTempFile;

outletWeightAreaTempFile.open("constant/outletWeightAreaTempFile");

if (!outletWeightAreaTempFile.is_open())

{

FatalErrorIn("codeExecute") << "Cannot open outletWeightAreaTempFile file for writing!" << exit(FatalError);

}

scalar Toutlet;

// MainTankはfloodとinteから

Toutlet = (gSum(mag(outletUList)*outletmagSf*outletTList) )

/ (gSum(mag(outletUList)*outletmagSf) );

// Write the header part of the file

outletWeightAreaTempFile << "Tout," << Toutlet <<";"<< std::endl;

outletWeightAreaTempFile << "Tprobe," << Tprobe <<";"<< std::endl;

outletWeightAreaTempFile << "Time,"<< mesh().time().timeName()<<";" << std::endl;

// Close the file

outletWeightAreaTempFile.close();

// ================================== End WeightAreaAverageTemp==================================

// ================================== Python同期処理 Start ==================================

Info << "===== Execute Python updateMOS_and_Run.py =====" << nl;

// ロックファイル作成

const fileName lockFile("constant/updateMOS.lock");

{

std::ofstream lock(lockFile.c_str());

if (lock.is_open())

{

lock << "LOCKED" << std::endl;

lock.close();

Info << "Lock file created: " << lockFile << endl;

}

else

{

FatalErrorIn("codeExecute") << "Cannot create lock file!" << exit(FatalError);

}

// Pythonスクリプトの実行（同期モードでも非同期でもOK、ここでは同期で実行）

const fileName pyScript("updateMOS_and_Run.py");

int ret = std::system(("python3 " + pyScript).c_str());

if (ret != 0)

{

Info << "Warning: Python script exited with nonzero status: " << ret << endl;

}

// ロックファイルが削除されるまで待機

while (Foam::isFile(lockFile))

{

Info << "Waiting for Python script to finish..." << endl;

Foam::sleep(0.01); // 1秒スリープ

}

Info << "Python script finished. Lock removed." << endl;

// ================================== Python同期処理 End ================================

#};

}

いろいろ書いていますが、やっているのは以下の事です。

流量で重みづけしたoutlet境界面の温度$\overline{T} = \frac{\sum_i \rho v_iA_i C_p T_i}{\sum_i \rho v_iA_i C_p }= \frac{\sum_i \rho Q_i C_p T_i}{\sum_i \rho Q_i C_p }$を求める。
座標$(5, 5, 2.5)$の温度を所得
それをconstant/outletWeightAreaTempFileファイルに記述
updateMOS_and_Run.pyのPythonファイルを実行
ただし、Pythonの実行が完全に終わるまで待機させる必要がある。事前にconstant/updateMOS.lockを作成しておき、Pythonファイルの実行が終わるとconstant/updateMOS.lockが削除される（この記述はupdateMOS_and_Run.pyに書かれている）ことでPythonの実行終了のフラグとしている。