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

計算実行：Allrun.run
グラフ化
まとめ

計算実行：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>&amp;1 &amp;

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

#!/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."

# ==========================================================

# ① 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 の再コンパイルをスキップします。"

# ==========================================================

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"

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

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

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

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

tail -f log.chtMultiRegionFoam

1	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

===== 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.100372

304.2260700472616

2.88

214363.63539380228

293.807

298.0

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);

// 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[1].plot(df["Time"], df["mout"], label="mout")
axes[1].set_xlabel("time [s]")
axes[1].set_ylabel("mdot [kg/s]")
axes[1].set_title("Mass Flux")
axes[1].grid(True)
axes[1].legend()

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

# --- ④ Heat Flow ---
axes[3].plot(df["Time"], df["heatFlowSensor.Q_flow"], label="Q_flow")
axes[3].axhline(y=100000, color='red', linestyle='--', linewidth=1.5, label='Y=1400')  # ★追加！
axes[3].set_xlabel("time [s]")
axes[3].set_ylabel("Q [W]")
axes[3].set_title("Heat Flow")
axes[3].grid(True)
axes[3].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()

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[1].plot(df["Time"], df["mout"], label="mout")

axes[1].set_xlabel("time [s]")

axes[1].set_ylabel("mdot [kg/s]")

axes[1].set_title("Mass Flux")

axes[1].grid(True)

axes[1].legend()

# --- ③ Inlet Temperature ---

axes[2].plot(df["Time"], df["Tout"], label="Tinlet")

axes[2].set_xlabel("time [s]")

axes[2].set_ylabel("T [K]")

axes[2].set_title("Tout (OpenFOAM Inlet Temperature)")

axes[2].grid(True)

axes[2].legend()

# --- ④ Heat Flow ---

axes[3].plot(df["Time"], df["heatFlowSensor.Q_flow"], label="Q_flow")

axes[3].axhline(y=100000, color='red', linestyle='--', linewidth=1.5, label='Y=1400') # ★追加！

axes[3].set_xlabel("time [s]")

axes[3].set_ylabel("Q [W]")

axes[3].set_title("Heat Flow")

axes[3].grid(True)

axes[3].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