OpenRadiossでくりまんじゅうのブロー成形にチャレンジしてみた
OpenRadiossのお勉強として公式の例題を参考にくりまんじゅう（ちいかわ©nagano）のブロー成形にトライしてみました。
#openradioss #ちいかわ #python #ブロー成形 #勉強記録
概要
OpenRadiossで遊んだ内容をまとめたチラシの裏です。
この前のお盆休みからOpenRadiossを触り始めたのですが（下記記事参照）これが意外と面白くてハマってしまいました＾ｐ＾；
なんか面白そうなのないかなー・・・と公式の例題集を眺めているとブロー成形でボトルを作る例題があり、
ｺﾚﾀﾞ！と思い、またしても以前FreeCADの練習で作ったくりまんじゅうを使って遊んでみました。
うまく表示されない場合は上記のAltair Community DocumentationでProduct=Radioss, Version=2024, GuideType=Tutrialにチェックを入れ、
Radioss 2024 DemoModelsを選択しダウンロードしてください。
radioss/exampleの中にあるRD-E-4400_Blow_molding_AMS.zipが今回参考にした例題です。
本記事では目次の通り、モデリングから計算および可視化までを一通りやってみた中での気づきや調べた内容のメモをまとめたものになります。
参考までに動作環境は下記の通りです。
Win11 + WSL2(Ubuntu 22.04 LTS)
LS-PrePost(R) 2024/R2(4.11.8) - 25Jul2024
FreeCAD 0.21.2
Paraview 5.12.1
※注意※
私はCAEもブロー成形もど素人なため、計算物理学や成形条件等についての知識ほぼゼロです lili orz lili
V&V（検証と妥当性確認）や精度は二の次でとりあえず動いたよー的な記事になります。ご容赦ください。。。
FreeCADでのモデル作成
今回もFreeCAD練習で作ったくりまんじゅうを使います。
（これを公開するのはさすがにまずそう？なのでご自分で試される際は球と円柱など単純な形状に置き換えて作成下さい）
例題と同様に、金型はZ軸方向の上下にそれぞれ配置し、パリソンと呼ばれる円筒形の樹脂は飲み口側と底側に大きなバリが残るのを想定したサイズとします。
これらはすべて2次元シェル要素としてメッシュ化されるものとします。
FreeCADで作成済のくりまんじゅうモデルファイルを開きます。
まずZ方向に直立しているくりまんじゅうをY軸方向に寝かせます。
このファイルのくりまんじゅうの身長は90mmなので例題集の寸法と大まかにそろえるために拡大します。
FreeCADのDraftワークベンチの「尺度」からサイズを変更したクローンを作成します。膨らむ方向の寸法をx,z方向でできるだけ揃えたかったので倍率を調整しました。このためちょっと縦長なくりまんじゅうになりました。。。
飲み口部分の円柱を作り頭頂部にめり込むように配置して拡大されたクローンと結合しておきます。
さらに直方体を二つ下図のようにくりまんじゅうの体の中心面（この場合はZ=0面）を対象に2つ合わさるように作成しそれぞれで差集合をとります。
出来上がったソリッドをダウングレードして面にして不要な面を削除します。
さらに位置を調整すればモデルは完成です。
最後に、上型、下型、パリソンをそれぞれigesファイルとして出力しておきます。
またここで、上型と下型の重心座標を調べてメモしておきます。
FreeCADのFCInfoマクロのCentre of Gravity X, Y, Zが重心座標です。
（インストールしていない方はツール -> Addon ManagerからFCInfoで検索しインストールしてください）
終わってから気づきましたがしっぽの部分が邪魔で型から抜けないですよねこれ。（結果しっぽの部分はうまく成形できなかったのでﾖｼｯ＾ｐ＾；）
メッシュ作成
今回もLS-PrePostを使用してメッシュを作成しました。
LS-PrePostを起動して先ほどFreeCADから出力したigesファイル3つをそれぞれインポートしたらMesh -> AutoMesherで2次元メッシュを作成します。
メッシュ化する部品（上型・下方・パリソン）の１つをAreaでドラッグして選択しMeshTypeをMixed, ElemSizeを0.6[mm]にしてMesh -> Accept -> Doneと順にクリックするとメッシュが形成されます。これを部品3つについて繰り返し行います。
くりまんじゅうがハーッ・・・ってしてる時の顔の部分が結構細かくて、メッシュサイズを小さくせざるを得ませんでした。
ちゃんと顔の部分もしっかりメッシュが切れているようです。
特に歪なメッシュも見つからなかったのでこのまま次の節点グループの取得に進みます。
節点グループの取得
例題にならい、パリソンの拘束点を設定するためにEntityCreationのSET_NODEから節点グループを設定します。
例題の拘束点を調べるためにradioss例題の節点・シェル情報をどうやってLS-PrePostに読ませればいいんだ・・・？
で詰まって途方に暮れていましたが、
もうこれ自分で変換したほうが早いんじゃね？
と思い、変換プログラムをPythonでサクッと作成して変換しました。
（/PART, /NODE, /SH3N, /SHELL, /GRNOD/NODEをそれぞれ変換するだけ）
先回作ったものをベースにちょこっと変えるだけで2時間くらいでできた＼(＾o＾)／
こんなことしなくてももっと楽な方法があると思います。
需要は無いと思うので公開はしませんが、この辺のプログラムはちゃんとまとめて整理しておきたいところ・・・
調べた結果、下図のようにパリソン端の円周エッジ部分（エッジ中の1点とそれ以外の点で別グループ）と
円柱の母線部分が節点グループとして設定されていたのでそれにならい節点グループを設定しておきます。
後程、成形中のパリソンがむやみに動かないよう拘束するためにこの情報を用います。
最後に.kファイルを出力してメッシュ化と節点グループの設定は終わりです。
Radioss形式への変換
先回作ったプログラムで変換。
LS-DYNAのキーワードがそのまま使えるっぽいので現状は*SET_NODE_LIST以外はおそらく変換不要のはず？
今回はとりあえずすべてradioss形式に変換しました。
LS-DYNAキーワードの自動変換を一度試したことがあるのですがマッピングに失敗？ or 自分の入力ミス？で
starter実行時にエラーや警告が出てしまったことがあり、それ以来可能な限りradioss形式で入力するようにしています。
この辺は私も知識不足なので詳しい方教えてくださいm(__)m
スターターファイルの作成
ようやくstarterファイルの作成にとりかかります。
冒頭で紹介した公式例題集のstarterをベースにしています。
Control Cards
解析全体の設定など。
starter実行時に警告がでたのでバージョンは140 -> 2022に変えて一部の不要なパラメータは削除しておきました。
#RADIOSS STARTER
# 参考: 公式例題集 - 44_blow_molding_amss
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/BEGIN
kurimanju_blow_molding
2022 0
Mg mm s
Mg mm s
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 1. CONTROL CARDS:
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/TITLE
kurimanju_blow_molding
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/ANALY
# N2D3D: 解析タイプを0=3次元, 1=軸対象, 2=平面ひずみに設定
# IPARITH: 並列オプション0=1=並列ON, 2=並列OFF
# ISUB: シェル要素のサブサイクリングフラグ0=無し, 2=有り
# N2D3D IPARITH ISUB
0 0 0
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/DEF_SOLID
# 全てのソリッドのPROPと厚肉シェル（本事例は厚肉シェルのみ使用）のデフォルト値を設定
# I_SOLID ソリッド要素定式化フラグ: 0->1 直行変形モードと剛体変形モードの補正を伴う粘性アワグラス定式化(Belytschko)
# ISMSTR 微小ひずみ定式化フラグ: 0->4 完全幾何非線形性(/DT/BRICK/CSTによる影響はありません)
# IFRAME 要素座標系定式化フラグ: 0->1 非共回転定式化
# I_SOLID ISMSTR IFRAME
0 0 0
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/IOFLAG
# 入出力設定
# 低減出力・STYファイル出力なし・回転自由度の計算を常に行わない
# IPRI IOUTP IOUTYY IROOTYY IDROT
0 0 0 0 0
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/SPMD
# https://2021.help.altair.com/2021/hwsolvers/ja_jp/rad/topics/solvers/rad/faq_rad_parallelism_r.htm
# DOMDEC Nproc Dkword Nthread
0 8 0 1 DAMPING
全節点に対して減衰係数を設定。
この辺の知識がないのでレイリー減衰について少し調べてみましたが、まだしっかり噛み砕いて理解できてませんorz
川重テクノロジー株式会社 技術レポート レーリー減衰について
α 質量に比例する減衰で低周波のゆっくりとした動きに対して大きくなる減衰に掛かる係数
β 剛性に比例する減衰で高周波の早い動きに対して大きくなる減衰に掛かる係数
αとβの値は例題の値をそのまま用いました。
β＝0なのでα減衰のみを考慮しているようです。
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 2. DAMPING:
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/DAMP/1
# 減衰係数の設定（レイリー減衰）
GLOBAL
# alpha beta grnod_ID skew_ID Tstart Tstop
.05 0 100000 0 0 0
/GRNOD/GENE/100000
# 1 ~ 412617のNODEをgrnod_ID=100000として登録
# ⇔全節点に上記の減衰設定を適用させる
ALL
1 412617MATERIALS
金型とパリソンの材料物性の設定。
パリソンの設定がちょっとだけ複雑です。
ひずみ速度による応力-ひずみ線図の変化を考慮しているようです。
FUNCTONの部分に記述されているものをまとめてプロットすると下記のようになっています。
ひずみ速度が大きくなればなるほど、応力（強度）が増しています。
NODES
node.incに重心情報を付け足したファイルをインクルード。
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 4. NODES:
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#include node_add_center_of_gravity.incインクルードファイルの中身です。
末尾に412616, 412617NODEを重心座標として追記してあります。
/NODE
1 -100.0 -200.0 27.0
2 -100.0 20.0 27.0
3 -100.0 -199.4005 27.0
4 -100.0 -198.8011 27.0
5 -100.0 -198.2016 27.0
...(省略)
412613 -25.0 -100.0 0.0134848
412614 -25.0 -100.5988 0.013443
412615 -25.0 -101.1976 0.0132821
412616 -0.0182 -92.7192 41.8850
412617 -0.0411 -93.0385 -42.2945BOUNDARY CONDITIONS
成形時の各種拘束条件の設定。
setnode.incから必要に応じコピペして作成しました。
金型はZ方向以外に動かないように固定します。
パリソンの両端をY方向に、母線をZ方向に動かないように拘束してあります。おそらくパリソンが型の中に吸い込まれないようにするためかと思います。
この辺の条件の詳細な意味は実際のブロー成形に携わったことがないのでわかりませんがひとまず例題の設定にならっておきます。
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 5. BOUNDARY CONDITIONS:
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/BCS/1
Tooling_12_const
# Tra rot skew_ID grnod_ID
110 111 0 100012
/GRNOD/NODE/100012
BCS_group_100012_of_NODE
412616 412617
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/BCS/11
PLASTIC_POINT
# Tra rot skew_ID grnod_ID
011 110 0 11
/GRNOD/NODE/11
PLASTIC_POINT
324511
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/BCS/12
PLASTIC_ENDS
# Tra rot skew_ID grnod_ID
010 000 0 12
/GRNOD/NODE/12
PLASTIC_ENDS
325702 325701 325700 325699 325698 325697 325696 325695 325694 325693
325692 325691 325690 325689 325688 325687 325686 325685 325684 325683
...(省略)
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/BCS/13
PLASTIC_LINE
# Tra rot skew_ID grnod_ID
001 110 0 13
/GRNOD/NODE/13
PLASTIC_LINE
324513 325440 324515 324516 325437 325436 324519 324520 324521 325432
325431 324524 325429 325428 325427 324528 325425 325424 325423 325422
...(省略) PARTS, GEOMETRICAL SETS
要素（シェル）情報およびPART情報の設定。
パリソンの厚みは2.0[mm]に設定します。あとは例題通りです。
例題ではパリソンのシェルごとに厚み情報が設定されていましたが今回は行いませんでした。
LS-PrePostのSET_ELEMENTで厚みを変えたい部分を選択して.kファイルに出力、それを変換してシェル情報と結合させれば簡単に部位ごとの厚みを変えることはできそうです。（めんどくさかったからやらなかった＼(＾o＾)／）
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 6. PARTS:
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#include shell.inc
/PART/1
TOOL1
1 1 0
/PART/2
TOOL2
2 2 0
/PART/3
PLASTIC
3 3 0#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 7. GEOMETRICAL SETS:
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/PROP/SHELL/1
Tool1
# Ishell Ismstr Ish3n Idrill
1 0 1 0
# hm hf hr dm dn
0 0 0 0 0
# N Istrain Thick Ashear Ithick Iplas
1 0 1 0 0 0
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/PROP/SHELL/2
Tool2
# Ishell Ismstr Ish3n Idrill
1 0 1 0
# hm hf hr dm dn
0 0 0 0 0
# N Istrain Thick Ashear Ithick Iplas
1 0 1 0 0 0
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/PROP/SHELL/3
Plastic
# Ishell Ismstr Ish3n Idrill
24 0 0 0
# hm hf hr dm dn
0 0 0 0 0
# Nはシェル厚み内の積分点の数で1~9を設定可能
# Ithick=1は厚みの変化を考慮するフラグ
# N Istrain Thick Ashear Ithick Iplas
5 0 2.0 0 1 1FUNCTIONS, PRESSURE LOADS
型の移動と圧力印加の設定。
（すでに上で説明したひずみ速度ごとの応力-ひずみ線図を示した関数は除く）
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 8. FUNCTIONS:
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/FUNCT/101
Tool Move
# X Y
0 0
0.40 1
0.60 1
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/FUNCT/111
Pressure on Plastic
# X Y
0 0
0.40 .020
0.60 .080
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 9. PRESSURE LOADS:
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/PLOAD/5
Pressure
# surfIDへの圧力設定
# 圧力の向きはSURF/SEGの n1n3 × n2n4 の方向であることに注意
# Ascale_x: 時間のスケール
# Fscale_y: 圧力のスケール
# 圧力(t) = Fscale_y * func * (t / Ascale_x)
# surf_ID functIDT sensor_ID Ascale_x Fscale_y
100013 111 0 1 -1.0
# /SURF/SEG/100013
#include surf_seg_parison_reindex.inc圧力の向きは n1n3 X n2n4 の向きであることに要注意です。
今回はなぜか内向きだったので圧力の向きにマイナスをかけて外向きにしました。(面の向きの設定方法がわからなかった・・・)
/SURF/SEGは例題では要素IDが1からの連番で再定義されていたのでそれにならいました。
/SURF/SEG/100013
Plastic
1 324511 324513 325703 325702
2 324513 325440 325704 325703
3 325440 324515 325705 325704
...(省略)金型の合わせ面のZ座標（シェル厚み考慮なしの中心値）とパリソンに印加する圧力の推移は下記のとおりです。
金型の動きは例題にならい、最も閉じた時点でやや隙間ができるようにしました。 また、圧力は計算時間短縮のため例題よりも短い時間で膨らむように調整しました。（結果、調整にちょっと失敗した＼(＾o＾)／）
INTERFACE
接触設定(/INTER/TYPE7)。
OpenRadioss初心者の私にはここが一番の鬼門でした。
今回のこのメモを通して少しは理解度が深まった気がします（´・ω・｀）
先に設定内容だけを記しておきます。
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 10. INTERFACES:
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
# 金型(Mast_id:メインサーフェスid)->パリソン(Slav_id:セカンダリ節点grnod)の接触設定
# パリソン->金型の接触設定
# Istf=4: インターフェース剛性はメイン及びセカンダリ剛性の最小値です
# Igap=2: 計算するギャップの可変ギャップ＋ギャップスケール補正
# Gapmin=0.01[mm]: 最小ギャップ
# Inacti=5: ギャップは時間とともに変化し初期ギャップは次のように調整されます
# gap0 = Gap - 初期貫通P0
/INTER/TYPE7/1
PLASTIC_INTER
# Slav_id Mast_id Istf Ithe Igap Ibag Idel Icurv Iadm
100006 100008 4 0 3 0 0 0 0
# Fscalegap GAP_MAX Fpenmax
.95 0 0
# Stmin Stmax %mesh_size dtmin Irem_gap
0 0 0 0 0
# Stfac Fric Gapmin Tstart Tstop
.1 .7 .01 0 0
# IBC Inacti VisS VisF Bumult
000 5 0 0 0
# Ifric Ifiltr Xfreq Iform sens_ID
0 0 0 2 0
/GRNOD/PART/100006
INTER_group_100006_of_PART
3
/SURF/PART/100008
INTER_group_100008_of_PART
1 2
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
# パリソン->金型の接触設定
/INTER/TYPE7/2
PLASTIC_INTER_sym
# Slav_id Mast_id Istf Ithe Igap Ibag Idel Icurv Iadm
100009 100011 4 0 3 0 0 0 0
# Fscalegap GAP_MAX Fpenmax
.95 0 0
# Stmin Stmax %mesh_size dtmin Irem_gap
0 0 0 0 0
# Stfac Fric Gapmin Tstart Tstop
.1 .7 .01 0 0
# IBC Inacti VisS VisF Bumult
000 5 0 0 0
# Ifric Ifiltr Xfreq Iform sens_ID
0 0 0 2 0
/GRNOD/PART/100009
INTER_group_100009_of_PART
1 2
/SURF/PART/100011
INTER_group_100011_of_PART
3
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
# パリソンの自己接触設定
# Igap=0なのでギャップはGapmin[mm]で固定
/INTER/TYPE7/3
PLASTIC_self
# Slav_id Mast_id Istf Ithe Igap Ibag Idel Icurv Iadm
100006 100011 4 0 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
10 .7 0.25 0 0
# IBC Inacti VisS VisF Bumult
000 1 0 0 0
# Ifric Ifiltr Xfreq Iform sens_ID
0 0 0 2 0ギャップ
接触の取り扱いには2種類あり、今回用いたTYPE7はペナルティ法と呼ばれるRadiossでは一般的な方法とのことです。
詳細は下記リンクのRadiossのドキュメントが非常にわかりやすいです。
Altair Radioss 2022 ユーザガイド 接触の取り扱い
各面にはギャップと呼ばれる接触判定領域があり、
その領域内にセカンダリ節点グループの節点が入ると面と節点が接触したと判定されます。
今回のシミュレーションでは金型とパリソンが両方動くので、金型面⇒パリソン節点の接触とパリソン面⇒金型節点の接触の2種類の設定が必要なようです。加えて、パリソン同士が自己接触するのでその設定も必要となり計3つの接触設定が起きています。
TYPE7のギャップ領域は下記図のような領域になっています。（下記リンクより引用）
Altair Radioss 2022 ユーザガイド 汎用インターフェース(/INTER/TYPE7)
Inacti, Igapの値を選ぶことで、時間ごとに可変ギャップを得ることも可能です。ここを理解するためには初期貫通について知っておく必要があります。
初期貫通
初期貫通と可変ギャップの関係については、
ユーザガイドの初期貫通のページの説明が丁寧でわかりやすかったので下記にリンクを記しておきます。
Altair Radioss 2022 ユーザガイド 初期貫通
例題では金型⇔パリソンはIgap=2で初期貫通を引いた値を初期ギャップにする設定でした。
それにならいIgap=2で計算開始しようとしたところIgap=3にしろと警告が出たので素直に3にしました＼(＾o＾)／
Igap=3の時は2の場合に加えてさらにメッシュのサイズも考慮してギャップを計算してくれるようで、より柔軟に可変ギャップを得ることができるようです。
パリソンの自己接触はIgap=0でGapMin固定です。例題にならい少し大きめの値にしておきました。
RIGID, IMMDISP, TH
剛体設定と型の移動設定など。
金型の移動についてはすでに説明済のため、ここで特記することはないです。
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 11. RIGID BODIES:
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/RBODY/1
TOOL1
# RBID ISENS NSKEW ISPHER MASS Gnod_id IKREM ICOG Surf_id
412616 0 0 0 0 100002 0 0 0
# Jxx Jyy Jzz
0 0 0
# Jxy Jyz Jxz
0 0 0
#Ioptoff
0
/GRNOD/PART/100002
RBODY_group_100002_of_PART
1
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/RBODY/2
TOOL2
# RBID ISENS NSKEW ISPHER MASS Gnod_id IKREM ICOG Surf_id
412617 0 0 0 0 100003 0 0 0
# Jxx Jyy Jzz
0 0 0
# Jxy Jyz Jxz
0 0 0
#Ioptoff
0
/GRNOD/PART/100003
RBODY_group_100003_of_PART
2
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 12. IMPOSED DISPLACEMENTS:
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/IMPDISP/1
TOOL1_MOV
# Ifunct DIR Iskew Isensor Gnod_id Frame Icoor
101 Z 0 0 100004 0 0
# Scale_x Scale_y Tstart Tstop
0 -25.0 0 0
/GRNOD/NODE/100004
IMPDISP_group_100005_of_NODE
412616
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/IMPDISP/2
TOOL2_MOV
# Ifunct DIR Iskew Isensor Gnod_id Frame Icoor
101 Z 0 0 100005 0 0
# Scale_x Scale_y Tstart Tstop
0 25.0 0 0
/GRNOD/NODE/100005
IMPDISP_group_100004_of_NODE
412617
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#- 13. TIME HISTORIES:
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
/TH/INTER/1
TH INTER
# var1 var2 var3 var4 var5 var6 var7 var8 var9 var10
DEF
# Obj1 Obj2 Obj3 Obj4 Obj5 Obj6 Obj7 Obj8 Obj9 Obj10
1 2
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|
#/AMS
#
/END
#---1----|----2----|----3----|----4----|----5----|----6----|----7----|----8----|----9----|---10----|エラーチェック
こちらで公開されているPythonのGUIアプリを使わせていただきました。
OpenRadioss Confluence - python/tk guis for job submission, anim-vtk conversion and T01-csv conversion
Run Starter Onlyにチェックを入れて実行するとStarterファイルの出力（~_0000.out）のみを得ることができます。
ファイル中には種々の情報に加え、エラーと警告のメッセージが表示されます。確認したところエラーも警告もありませんでした。
これでスターターファイルの作成は終わりです。
エンジンファイルの作成
engineファイルを作成します。
ほぼ例題と同じですがAMSはよくわからなかったので今回は使用していません。
試しに使ってみましたが、計算途中でエラーが発生したり結果がおかしかったりで使いこなせませんでした。
0.0[s]~0.6[s]まで0.01[s]刻みで、
ミーゼス応力やシェル厚みを含んだアニメーションファイルを出力させます。
/TITLE
kurimanju_blow_molding
/VERS/140
/DT/INTER/DEL
6.700000e-001 1.000000e-007
/DT/NODA/CST
6.700000e-001 7.700000e-007
/ANIM/DT
0.000000e+000 1.000000e-002
/TFILE
1.000000e-003
/ABF
1.000000e-003 1e30
/RFILE
10000 0 0
/PRINT/-100/55
/RUN/kurimanju_blow_molding/1
0.6000000e-000
/ANIM/ELEM/EPSP
/ANIM/ELEM/VONM
/ANIM/VECT/VEL
/ANIM/SHELL/THIC
/MON/ONおそらくもっとも重要なのは、/DT/NODA/CSTの設定だと思います。
公式ユーザガイドおよびstarterファイルの出力の最後のほうに下記のような図と表があるのでこれを参考に設定しておけばよいようです。
Altair Radioss ユーザガイド 節点時間ステップコントロール
スケールファクターは0.67、初期DM/Mは0.03が推奨値とのことです。
今回のstarterファイル出力ではそれに対応する値は7.70e-7[s]でした。
/DT/NODA/DELは時間ステップが設定値以下になるような節点を除去する設定です。
DT/NODA/CSTの設定値の1/100位を目安にするのが良いとユーザガイドに書いてありましたが気づかずに1.0e-7にして計算実行してしまいました。
（・・・計算した結果、除去節点は発生しておらず結果オーライでした。次から気を付けます＾ｐ＾；）
計算実行
ようやく計算開始！
上記エラーチェックの項で使用したGUIアプリを使用させていただきました。VTKファイルへの変換をONにするチェックを入れておくと後々楽です。
8コア並列で約5時間かかりました・・・
事前に時間ステップを大きな値にして動作の確認を可能な限り行いましたが、当然その場合は計算がエラー停止することが多く、条件のあたり付けに苦労しました。。。
結果
計算結果のチェックをして、paraviewで可視化します。
計算結果のチェック
Altair Radioss ユーザガイド 計算チェック
を参考に、結果のチェックを行います。
ENERGY_ERRORの値はマイナスかつ変化量が15%より小さい必要があるとのことで、今回の結果は問題なさそうです。
TIME_STEPに関しても極端に小さくなっている箇所はなく、これも同じく問題なさそうです。
また、MASS.ERR, MASS_ADDEDは全ステップを通して0でした。
計算結果には特に大きな問題はなさそうで、正しく計算されたのではないかと思います。
CYCLE,TIME,TIME-STEP,ELEMENT,ID,ERROR,I-ENERGY,K-ENERGY_T,K-ENERGY_R,EXT-WORK,MAS.ERR,TOTALMASS,MASSADDED
100,0.4793E-03,0.4793E-05,NODE,412063,-0.0%,0.4519E-04,1.755,0.2651E-16,1.755,0.000,0.9614E-03,0.000
200,0.9585E-03,0.4793E-05,NODE,412063,-0.0%,0.1627E-02,1.756,0.2215E-14,1.758,0.000,0.9614E-03,0.000
300,0.1438E-02,0.4793E-05,NODE,412063,-0.0%,0.7101E-02,1.756,0.2910E-13,1.763,0.000,0.9614E-03,0.000
...
188900,0.5994,0.2575E-05,NODE,357613,-1.4%,0.1670E+05,0.1678E-01,0.8538,0.1693E+05,0.000,0.9614E-03,0.000
189000,0.5996,0.2573E-05,NODE,357613,-1.4%,0.1670E+05,0.1539E-01,0.8696,0.1693E+05,0.000,0.9614E-03,0.000
189100,0.5999,0.2574E-05,NODE,357613,-1.4%,0.1670E+05,0.1642E-01,0.8862,0.1693E+05,0.000,0.9614E-03,0.000また、内部エネルギーの推移は下記の通りでした。（例題集で同様にプロットしてたので載せてみました）
paraviewでの可視化
金型を横から見たときと下から見上げた時の動画です。
金型のみ50%透過、色はシェルの厚みです。
ちょっと型を閉じる前に印加する圧力が高かったようで、型が閉じきる前（0.38[s]付近）にブワっと膨らんでしまっています。もうちょっとこの辺は調整が必要かも？
型が閉じてからは強めの圧力が印加され、くりまんじゅうの顔部分に沿ってパリソンがじんわり変形しているのが確認できます。
とりあえず例題の結果からは大きく外れることなくシミュレーションを行えたようで良かったです。
下記はt=0.60[s]での成形後のくりまんじゅうです。
一番気にしていた顔の部分はなんとかくりまんじゅうがハーッしてることが認識できるレベルで成形できました！
でも耳・手・しっぽは膨らみが足りませんでした・・・
（しっぽに関しては前述した通りこの形状のままだと型から抜けなくなっちゃうので要修正）
この辺りはもう少し条件を振ったらどうなるのか？試してみたかったのですが今回メッシュを無駄に細かく切ってしまったために計算時間が延びてしまい色々試すことができませんでした。
上手にメッシュ作成できるようになりたいです。
感想
ほぼ例題集の内容をなぞっただけですが、自分でモデル・メッシュの作成からトライしてみると色々な気づきがありました。
特に接触設定は少しずつ解るようになってきました。
ブロー成形についても興味がわきました。
実際は今回考慮しなかった温度によっても大きく出来栄えが変わるみたいです。また、3次元ブロー成形のような複雑な形状をブロー成形する手法もあるようで面白そうです。
とにかく楽しかったです。
今後もなにか面白そうな題材を見つけて遊んでいきたいです。
参考資料
この記事が参加している募集
この記事が気に入ったらサポートをしてみませんか？