計算実行:Allrun.run

冒頭でも説明しましたが再度Allrun.runスクリプトを記載しておきます。

#!/bin/bash
# Allrun.run

set -e
cd "$(dirname "$0")"

echo "[INFO] Remove old output.dat if exists..."
if [ -f output.dat ]; then
    rm -f output.dat
    echo "[INFO] output.dat removed."
else
    echo "[INFO] No existing output.dat."
fi

# ==========================================================
# ① preload.mos の再コンパイル確認
# ==========================================================
echo -n "preload.mos を再コンパイルしますか? [y/N]: "
read answer

if [[ "$answer" == "y" || "$answer" == "Y" ]]; then
    echo "[INFO] Running preload.mos (initial Modelica compile)..."
    ./runMos.sh preload.mos
    echo "[INFO] Preload complete."
else
    echo "[INFO] preload.mos の再コンパイルをスキップします。"
fi
# ==========================================================

echo "boundarySyncFile をコピーします..."
cp -v constant/fluid/boundarySyncFile_base constant/fluid/boundarySyncFile

echo "chtMultiRegionFoam をバックグラウンドで実行します..."
chtMultiRegionFoam > log.chtMultiRegionFoam 2>&1 &

SOLVER_PID=$!
echo "chtMultiRegionFoam started (PID: ${SOLVER_PID})"
echo "ログ: log.chtMultiRegionFoam"

やっているのは以下です。

  1. output.dat の削除
    古い結果ファイル output.dat が存在する場合は削除する。
    存在しなければ「無い」とメッセージを表示する。
  2. preload.mos の再コンパイルの確認
    「preload.mos を再コンパイルしますか?」と質問する(y/Y → 実行)
    ./runMos.sh preload.mos によりModelica の初期コンパイルを行う
  3. boundarySyncFile のコピー
    constant/fluid/boundarySyncFile_baseをconstant/fluid/boundarySyncFile へ上書きコピーする。
  4. chtMultiRegionFoam をバックグラウンド起動
    chtMultiRegionFoam > log.chtMultiRegionFoam 2>&1 &
    標準出力とエラーを log.chtMultiRegionFoam に保存
  5. 実行しているプロセス ID を取得
    SOLVER_PID=$!
    PID とログファイル名を表示する

計算実行するとlog.chtMultiRegionFoamにOpenFOAMの実行とともにOpenModelicaも実行されているのを確認できます。

ターミナルでlogを見たい場合は、

tail -f log.chtMultiRegionFoam

とするとリアルタイムでログが確認できます。
「ctrl + c」で閉じることができます。

===== codeTempFieldValue (coded) start =====
 === probe Temperature === 
41177 Tprobe = 293.807431638
 === probe Temperature end === 
===== Execute Python updateMOS_and_Run.py =====
Lock file created: "constant/updateMOS.lock"
{"", ""}
LOG_SUCCESS       | info    | The initialization finished successfully without homotopy method.
LOG_SUCCESS       | info    | The simulation finished successfully.
0
0.100372
304.2260700472616
2.88
214363.63539380228
293.807
298.0
true
true
true
true
true
true
true

毎回OpenModelicaの計算を挟んでいるので、少し計算の遅さを感じます。

OpenModelicaのコンパイル用のmosファイルはこちらです。

preload.mos

// preload.mos
// === 元の動く形を維持しつつ、パスだけ変数にする ===

// Modelica ライブラリ読み込み
loadModel(Modelica);

// ★ run_OM.sh から渡された現在の作業フォルダ
caseDir := getEnvironmentVar("OM_WORKDIR");

// ★ モデルファイルパス
modelPath := caseDir + "/HVACSystem001.mo";

// モデルファイルを読み込む(元の形)
loadFile(modelPath);

// === 実行ファイル名 ===
filePrefix := "HVACSystem";

// === build(元の形を維持)===
buildModel(HVACSystem001,
  startTime         = 0,
  stopTime          = 300,
  numberOfIntervals = 1000,
  tolerance         = 1e-6,
  method            = "dassl",
  outputFormat      = "mat",
  fileNamePrefix    = filePrefix
);

// ★ 実行ファイルパス(/tmp/OpenModelica は固定)
exePath := "/tmp/OpenModelica/" + filePrefix;

// ★ 結果ファイルパス(現在のフォルダに出す)
resultPath := caseDir + "/HVACSystemResult.mat";

// === ★ 元の形を維持した system() ===
system(exePath + " -override Tinlet=293.153,TSensor=297.223 -r=" + resultPath);

グラフ化

output.datがOpenModelicaの結果なのでそれをグラフ化して、思っている挙動を示しているかを確認します。

graph.ipynb

import pandas as pd
import matplotlib.pyplot as plt
from pathlib import Path
import numpy as np

# ====== ファイルパス設定 ======
file_path = Path("output.dat")

# ====== データ読み込み ======
df = pd.read_csv(file_path, sep=r"[\t,]+", engine="python")

# ====== 図と軸の設定 ======
fig, axes = plt.subplots(1, 4, figsize=(20, 4))  # 横4つのグラフ

# X軸の範囲と刻み幅(30秒刻みに変更)
x_min, x_max = df["Time"].min(), df["Time"].max()
# x_ticks = np.arange(x_min, x_max + 30, 30)  # ★←ここで間引き
# x_ticks = np.arange(0,1, 0.1)  # ★←ここで間引き
# y_ticks = np.arange(273,310, 0.1)  # ★←ここで間引き

# --- ① Target and Measured T ---
axes[0].plot(df["Time"], df["TSensor"], label="TSensor")
axes[0].plot(df["Time"], df["ramp.y"], label="Target")
axes[0].set_xlabel("time [s]")
axes[0].set_ylabel("T [K]")
axes[0].set_title("Target and Measured T")
axes[0].grid(True)
axes[0].set_ylim(280, 330)  # 縦軸範囲設定
axes[0].legend()

# --- ② Mass Flux ---
axes.plot(df["Time"], df["mout"], label="mout")
axes.set_xlabel("time [s]")
axes.set_ylabel("mdot [kg/s]")
axes.set_title("Mass Flux")
axes.grid(True)
axes.legend()

# --- ③ Inlet Temperature ---
axes.plot(df["Time"], df["Tout"], label="Tinlet")
axes.set_xlabel("time [s]")
axes.set_ylabel("T [K]")
axes.set_title("Tout (OpenFOAM Inlet Temperature)")
axes.grid(True)
axes.legend()

# --- ④ Heat Flow ---
axes.plot(df["Time"], df["heatFlowSensor.Q_flow"], label="Q_flow")
axes.axhline(y=100000, color='red', linestyle='--', linewidth=1.5, label='Y=1400')  # ★追加!
axes.set_xlabel("time [s]")
axes.set_ylabel("Q [W]")
axes.set_title("Heat Flow")
axes.grid(True)
axes.legend()

# ====== すべてのグラフに横軸tick設定&回転 ======
for ax in axes:
    ax.locator_params(axis='x', nbins=5)
    ax.tick_params(axis='x', rotation=45)

plt.tight_layout()
plt.show()

以下のようになっていたらOK。

ただし、現状温度上限が激しくて途中で発散してしまいましたが、動作はしていることを確認できたので、今後修正していきます。

まとめ

OpenFOAMとOpenModelicaを、FMU4FOAMを使わずに連携させる仕組みを作りました。
OpenFOAMで取得した温度や流量をPython経由でOpenModelicaへ渡し、PID制御結果を再びOpenFOAMに反映する流れです。
毎ステップでMOSファイルを自動生成して計算を回すことで、モデル間のやりとりを完全に自動化できました。

今回は、FMU4FOAMが対応していないマルチリージョン環境でも動作するプログラムを作成できたことになります。

もし、OpenFOAMのマルチリージョンソルバでOpenModelicaと連成させたい場合に参考にしてください。

1 2 3 4 5