【解析設定の詳細】OpenFOAMの設定ファイルとFrontISTRへの温度の受け渡しプログラム
こんにちは（@t_kun_kamakiri）
この記事では、OpenFOAMで計算した時刻歴温度分布をFrontISTRへ渡し、各時刻の熱膨張を計算する仕組みを、フォルダ構成、メッシュ生成、座標変換、温度補間、FrontISTR入力の順に説明します。
chtMultiRegionFoamの設定ファイルを一つずつ読む（輻射・発熱面つき）」で説明しています。本記事では重複を避け、OpenFOAMの計算が完了した後からFrontISTRへ渡す部分を中心にします。
- OpenFOAMのメッシュファイルをFrontISTR形式へ直接変換しているわけではない
- 同じ寸法・同じ分割数のFrontISTRメッシュをPythonで新しく作る
- OpenFOAMのセル中心温度をFrontISTRの節点温度へ補間する
- FrontISTR自身はOpenFOAMファイルを読まず、変換後の
!TEMPERATUREを読む - OpenFOAMの1保存時刻につき、FrontISTRの線形静解析を1回行う
この連成で行っていること
今回の連成はOpenFOAMからFrontISTRへの一方向・準静的連成です。熱流体計算で得た温度だけを構造解析へ渡します。FrontISTRの変形をOpenFOAMへ戻して流路形状を更新する反復計算は行いません。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
OpenFOAM（101_0） 各時刻の solid/C, solid/T 各時刻の heaterMat/C, heaterMat/T │ │ Pythonで読み込み ▼ セル中心座標とセル温度を結合 │ │ 座標を平行移動 │ 8近傍の逆距離加重補間 ▼ FrontISTRの全節点温度 │ │ .msh / .cntを生成 ▼ FrontISTR（101_1） t=0, 5, 10, ... 600秒を個別に線形静解析 │ ▼ 変位・応力・温度をPVDで時刻歴表示 |
ここでいう「準静的」とは、温度はOpenFOAMの時間変化を使いますが、構造解析では慣性や前時刻の応力履歴を引き継がず、その瞬間の温度で釣り合う静的な変形を時刻ごとに求める、という意味です。
101_1のフォルダ構成
101_1_frontistr_cht_box_thermal_expansion直下は次の構成です。
|
1 2 3 4 5 6 7 8 9 10 |
101_1_frontistr_cht_box_thermal_expansion/ ├── README.md ├── ani/ ├── case/ ├── config/ ├── data/ ├── docs/ ├── python/ ├── requirements.txt └── tests/ |
| 場所 | 何をする場所か | 通常編集するか |
|---|---|---|
README.md | ケースの目的、実行方法、主要結果、GitHubへ含める範囲を説明する入口 | 最初に読む |
config/ | ヤング率、ポアソン比、密度、線膨張係数、基準温度をYAMLで管理する | 実材料に合わせて編集 |
python/ | FEMメッシュ生成、OpenFOAM温度読込み、温度写像、FrontISTR入力生成、実行、後処理を行う | 連成方法を変更するときに編集 |
case/ | 時刻ごとのFrontISTR入力と計算結果を保存する | 自動生成。直接編集しない |
data/ | 時刻歴CSV、グラフ、ParaView用PVD、要約YAML、公開用GIFを保存する | 自動生成された結果を確認 |
ani/ | ParaViewから出したPNG連番と元GIFを保存する | 可視化を作り直すときに使用 |
docs/ | 連成手順、FrontISTRキーワード、本記事などの説明資料 | 通常は参照 |
tests/ | OpenFOAMフィールド読込み、座標変換、温度補間、材料YAMLを検証する単体テスト | コード変更後に実行 |
requirements.txt | NumPy、Matplotlib、PyYAML、PillowなどPython依存関係 | 初回環境構築で使用 |
python/の中は役割ごとに分けています。
| ファイル | 役割 |
|---|---|
box_mesh.py | 直方体のFrontISTR六面体メッシュ、節点群、要素群を作る |
openfoam_temperature.py | OpenFOAMのCとTを読み、座標整合と節点補間を行う |
fistr_case.py | .msh、.cnt、hecmw_ctrl.datを作り、fistr1を実行して結果を読む |
run_thermal_expansion.py | 1時刻分の一連の処理をまとめる |
run_thermal_expansion_timehistory.py | OpenFOAMの全保存時刻を並べ、1時刻処理を繰り返してPVDとCSVを作る |
compress_preview_gif.py | 公開用GIFを縮小・減色する |
どれがmainプログラムか
通常の実行入口（main）はpython/run_thermal_expansion_timehistory.pyです。OpenFOAMの全保存時刻を処理し、FrontISTR解析、CSV、グラフ、PVDまで作ります。
|
1 2 3 |
# 通常はこちらを実行する python3 python/run_thermal_expansion_timehistory.py \ --of-case ../101_0_openfoam_cht_radiation_box |
python/run_thermal_expansion.pyもmain()を持ちますが、これは1時刻だけを確認するための入口です。変換方法を変更したときのデバッグや、まず300秒だけ試したい場合に使います。
|
1 2 3 4 |
# 300秒だけ試す python3 python/run_thermal_expansion.py \ --of-case ../101_0_openfoam_cht_radiation_box \ --time 300 |
box_mesh.py、openfoam_temperature.py、fistr_case.pyはmainから呼ばれる部品です。単独実行するものではありません。compress_preview_gif.pyだけは解析後に必要に応じて単独実行する補助ツールです。
各プログラムの関数と意味
| プログラム | 主な関数 | 何をしているか |
|---|---|---|
box_mesh.py | build_box_mesh() | 節点座標、六面体要素の接続、NALL・BOTTOM・TOP・EALLのID一覧を作る |
cell_centers() | 理想的な一様格子のセル中心を作る検証用関数。実際の温度写像ではOpenFOAMのCを直接読む | |
openfoam_temperature.py | _read_of_scalar_field() | OpenFOAMのTを読み、uniform/nonuniformの両形式を数値配列へ変える |
_read_of_vector_field() | OpenFOAMのセル中心CをN×3の座標配列へ変える | |
load_solid_cell_temperatures() | solidとheaterMatのC・Tを読み、1つの128セル点群へ結合する | |
align_cell_centers_to_node_mesh() | OpenFOAM点群とFrontISTR点群の外接直方体中心を一致させる | |
interpolate_to_nodes() | 8近傍IDWでセル中心温度からFrontISTR節点温度を作る | |
fistr_case.py | write_mesh() | FrontISTRの.mshへ節点、要素、グループ、材料、断面を書く |
write_hecmw_ctrl() | メッシュ・制御・結果ファイル名をhecmw_ctrl.datへ書く | |
write_cnt() | 解析種別、底面拘束、225節点温度、材料、ステップ、VTK出力を.cntへ書く | |
run_fistr() | 対象時刻フォルダでfistr1を実行し、ログを保存する | |
add_temperature_to_visualization() | FrontISTR標準VTUへTEMPERATUREとOpenFOAM実時刻を追加する | |
read_displacement() | 荷重適用後の.res.0.1から節点変位を読む | |
run_thermal_expansion.py | ensure_cell_centres() | 1時刻のCがなければOpenFOAMのpostProcessで生成する |
run_one_time() | 読込み、メッシュ生成、補間、入力生成、実行、結果集計を1本につなぐ中核関数 | |
main() | コマンドライン引数とYAMLを読み、指定した1時刻を処理する | |
run_thermal_expansion_timehistory.py | list_time_dirs() | OpenFOAMの数値名フォルダを時間順に並べる |
ensure_cell_centres_for_times() | 不足している全時刻のCをリージョン単位で一括生成する | |
write_pvd_collection() | 時刻別PVTUをOpenFOAM時刻付きPVDへまとめる | |
main() | 全時刻でrun_one_time()を繰り返し、summary.yaml、CSV、PNG、PVDを作る通常の入口 |
「メッシュ変換」ではなく「同等メッシュの再生成」
今回の実装では、OpenFOAMのpoints、faces、owner、neighbourを読み、FrontISTRの節点・要素へ一般的に変換しているわけではありません。OpenFOAM側の固体が単純な直方体・一様格子なので、同じ外形寸法と同じ分割数の構造格子をFrontISTR側で再生成しています。
| 項目 | OpenFOAM | FrontISTR |
|---|---|---|
| 固体外形 | 200 × 200 × 400 mm | 200 × 200 × 400 mm |
| 分割数 | 4 × 4 × 8セル | 4 × 4 × 8要素 |
| 温度の位置 | 128セルの中心 | 225節点 |
| 原点 | 流体領域中央のためx,y=0.4 m付近 | (0, 0, 0) |
| 要素 | 有限体積セル | 1次六面体要素 TYPE=361 |
OpenFOAMではsolidが120セル、heaterMatが8セルで、合わせて128セルです。FrontISTRでは両者を分けず、128個の六面体要素からなる1つの連続体として扱います。既定値では節点数は
N_{\mathrm{node}}
= (n_x+1)(n_y+1)(n_z+1)
= 5\times5\times9
=225
\end{align*}
要素数は
N_{\mathrm{elem}}
= n_x n_y n_z
=4\times4\times8
=128
\end{align*}
です。box_mesh.pyは次の式で節点座標を作ります。
x_i=\frac{iL_x}{n_x},\qquad
y_j=\frac{jL_y}{n_y},\qquad
z_k=\frac{kL_z}{n_z}
\end{align*}
|
1 2 3 4 5 6 7 |
for k in range(nz + 1): z = lz * k / nz for j in range(ny + 1): y = ly * j / ny for i in range(nx + 1): x = lx * i / nx node_ids[i, j, k] = node_id |
8節点を[n0,n1,n2,n3,n4,n5,n6,n7]の順に結び、FrontISTRの1次六面体要素TYPE=361を作ります。さらに、全節点NALL、底面節点BOTTOM、上面節点TOP、全要素EALLを自動作成します。底面は5×5=25節点で、後ほど全固定に使います。
OpenFOAMとFrontISTRで、外形寸法、軸方向、分割数が一致している必要があります。回転、拡大縮小、曲面、非構造格子にはそのまま使えません。一般形状では、メッシュ変換ツールまたは空間探索・面投影を含む専用マッピングが必要です。
OpenFOAMから読み出す値
OpenFOAMから必要なのは、各保存時刻・各固体セルの中心座標Cと温度Tです。たとえば300秒では次の4ファイルを使います。
|
1 2 3 4 5 6 7 8 |
101_0_openfoam_cht_radiation_box/ └── 300/ ├── solid/ │ ├── C # solidセル中心座標 │ └── T # solidセル温度 └── heaterMat/ ├── C # heaterMatセル中心座標 └── T # heaterMatセル温度 |
Tはソルバー結果に含まれますが、Cがない場合は次のOpenFOAM標準処理で生成します。
|
1 2 |
postProcess -func writeCellCentres -region solid -time 300 postProcess -func writeCellCentres -region heaterMat -time 300 |
全時刻処理ではensure_cell_centres_for_times()が不足している時刻だけをまとめて生成します。CとTは同じリージョン・同じメッシュのinternalFieldなので、同じセル順序です。プログラムはuniformとnonuniform List<scalar>の両形式に対応し、宣言セル数と実データ数が一致するかも確認します。
|
1 2 3 4 5 |
centers = _read_of_vector_field(c_path) # shape = (N, 3) temps = _read_of_scalar_field(t_path) # shape = (N,) # solid 120セルとheaterMat 8セルを1つの点群にする return np.concatenate(all_centers), np.concatenate(all_temps) |
温度単位はOpenFOAMもFrontISTRもKを使うため、ここでは℃への変換は行いません。
座標系を合わせる
OpenFOAMの固体は1 m角の流体領域中央にあり、外形はx,y=0.4～0.6 m、z=0～0.4 mです。FrontISTRメッシュは原点始まりで、x,y=0～0.2 m、z=0～0.4 mです。このまま距離を計算すると別の場所にある点群と判断されるため、補間前に原点を合わせます。
OpenFOAMセル中心群の外接直方体中心を$\mathbf{c}_{\mathrm{OF}}$、FrontISTR節点群の外接直方体中心を$\mathbf{c}_{\mathrm{FEM}}$とすると、平行移動量は
\mathbf{t}
=\mathbf{c}_{\mathrm{FEM}}-\mathbf{c}_{\mathrm{OF}},\qquad
\mathbf{x}’_{\mathrm{cell}}
=\mathbf{x}_{\mathrm{cell}}+\mathbf{t}
\end{align*}
です。本ケースでは
\mathbf{t}=(-0.4,-0.4,0)\ \mathrm{m}
\end{align*}
|
1 2 3 4 |
source_center = 0.5 * (cell_centers.min(axis=0) + cell_centers.max(axis=0)) target_center = 0.5 * (node_coords.min(axis=0) + node_coords.max(axis=0)) translation = target_center - source_center aligned_centers = cell_centers + translation |
この処理は平行移動だけです。回転やスケーリングはしません。したがって、移動後の外形寸法と軸が一致していることを確認する必要があります。適用した移動量は時刻ごとのsummary.yamlにも保存します。
セル中心温度を節点温度へ変換する
OpenFOAMは有限体積法なので温度はセル中心にあります。一方、今回FrontISTRへ与える!TEMPERATUREは節点値です。4×4×8個のセル値128点を、その格子境界にある225節点へ移す必要があります。
各FrontISTR節点$\mathbf{x}_n$について、移動後のOpenFOAMセル中心から近い順に$k=8$個を選びます。距離$d_i$、重み$w_i$、節点温度$T_n$を次式で求めます。
d_i&=\left\|\mathbf{x}_n-\mathbf{x}_{c,i}\right\|,\\
w_i&=\frac{1}{d_i},\\
T_n&=\frac{\displaystyle\sum_{i=1}^{k}w_iT_i}{\displaystyle\sum_{i=1}^{k}w_i},\qquad k=8
\end{align*}
近いセルほど強く、遠いセルほど弱く反映する平均です。同じ一様格子の内部節点では、周囲8セルまでの距離がほぼ等しいため、8セルの平均に近くなります。節点とセル中心がほぼ一致して$d_i<10^{-9}$ mとなった場合は、ゼロ除算を避けてそのセル温度を直接使います。
|
1 2 3 4 5 6 7 8 9 10 11 |
for idx in range(n_nodes): d = np.linalg.norm(cell_centers - node_coords[idx], axis=1) nearest = np.argsort(d)[:k_eff] dn = d[nearest] if dn[0] < 1.0e-9: result[idx] = cell_temperatures[nearest[0]] continue w = 1.0 / dn result[idx] = np.sum(w * cell_temperatures[nearest]) / np.sum(w) |
IDWは実装が単純で、今回の整った格子には使いやすい一方、エネルギー保存を保証する補間ではありません。局所的な最高温度は周囲との平均で低くなります。300秒ではOpenFOAMセル最高温度293.7681 Kに対し、補間後の節点最高温度は293.6414 Kでした。高温部を厳密に保持したい場合は、近傍数、形状関数補間、保存型マッピングを検討します。
1時刻分の処理をコードで追う
run_thermal_expansion.pyのrun_one_time()が、1時刻分の変換と解析をまとめています。処理順は次の通りです。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
# 1. CがなければOpenFOAMのpostProcessで作る ensure_cell_centres(of_case, time_dir) # 2. solid + heaterMatのセル中心座標と温度を読む centers, temps = load_solid_cell_temperatures(of_case, time_dir) # 3. FrontISTRの節点・要素メッシュを作る mesh = fistr_case.write_mesh(case_dir, nx, ny, nz, ...) node_coords = np.array([xyz for _nid, xyz in mesh["nodes"]]) # 4. 座標を合わせ、節点温度へ補間する aligned_centers, translation = align_cell_centers_to_node_mesh(centers, node_coords) node_temps = interpolate_to_nodes(node_coords, aligned_centers, temps, k=8) # 5. FrontISTR制御ファイルへ節点温度を書く fistr_case.write_cnt(case_dir, node_ids, node_temps, ...) # 6. FrontISTRを実行し、VTKへ温度と実時刻を追加する fistr_case.run_fistr(case_dir) fistr_case.add_temperature_to_visualization(case_dir, node_coords, node_temps, time_value) |
300秒を処理すると、次のようなフォルダができます。
|
1 2 3 4 5 6 7 8 9 10 |
case/t_300/ ├── box_thermal_expansion.msh # 節点・要素・節点群・材料・断面 ├── box_thermal_expansion.cnt # 解析種別・拘束・節点温度・材料・出力 ├── hecmw_ctrl.dat # 入出力ファイルの対応 ├── box_thermal_expansion.res.0.0 # 荷重適用前の結果 ├── box_thermal_expansion.res.0.1 # 荷重適用後の結果 ├── box_thermal_expansion_vis_psf.* # ParaView用PVTU/VTU ├── log.fistr1 # 標準出力・エラー ├── FSTR.msg / FSTR.sta # 詳細メッセージ・進行状況 └── summary.yaml # 温度範囲、移動量、変位要約 |
FrontISTRの材料設定YAML
人が変更する材料値はconfig/material_properties_steel.yamlへ集約しています。
|
1 2 3 4 5 |
young_modulus_Pa: 2.05e+11 poisson_ratio: 0.3 density_kg_m3: 7850.0 thermal_expansion_coeff_per_K: 1.2e-5 reference_temperature_K: 293.15 |
| 項目 | 意味 | 今回の値 |
|---|---|---|
young_modulus_Pa | 弾性変形の硬さ。熱変形を拘束したときの応力へ影響 | 205 GPa |
poisson_ratio | 軸方向ひずみに対する横方向ひずみ | 0.3 |
density_kg_m3 | 密度。今回の慣性なし静解析では結果への直接影響はない | 7850 kg/m³ |
thermal_expansion_coeff_per_K | 1 K上昇したときの自由熱ひずみ | 12×10-6/K |
reference_temperature_K | 熱ひずみを0とする基準温度 | 293.15 K |
YAMLがFrontISTR入力になるまで
FrontISTRはYAMLを直接読みません。mainプログラムがPyYAMLのsafe_load()で辞書へ変換し、その値をwrite_mesh()とwrite_cnt()へ渡します。
|
1 2 3 4 5 6 |
# run_thermal_expansion_timehistory.py mat = yaml.safe_load(Path(args.material).read_text(encoding="utf-8")) summary = run_one_time( of_case, time_dir, case_dir, mat, args.nx, args.ny, args.nz ) |
run_one_time()では、同じ辞書から必要なキーを明示して渡します。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
# .mshへ渡す値 fistr_case.write_mesh( case_dir, nx, ny, nz, lx, ly, lz, young_modulus=mat["young_modulus_Pa"], poisson_ratio=mat["poisson_ratio"], density=mat["density_kg_m3"], thermal_expansion_coeff=mat["thermal_expansion_coeff_per_K"], ) # .cntへ渡す値 fistr_case.write_cnt( case_dir, node_ids, node_temps, reference_temperature=mat["reference_temperature_K"], young_modulus=mat["young_modulus_Pa"], poisson_ratio=mat["poisson_ratio"], thermal_expansion_coeff=mat["thermal_expansion_coeff_per_K"], ) |
| YAMLキー | .mshで使う場所 | .cntで使う場所 |
|---|---|---|
young_modulus_Pa | !MATERIAL / !ITEM=1 | !ELASTICの1番目 |
poisson_ratio | !MATERIAL / !ITEM=1 | !ELASTICの2番目 |
density_kg_m3 | !MATERIAL / !ITEM=2 | 今回の.cntには書かない |
thermal_expansion_coeff_per_K | !MATERIAL / !ITEM=3 | !EXPANSION_COEFF |
reference_temperature_K | 使わない | !REFTEMPと!INITIAL_CONDITION |
sources、note | 人が根拠と注意点を読むための情報。計算には使わない | |
つまり、YAMLは人が編集しやすい材料設定の原本で、PythonがFrontISTRの2種類の記法へ展開します。材料値を変更した後は既存のcase/t_*を直接直すのではなく、YAMLを変更してmainプログラムを再実行します。
OpenFOAMではsolidとheaterMatを別リージョンとして扱いますが、今回のFrontISTRモデルは両者を1つの鋼材相当ブロックとして扱います。実際のヒートマットと母材で弾性率や線膨張係数が異なる場合は、要素群と材料を分け、接着層や接触条件も追加する必要があります。
.mshを一つずつ読む
!NODE：節点番号と座標
|
1 2 3 4 5 |
!NODE 1, 0, 0, 0 2, 0.05, 0, 0 3, 0.10, 0, 0 ... |
形式は節点ID, x, y, zです。座標単位はmです。節点IDは1から始まり、x方向、y方向、z方向の順に増えます。
!ELEMENT, TYPE=361：六面体要素
|
1 2 |
!ELEMENT, TYPE=361 1, n0, n1, n2, n3, n4, n5, n6, n7 |
TYPE=361は8節点1次六面体ソリッド要素です。節点順序が不正だと負体積や反転要素になるため、box_mesh.pyでは下面4点、対応する上面4点の順に一定の接続規則で作ります。
!NGROUPと!EGROUP：名前付き集合
|
1 2 3 4 5 6 7 |
!NGROUP, NGRP=NALL # 全225節点 ... !NGROUP, NGRP=BOTTOM # z=0の25節点 ... !NGROUP, NGRP=TOP # z=0.4の25節点 ... !EGROUP, EGRP=EALL # 全128要素 |
BOTTOMは拘束、TOPは上面変位の集計、NALLは初期温度、EALLは材料割当てに使います。
BOTTOM
TOP
このようにNOTEセットグループを作っておくことで境界条件の設定に紐づけることができます。
今回は後ほどBOTTOMに固定条件を与えます。
!MATERIALと!SECTION
|
1 2 3 4 5 6 7 8 9 10 |
!MATERIAL, NAME=STEEL, ITEM=3 !ITEM=1, SUBITEM=2 2.05e+11, 0.3 !ITEM=2, SUBITEM=1 7850 !ITEM=3, SUBITEM=1 1.2e-5 !SECTION, TYPE=SOLID, EGRP=EALL, MATERIAL=STEEL 1.0 |
ITEM=1はヤング率とポアソン比、ITEM=2は密度、ITEM=3は線膨張係数です。!SECTIONは全要素EALLを3次元ソリッドとし、材料STEELを割り当てます。
FrontISTRは!SECTIONが参照する材料名をメッシュ読込み時に解決するため、.mshにも材料定義が必要でした。.mshでは!ELASTICではなくHEC-MWの!ITEM形式を使います。同じ値を.cntにも書くので、両方が食い違わないようYAMLから生成します。
hecmw_ctrl.datを読む
|
1 2 3 4 5 6 7 8 9 10 11 |
!MESH, NAME=fstrMSH, TYPE=HECMW-ENTIRE box_thermal_expansion.msh !CONTROL, NAME=fstrCNT box_thermal_expansion.cnt !RESULT, NAME=fstrRES, IO=OUT box_thermal_expansion.res !RESULT, NAME=vis_out, IO=OUT box_thermal_expansion_vis |
FrontISTR本体fistr1が、どのファイルをメッシュ、解析制御、数値結果、可視化結果として使うかを対応付けるファイルです。TYPE=HECMW-ENTIREは単一のHEC-MWメッシュファイルを読む指定です。
.cntを一つずつ読む
!SOLUTION, TYPE=STATIC
|
1 2 3 |
!VERSION 3 !SOLUTION, TYPE=STATIC |
入力形式バージョン3、線形静解析を指定します。OpenFOAMの300秒をFrontISTR内部で300秒間積分する指定ではありません。300秒の温度場を1つの静的熱荷重として解きます。
!WRITE：出力指定
|
1 2 |
!WRITE,RESULT,FREQUENCY=1 !WRITE,VISUAL |
節点変位・応力などの結果ファイルと、ParaView向け可視化ファイルを出力します。
!SOLVER：連立方程式の解法
|
1 2 3 |
!SOLVER,METHOD=CG,PRECOND=1,NSET=0,ITERLOG=NO,TIMELOG=NO 5000, 1 1.0e-08, 1.00, 0.0 |
METHOD=CG：共役勾配法PRECOND=1：前処理を使用5000：最大反復回数1.0e-08：収束判定値ITERLOG=NO、TIMELOG=NO：詳細ログを抑制
このソルバーが最終的に解く基本形は、熱ひずみから作られる等価節点荷重$\mathbf{f}_{\mathrm{th}}$を使った
\mathbf{K}\mathbf{u}=\mathbf{f}_{\mathrm{th}}
\end{align*}
です。$\mathbf{K}$は剛性行列、$\mathbf{u}$は節点変位です。
!REFTEMPと初期温度
|
1 2 3 4 |
!REFTEMP 293.15 !INITIAL_CONDITION, TYPE=TEMPERATURE NALL, 293.15 |
!REFTEMPは熱ひずみが0になる温度です。OpenFOAMの初期温度と同じ293.15 Kにします。等方材料の自由熱ひずみは
\boldsymbol{\varepsilon}_{\mathrm{th}}
=\alpha\left(T-T_{\mathrm{ref}}\right)\mathbf{I}
\end{align*}
です。$\alpha$は線膨張係数、$\mathbf{I}$は単位テンソルです。応力は全ひずみから熱ひずみを差し引いて
\boldsymbol{\sigma}
=\mathbf{D}\left(\boldsymbol{\varepsilon}(\mathbf{u})-\boldsymbol{\varepsilon}_{\mathrm{th}}\right)
\end{align*}
と評価されます。!INITIAL_CONDITIONは全節点の初期値を設定し、次の!TEMPERATUREが各節点の実際の熱荷重を与えます。
!BOUNDARY：底面固定
|
1 2 3 4 |
!BOUNDARY, GRPID=1 BOTTOM,1,1 BOTTOM,2,2 BOTTOM,3,3 |
形式は節点群, 最初の自由度, 最後の自由度です。自由度1、2、3はX、Y、Z変位なので、z=0の底面25節点をXYZ全方向に固定します。
この拘束は剛体移動を防ぎますが、実物の支持が滑り、ボルト締結、接触などの場合は結果が変わります。特に熱応力は拘束条件に敏感なので、実機評価では支持方法をモデル化し直す必要があります。
!TEMPERATURE：OpenFOAMから渡された値
|
1 2 3 4 5 6 |
!TEMPERATURE, GRPID=1 1, 293.1607598 2, 293.1631103 3, 293.1752641 ... 225, 293.2206677 |
左がFrontISTR節点ID、右がIDW補間後の温度[K]です。FrontISTRが認識するのはこの値であり、OpenFOAMのTファイルではありません。PythonがOpenFOAM形式からFrontISTR形式への橋渡しをしています。
!MATERIAL：弾性と線膨張
|
1 2 3 4 5 |
!MATERIAL, NAME=STEEL !ELASTIC 2.05e+11, 0.3 !EXPANSION_COEFF 1.2e-5 |
!ELASTICはヤング率[Pa]とポアソン比、!EXPANSION_COEFFは線膨張係数
です。温度依存性は入れておらず、全温度範囲で一定です。
!STEP：拘束と温度荷重を有効化
|
1 2 3 |
!STEP, SUBSTEPS=1, CONVERG=1.000E-07 BOUNDARY,1 LOAD,1 |
SUBSTEPS=1で温度荷重を1段階で与えます。BOUNDARY,1はGRPID=1の底面拘束、LOAD,1はGRPID=1の節点温度を有効にします。
!VISUAL：ParaView出力
|
1 2 3 4 |
!VISUAL, method=PSR !surface_num=1 !surface 1 !output_type = VTK |
表面結果をVTK形式で出力します。FrontISTR標準出力にはDISPLACEMENT、NodalSTRESS、NodalMISESが含まれますが、熱荷重として使ったTEMPERATUREは標準VTKに含まれません。
そこでadd_temperature_to_visualization()が解析後のVTU節点座標を入力メッシュへ照合し、TEMPERATURE配列を追加します。同時にTimeValueをOpenFOAMの実時刻へ置き換えます。これによりParaView上で温度、変位、応力を同じ時刻で切り替えられます。
全時刻をFrontISTRへ渡す仕組み
run_thermal_expansion_timehistory.pyはOpenFOAMケース直下の数値名フォルダを探し、数値として並べ替えます。
|
1 |
0, 5, 10, 15, ... 595, 600 |
そして各時刻についてcase/t_<time>を作り、run_one_time()を呼びます。既定条件では121回のFrontISTR静解析です。
|
1 2 3 |
for time_dir in times: case_dir = case_root / f"t_{time_dir}" summary = run_one_time(of_case, time_dir, case_dir, material, nx, ny, nz) |
最後に、各時刻のPVTUをdata/thermal_expansion_timehistory.pvdへまとめます。
|
1 2 3 4 5 |
<DataSet timestep="0" file="../case/t_0/...pvtu"/> <DataSet timestep="5" file="../case/t_5/...pvtu"/> <DataSet timestep="10" file="../case/t_10/...pvtu"/> ... <DataSet timestep="600" file="../case/t_600/...pvtu"/> |
このPVDをParaViewで開くと、OpenFOAMと同じ時刻値でFrontISTR結果を連番表示できます。ただし、PVDが時系列に見せているだけで、FrontISTR解析同士が内部状態を受け渡しているわけではありません。
実行コマンド
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
# 初回だけPython依存関係を導入 cd sample/101_1_frontistr_cht_box_thermal_expansion python3 -m pip install -r requirements.txt # OpenFOAM環境を読み込む（postProcessを使うため） source /usr/lib/openfoam/openfoam2512/etc/bashrc # FrontISTR実行ファイルを指定 export FISTR1=/home/kamakiri/local/frontistr/bin/fistr1 # 全保存時刻を計算 cd python python3 run_thermal_expansion_timehistory.py \ --of-case ../../101_0_openfoam_cht_radiation_box |
単一時刻だけ確認する場合は次のコマンドです。
|
1 2 3 |
python3 run_thermal_expansion.py \ --of-case ../../101_0_openfoam_cht_radiation_box \ --time 300 |
時刻を間引く場合は--every 2、範囲を限定する場合は--start-timeと--end-timeを使います。
FrontISTRで使っている単位
FrontISTRの入力ファイルには「このモデルはm単位」といった単位情報がありません。すべての数値を一貫した単位系で入力する必要があります。本ケースはSI単位系です。
| 量 | 本ケースの単位 | 例 |
|---|---|---|
| 座標 | m | ブロックは0.2 × 0.2 × 0.4 m |
| 変位 | m | .resとVTUのDISPLACEMENTはm |
| 力 | N | 今回の構造解析では外力なし |
| 応力・ヤング率 | Pa = N/m² | E=2.05×1011 Pa |
| 密度 | kg/m³ | 7850 kg/m³ |
| 温度 | K | 基準温度293.15 K |
| 線膨張係数 | 1/K | 1.2×10-5/K |
| OpenFOAM時刻 | s | PVDの0～600 s。FrontISTR静解析の時間積分値ではない |
data/timehistory.csvでは読みやすくするため、変位だけmからmmへ変換しています。一方、PVTU/VTUの変位はmのままです。ParaViewで変形倍率を設定するときはこの点に注意します。
EasyISTRでよく使われるmm・N・MPa系へ数値を変えて作り直す場合は、座標を200、200、400 mm、ヤング率を205000 MPa、変位をmmとして全体を統一します。線膨張係数と温度差の数値は変わりません。密度は採用する質量・時間単位まで含めて整合させる必要があるため、動解析では特に注意が必要です。今回の静的熱膨張では密度は結果へ直接影響しません。
caseフォルダをEasyISTRで読めるか
結論からいうと、FrontISTR入力としては正しいのでfistr1で直接実行できますが、既存の.mshと.cntをEasyISTRで完全に読み戻して編集できるかはEasyISTRの版と操作経路に依存します。FrontISTR公式資料ではEasyISTRは対応プリポストとして案内され、EasyISTRが.mshと.cntを生成する流れが示されています。しかし、既存の任意の.cntを全キーワード込みでプロジェクトへ逆読込みすることまでは保証されていません。
本ケースにはEasyISTRで扱う際に特に注意する点があります。
case直下ではなく、case/t_300など1時刻のフォルダを対象にするbox_thermal_expansion.mshはHEC-MW単一領域メッシュで、FrontISTR自身は読めるbox_thermal_expansion.cntには225節点分の!TEMPERATUREがあり、EasyISTRで設定を保存し直すと失われる可能性がある.mshと.cntの材料定義が重複しているため、GUI側の再生成で記法が変わる可能性がある- 単位メタデータがないため、EasyISTR上でも0.2をmとして扱う必要がある
試す場合は元のcase/t_300を変更せず、別フォルダへ複製します。
|
1 |
cp -a case/t_300 case/t_300_easyistr_test |
EasyISTRへメッシュを読み込めた場合も、保存前後で少なくとも次を比較してください。
|
1 2 3 |
diff -u \ case/t_300/box_thermal_expansion.cnt \ case/t_300_easyistr_test/box_thermal_expansion.cnt |
!TEMPERATURE、!REFTEMP、!EXPANSION_COEFF、BOTTOM拘束が残っていることを確認します。単に形状と結果を見たい場合は、EasyISTRへ戻すよりdata/thermal_expansion_timehistory.pvdまたは時刻別PVTUをParaViewで開く方が確実です。
参考：FrontISTR公式FAQではEasyISTRを対応プリポストとして紹介しています。また、FrontISTR公式の逐次実行手順では、本ケースと同じhecmw_ctrl.dat、単一領域.msh、解析制御.cntの構成が説明されています。
変換が正しいか確認する
| 確認項目 | 既定値・確認方法 |
|---|---|
| OpenFOAM固体セル数 | solid 120 + heaterMat 8 = 128 |
| FrontISTR節点・要素数 | 225節点、128要素 |
| 座標移動量 | (-0.4, -0.4, 0) m |
| 温度単位 | OpenFOAM、YAML、!TEMPERATUREの全てK |
| 基準温度 | OpenFOAM初期温度と!REFTEMPが293.15 K |
| 底面拘束 | BOTTOMの25節点、自由度1～3 |
| 300秒の温度範囲 | セル293.1531～293.7681 K、節点293.1575～293.6414 K |
| ParaView | PVDに121個のPVTU、各VTUにTEMPERATUREとDISPLACEMENT |
単体テストは次のコマンドで実行します。
|
1 2 |
cd sample/101_1_frontistr_cht_box_thermal_expansion python3 -m unittest discover -s tests -v |
解析エラーではcase/t_<time>/log.fistr1、FSTR.msg、FSTR.staを確認します。値がおかしい場合はsummary.yamlのセル温度範囲、節点温度範囲、座標移動量、基準温度を確認します。
まとめ
- OpenFOAMメッシュを直接変換せず、同じ寸法・分割数のFrontISTRメッシュを再生成する
- OpenFOAMの
CとTを読み、solidとheaterMatを128セルの点群として結合する - 座標を(-0.4, -0.4, 0) m平行移動し、8近傍IDWで225節点の温度を作る
- 節点温度を
.cntの!TEMPERATUREへ書き、底面全固定の線形静解析を行う - OpenFOAMの各保存時刻についてFrontISTRを独立に実行し、PVDで1つの時刻歴として表示する
- 材料値はYAML、処理はPython、生成結果は
caseとdataに分けて管理する
OpenFOAM側の設定内容は別記事、計算結果と熱変形の読み方は概要記事と分けることで、本記事は「どの値を、どのコードで、どのFrontISTR入力へ渡しているか」を調べるための技術資料として使える構成にしています。
Follow @t_kun_kamakiri
