FrontISTR(EasyISTR)：非定常伝熱解析（円柱の冷却）
1.概要
この記事はFrontISTR（EasyISTR）で図1-1のような円柱の冷却解析を行ったときの記録です．解析のモチベーションとしては，300℃に加熱された十分に長い鉄の棒を水に漬けて冷却したとき，何秒くらいで素手で安全に触れるくらいまで冷却されるのか？というものです（図1-2）．
図1-1. 円柱の冷却シミュレーション：非定常伝熱解析
図1-2. 解析モチベーション
今回は主にFrontISTRとEasyISTRを使用していますが，FrontISTRはFEM解析を実施する際のソルバーで(参考[1])，EasyISTRはFrontISTRのプリソフトになります(参考[2])．EasyISTRではモデル＆メッシュ作成が出来ないので，その他のソフトで別途用意する必要があります．そのため今回は，以前書いた記事(参考[3])で作成したモデルとメッシュを使用しています（Salome-Mecaで作成）．また，ポスト処理（結果確認）にはParaView(参考[4])を使用しました．使用した各ソフトの機能をまとめると，図1-3のようになります．必須ではありませんが，LibreOfficeやgnuplotなどのソフトを入れておくと，EasyISTRの機能をフルに使えるので私はこれらのソフトも入れています．
図1-3. 各ソフトの機能と役割
今回の解析はmm, ton, secの単位系で実施していますが，物性値や解析条件を記載する際はMKS単位系(m, kg, sec)(SI単位系)での値も併記しました．以下でEasyISTRを利用したFrontISTRでの解析実施方法を説明します．
2.作業ディレクトリの設定とメッシュの取り込み
2-1.作業ディレクトリの設定
まずは解析を実施するディレクトリを設定します．ここで設定した作業ディレクトリの中に解析インプットファイルや解析結果ファイルが出力されます．既に作成済みのディレクトリで解析を行っても良いですが，今回はEasyISTR上から新たに作業用のディレクトリを作成し，その作成したディレクトリを作業用のディレクトリとして設定することにします．その手順は下記になります（図2-1）．
- workFolder移動をクリック
- 作業用ディレクトリを作るディレクトリを選択：今回は「1_work」とした
- 右クリックから「新しいfolderを追加」を選択
- folder名を入力：今回は「1_cylinder_heat」とした
- 新たに作ったディレクトリ「1_cylinder_heat」を選択
- 「選択dirを作業folderに設定」をクリック
図2-1. 作業ディレクトリの設定
2-2.メッシュファイルの取り込みとメッシュ変換
次はメッシュファイルを取り込んで，FrontISTR形式のメッシュファイルに変換します．今回はSalome-Mecaで作成したunv形式のメッシュファイルをEasyISTRに取り込んで，FrontISTR形式のメッシュファイルに変換します．メッシュファイルを取り込む前は，ダミーの解析モデルが表示されていますが（図2-2），メッシュファイルを取り込むと取り込んだメッシュファイルのモデルが表示されます．
まずは，取り込むメッシュファイルを指定します．その手順は下記です（図2-2, 2-3）．
- 設定項目から「FistrModel.msh」を選択（このmshファイルがダミーのモデルから取り込むメッシュファイルに置き換わる）
- 「unv2fistr」を選択：今回はunv形式のメッシュフォーマットからFrontISTR形式のメッシュフォーマットに変換する
- 「参照」をクリックし，取り込みたいunvファイルを指定する．
図2-2. 取り込むメッシュファイルの指定
- 取り込むunvファイルを選択：今回はcylinder_ununiform.unvを選択
- 「Open」をクリックし，選択完了
- ファイル変換をクックし，unv形式からFrontISTR形式へメッシュフォーマットを変換
- 「mm ⇔ m」などの単位変換をしたい場合はスケール変換もできる（今回はmmでモデル作成しているので実施せず）
図2-3. unvファイルの選択とメッシュフォーマット変換
メッシュ変換とモデルの読み込みに成功すると，モデルの表示画面に取り込んだモデルが表示されます．メッシュ内容のところにメッシュのサマリーが表示されるので，ちゃんとモデルが読み込めているか簡単なチェックができます．（図2-4）．ちなみに今回は下記のような円柱のモデルを作成しており，きちんとメッシュ変換とモデルの読み込みができていることが分かります．
- 円柱の寸法
- 半径：50 [mm] (=0.050 [m])
- 高さ：300 [mm] (=0.300 [m])
| 個数 | |
|---|---|
| 6面体2次要素 | 9,600 |
| プリズム2次要素 | 400 |
| 全要素数 | 10,000 |
| 全節点数 | 43,021 |
図2-4. モデルの読み込みとメッシュサマリー
2-3.解析に使用する節点・面・要素集合
次に，解析に使用する節点集合，面集合，要素集合を作成します．以前書いた記事ではその手順を省略しましたが，Salome-Mecaでこれらの集合を予め作っておくこともできます．今回はSalome-Mecaで解析に使用する集合を予め作った上で，EasyISTRに読み込ませたため，その集合がここに引き継がれています．今回は節点・面・要素集合の作り方の手順は省略しますが，EasyISTR上からもこれらの集合はマウスピックすることで簡単に作成できます（図2-5）．
図2-5. 解析に使用する節点・面・要素集合
今回解析に使用するのは下記3つの集合です．これらの集合名は自分で任意に付けたものです．全節点集合は初期温度を設定するときに使用し，円柱側面の面集合は冷却面を指定するときに使用し，全要素集合は材料特性を割り当てるときに使用しています．
- 全節点集合：node_all → 初期温度の指定に使用
- 円柱側面の面集合：side_face → 冷却面の指定に使用
- 全要素集合：elem_all → 材料特性の割り当てに使用
それぞれの集合がきちんと意図した通りに作成できているかを確認していきます．各集合は設定項目の「GROUP編集」を選択することで確認できます．また，この「GROUP編集」では各集合の名前を変更したり，不要な集合を削除することもできます．各集合を確認した結果を図2-6，2-7，2-8に示します．各集合が意図した通りに作成できていることが確認できます．
図2-6. 全節点集合（node_all）の確認
図2-7. 円柱側面の面集合（side_face）の確認
図2-8. 全要素集合（elem_all）の確認
3.解析の設定
3節では下記3点を実施します．
- 解析タイプの選択（3-1節）
- 材料特性の割り当て（3-2節）
- 境界条件の設定（3-3節）
3-1.解析タイプの選択
今回は加熱された円柱の温度履歴を求めたいので，解析タイプは「熱伝導（動解析）」を選択します．この解析タイプを選択すると非定常伝熱解析を実施できます．これを設定するための手順は下記になります（図3-1）．
- 「解析の種類」を選択
- 「解析の種類」から「熱伝導（動解析）」を選択
- 「設定」をクリックし，設定を保存
図3-1. 解析タイプの選択
今回は解析タイプに「熱伝導（動解析）」を選択したので，時間増分や計算時間などを指定する必要があります．ここで言う「計算時間」とは何秒後の現象まで解析を実施するのかを意味します．また，他の多くのCAEソフトと同じようにFrontISTRでは時間刻みを固定するのか，解析時間を短縮するために徐々に時間増分を増加させていくのかを選択できます．今回は解析時間を短縮するために，時間増分を固定せずに解析を実施します．時間増分を固定する際は，単位時間当たりの伝熱量がメッシュ刻みに比べて大きすぎる場合などに計算が発散し，エラーを起こしてしまう場合があるので注意しましょう．
今回は下記のように時間増分などを設定しています（図3-2）．また，初期時間増分＞最小時間増分となるように設定しないと解析エラーで最初から動かないので注意しましょう．慌てて解析設定をしていると良く間違えます．許容変化温度は1時間増分あたり（1計算stepあたり）に許容する温度変化量です．ここで設定した許容変化温度を超えない範囲で，FrontISTRが時間増分を解析進行に伴い，徐々に大きく取ってくれます．時間増分を固定しない場合に有効な設定となります．
- 初期時間増分：0.01 [sec]
- 計算時間：200 [sec]
- 最小時間増分：0.001 [sec]
- 許容変化温度：20.0 [℃]
図3-2. 時間増分などの設定
また，今回は非線形最大反復数と収束判定値はデフォルトの設定から変更していません．基本的にこの2つの項目は変更する必要がありません．解析を無理やり収束させたい場合などに，非線形最大反復数と収束判定値を大きくすることで解析を収束させることができます．非線形最大反復数は大きくしても解析時間が伸びるだけで解析精度は落ちませんが，収束判定値は解析精度の観点からあまり大きく取るべきではありません．
解析が収束しないときは，この2つのパラメータを変更することを考えるのではなく，まずはモデリングやその他の解析条件設定を見直した方が良いです．伝熱解析で解析が収束しない原因の殆どは，メッシュ刻み対し単位時間当たりの伝熱量が大きいことにあります．例えば下記のような状況が考えられます．
- 時間刻みが大きすぎる
- メッシュが粗すぎる
- 発熱境界条件で指定している発熱量が大きすぎる
- 熱伝達境界，熱流束境界で指定している値が大きすぎる
- 温度境界条件があまりにも非合理的過ぎる（非物理的過ぎる）
- 物性値の非線形性が強すぎる
物性値の非線形性以外は，多くの場合，解析者のモデリング（メッシング）や解析条件設定に問題があります．物性値が大きく異なる異種材料界面や非常に大きなピーク値（または非常に小さなピーク値）を持つ物性値を取り扱う際などは，どう頑張ってもモデリングが上手くできない場合もあります．こういうどうしようも無い問題についてどうしても計算したいときに，非線形最大反復数や収束判定値の見直しをするべきです．
メッシングをやり直すのは結構な手間がかかる場合があるので，考えている問題設定（物理現象）に応じてメッシュの細かさを予め決めておくのが適当です．今回の解析は，十分に長い高温の円柱を水に漬けて円柱側面から冷却する，という結構強い冷却現象を考えているので，冷却面に相当する円柱側面近辺はメッシュを細かくしています．また，どのようにして円柱を水に漬けるのかにもよりますが，円柱側面からほぼ均等に冷却されることが予想されるため，対称性の高いメッシングをしています．
3-2.材料特性の割り当て
次は，材料特性の割り当てを行います．今回はモデル全体が単一の材料（鉄）であるとします．EasyISTRでFrontISTR用のインプットファイル作る際は，予め用意した材料データベース（csvファイル）を指定し，その中から材料種を選択するというやり方で材料特性を割り当てます．また，デフォルトで用意されている材料データベースから変更しない場合は図3-3，3-4の手順は不要です．今回の解析はmm，ton，secの単位系で解析を実施しています．一方，デフォルトで用意されている材料データベースはMKS単位系の値で作成されているので，自分で用意したデータベースファイルを解析に使用します．使用する材料データベースを変更するための手順は下記です（図3-3，3-4）．
-
設定項目の「材料物性値」を選択
-
「参照」から使用する材料DBを選択（デフォルトではmat.csvが選択されている）
図3-3. 使用する材料DBの選択（変更） -
任意のディレクトリに置いた自分で作成した材料データベースファイル（csvファイル）を選択：今回の場合は「my_material.csv」というファイル名
-
「Open」をクリックし，選択完了
図3-4. 自分で用意した材料DBの選択
次は材料特性を割り当てるための手順です．この設定は使用する材料データベースをデフォルトのデータベースから変更していなくても必要な設定です．材料特性の割り当てをするための手順は下記になります（図3-5，3-6）．
- 図3-5の手順は材料特性の割り当てをする要素集合を選択するための手順
- 材料特性の割り当てる要素集合を選択：今回の場合はelem_all（全要素集合）を選択
- 材料特性の定義：「選択」をクリックし選択した要素集合を右に移動させる
- 「設定」をクリックし，設定を保存
図3-5. 材料特性を割り当てる要素集合の選択
- 図3-6は図3-5の手順で選択した要素集合に対して，要素集合を実際に割り当てる手順
- 図3-5の手順で選択した要素集合を選択：今回の場合はelem_all
- 材料名から「Steel」を選択し，材料特性を割り当てる
- 「設定」をクリックした後に「物性値の確認」をクリックすると割り当てた物性値を確認できる
図3-6. 材料特性の割り当てと割り当てた材料特性の確認
今回は下記のような物性値を採用しました．この物性値はEasyISTRがデフォルトで用意している材料データに単位変換を施しただけのものになります．常温における鉄の一般的な物性値とそう変わらないものだったので，簡単のためにそのまま流用しています．解析にはmm，ton，sec単位系のものを使用しています．
| mm, ton, sec単位系 | MKS単位系 | |
|---|---|---|
| 密度 |
|
|
| 熱伝導率 |
|
|
| 比熱 |
|
|
3-3.境界条件の設定
3-3-1.初期温度の設定
ここでは解析モデルに初期温度を設定する手順を説明します．多くのCAEソフトがそうであるように，FrontISTRでも温度データは要素の積分点ではなく節点で持ちます．そのため，温度の条件は節点集合に対して与えます．今回は，鉄の棒が初期状態において300℃の均一温度状態にあるものとします．全節点集合（node_all）に対して，初期温度を設定するための手順は下記になります．
- 初期温度を設定する節点集合を選択：今回の場合は全節点集合(node_all)を選択
- 初期温度：300 [℃]
- 「設定」をクリックし，設定を保存
図3-7. 初期温度の設定
3-3-2.熱伝達境界の設定
今回は円柱（鉄の棒）の側面のみから冷却される解析条件を考え，円柱の上面と底面からは冷却されないものとします（円柱の上面と底面は断熱）．この解析条件は鉄の棒が十分に長いときに良い近似となります．また，今回は鉄の棒を水に漬けて撹拌しながら冷却するときの熱伝達係数を
熱伝達係数は，冷媒（水）との接触状況（撹拌の有無）や沸騰状態，鉄の表面温度，鉄の表面粗さ，水へ混入している不純物の割合などによって変化しますが，これはそれなりに妥当な値です(参考[5])．熱伝達係数に影響を及ぼすパラメータは多岐に渡りますが，一番効くのは撹拌状況（沸騰状態）と鉄の表面温度です．水の中での撹拌状況（沸騰状態）までモデル化したり，実験による測定をしだすとキリが無いので今回は簡単のために熱伝達係数をエイやで決め打ちしています．また，熱伝達係数の鉄の表面温度への依存性も無視しています．
円柱側面へ熱伝達境界を設定するための手順は下記になります（図3-8, 3-9）．手順としては熱伝達境界を与える面集合を指定し（図3-8），その面集合に対し熱伝達係数と雰囲気温度（水の温度）を指定します（図3-9）．
- 「境界条件」から「SFILM（熱伝達率）」を選択
- 熱伝達境界を与える面集合を選択：今回は円柱側面の面集合side_faceを選択
- 「選択」ボタンを押し，右側の「設定するgroup」へ移動させる
- 最後に「設定」ボタンをクリックし，設定を保存
図3-8. 熱伝達境界を与える面集合の指定
- 「SFILM（熱伝達率）」から熱伝達境界を与える面集合を選択：今回はside_faceを指定
- 雰囲気温度：20 [
\rm {}^\circ\mathrm{C} - 熱伝達率：5.0 [
\rm mW/(mm^2 \cdot {}^\circ\mathrm{C}) - 最後に「設定」をクリックし，設定を保存
図3-9. 雰囲気温度と熱伝達係数の指定
3-3-3.熱流束境界の設定
今回は鉄の棒（円柱）の上面と底面に断熱境界を与えます．このような解析設定は鉄の棒が十分に長いときに良い近似となります．断熱境界を設定するには，熱流束が時間によらず常にゼロになるように設定すれば良いです．しかしながら，多くのCAEソフトがそうであるように，FrontISTRでも境界条件が与えられていないモデルの外面は自動的に断熱境界（熱流束が常にゼロ）として取り扱われます．FrontISTRで解析するのは殆ど初めてだったので，断熱境界を明示的に設定した場合としなかった場合とで計算結果が変わらないか念のため確認しておきました．その結果，両者ともに計算結果は変わらなかったので，断熱境界の取り扱いは他のCAEソフトと同じと考えて間違いないと思います．
ここでは，断熱境界を明示的に与える場合の手順を省略しつつ簡単に説明します．本来なら必要ない手順です．その手順は下記になります（図3-10）．
- SFLUX（面熱流束）を選択
- 面集合の中から熱流束境界を指定する面集合を選択：今回は円柱上面のtop_faceと円柱底面のbottom_face
- その後，top_faceとbottom_faceに対し，熱流束をゼロと指定する（設定画面のスクリーンショットは省略）
図3-10. 断熱境界の明示的な指定方法
4.ジョブ投入と結果の確認
4-1.ジョブ投入（解析実施）
3節までで解析の設定は終わったので，FrontISTRに投入します．コンソール画面からコマンド入力してEasyISTRを経由せずにジョブを投入することもできますが，今回はEasyISTRからジョブを投入する方法を説明します．
私の環境の場合は，EasyISTRを経由せずにコンソール画面からコマンド入力して直接ジョブ投入した方が30%程度計算時間が短縮されました．計算時間が気になる方はコンソール画面からジョブ投入すると若干計算時間が短縮されるかもしれません．また，今回は並列計算ではない普通のジョブの投入方法について説明します．ジョブの投入方法は下記になります（図4-1）．
- 「solver」を選択
- 結果出力頻度（何計算ステップごとに解析結果を出力するか？）を指定：今回は「1」を指定（全ステップの計算結果が出力される）
- 「設定」をクリックし，「solver」の設定を保存
- 「FrontISTR実行」をクリックして解析実行
図4-1. ジョブ投入手順
4-2.ジョブ進行状況の確認
EasyISTRには，明らかにおかしな計算結果を行っていないかを確認するために，ジョブの進行状況（最大温度，最小温度，解析ステップ数，時刻）を確認するための機能が備わっています．解析の設定ミスを計算途中に簡易的に確認できるのは非常に便利です．解析実施のための必須の機能ではないため，gnuplotのインストールとパスの設定が必要です．ここでは，ジョブの進行状況を確認するための手順を説明します．その手順は下記です（図4-2）．
- 「step状況表示」をクリック
- ある程度解析が進むとgnuplotの画面が立ち上がりグラフが表示される
- 横軸：解析ステップ数
- 縦軸（第一軸）：時刻
- 縦軸（第二軸）：温度
図4-2. ジョブの進行状況（最大温度，最小温度，解析ステップ数，時刻）
4-3.解析結果の確認
ParaViewがインストール済みで，EasyISTRの設定ファイルにParaViewの実行ファイルへのパスをきちんと記述している場合は，EasyISTRからParaViewを起動させることができます．EasyISTRからParaViewを起動する手順は下記になります（図4-4）．
- 「post」を選択
- 「ParaView起動」をクリック
- 解析結果を読み込み済みの状態でParaViewが立ち上がる
図4-4. ParaView起動手順
今回はParaViewの操作方法自体は省略しますが，ParaViewでは3Dモデルのコンター図など表示させることができます．そのコンター図をアニメーションにしたものが図4-5になります．図4-5には温度を表示しています．
図4-5のアニメーションは各計算ステップで出力された温度データを，等間隔に繋いでアニメーション化したものです．各計算ステップの時間増分は固定では無く，計算ステップの進行に伴い，時間増分は増加しています（図4-6）．そのため，図4-5に示したアニメーションは，実際の時間間隔と異なる間隔でアニメーションが進行しています．一方で，円柱の表面から温度が低下する様子は確認でき，明らかにおかしな解析になっていないことは確認できます．
図4-5. 温度履歴のアニメーション
図4-6. 計算ステップと時間経過（時刻）の関係
図4-7に円柱の中心と円柱の表面の2点における温度の時刻歴データを示します．図4-7の横軸の時間は計算ステップではなく，実際の時間スケールに合わせています．このグラフから円柱全体（鉄の棒）の温度が人肌程度（約40℃）になるまでに約100sec程度かかることが，少なくとも計算上，分かりました．つまり，300℃に加熱された長い鉄の棒は，水に漬けて2分程度冷却すれば，素手で触っても安全であることが分かります（安全にみるなら5分くらいは冷却する）．
図4-7. 円柱の中心と円柱の表面における時刻歴データ
感想
今回は直径が100mm（半径50mm）の300℃に加熱された鉄の棒を水に漬けて冷却するというものでした．100mmというと成人女性の手のひらサイズくらいの大きさです．これくらいのサイズの鉄の棒に対する冷却実験を何回かやったことがありますが，常温になるまでに数分掛かります．
私の数少ない経験を以てして今回の計算結果が妥当であると言い切ることはできませんが，少なくとも致命的な解析ミスや出鱈目な解析条件を入力していることはなさそうです．今のところ，鉄の棒に対する冷却実験は実際に行う予定はありませんが，今後の記事で意図した解析が本当にできていたのか，もう少し詳しく見ていこうと思います．
数年前にチュートリアルだけ動作させたことはあったものの，今回初めて割と真面目にEasyISTRとFrontISTRを使ってみました．EasyISTRもFrontISTRも国産のCAEソフトなので応援する意味も込めて，今後使っていこうと思います．EasyISTRのインターフェースはシンプルで分かりやすく使いやすかったです．構造系のCAEソフトを触ったことがある人であれば，恐らくマニュアルも殆ど見ずに設定できるのではないかと思います．
環境
・FrontISTR: ver.5.4
・EasyISTR: ver.3.34
・gnuplot: ver.5.2.4
・ParaView: ver.5.11.0
・Windows 10
参考
[1]. FrontISTR：公式HP
[2]. EasyISTR：公式HP
[3]. Salome：任意の密度関数で不均一なメッシュを切る
[4]. ParaView：公式HP
[5]. 三塚, "高温鋼材の冷却技術", 鉄と鋼, 79(1993), p.N405-N416. (DOI)
[6]. gnuplotの各バージョンのインストーラーへのリンク
Discussion