【OpenFOAMとOpenMoelica連携（３）】chtMultiRegionFoamによるマルチリージョンソルバ

OpenFOAMの境界条件
output.dat
1. 注意

OpenFOAMの境界条件

0/fluid/Uと0/fluid/Tの全文を載せときます。

0/UのinletにOpenModelicaで計算した質量流moutを設定
0/TのinletにOpenModelicaで計算したToutを設定

マルチリージョンソルバは各領域に対して境界条件を設定する必要があることに注意してください。

0/fluid/U

#include "<case>/constant/inletTempFile";
dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
     #includeEtc "caseDicts/setConstraintTypes"

    "f-inlet(XMin|XMax|YMin|YMax)"
    {
        type            codedFixedValue;
        name            pythonSyncFixedValue;
        value           uniform (0 0 0);  // ダミー

        codeInclude
        #{
            #include "fvCFD.H"
        #};

        code
        #{
            // boundarySyncFile を登録（初回のみ）
            if (!db().foundObject<IOdictionary>("boundarySyncFile"))
            {
                autoPtr<IOdictionary> pythonSyncPtr
                (
                    new IOdictionary
                    (
                        IOobject
                        (
                            "boundarySyncFile",
                            db().time().constant(),
                            db(),
                            IOobject::MUST_READ,
                            IOobject::NO_WRITE
                        )
                    )
                );
                db().store<IOdictionary>(pythonSyncPtr);
            }

            IOdictionary& pythonSync =
                db().lookupObjectRef<IOdictionary>("boundarySyncFile");
            pythonSync.regIOobject::read();

            // mout [kg/s]
            const scalar moutSync = readScalar(pythonSync.lookup("mout"));

            const fvPatch& p = patch();
            const scalar A   = gSum(p.magSf());      // パッチ全体面積 [m2]

            const scalar rho = 1.2;                 // 空気密度（簡易一定）

            // 法線方向速度 [m/s] （法線が外向きなので、流入にしたいならマイナス）
            const scalar Un = -moutSync / (rho * A);

            // 法線ベクトル
            const vectorField nHat = p.Sf()/p.magSf();
            vectorField Ubc(p.size());

            forAll(Ubc, faceI)
            {
                Ubc[faceI] = Un * nHat[faceI];
            }

            Info<< "pythonSyncFixedValue: time = " << db().time().value()
                << " patch = " << p.name()
                << " A = " << A << " [m2]"
                << " mout = " << moutSync << " [kg/s]"
                << " Un = " << Un << " [m/s]"
                << nl << endl;

            operator==(Ubc);
        #};
    }

    // f-inletXMin
    // {
    //     type            fixedValue;
    //     value           uniform (0 0 -1.0);
    // }
    f-outlet
    {
        type pressureInletOutletVelocity;
        value uniform (0 0 0);
    }
    // f-inletXMax
    // {
    //     type            fixedValue;
    //     value           uniform (0 0 -1.0);
    // }
    // f-inletYMax
    // {
    //     type            fixedValue;
    //     value           uniform (0 0 -1.0);
    // }
    f-walls
    {
        type            noSlip;
    }
    f-inwalls
    {
        type            noSlip;
    }

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

OpenModelicaでは質量流を計算しているため、OpenFOAMの境界条件で流速条件に変えています。

空気の密度$1.2\, [\text{kg}/\text{m}^3]$、流入面積$0.06\, [\text{m}^2]$となっているかを計算中のログで出力するようにしています。

0/fluid/T

#include "<case>/constant/inletTempFile";
dimensions      [0 0 0 1 0 0 0];

internalField   uniform 293.15;//初期値

boundaryField
{
     #includeEtc "caseDicts/setConstraintTypes"

    "f-inlet(XMin|XMax|YMin|YMax)"
    {
        // type            fixedValue;
        // value           uniform 303.15; //流入温度(30deg.)
        type            codedFixedValue;
        name            pythonSyncFixedValue;
        value           uniform $Tout;
        code
        #{
            // pythonSyncファイルをobjectRegistryへ登録（初回のみ）
            if(!db().foundObject<IOdictionary>("boundarySyncFile"))
            {
                autoPtr<IOdictionary> pythonSyncPtr
                (
                    new IOdictionary
                    (
                        IOobject
                        (
                            "boundarySyncFile",
                            db().time().constant(),
                            db(),
                            IOobject::MUST_READ,
                            IOobject::NO_WRITE
                        )
                    )
                );
                db().store<IOdictionary>(pythonSyncPtr);
            }

            // objectRegistry上のpythonSyncを呼び出す
            IOdictionary& pythonSync = db().lookupObjectRef<IOdictionary>("boundarySyncFile");
            // 内容の更新
            pythonSync.regIOobject::read();
            // pythonSyncファイルからUColdSyncの値を取得
            const scalar& ToutSync = pythonSync.get<scalar>("Tout");
            // fixedValueの値を渡す
            operator==(ToutSync);
        #};
    }

    // f-inletXMin
    // {
    //     type            fixedValue;
    //     value           uniform 303.15; //流入温度(30deg.)
    // }
    f-outlet
    {
        type            inletOutlet;
        inletValue      $internalField;
        value           $internalField;
    }

    // f-inletYMax
    // {
    //     type            fixedValue;
    //     value           uniform 303.15; //流入温度(30deg.)
    // }
    f-walls
    {
        type            zeroGradient; //断熱条件
    }
    f-inwalls
    {
        type            zeroGradient; //断熱条件
    }
    fluid_to_solid
    {
        type            compressible::turbulentTemperatureRadCoupledMixed;
        Tnbr            T;
        kappaMethod     fluidThermo;
        value           uniform 293.15;
    }
}

consant/boundarySyncFileファイルに記述されたmoutとToutを計算毎に読み込む記述が書かれています。

あまり重要ではないですが初期条件は#include "<case>/constant/inletTempFile";として別ファイルにしています。

constant/inletTempFile

FoamFile
{
 version 2.0;
 format ascii;
 class dictionary;
 location constant;
 object inletTempFile;
}

mout	0;
Tout	293.15
Time	0.0;

output.dat

output.datはupdateMOS_and_Run.pyの中に書かれたOpenModelicaの結果ファイルです。以下のように書き出されます。

output.dat

Time	mout	Tout	heatFlowSensor.Q_flow	TSensor	ramp.y
0.0011999	2.88	293.266	244256	293.15	298
0.00263965	2.88	293.405	243839	293.15	298
0.00436705	2.88	293.572	243327	293.156	298
0.00643992	2.88	293.777	242644	293.196	298
0.00892676	2.88	294.031	241771	293.269	298
0.0118421	2.88	294.327	240852	293.305	298
0.0151831	2.88	294.658	239894	293.302	298
0.0189567	2.88	295.044	238741	293.332	298
0.0231601	2.88	295.496	237376	293.403	298
0.027432	2.88	295.966	235999	293.47	298
0.0317261	2.88	296.442	234644	293.524	298
...（省略）...

output.datの出力はupdateMOS_and_Run.pyからimportしているsync_to_table.pyに書かれています。

sync_to_table.py

from pathlib import Path
from typing import Dict, List, Tuple

def _strip_semicolon(s: str) -> str:
    return s[:-1] if s.endswith(";") else s

def parse_boundary_sync(sync_path: Path) -> Dict[str, float]:
    """
    constant/boundarySyncFile を読み込み、{key: value} を返す。
    FoamFile ブロックを無視し、 'Key  Value;' 形式に対応。
    """
    vals: Dict[str, float] = {}
    in_foam_block = False

    with sync_path.open(encoding="utf-8") as f:
        for raw in f:
            line = raw.strip()
            if not line:
                continue

            if line.startswith("FoamFile"):
                in_foam_block = True
                continue
            if in_foam_block:
                if "}" in line:
                    in_foam_block = False
                continue

            line = _strip_semicolon(line)
            parts = line.split()
            if len(parts) < 2:
                continue

            key, val_str = parts[0], parts
            try:
                val: float | str = float(val_str)
            except ValueError:
                val = val_str
            vals[key] = val

    return vals


def _fmt(v):
    if isinstance(v, (float, int)):
        return f"{v:g}"
    return "" if v is None else str(v)


def _read_existing_header(path: Path, sep: str) -> List[str] | None:
    if not path.exists() or path.stat().st_size == 0:
        return None
    with path.open("r", encoding="utf-8") as f:
        first = f.readline().rstrip("\n")
    return first.split(sep)


def _append_missing_columns_in_file(path: Path, sep: str, existing: List[str], target: List[str]) -> None:
    """
    既存ヘッダ existing に対し、target にあるが existing に無い列を末尾追加。
    既存データ行には空欄を付加して列数を合わせる。
    """
    missing = [c for c in target if c not in existing]
    if not missing:
        return

    new_header = existing + missing
    # 全行を読み込み
    with path.open("r", encoding="utf-8") as f:
        lines = f.readlines()

    # 先頭行を新ヘッダへ
    lines[0] = sep.join(new_header) + "\n"

    # データ行を末尾に空欄を追加
    for i in range(1, len(lines)):
        lines[i] = lines[i].rstrip("\n") + (sep + sep.join([""] * len(missing)) if missing else "") + "\n"

    with path.open("w", encoding="utf-8") as f:
        f.writelines(lines)


def _ensure_header(out_path: Path, headers: List[str], sep: str) -> List[str]:
    """
    output.dat のヘッダ整合性を保証:
    - ファイルが無ければ headers を書く
    - あれば既存ヘッダを読み、不足列を末尾に追加して既存行を補完
    - 返り値: 実際のヘッダ（既存＋不足列）
    """
    existing = _read_existing_header(out_path, sep)
    if existing is None:
        out_path.write_text(sep.join(headers) + "\n", encoding="utf-8")
        return headers

    # 既存に不足があれば追加
    _append_missing_columns_in_file(out_path, sep, existing, headers)
    return _read_existing_header(out_path, sep) or headers


def append_output_row(
    base_dir: Path,
    vals: Dict[str, float],
    headers: List[str] | None = None,
    out_name: str = "output.dat",
    sep: str = "\t",
) -> Tuple[Path, str]:
    """
    vals を out_name に追記。列は headers の順序で出力。
    既存ファイルがあれば自動的に不足列を追加し整合させる。
    """
    if headers is None:
        headers = ["Time", "mout", "Tout", "pipe.Q_flow"]

    out_path = base_dir / out_name
    actual_headers = _ensure_header(out_path, headers, sep)

    row_vals = [_fmt(vals.get(h, None)) for h in actual_headers]
    row = sep.join(row_vals)

    with out_path.open("a", encoding="utf-8") as g:
        g.write(row + "\n")

    return out_path, row


# ===== 使用例 =====
if __name__ == "__main__":
    base = Path("/mnt/d/work/openfoam/20251017_OpenCAE2025/heatedRoom")
    sync_file = base / "constant/boundarySyncFile"

    vals = parse_boundary_sync(sync_file)

    # ここで列が増えても OK（自動で不足列を追加し、既存行は空欄で補完）
    out_path, row = append_output_row(
        base,
        vals,
        headers=["Time", "mout", "Tout", "pipe.Q_flow", "TSensor", "TTarget"],
        sep="\t"
    )
    print(f"[INFO] Appended to {out_path}: {row}")