【OpenRadioss入門】熱源移動による熱解析の設定方法

こんにちは（@t_kun_kamakiri）

Altairのサイトにこれまで商用利用としていたRadiossを2022年9月からオープンソース化するという内容の記事が投稿されていました。RadiossといえばLS-DYNAやPAM-Crashなどと並ぶAltairの衝突系（衝撃解析）の解析ソフトのイメージが強いですが、それが無償で使えるようになったというのはなかなかの衝撃でした。

今回解説するのはこちらの熱源が移動する解析です。
設定としては適当で「とりあえず計算できればいいや！」くらいの感覚ですので実現象を再現しうる精度にはなっていないでしょうが、設定方法は参考になると思います。

あわせて読みたい

Altair社のマニュアルに熱解析の例題があり、そちらを参考にしています。
RD-T： 5300 熱解析

OpenRadiossの公式

Altairの例題集は大変参考になります。
しかし、モデル作成からRadiossは調べてキーワードというカードで設定を記述する必要があります。
そのため、GUI操作で順番にボタンをポチポチ押して行って設定するという感じではなく、適切な解析条件を自分で選択して解析設定を行う必要があります。
無償のOpenRadiossにはAltair社のHyper WorksのようなGUIが無いため、設定に苦労します。本記事では、できるだけテキスト編集をメインにOpenRadiossの基本的な設定を解説します。

WSL2

あわせて読みたい

インストール方法

モデルはこちらです。

github

解析設定の手順
モデル作成
メッシュ作成
節点と要素
材料特性の設定
プロパティの設定
PARTの作成
熱源の拘束条件
熱源の強制速度
100×100の正方形の面の拘束
接触条件の設定
熱伝達の設定
節点の履歴
Engine設定
解析実行
結果処理
まとめ

解析設定の手順

解析設定にはStarterファイルとEngineファイルの作成が必要です。
今回は以下のファイル名で設定ファイルを作成します。

001thermal_0000.rad：Starter ファイル（解析設定ファイル）
001thermal_0001.rad：Engine ファイル（解析制御と出力設定ファイル）

あわせて読みたい

ユーザーマニュアル

モデル作成

モデル作成はCADソフトを使います。
無償で使えるCADソフトはぱっと思いつくのは以下があります。

FreeCADはその名の通りCADソフトですが、メッシャーもついており出力形式もSTEP、igs、stlなど他のソフトでも読み込める形式で出力が可能です。

Salome-mecaは、shaperというCAD機能もありますしメッシャーもあり、メッシュ生成後に節点を接合したり修正したり非常に多くの機能があります。FreeCADと同じように様々な出力形式に対応しています。

Fusion360は、個人的に使ったことないのでホームページを見てください_(._.)_

今回はFreeCADを採用してモデルだけを作成してSTEPファイルをとして出力することを考えます。
※FreeCADの詳しい解説はここではしません。
メッシャーは別のソフトを使います。

作成するモデルは以下の寸法の100×100の正方形の面と、5×5の正方形の面です。
解析設定で両者の厚みを1mmに設定しますが、モデルは厚みのない面（サーフェス）として作成します。
ですので、5×5の正方形の面は1mm離してモデル作成をします。

FreeCADを起動してX-Y平面にスケッチを開始します。2つの面を作成します。作成した2つ面を選択してSTEPファイルとして出力します。

今回モデル作成はFreeCADで行いましたが、STEPファイルが出力できるソフトであれば何でもよいです。

メッシュ作成

続いてメッシュ作成に移ります。
メッシュ作成もFreeCADや別のソフトでもできるのですが、今回はGmshを使います。

Gmshを使う理由はOpenRadiossの形式として出力に対応しているからです。

あわせて読みたい

先ほど出力したSTEPファイルをGmshで読み込んで、グループを作ります。

メッシュ作成を行います。
メッシュのサイズ感は適当に「Min = 1mm」「Max = 10mm」とします。
以下のようにRadiossの出力形式（.rad)で出力します。
先ほど作成したグループを出力してくとOpenRadiossでの節点グループに使えるので出力するようにしておきます。ただし、節点グループはLS-DYNA?の形式で「*SET_NODE_LIST」のキーワードでカンマ区切りで出力されます。
これをOpenRadiossの10カラム区切りに編集する必要があります。
以下の流れで10カラム区切りにします。以上の作成した節点グループは使うかわからないですが、「/GRNOD/NODE/」として書いておくと良いでしょう。

Gmshから出力したファイルを確認すると以下のようになっています。

#RADIOSS STARTER
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
# Created by Gmsh, Radioss Mesh Interface by PaulAltair sharp@altair.com
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/BEGIN
D:\OpenRadioss\work\001_example\radioss\100_thermal\model
      2022         0
                  Mg                  mm                   s
                  Mg                  mm                   s
/NODE
         1                 -50                  50                   0
         2                 -50                 -50                   0
 （省略）
/SHELL/2000001
#SET_ELEMENT=PhysicalSurface11
         1         1        26        89         5
         2        26        37        90        89
         3        37        36        91        90
（省略）
/SHELL/2000002
#SET_ELEMENT=PhysicalSurface22
       605       530       539       554       534
       606       539       542       555       554
（省略）
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/GRNOD/NODE/1
#title                                                                                              
under_plate                                                                                       
# item_ID1  item_ID2  item_ID3  item_ID4  item_ID5  item_ID6  item_ID7  item_ID8  item_ID9 item_ID10
         2         1         3         4         9         8        15        10
         5        47        48        49        50        51         6        53
（省略）
/GRNOD/NODE/2
#title                                                                                              
thermal_plate                                                                                       
# item_ID1  item_ID2  item_ID3  item_ID4  item_ID5  item_ID6  item_ID7  item_ID8  item_ID9 item_ID10
       531       530       533       532       540       541       536       535
       534       539       544       545       546       549       550       551
（省略）
/END

#RADIOSS STARTER

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

# Created by Gmsh, Radioss Mesh Interface by PaulAltair sharp@altair.com

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

/BEGIN

D:\OpenRadioss\work\001_example\radioss\100_thermal\model

2022 0

Mg mm s

/NODE

1 -50 50 0

2 -50 -50 0

（省略）

/SHELL/2000001

#SET_ELEMENT=PhysicalSurface11

1 1 26 89 5

2 26 37 90 89

3 37 36 91 90

（省略）

/SHELL/2000002

#SET_ELEMENT=PhysicalSurface22

605 530 539 554 534

606 539 542 555 554

（省略）

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

/GRNOD/NODE/1

#title

under_plate

# item_ID1 item_ID2 item_ID3 item_ID4 item_ID5 item_ID6 item_ID7 item_ID8 item_ID9 item_ID10

2 1 3 4 9 8 15 10

5 47 48 49 50 51 6 53

（省略）

/GRNOD/NODE/2

#title

thermal_plate

# item_ID1 item_ID2 item_ID3 item_ID4 item_ID5 item_ID6 item_ID7 item_ID8 item_ID9 item_ID10

531 530 533 532 540 541 536 535

534 539 544 545 546 549 550 551

（省略）

/END

「/GRNOD/NODE」の節点グループは後ほど拘束条件や接触条件に使います。

節点と要素

モデルはメッシュという要素に分割されて作られています。

今回は扱うモデルは厚みのないシェル要素で、下の絵のように要素は3つもしくは4つの節点から構成されています。

節点は節点のIDと$x,y,z$ $x, y, z$ の座標を指定
要素は3つもしくは4つの節点IDを指定

節点が3つで構成される要素は「/SH3N」のキーワードを使い、節点が4つの場合は「/SHELL」というキーワードを使います。

#RADIOSS STARTER
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
# Created by Gmsh, Radioss Mesh Interface by PaulAltair sharp@altair.com
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/BEGIN
D:\OpenRadioss\work\001_example\radioss\100_thermal\model
      2022         0
                  Mg                  mm                   s
                  Mg                  mm                   s
/NODE
         1                 -50                  50                   0
         2                 -50                 -50                   0
 （省略）
/SHELL/2000001
#SET_ELEMENT=PhysicalSurface11
         1         1        26        89         5
         2        26        37        90        89
         3        37        36        91        90
（省略）
/SHELL/2000002
#SET_ELEMENT=PhysicalSurface22
       605       530       539       554       534
       606       539       542       555       554
（省略）
/END

#RADIOSS STARTER

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

# Created by Gmsh, Radioss Mesh Interface by PaulAltair sharp@altair.com

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

/BEGIN

D:\OpenRadioss\work\001_example\radioss\100_thermal\model

2022 0

Mg mm s

/NODE

1 -50 50 0

2 -50 -50 0

（省略）

/SHELL/2000001

#SET_ELEMENT=PhysicalSurface11

1 1 26 89 5

2 26 37 90 89

3 37 36 91 90

（省略）

/SHELL/2000002

#SET_ELEMENT=PhysicalSurface22

605 530 539 554 534

606 539 542 555 554

（省略）

/END

また、PART IDが以下のように対応していることを確認します。

001thermal_0000.rad

Gmshでメッシュ作成からRadioss形式(.rad)に出力した際にシェル要素には自動的にPART IDと関連付けされています。
100×100の正方形の面にはPART ID = 2000001が、5×5の正方形の面にはPART ID = 2000002がそれぞれ自動で付けられています。
こちらはこのままのIDでも良いですし、本記事のように後ほど開設するように変更しても良いです。

節点、シェル要素（ソリッド要素）がPARTと紐づいていることを確認したら、次に行うのはPARTに対するプロパティと材料特性の設定です。

材料特性の設定

OpenRadiossには多くの材料特性の設定が用意されており、こちらの材料一覧から確認することができます。

今回使うのは以下の材料設定です。

則名	グループ	タイプ	モデル概要	適用分野	粘性	直交異方性	組み込み熱オプション**	ひずみ速度	定式化
/MAT/LAW2 (PLAS_JOHNS)*	弾塑性材料	金属	Johnson-Cook材料モデルを使用した等方性弾塑性材料	金属、弾塑性材料、軟鋼、ひずみ速度依存、温度依存、等方および移動硬化、熱硬化性、低合金鋼、ST52、およびDP600	×	×	○	○	Lagrange

この材料則は、Johnson-Cook材料モデルを使用して等方性弾塑性材料を表します。

マニュアルを見ながら必要な項目を埋めていきます。
今回は使う材料カードは温度特性も解くことができるものですので、併せて「/HEAT/MAT」のキーワードも使います。

/HEAT/MAT
熱解析用の熱パラメータを記述します。
/THERM_STRESS/MAT
Radioss材料（シェルおよびソリッド）の熱膨張プロパティを追加するために使用します。

他にも/THERM_STRESS/MATという熱膨張を扱うキーワードも使えるのですが、ここでは使用しません。

MAT ID = 1は100×100の正方形の面の材料特性に使い、MAT ID = 2は5×5の正方形の面の材料特性に使うようにします。

##HWCOLOR materials 1 5
/MAT/PLAS_JOHNS/1
PLATE -Alumaterial1                                                                                 
#              RHO_I
              0.0028
#                  E                  Nu     Iflag
              7000.0                0.33         0
#                  a                   b                   n           EPS_p_max            SIG_max0
               260.0               450.0                 0.5                                        
#                  c           EPS_DOT_0       ICC   Fsmooth               F_cut               Chard
                                                                                                    
#                  m              T_melt              rhoC_p                 T_r
                                                        2.51               298.0
/HEAT/MAT/1
#                 T0             RHO0_CP                  AS                  BS     IfORM
               398.0                2.51                0.23                             1
#                 T1                  AL                  BL               EFRAC
                                                                                
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
##HWCOLOR materials 2 6
/MAT/PLAS_JOHNS/2
Source -Alumaterial1                                                                                
#              RHO_I
              0.0078
#                  E                  Nu     Iflag
            210000.0                0.33         0
#                  a                   b                   n           EPS_p_max            SIG_max0
               206.0               450.0                 0.5                                        
#                  c           EPS_DOT_0       ICC   Fsmooth               F_cut               Chard
                                                                                                    
#                  m              T_melt              rhoC_p                 T_r
                                                        2.51               298.0
/HEAT/MAT/2
#                 T0             RHO0_CP                  AS                  BS     IfORM
               800.0                2.51                0.23                             1
#                 T1                  AL                  BL               EFRAC

##HWCOLOR materials 1 5

/MAT/PLAS_JOHNS/1

PLATE -Alumaterial1

# RHO_I

0.0028

# E Nu Iflag

7000.0 0.33 0

# a b n EPS_p_max SIG_max0

260.0 450.0 0.5

# c EPS_DOT_0 ICC Fsmooth F_cut Chard

# m T_melt rhoC_p T_r

2.51 298.0

/HEAT/MAT/1

# T0 RHO0_CP AS BS IfORM

398.0 2.51 0.23 1

# T1 AL BL EFRAC

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

##HWCOLOR materials 2 6

/MAT/PLAS_JOHNS/2

Source -Alumaterial1

# RHO_I

0.0078

# E Nu Iflag

210000.0 0.33 0

# a b n EPS_p_max SIG_max0

206.0 450.0 0.5

# c EPS_DOT_0 ICC Fsmooth F_cut Chard

# m T_melt rhoC_p T_r

2.51 298.0

/HEAT/MAT/2

# T0 RHO0_CP AS BS IfORM

800.0 2.51 0.23 1

# T1 AL BL EFRAC

5×5の正方形が熱源なので800[K]としています。
100×100の正方形は初期は398[K]という設定です。

以下、こちらからの引用です。

１．これは、真の応力およびひずみ出力を伴うひずみ速度および温度効果を含む弾塑性材料モデルです。

２．このモデルでは、相当応力が塑性降伏応力よりも低い時、材料は線形弾性材料として挙動します。もっと高い応力値では、材料は塑性挙動で、応力は以下のように計算されます：

\begin{align*}
\sigma = \big(a+b\varepsilon_{p}^{n}\big)\bigg(1+c\ln\frac{\dot{\varepsilon}}{\dot{\varepsilon}_{0}}\bigg)\big(1-(T^{*})^{m}\big)
\end{align*}

$T^{*}=\frac{T-T_{r}}{T_{melt}-T_{r}}$

$\varepsilon_{p}$：塑性ひずみ
$\dot{\varepsilon}$：ひずみ速度
$T$：温度
$T_{r}$：雰囲気温度
$T_{melt}$：溶融温度

３．I_flag=1の場合、σy, UTSおよび $ε_{U T S}$ に実験工学応力およびひずみデータを入力でき、パラメータa、bおよびnが計算されてStarter出力ファイルに出力されます。

$ε_{p}$ が1つの積分点で $ε_{p}^{m a x}$ の値に到達してから、要素タイプに基づく場合：
シェル要素：対応するシェル要素が削除されます。

今回は以下のように設定しています。

SIG_Y = 0.09025：降伏応力
UTS = 0.175：工学応力UTS > σy
EUTS = 0.24：UTSにおける工学ひずみ
EPS_p_max = 0.75：破壊塑性ひずみ
ICC（設定なし）：デフォルトは1なのでひずみ速度効果あり

つまり、蘇生ひずみεpが0.75を超えるとシェル要素が削除されるというモデル設定です。
さらに詳しく知りたい方はマニュアルを参照ください。

テキスト設定の入手方法
GUIが無い場合はテキストによる設定の記述が必要となり、マニュアルを見ながら正確に設定を書くのは正直困難です。
設定のテンプレートがあれば楽ですよね。

その場合は、マニュアルの材料設定のページにいくと、以下のように材料設定が使われている例題の紹介があるので、どれかひとつ参考になるモデルからテキスト設定をコピペすることを考えます。

RD-E：0300 S字型ビームの圧壊

RD-E：0500 ビームフレーム

RD-E：0601 流体-構造連成

RD-E：0800 Hopkinsonバー

RD-E：4802 クラッシュボックス内のソリッドおよびスプリングスポット溶接

RD-E：5300 熱解析

例題集のモデルはモデルファイルへのアクセスから入手できるので、必要なキーワードはここから入手していけばよいでしょう。

以上で材料設定が終わりです。

※ここでは材料特性の設定を記述しただけで、PARTとは紐づいていません。

プロパティの設定

シェル要素に対して要素の定義と厚みの設定を行います。

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/PROP/SHELL/1
PID 1                                                                                               
#   Ishell    Ismstr     Ish3n    Idrill                            P_Thick_Fail
         1         1                                                            
#                 Hm                  Hf                  Hr                  Dm                  Dn
                                                                                                    
#        N   Istrain               Thick              Ashear              Ithick     Iplas
         1                           1.0                                                  
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/PROP/SHELL/2
PID 2                                                                                               
#   Ishell    Ismstr     Ish3n    Idrill                            P_Thick_Fail
         1         0                                                            
#                 Hm                  Hf                  Hr                  Dm                  Dn
                                                                                                    
#        N   Istrain               Thick              Ashear              Ithick     Iplas
         1                           1.0

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

/PROP/SHELL/1

PID 1

# Ishell Ismstr Ish3n Idrill P_Thick_Fail

1 1

# Hm Hf Hr Dm Dn

# N Istrain Thick Ashear Ithick Iplas

1 1.0

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

/PROP/SHELL/2

PID 2

# Ishell Ismstr Ish3n Idrill P_Thick_Fail

1 0

# Hm Hf Hr Dm Dn

# N Istrain Thick Ashear Ithick Iplas

1 1.0

ほとんどがデフォルト設定でよく、すべて設定しなければならないというわけではありません。
大事なところは以下です。

Ishell = 1の場合：/DEF_SHELLが定義されていない場合のデフォルト。Q4、変形モードと剛体モードに直交する粘弾性アワグラスモード（Belytschko）
Thick = 1.0：厚み設定

その他の設定は必要に応じてマニュアルを見て設定します。

※ここではプロパティの設定を記述しただけで、PARTとは紐づいていません。

PARTの作成

ここで先ほど記述したプロパティと材料特性をPARTに紐づけます。

PARTは要素（ビーム要素、シェル要素、ソリッド要素）で構成されており、PARTへはプロパティのIDと材料のIDを指定します。

設定は10カラム（10行）ごとに区切られているので、記述する位置を間違えないようにしましょう。「#」で始まる行はコメント行になるので、設定のメモを残しておくとわかりやすいです。

先ほどGmshでメッシュ作成からRadioss形式(.rad)に出力した際にシェル要素には自動的にPART IDと関連付けされていましたが、以下のように変更しておきましょう。

するとPART IDはそれぞれ100×100の正方形の面にはPART ID = 1が、5×5の正方形の面にはPART ID =2となります。
ここにプロパティのIDと材料のIDを書いてPARTと関連付けます。

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/PART/1
#part_title                                                                                         
Plate                                                                                     
#  prop_ID    mat_ID                         Thick
         1         1
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/PART/2
#part_title                                                                                         
Source                                                                                   
#  prop_ID    mat_ID                         Thick
         2         2

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

/PART/1

#part_title

Plate

# prop_ID mat_ID Thick

1 1

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

/PART/2

#part_title

Source

# prop_ID mat_ID Thick

2 2

これでPARTにプロパティと材料が紐づけられました。

解析設定としてPARTにプロパティと材料を指定すると使用されるため、例えば使う可能性があるプロパティや材料を設定の記述だけをテキストに書いておくこともできます。
その際はプロパティや材料はそれぞれIDが被らないようにしておきましょう。

熱源の拘束条件

今回の解析は5×5の正方形が熱源となりy軸方向には自由に動けるようにしますが、その他の軸と回転は拘束するようにします。
拘束条件を与えるキーワードして以下があります。

/BCS
並進および回転移動に対する節点グループ上の境界条件を定義します。

こちらは節点グループに強制速度を与えるため節点グループに「/GRNOD」を使う必要があります。しかし、節点グループはGmshで出力後に編集したように既に作成済みです。
「/GRNOD/NODE/1」の1がIDになっているので、それを使えば良いでしょう。

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/GRNOD/NODE/2
#title                                                                                              
thermal_plate                                                                                       
# item_ID1  item_ID2  item_ID3  item_ID4  item_ID5  item_ID6  item_ID7  item_ID8  item_ID9 item_ID10
       531       530       533       532       540       541       536       535
       534       539       544       545       546       549       550       551
（省略）

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

/GRNOD/NODE/2

#title

thermal_plate

# item_ID1 item_ID2 item_ID3 item_ID4 item_ID5 item_ID6 item_ID7 item_ID8 item_ID9 item_ID10

531 530 533 532 540 541 536 535

534 539 544 545 546 549 550 551

（省略）

拘束条件のキーワード「/BC」の「grnod_ID = 2」と節点グループのIDを与えることで拘束条件が設定されます。

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/BCS/2
PID 2 NODE GR                                                                                       
#  Tra rot   skew_ID  grnod_ID
   101 111         0         2

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

/BCS/2

PID 2 NODE GR

# Tra rot skew_ID grnod_ID

101 111 0 2

Tra：x,y,zの拘束設定（0がフリー、1が拘束）
rot：ωx、ωy、ωzの拘束設定（0がフリー、1が拘束）

なので節点グループID = 2はy方向のみ自由に動くことができる設定になっています。

熱源の強制速度

今回の解析は5×5の正方形が熱源となりy軸方向に0.5[m/s]で移動するものとします。
強制移動をさせる場合のキーワードして以下があります。

/IMPACC
節点のグループに対する強制加速度を定義します。
/IMPDISP
節点のグループに対する強制変位を定義します。
/IMPVEL
節点のグループに対する強制速度を定義します。

今回は/IMPVELを使います。

こちらも「/IMPVEL」の「grnod_ID = 2」と節点グループのIDを与えることで拘束条件が設定されます。

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/IMPVEL/3
Vy                                                                                                  
#funct_IDT       Dir   skew_ID sensor_ID  grnod_ID  frame_ID     Icoor
         1         Y         0         0         2                   0
#            Scale_x             Scale_y              Tstart               Tstop
                 1.0                 0.5                 0.01.00000000000000E+31

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

/IMPVEL/3

#funct_IDT Dir skew_ID sensor_ID grnod_ID frame_ID Icoor

1 Y 0 0 2 0

# Scale_x Scale_y Tstart Tstop

1.0 0.5 0.01.00000000000000E+31

funct_IDTには強制速度の時間変化を与える関数テーブルをしているする必要があります。

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/FUNCT/1
#title                                                                                              
curve1                                                                                              
#                  X                   Y
                 0.0                 1.0
1.00000000000000E+30                 1.0

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

/FUNCT/1

#title

curve1

# X Y

0.0 1.0

1.00000000000000E+30 1.0

これにより関数は常に値を1.0を取るようになっていますが、「/IMPVEL」キーワードの「Scale_y = 0.5」とスケールしているため、熱源移動のy方向の速度は常に0.5[m/s]にしています。

補足ですが、「/GRNOD/NODE/」のように節点グループの節点を全て記述するという方法もありますが、PARTに属するすべての節点を指定していることと同じなので以下のように書いても良いです。

/GRNOD/PART

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/GRNOD/PART/2
#title                                                                                              
thermal_plate                                                                                       
# item_ID1  item_ID2  item_ID3  item_ID4  item_ID5  item_ID6  item_ID7  item_ID8  item_ID9 item_ID10
        2

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

/GRNOD/PART/2

#title

thermal_plate

# item_ID1 item_ID2 item_ID3 item_ID4 item_ID5 item_ID6 item_ID7 item_ID8 item_ID9 item_ID10

「item_ID1 = 2」とPART ID を指定すると、PART ID = 2の全ての節点がグループ化されます。こちらの方が設定自体は楽ですよね。

PARTの1つの節点群をグループにしたい場合は、「/GRNOD/NODE/」を使うなどモデルによって使い分けると良いでしょう。

100×100の正方形の面の拘束

100×100の正方形の面は完全に拘束します。

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/BCS/1
PID 1 NODE GR                                                                                       
#  Tra rot   skew_ID  grnod_ID
   111 111         0         1

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

/BCS/1

PID 1 NODE GR

# Tra rot skew_ID grnod_ID

111 111 0 1

こちらも「grnod_ID」には節点グループの指定が必要ですが既に以下のように節点グループは作っています。

Tra：x,y,zの拘束設定（0がフリー、1が拘束）
rot：ωx、ωy、ωzの拘束設定（0がフリー、1が拘束）

なので節点グループID = 1は完全に拘束の設定です。

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/GRNOD/NODE/1
#title                                                                                              
under_plate                                                                                       
# item_ID1  item_ID2  item_ID3  item_ID4  item_ID5  item_ID6  item_ID7  item_ID8  item_ID9 item_ID10
         2         1         3         4         9         8        15        10
         5        47        48        49        50        51         6        53
（省略）

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

/GRNOD/NODE/1

#title

under_plate

# item_ID1 item_ID2 item_ID3 item_ID4 item_ID5 item_ID6 item_ID7 item_ID8 item_ID9 item_ID10

2 1 3 4 9 8 15 10

5 47 48 49 50 51 6 53

（省略）

接触条件の設定

100×100の正方形と5×5の正方形の熱源(800K)との接触を定義します。

接触条件の設定は解析で理解が難しいものの1つなので、自分も理解が乏しいです。
接触の定義は大別するとラグランジュ法とペナルティ法があります。

今回はペナルティ法を使うため、OpenRadiossで設定できるペナルティ法をマニュアルから引用します。

ペナルティ法を伴うインターフェースタイプ※Altairマニュアル引用

熱伝導などを扱える接触は/INTER/TYPE7なのでこちらを使います。

こちらはサーフェスと節点群との接触を定義する必要があります。

スレーブ：Slav_idにはサーフェス
マスター：Mast_idには節点グループ

それぞれPART IDを指定することで、グループは簡単に作っています。
基本的な設定は以下のようになっています。

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/INTER/TYPE7/2
New INTER 2
#  Slav_id   Mast_id      Istf      Ithe      Igap                Ibag      Idel     Icurv      Iadm
        21        12         0         1         0                   0         0         0         0
#          Fscalegap             GAP_MAX             Fpenmax
                   0                   0                   0
#              Stmin               Stmax          %mesh_size               dtmin  Irem_gap
                   0                   0                   0                   0         0
#              Stfac                Fric              Gapmin              Tstart               Tstop
                1E-7                 .15                1.01                   0                   0
#      IBC                        Inacti                VisS                VisF              Bumult
       000                             0                   0                   0                   0
#    Ifric    Ifiltr               Xfreq     Iform   sens_ID
         0         0                   0         1         0
#               Rthe                                    Tint Ithe_form
                .015                                     398         1
#               Frad                Drad              Fheats              Fheatm
                   0                   0                   0                   0
/GRNOD/PART/21
PID 1 NODE GROUP
         1
/SURF/PART/12
PID 2 SURF GROUP
# item_ID1  item_ID2  item_ID3  item_ID4  item_ID5  item_ID6  item_ID7  item_ID8  item_ID9 item_ID10
         2

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

/INTER/TYPE7/2

New INTER 2

# Slav_id Mast_id Istf Ithe Igap Ibag Idel Icurv Iadm

21 12 0 1 0 0 0 0 0

# Fscalegap GAP_MAX Fpenmax

0 0 0

# Stmin Stmax %mesh_size dtmin Irem_gap

0 0 0 0 0

# Stfac Fric Gapmin Tstart Tstop

1E-7 .15 1.01 0 0

# IBC Inacti VisS VisF Bumult

000 0 0 0 0

# Ifric Ifiltr Xfreq Iform sens_ID

0 0 0 1 0

# Rthe Tint Ithe_form

.015 398 1

# Frad Drad Fheats Fheatm

0 0 0 0

/GRNOD/PART/21

PID 1 NODE GROUP

/SURF/PART/12

PID 2 SURF GROUP

# item_ID1 item_ID2 item_ID3 item_ID4 item_ID5 item_ID6 item_ID7 item_ID8 item_ID9 item_ID10

接触剛性$K$、温度依存する摩擦係数$\mu$、熱摩擦、熱交換、放射など、非常に多くのパターンの設定が可能ですので、詳しくはマニュアルを見ることをお勧めします。

熱伝達の設定

オプションとして外気との熱伝達（自然対流、強制対流）の考慮も行うことができます。
それには「/CONVEC」キーワードを使います。

SURF_IDでサーフェスを指定。
「/SURF/PART/」でPARTを指定。
FUNCT_IDで参照温度の時間変化の関数IDを指定
「/FUNCT」で関数テーブルを作成
Hで熱伝達率を設定

/CONVEC/1
CONVEC with water
#  SURF_ID  FUNCT_ID SENSOR_ID
        11         6         0
#             ASCALE              FSCALE              TSTART               TSTOP                   H
                   0                   0                   0                   0                0E-4
/SURF/PART/11
New SURF
# item_ID1  item_ID2  item_ID3  item_ID4  item_ID5  item_ID6  item_ID7  item_ID8  item_ID9 item_ID10
         1
/FUNCT/6
CONVEC, temp vs time
#                  X                   Y
                   0                 298                                                            
                  60                 298

/CONVEC/1

CONVEC with water

# SURF_ID FUNCT_ID SENSOR_ID

11 6 0

# ASCALE FSCALE TSTART TSTOP H

0 0 0 0 0E-4

/SURF/PART/11

New SURF

# item_ID1 item_ID2 item_ID3 item_ID4 item_ID5 item_ID6 item_ID7 item_ID8 item_ID9 item_ID10

/FUNCT/6

CONVEC, temp vs time

# X Y

0 298

60 298

熱伝達は以下のように計算されます。

\begin{align*}
\dot{Q}=H\big(T-T_{inf}(T)\big)
\end{align*}

設定は記述しましたが、熱伝達率が大きすぎるとあまり熱伝達の様子が見れなかったので今回はH=0としています。つまり、熱伝達は起こっていません。

しかし、現実には外気との熱伝達があるため設定するほうが良いです。

節点の履歴

必要に応じて出力するようにしておくと便利です。
時刻歴応答「/TH」からはじまるキーワードを使います。

/TH/ACCEL：加速度
/TH/INTER：接触力など
/TH/NODE：節点の変位、速度、加速度
/TH/MONVOL：エアバッグの体積
/TH/PART：重心位置、慣性モーメントなど
/TH/RBODY：剛体の荷重
/TH/RWALL：剛体壁の接触力
/TH/SH3N：3節点のシェル要素の応力、ひずみ
/TH/SHEL：4節点のシェル要素の応力、ひずみ

以下は節点座標の時刻歴を出力する設定です。

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/TH/NODE/1
#thgroup_name
TH_Measurting_nodes                                                                                 
#  var_ID1   var_ID2   var_ID3   var_ID4   var_ID5   var_ID6   var_ID7   var_ID8   var_ID9  var_ID10
DEF       DEF       
#  node_ID   skew_ID                                         node_name
       103         0                                                                                
       620         0

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

/TH/NODE/1

#thgroup_name

TH_Measurting_nodes

# var_ID1 var_ID2 var_ID3 var_ID4 var_ID5 var_ID6 var_ID7 var_ID8 var_ID9 var_ID10

DEF DEF

# node_ID skew_ID node_name

103 0

620 0

以下は剛体の荷重の出力設定です。

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/TH/RBODY/2
#thgroup_name
T_RBODY                                                                                             
#  var_ID1   var_ID2   var_ID3   var_ID4   var_ID5   var_ID6   var_ID7   var_ID8   var_ID9  var_ID10
DEF       
#  Obj_ID1   Obj_ID2   Obj_ID3   Obj_ID4   Obj_ID5   Obj_ID6   Obj_ID7   Obj_ID8   Obj_ID9  Obj_ID10
         1

#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|

/TH/RBODY/2

#thgroup_name

T_RBODY

# var_ID1 var_ID2 var_ID3 var_ID4 var_ID5 var_ID6 var_ID7 var_ID8 var_ID9 var_ID10

DEF

# Obj_ID1 Obj_ID2 Obj_ID3 Obj_ID4 Obj_ID5 Obj_ID6 Obj_ID7 Obj_ID8 Obj_ID9 Obj_ID10

必要に応じて設定を追加すれば良いでしょう。

以上までの設定を001thermal_0000.radに保存します。
これがRadiossのstarter inputファイルになります。

Engine設定

以上まではRadiossのstarter inputファイルの設定でしたが、Radioss Engineファイルはシミュレーションの解析制御と出力を記述するファイルです。

ここでは最低限の設定だけを記述します。

/RUN：計算実行に関するキーワード。60msまでの解析
/DT/NODA/CST：モデルの時間ステップを増やすためのマススケーリング
/TFILE：時刻歴ファイルであるT-ファイルの書き出し時間間隔
/ANIM/DT：アニメーションファイル(A-ファイル)を出力する開始時間と出力間隔(0.5ms）

001thermal_0001.rad

/ANIM/DT
#   TSTART     TFREQ
0.000000 3.000000
/PRINT/-100
# Proc  nMIMD   minMIMD
/PROC/4/0
/RUN/001thermal/1
             60.0001
/TFILE/4
#            dT_HIS
0.100000
/TITLE/3
                                                                                                    
/VERS/140
/ANIM/VECT/CONT
/ANIM/NODA/TEMP
/ANIM/SHEL/TEMP
#/ANIM/GZIP
/MON/ON
/PARITH/ON
/DT/NODA/CST
#0.1 1e-009
0.67 5e-4

/ANIM/DT

# TSTART TFREQ

0.000000 3.000000

/PRINT/-100

# Proc nMIMD minMIMD

/PROC/4/0

/RUN/001thermal/1

60.0001

/TFILE/4

# dT_HIS

0.100000

/TITLE/3

/VERS/140

/ANIM/VECT/CONT

/ANIM/NODA/TEMP

/ANIM/SHEL/TEMP

#/ANIM/GZIP

/MON/ON

/PARITH/ON

/DT/NODA/CST

#0.1 1e-009

0.67 5e-4

さらに詳しく知りたい方はマニュアルを参考にしましょう。
今回は歪みや応力もアニメーションにしたかったので以下を追加しました。

/ANIM/VECT/CONT：接触力
/ANIM/NODA/TEMP：節点の温度
/ANIM/SHEL/TEMP：シェル要素の温度

解析実行

以上により解析の設定が終わりました。

解析に必要なファイルは以下の2つです。

001thermal_0000.rad：Starter ファイル（解析設定ファイル）
001thermal_0001.rad：Engine ファイル（解析制御と出力設定ファイル）

OpenRadiossの計算実行は、こちらで公開されているPythonで作られたGUIを使うと簡単に行えます。

WSLを使っているのでLinux環境での実行方法を参考にします。

これで計算投入用のGUIが立ち上がるので、アニメーションとデータ出力の項目にチェック（②)を入れて計算を実行します。
ここにチェックを入れておくことでA-ファイル（アニメーション用ファイル）はParaViewでも閲覧できるvtk形式に変換し、T-ファイル（データ用ファイル）はcsvファイルに変換してくれます。

OpenRadiossのインストール先は下記のスクリプトで自身の環境に合わせる必要があります。

openradioss_run_script_ps.sh

#!/bin/bash
set -e

pushd "$(dirname "$0")" >/dev/null
REL_PATH="../"
pushd "$REL_PATH" >/dev/null
ABS_PATH=/mnt/d/OpenRadioss/OpenRadioss #$PWD
popd >/dev/null

export OPENRADIOSS_PATH=$ABS_PATH
export RAD_CFG_PATH=$OPENRADIOSS_PATH/hm_cfg_files
export RAD_H3D_PATH=$OPENRADIOSS_PATH/extlib/h3d/lib/linux64
export LD_LIBRARY_PATH=$OPENRADIOSS_PATH/extlib/h3d/lib/linux64:$OPENRADIOSS_PATH/extlib/hm_reader/linux64
export OMP_STACKSIZE=400m
export PATH=$OPENRADIOSS_PATH/extlib/hm_reader/linux64:$PATH
export PATH=$OPENRADIOSS_PATH/extlib/intelOneAPI_runtime/linux64:$PATH
（省略）