【OpenFOAM×FrontISTR】ヒートマットで温めた金属ブロックの熱膨張。片方向連成で時刻歴計算

工作機械や電子機器の設計では、次のような場面によく出会います。

どこか一部が発熱し、そこから周囲へ熱が伝わりながら、全体がゆっくり変形していく
熱解析（温度分布）はできても、そこから構造解析（熱変形）へつなぐのが手間
ヒーターのON/OFFのような時間変化を含めた「変形の時刻歴」まで見たい

この記事では、OpenFOAM（共役熱伝達＋輻射）で温度分布を計算し、その結果をFrontISTR（構造解析）へ時刻ごとに引き継いで熱膨張の時刻歴を求めるという、2つのオープンソースCAEをつなぐ流れを、実際に動かしながら説明します。

この記事でわかること

OpenFOAMのchtMultiRegionFoamで、固体・流体・発熱体を3つの独立リージョンに分ける方法
輻射（fvDOM）を含めたまま、発熱源を「セルへの体積発熱」ではなく「面への熱量境界条件」として与える方法
ヒーターのON/OFFのような時間変化する発熱を、OpenFOAM側で表現する方法
OpenFOAMの温度計算結果を、FrontISTRの構造解析へ時刻歴で引き継ぐ方法
実際につまずいたエラーとその対処（領域分割の衝突、境界条件の変数展開、FrontISTR特有のファイル分割のクセ）

解析設定の詳細は別途記事にしようと思います。

対象とする課題
使った手法とその位置づけ
全体構成とジオメトリ
発熱源を「セルではなく面」で与える
輻射（fvDOM）を含める理由
FrontISTRへ温度を引き継ぐ
得られた結果
応用できること
まとめ

対象とする課題

工作機械のベッドや電子機器の筐体など、多くの構造物は「どこか一部が発熱し、そこから熱が伝わりながら全体がゆっくり変形していく」という状況にさらされています。これを素直に再現しようとすると、次の3つを同時に扱う必要があります。

発熱体まわりの熱の伝わり方（伝導・対流・輻射）
発熱パターンの時間変化（ヒーターのON/OFFなど）
温度分布から生じる構造の熱変形

今回はこれを、次のようなシンプルなモデルで確認しました。1000mm角の空気の中に、200×200×400mmの金属ブロックが床の上に立っていて、その側面の一部に「ヒートマット」を貼り付けたモデルです。

計算開始（0秒）の様子です。固体・流体とも一様に20℃の状態からスタートし、固体側面の一部に貼ったヒートマットから、これから加熱が始まります。

今回の問いはシンプルです。

今回の問い

ヒーターを0〜5分ON（20W）、5〜10分OFFにしたとき、金属ブロックの温度分布と変形はどのように時間変化するか。

使った手法とその位置づけ

使った手法は、OpenFOAM（chtMultiRegionFoam、共役熱伝達＋輻射）と、FrontISTR（線形静解析）の片方向連成です。

OpenFOAM側：流体（空気）・固体・ヒートマットの3リージョンで、熱伝導・自然対流・輻射（fvDOM）・時間変化する発熱を同時に解く
FrontISTR側：OpenFOAMが保存した各時刻の固体温度分布を読み込み、その時刻ごとに「今の温度分布なら、変形はいくらか」という線形静解析を行う

他の方法と比べると、次のような位置づけになります。

方法	特徴
温度分布を仮定した簡易式（$\Delta L=\alpha L \Delta T$）	手軽だが、局所加熱による温度分布の偏りや反りは表現できない
OpenFOAMの熱解析のみ	温度分布は分かるが、変形量は分からない
FrontISTRの構造解析のみ（温度分布を手入力）	変形は分かるが、その温度分布が物理的に妥当かは別途確認が要る
完全連成（熱↔構造を毎時刻双方向でやり取り）	もっとも正確だが、実装・計算コストが大きい
今回：片方向・準静的連成	熱計算はOpenFOAMで物理的に求め、各保存時刻の温度分布を「その瞬間の静的な状態」とみなしてFrontISTRで変形を計算する

今回のように、変形量が非常に小さく（後述しますが数μmオーダー）、変形が流れ場や温度分布にほぼ影響しない場合は、片方向連成で十分な精度が得られます。

全体構成とジオメトリ

流体・固体・ヒートマットの寸法や位置は、すべて1つの設定ファイル（system/include/caseSettings）にまとめてあります。ここの数値を変えるだけで、メッシュ生成（blockMesh）と領域分割（topoSet）の両方に反映される作りです。

固体・流体・ヒートマットの3領域は、単一のブロックメッシュをtopoSetで3つのcellZoneに切り分け、splitMeshRegionsで実際の別リージョンへ分割するという、OpenFOAM公式チュートリアル（multiRegionHeater）と同じ考え方で作っています。

実際のcaseSettingsは次のようになっています。寸法や発熱条件はここだけを見れば分かるようにしてあり、格子数nx,ny,nzや各領域の切り出し座標は#evalで自動計算させています（詳しい設定は別記事にする予定です）。

// system/include/caseSettings（寸法・発熱条件を一元管理）
domainLx  1.0;    // 流体ドメイン X [m]（1000mm）
solidLx   0.2;    // 固体 X [m]（200mm）
solidLz   0.4;    // 固体 Z（高さ）[m]（400mm）
cellSize  0.05;   // メッシュ1辺の目安 [m]

heaterWattage    20.0;  // ヒーター出力 [W]
heaterOnTime     300;   // ON継続 [s]（0〜5分）
heaterOffTime    600;   // OFFのまま保持 [s]（5〜10分）
heaterMatWidth   0.10;  // ヒートマットの幅 [m]
heaterMatHeight  0.20;  // ヒートマットの高さ [m]

// 格子数や各cellZoneのbox座標は #eval で自動計算
nx  #eval #{ round($domainLx / $cellSize) #};

// system/include/caseSettings（寸法・発熱条件を一元管理）

domainLx 1.0; // 流体ドメイン X [m]（1000mm）

solidLx 0.2; // 固体 X [m]（200mm）

solidLz 0.4; // 固体 Z（高さ）[m]（400mm）

cellSize 0.05; // メッシュ1辺の目安 [m]

heaterWattage 20.0; // ヒーター出力 [W]

heaterOnTime 300; // ON継続 [s]（0〜5分）

heaterOffTime 600; // OFFのまま保持 [s]（5〜10分）

heaterMatWidth 0.10; // ヒートマットの幅 [m]

heaterMatHeight 0.20; // ヒートマットの高さ [m]

// 格子数や各cellZoneのbox座標は #eval で自動計算

nx #eval #{ round($domainLx / $cellSize) #};

発熱源を「セルではなく面」で与える

最初は、ヒートマットを固体のサブゾーン（体積）として扱い、fvOptionsの体積発熱源（scalarSemiImplicitSource）で20Wを与えていました。しかし、これは物理的には正しくありません。ヒートマットは本来「面」に貼り付いた発熱体であり、セル（体積）に均等にばらまかれる発熱ではないからです。

そこで、ヒートマットを固体・流体とは別の独立したリージョンにしました。厚さ1セル分の薄い直方体で、2つの面を持ちます。

外側面（heaterMat_to_fluid）：熱量の境界条件（externalWallHeatFluxTemperature, mode=power）
内側面（heaterMat_to_solid）：通常の熱伝導結合（固体へ伝わる面）

外側面はケーシングで室内空気とは断熱されている想定とし、投入した熱量はすべて内側面から固体へ伝導します。これにより、「面に熱量を与える」という物理的に自然な表現になりました。時間変化する発熱は、この熱量境界条件のQ（パワー）を時刻歴のテーブルとして与えることで表現します。実際の設定（system/heaterMat/changeDictionaryDictの温度Tの境界条件）は次の通りです。

heaterMat_to_fluid          // ヒートマットの外側面
{
    type   externalWallHeatFluxTemperature;
    mode   power;           // 熱量[W]で与える
    Q      table            // 時間変化テーブル
    (
        (0     20)          //   0秒: 20W
        (300   20)          // 300秒: 20W（ここまでON）
        (302   0)           // 302秒: 0W（切替）
        (600   0)           // 600秒: 0W（OFF保持）
    );
    kappaMethod  solidThermo;
    value        uniform 293.15;   // 初期20℃
}

heaterMat_to_fluid // ヒートマットの外側面

{

type externalWallHeatFluxTemperature;

mode power; // 熱量[W]で与える

Q table // 時間変化テーブル

(

(0 20) // 0秒: 20W

(300 20) // 300秒: 20W（ここまでON）

(302 0) // 302秒: 0W（切替）

(600 0) // 600秒: 0W（OFF保持）

);

kappaMethod solidThermo;

value uniform 293.15; // 初期20℃

}

Q（W）を時刻とセットの表で与えるだけで、「0〜300秒は20W、300秒で切って以降は0W」という矩形波状のON/OFFが表現できます。

輻射（fvDOM）を含める理由

固体表面から周囲への熱の逃げ方は、対流だけでなく輻射（ふく射）も無視できません。特に、局所的に温度が上がった面からの熱損失を正しく評価するには、輻射を含めた方が物理的に妥当です。

今回はfvDOM（Discrete Ordinates Method）という輻射モデルを使い、固体表面と周囲壁面をグレー拡散面（輻射率0.85〜0.9）として扱いました。空気自体はほぼ輻射を吸収・放射しないため、吸収率・放射率は小さい値（0.01）にしています。

FrontISTRへ温度を引き継ぐ

OpenFOAM側で保存された各時刻の固体温度分布（セル中心の値）を、FrontISTR側の節点（構造解析用メッシュの点）へ写像します。手順は次の通りです。

OpenFOAMのpostProcess -func writeCellCentresで、固体・ヒートマットの各セルの中心座標を書き出す
各FrontISTR節点について、最近傍のOpenFOAMセルの温度を距離で重み付けして平均する（逆距離加重平均）
得られた節点温度を、FrontISTRの!TEMPERATUREカードとして書き込む
熱ひずみは $\varepsilon_{th}=\alpha\,(T-T_{\mathrm{ref}})$ で計算される（FrontISTRの!REFTEMP機構）。基準温度 $T_{\mathrm{ref}}$ はOpenFOAM側の初期温度（20℃）に合わせている

底面（floor）は、X・Y・Z方向すべて固定しています。これは、OpenFOAM側でも同じ面を熱的な設置面として扱っていることと対応させたものです。

この写像とFrontISTR実行は、次のPythonスクリプト群で自動化しています（それぞれの中身は連成の詳細記事で解説します）。

box_mesh.py：OpenFOAMの固体と同じ寸法・同じ分割のFEMメッシュを作る
openfoam_temperature.py：各時刻のセル中心温度を読み、FEM節点へ逆距離補間する
fistr_case.py：FrontISTRの.msh/.cntを書き出し、fistr1を実行して変位を読む
run_thermal_expansion_timehistory.py：全時刻を順に回して、変位の時刻歴とグラフを作る

実行は次の1コマンドだけです。

python3 run_thermal_expansion_timehistory.py \
    --of-case ../../101_0_openfoam_cht_radiation_box

1 2	python3 run_thermal_expansion_timehistory.py \ --of-case ../../101_0_openfoam_cht_radiation_box

得られた結果

0〜600秒（0〜5分がヒーターON、5〜10分がOFF）の計算結果です。まず、加熱のピークである300秒（ヒーターOFFになる瞬間）を見てみます。

ヒートマット付近が局所的に温まり、上昇気流（自然対流のプルーム）が発生しています。右のFrontISTRモデルでも、ヒーター側だけが赤く高温になり、反りが生じているのが分かります。

次に、冷却が進んだ600秒（計算終了時）です。

ピーク時（300秒）より温度は下がっていますが、熱が固体内部に広がり、分布はより均一に近づいています。

0〜600秒の温度・変形の時刻歴をアニメーションにすると、次のようになります（左がOpenFOAMの温度場、右がその温度を引き継いだFrontISTRモデル）。

変位と温度の時刻歴をグラフにすると、次の通りです。

主な数値は次の通りです。

時刻	固体の最高温度	上面の平均鉛直変位	上面の変位のばらつき（反りの指標）
0秒（開始）	20.00℃	0 mm	0 mm
300秒（ヒーターOFFの瞬間）	20.62℃	0.00053 mm	0.00010 mm
400秒（冷却中）	20.35℃	0.00053 mm	0.00016 mm
600秒（終了）	20.19℃	0.00052 mm	0.00021 mm

ここで面白いのは、上面の平均的な鉛直変位（全体としての膨張量）はヒーターが切れた直後からほぼ横ばいなのに対し、「反り」の指標である変位のばらつきは、ヒーターを切った後もじわじわ増え続けていることです。

これは、ヒーターを切った瞬間に温度分布が均一に戻るわけではなく、ヒートマット側に溜まった熱が固体内部を伝わって反対側へ広がっていく途中だからだと考えられます。全体の膨張量は温度の「平均」にほぼ追従しますが、反り（曲げ）は温度分布の「偏り」に敏感なため、平均が落ち着いた後もしばらく変化し続ける、という熱応答の時間差が見えた形です。600秒時点の画像で温度分布がピーク時より均一に近づいているのも、この「熱が広がっていく途中」の様子と対応しています。

なお、20W・10分間という条件では、変位量は数百nm〜1μm程度とごく小さい値になりました。これは物理的に妥当な結果で（鋼材の塊を20Wで温めてもすぐには大きく膨張しません）、より大きな変形を見たい場合はヒーターの出力や加熱時間を増やせば確認できます。