ポスト処理

2013年1月 5日 (土)

ParaVISで応力やひずみを部分的に表示する

SalomeのParaVIS(ParaVIEW)で計算結果の応力やひずみ等を表示する際、モデル全体ではなく、ある一部分だけ表示し、残りは単一の色または非表示にしたい場合があります。

これについては、メッシュ作成時、またはCode_Asterの中でグループを作成しておけば、そのグループごとに表示をON/OFFして、表示方法の設定を変えることが出来ます。

下記は、Salome-Meca2012.2(SALOME 6.5 + Code_Astr 11.2.0)での例です。OSはdebian6です。

これは、先日実施した「ゴムバンドの摩擦力による動力伝達(前編後編問題点)」の結果表示例です。

Group01


これは「WarpByVector」フィルターを用いて変形状態でMises応力を表示したものですが、ローラー部分についても応力分布が表示されているため、見づらくなっています。このローラー部分は評価対象では無いので、応力を消したいところです。

そこでまず、応力を表示したいゴムバンドである「bandグループ」のみに絞った表示をします。

Pipeline Browser上で「WarpByVector1」を選択した状態にしておき、そこでメニューから「Filters」-「Alphabetical」で「Extract Group」をクリックして追加します。そして追加された「ExtractGroup1」のPropertiesに「Group」のタブがあるので、その中で表示したいグループのみ(ここでは「band」)にチェックを入れ、その他はチェックを外し、Applyを押すと、チェックしたグループのみが表示されます。

Group02


表示内容については、ツールバーで設定します。今回はMises応力分布を「Surface With Edges」で表示しています。

ここでもしPipeline Browser上で「WarpByVector1」の表示がONになっているようなら、OFFにしてください(左端の眼のマークをクリックしてグレーにする)。また、「ExtractGroup1」が表示OFFになっていたらONにしてください(左端の眼のマークをクリックして黒にする)。

次に、応力を表示したくないローラーである「drグループ」、「dnグループ」です。
再びPipeline Browser上で「WarpByVector1」を選択し、もう一度「Extract Group」を追加します。
ここで、必ず「WarpByVector1」を選択した上で「Extract Group」を追加する事が重要です。
上で作成した「ExtractGroup1」を選択したままで追加すると、「ExtractGroup1」から更に絞り込むという設定になってしまいます。

「ExtractGroup2」が出来るので、PropertiesのGroup設定で、drとdnのみにチェックが入った状態にしてApplyを押します。
これで、ローラー(dr, dn)が表示されますので、ツールバーで表示内容を「Solid Color」にします。これで単色で塗られます。

Group03


Solid Colorの色の設定は、「ExtractGroup2」の「Display」タブに「Set Solid Color」のボタンがあるので、それを押して変更できます。

Group04

なお、今回は「Extract Group」フィルターを用いましたが、同じような事を「Exract Block」や「Extract Selection」でも行えます。後者は、グループを作成していなくても、ParaVISの画面上で選択さえ出来れば、その要素(節点)のみを抽出できるので便利です。

Group05


今回は以上です。

2012年9月 2日 (日)

非線形解析の設定と、出力される応力・ひずみの種類について(弾塑性編)

前回、「非線形解析の設定と、出力される応力・ひずみの種類について」を線形弾性材料(ELAS)と超弾性材料(ELAS_HYPER)にて検証しましたが、今回は同じ検証を弾塑性材料で実施しました。

検証の趣旨や方法、理論式等については前回と同じなので、そちらをご参照ください。
ソフトウェアも同じく、DEXCS2011-Salome-A3(64)(SALOME6.3 + Code_Aster STA10.5)の利用です。

材料物性は下記の設定としました。
(「公称応力-工学ひずみ」での定義)

ヤング率:1000 [MPa]
ポアソン比: 0.499
初期降伏応力:1 [MPa] (ひずみ:0.001)
ひずみ硬化率:10 [MPa]

理論式が非圧縮を前提としているため、弾性域についても、ほぼ非圧縮(ポアソン比0.499)としています)。
上記の物性は「公称応力-工学ひずみ」による定義ですが、これを「コーシー応力(真応力)-対数ひずみ(真ひずみ)」に換算すると、下表の通りになります。

Nl0401

また、両曲線をグラフにすると下図のようになります。

Nl0402

これを前回同様に、検証用解析モデルの形状である、初期断面積=1[mm^2]、初期長さ1[mm]の条件を理論式に入れて、横軸を変位(0~1mm)にしてプロットしたものが下図です。

Nl0403

次に、Code_Asterによる計算の設定です。

塑性域における応力ーひずみの関係は、DEFI_FONCTIONにて定義します。
データの最初の点を降伏点として、「応力ー全ひずみ」の形式でデータを入れていくのですが、この「応力ー全ひずみ」の関係を、どの応力とひずみの組み合わせで入力するのか、選択肢が生じます。

基本的には共役な応力とひずみの関係である、「公称応力 - 工学ひずみ」、「コーシー応力 - 対数ひずみ」、「第2ピオラキルヒホッフ応力 - グリーンひずみ」の3通りが考えられますが、第2ピオラキルヒホッフ応力はあまり一般的でないと思いますので、今回は前者2つ(工学ひずみと対数ひずみ)だけで検証することにしました。

Eficasによる設定画面は下記です。

公称応力-工学ひずみ

Nl0404


コーシー応力-対数ひずみ

Nl0405


それから、非線形解析の設定(STAT_NON_LINE)ですが、今回は「COMP_ELAS」は試さず、「COMP_INCR」のみにしました。
この両者の違いについては、Code_Asterの公式サイトの「Support」-「Code_Aster trainings」に分かりやすい解説があります(下記リンク)。

Non linear capabilities

「COMP_ELAS」は弾性挙動(線形または非線形)となり、負荷過程と除荷過程で同じS-S線上をたどるので、永久ひずみが生じません(これについては、後日、別に検証記事を書きます)。

「RELATION」は「VMIS_ISOT_TRAC」です。Von-Misesの降伏条件による、等方硬化則の設定となります。

Nl0406

「DEFORMATION」については、下記の5種としました。

・PETIT
・PETIT_REAC
・GROT_GDEP
・SIMO_MIEHE
・GDEF_LOG

結局、入力データ2種類と、DEFORMATION5種類の10ケースを試すことになりました。

結果ですが、まず計算自体が出来たかできたかどうかについては下表の通りです。

Nl0407

“S-S”は入力データとして用いた「応力-ひずみ曲線」、“Eng”は「公称応力-工学ひずみ」、“Log”は「コーシー応力-対数ひずみ」の意味です。
GROT_GDEPにおいて、変位が0.7mmを超えると計算が止まってしまうことを除けば、全ての組み合わせで計算することが出来ました。

そして、これも前回同様、理論による応力・ひずみと、Code_Asterの計算結果を重ねてプロットしました。

ただし今回は、ひずみの比較は載せていません。
ひずみの出力は、どのような組み合わせにおいても

・EPSI_NOEU:工学ひずみ
・EPSG_NOEU:グリーンひずみ

となったためです。全て下図のようになりました。

Nl0408

したがって、結果は応力の比較のみです。図が大きいので前回同様にPDF化しています。
また、変位0~1mmのグラフでは、降伏点付近が見にくいので、2ページ目に変位0~0.05mmの拡大グラフも載せています。

nl04-1.pdf

●入力データに「公称応力-工学ひずみ」を用いた場合。

PETITを除き、変位(=工学ひずみ)が微小領域(およそ4~5%程度)を超えると、SIEF_NOEUはどの応力にも一致せず、正しい評価が出来ません。反力も同様です。
PETITでは大ひずみ領域でも公称応力が得られるのですが、変形がおかしくなっています。

Nl0409


これはX方向への変形1mm(X方向伸張比=2)の状態ですが、Y方向、Z方向の長さがそれぞれ半分になっています。これでは「非圧縮」の条件を満たしていません。
2*0.5*0.5=0.5で、元の体積=1から半分に減っています。

したがって、入力データに「公称応力-工学ひずみ」を用いた場合、工学ひずみ4~5%程度の微小変形領域までの利用に留めるべきと思われます。

●入力データに「コーシー応力-対数ひずみ」を用いた場合

GDEF_LOGでは、SIEF_NOEUがコーシー応力に一致し、反力も正確です。
PETIT_REACとSIMO_MIEHEでは、大ひずみ域で若干のズレがあるものの、SIEF_NOEUがコーシー応力にほぼ一致し、反力もほぼ正確です。

PETITは微小ひずみ領域のみコーシー応力に一致しています。
GROT_GDEPは微小ひずみ領域でも、ズレが大きくて使えません。GROT_GDEPを使うのであれば、「公称応力ー工学ひずみ」の入力のほうが良いです。

-------------------------------------------------------------

応力については、以上のような結果となりました。
しかし、応力以外にも、変形や回転が正確なのかどうかが気になるところです。
少し長くなりますが、今回はついでに変形についてもチェックしてみました。

上で少し触れましたが、PETITの場合は非圧縮の条件から外れてしまう結果となりました。
非圧縮の単軸変形の場合、引っ張り方向の伸張比をλ1、他の2方向の伸張比をそれぞれλ2、λ3とすると、下記の関係式が成り立ちます。

λ1×λ2×λ3 = 1

ここで、λ2=λ3 なので、

λ2 = λ3 = 1 / SQRT(λ1)

となります(SQRTは平方根)。

この理論値と、各計算結果を重ねてプロットしました。
横軸はX方向(引っ張り方向)の変位、縦軸はZ方向の伸張比です。


Nl0410


入力データに関わらず、PETITおよびGROT_GDEP以外は理論値と重なっています。つまり正確な変形をしています。

ここまでの結果をまとめると、下表の通りになります。


Nl0411_2



これを見ると、上のほうで紹介した公式サイトのドキュメント「Non linear capabilities」で説明されている内容とほぼ整合するようです(○△×の判定は、私の主観によるものです)。

ただし、今回は大回転については未検証です。
荷重を与えながら剛体回転させる、または変形によって応力を生じさせてから剛体回転させるという設定は、ちょっとややこしそうです。そこで公式サイトの例題を探ってみたのですが、

Documentation version 11 v6.08.106.pdf (SSND106)

が、それに該当するようです。commファイルの内容も複雑なので全ては理解していないのですが、SIMO_MIEHEとGDEF_LOG(およびGDEF_HYPO_ELAS)は、大回転も可能のようです。

一方、PETIT_REACは、上述の「Non linear capabilities」や、「R5.03.24 (4.2.2項)」等を見ると、大回転には向いていないようです。

結局、弾塑性解析で「大ひずみ」・「大変形」・「大回転」の全てに対応するには、下記のようにするのが良さそうです。

・永久ひずみを考慮する場合、「COMP_INCR」を用いる。
・DEFORMATIONは「SIMO_MIEHE」か「GDEF_LOG」。
・入力データは「コーシー応力 - 対数ひずみ」。
・計算結果の応力(SIEF_xxxx)は「コーシー応力」として評価する。

なお、今回のメッシュファイルとcommファイルは下記です。
commファイルは、入力データがEngとLogの2種類です(DEFORMATIONはGDEF_LOGにしています)。

nl04.zip

2012年8月25日 (土)

非線形解析の設定と、出力される応力・ひずみの種類について

大変形・大ひずみをともなう非線形解析を行う場合、使用するFEM解析ソフトの設定によって、出力される応力とひずみの種類が変わってきます。Code_Asterについても、自分が良く使う解析の設定で、どのような応力・ひずみが出力されているのか把握しておく必要があると思いました。

そこで今回、私が良く使う等方性線形材料(ELAS)と超弾性材料(ELAS_HYPER)について、非線形解析(STAT_NON_LINE)の設定によって出力がどのようになるのか検証しました。

当初は弾塑性材料も検証対象に入れるべきと思ったのですが、塑性変形については色々と設定することが多く、私自身もあまり塑性については詳しくないため、ややこしいので今回は見送りました。また機会があれば塑性についても試してみます。

ソフトウェアは、DEXCS2011-Salome-A3(64)(SALOME6.3 + Code_Aster STA10.5)の利用です。

まず検証用の解析モデルですが、下図のように、1辺1mmの正方形の六面体1次要素1個のみです(要素数1、節点数8)。1要素だけの挙動を検証することにより、メッシュの粗密などの影響を排除しています。

Nl0201

この要素の面のうち、X,Y,Z軸に垂直な面のそれぞれ1面ずつを面内拘束しました(fixX, fixY, fixZ)。そして、X軸に垂直な面の、拘束した側と反対の面(dispX)を、強制変位で1mm引っ張ることにしました。これで、X方向に単軸引っ張り(工学ひずみ 0 -> 1)の状態になります。

Nl0202


材料物性は、下記の通りです。

・線形材料:ヤング率=1.0 [MPa], ポアソン比=0.499

・超弾性材料:C10=0.4057582, C01=-0.254934, C20=0.0019713
          (Cxxの単位はMPa)
         ポアソン比=0.499 

線形材料は、工学ひずみ=1のときに反力=1[N], 公称応力=1[MPa]となるよう、ヤング率=1[MPa]に設定しました。ほぼ非圧縮を仮定し、ポアソン比は0.499です。

超弾性材料は、「ゴム(超弾性)の物性値を同定する方法」にて作成したエクセルによる同定ツールを用い、単軸引っ張りの工学ひずみ=0~2の範囲で、ヤング率=1[MPa]の線形材料と同じS-Sとなるよう、定数を調整しました。

Nl0203


超弾性材料でありながら、特定の範囲で線形挙動となる、変な材料になりました。
ポアソン比はほぼ非圧縮の0.499です。

次に、単軸変形と非圧縮を仮定して、理論式から各種のひずみと応力を求めました。
一般的に、FEM解析ソフトウェアで扱われる(出力される)のは、下記のひずみと応力になります。

■ひずみ
(元の長さ = L、変位 = u とする)

・工学ひずみ: ε= u / L
・グリーン・ラグランジュひずみ:  E = ε + 0.5*ε*ε
  (以下、“グリーンひずみ”と略します)
・対数ひずみ: e = ln (1+ε)

■応力
(反力 = f、 元の断面積=A とする)

・公称応力: P = f/A
         線形材料の場合、P = ヤング率 * ε
・第2ピオラ・キルヒホッフ応力: S = P / (1+ε)
・コーシー応力(真応力): T = P * (1+ε)

(これらは単軸変形と非圧縮の仮定下でのみ成り立ちますので、ご注意ください。)

理論式については、下記書籍を参考にしました。
「解析塾秘伝 非線形構造解析の学び方!」 著:石川覚志、日刊工業新聞社

各種の応力とひずみの説明についても、上記書籍が参考になります。
(初版本は、第2ピオラキルヒホッフ応力の記述に誤りがありました。正誤表を合わせて見てください)。

これらの式に、ヤング率=1[MPa]、初期断面積=1[mm^2]、初期長さ1[mm]の条件を入れて、横軸を変位(0~1mm)にしてプロットしたものが下図です。

Nl0204


今回の形状・条件では、変位=工学ひずみ=公称応力=反力 になります(全て0~1になります)。

さて、次はCode_Asterでの計算です。

細かい説明は省略しますが、上述の解析モデル・条件で、dispXへのX方向変位1mmを、100分割で与えました。出力は20分割にしています。
解析に用いたファイル(の一例)は下記にまとめて圧縮しています。

nl02.zip

nl02.mmed (メッシュファイル)
nl02_elas.comm (線形材料のcommファイルの一例)
nl02_hyper.comm (超弾性材料のcommファイルの一例)

線形材料(RELATION:ELAS)および超弾性材料(RELATION:ELAS_HYPER)の両材料について、「COMP_INCR or COMP_ELAS」、「DEFORMATION」を変えて計算を実行しました。その結果、まずは計算自体が出来たかどうかの結果が下表です。

Nl0205

○は最後まで計算できたもの、△は変位1mmに達する前に途中で停止してしまったもの、×は、「設定が適切でない」とのエラーメッセージが出て、計算が実行できなかったものです。
超弾性材料は「GROT_GDEP」のみでしか計算できないことが分かりました。
線形材料も、「PETIT」、「PETIT_REAC」、「GROT_GDEP」のみです。

計算できた7ケースについては、それらの各種出力(X方向成分)を、理論式によるひずみと応力のグラフに重ねてプロットしました。グラフ数が多いので、別途PDFファイル化しました。

nl02-1.pdf

●ひずみについて

7ケース全てについて、

・EPSI_NOEU:工学ひずみ
・EPSG_NOEU:グリーンひずみ

となりました。これは、公式ドキュメント「U2.01.05」に書かれている通りです。

●応力・反力について

REAC_NODAを集約したものが反力ですが、今回は断面積の初期値が1[mm^2]なので、これがそのまま公称応力になります(反力の集約方法は「節点反力を集約する方法」と同じ方法を用いました)。

(1)と(2)(線形材料でPETIT)については、反力は理論通りで、SIEF_NOEUは公称応力となりました。

(3)と(4)(線形材料でGROT_GDEP)については、計算が変位0.7mmを超えたところで停止しています。反力の値は理論からズレており、明らかにおかしな挙動となっています。SIEF_NOEUも、3つの応力理論値のいずれとも合致しませんでした。

しかし、微小変形領域(変形が数%程度)では、これら3種の応力の区別は実用上ほとんど無く、SIEF_NOEUの値もそれに近いものとなっています。線形材料で大回転を扱う場合、微小ひずみ領域に限ってはこのオプションで良さそうです。
(*「異なる「DEFORMATION」設定の混在」参照)

(5)についても、ズレ方は(3)・(4)と違いますが、同様のことが言えます(大回転については未確認です)。

(6)と(7)は超弾性材料のGROT_GDEPですが、反力は理論どおりで、SIEF_NOEUはコーシー応力(真応力)として出力されています。これも問題なさそうです。

なお今回の結果では、COMP_ELASとCOMP_INCRの違いがありませんでした。

以上の結果をまとめると下表の通りとなります。

Nl0206

結論として、

やはり線形材料(ELAS)では数%以上の大ひずみ変形は避けた方が良いようです。微小ひずみで大回転をする場合は「GROT_GDEP」で良さそうです。

超弾性材料(ELAS_HYPER)の場合は、GROT_GDEPのみです。応力はコーシー応力(真応力)になります。

今回は以上です。

2012年8月 7日 (火)

計算途中にエラーで止まってしまった非線形解析の結果を見る方法

先日実施した「ゴムと金属の接触・摩擦解析(平面ひずみ)」を題材にして、多数ステップの非線形解析で計算途中にエラーで止まってしまったJOBについて、計算が出来ているところまでの結果を見る方法について考えます。

ソフトウェアは、DEXCS2011-Salome-A3(64)(SALOME6.3 + Code_Aster STA10.5)の利用です。

ステップ数が多く時間がかかる非線形解析をする場合、最終インクリメント(最終ステップ)までもう少しという所で計算が止まってしまい、ガッカリすることが良くあります。そのような時は、なぜそこから先に計算が進まないのか、変形や応力の様子を見て検証したいものです。しかし、Code_Asterで普通に計算していると、「STAT_NON_LINE(非線形静解析)」でエラーになった場合、その時点でJOBが止まってしまうので、その下の結果出力の部分を実行してくれず、結果が出力されません。

計算が出来ている所までの結果を出すのに、一番真っ当な方法は、ASTKで「.base」フォルダを指定して、そこに計算結果のデータが保存されるようにしておいて、計算が完了するかエラーで止まったら、その「.base」フォルダからデータを読んで結果を出力するJOBを新たに実行することだと思います。
(むしろ、本来は計算JOBと結果出力JOBを分けるべきなのかな? と思います)。

その結果出力用のJOBのコマンドファイルは、「POURSUITE();」で始まり、結果出力する部分のみを記述します。

この方法については、「Open CAE Users Wiki」-「オープンCAE初心者勉強会」-「 SALOME-Mecaの活用技術:MK2氏提供」に、「ポスト処理 MECA_STATIQUEの場合の一例 資料:Post_Lin_ver10.pdf 」として情報公開して頂いております(これは線形解析の例ですが、非線形解析でも同じです)。大変分かりやすく説明されているので、私がここで解説するよりは、こちらの資料をお読み頂いたほうが良いかと思います。

-----------------------

さて、ここまでが前置きで、今回は同じ事を別の方法で行ってみました。

公開されている公式フォーラムをさっと検索するだけで、下記のような情報見つかりました。

http://www.code-aster.org/forum2/viewtopic.php?id=12992

http://www.code-aster.org/forum2/viewtopic.php?id=13650

http://www.code-aster.org/forum2/viewtopic.php?id=14104

同じようなことを考える人が結構多くいるようです。

Pythonを使うようで、必要に応じて色々な記述が出来るようですが、下記の通りにするのがもっともシンプルなようです(題材は上述の通り、先日実施した「ゴムと金属の接触・摩擦解析(平面ひずみ)」です)。

(1)DEBUTにPAR_LOT='NON'を設定する(Eficasで設定できます)。
(2)STAT_NON_LINEを、下記のような記述で挟む。Eficasでは編集できないので、Eficasからcommファイルの出力後に、テキストエディタ等で編集する。

Wedge0201_2

ここで重要なのは、「RESU=・・・」の前と、「pass」の前にそれぞれインデント(字下げ)が必要なことです。これはPythonの文法によるものと思います。Eficasから出力された時点では、「RESU=・・・」の前にはインデントは無いので、注意が必要です。

元の題材は最終ステップまで正常に計算できる設定になっていましたので、ここでゴムの押し込み量を元の「-1.6mm」から「-3.0mm」に変更して、わざと途中で計算が止まるようにしました(自動ステップ調整機能があっても、最も細かい時間幅まで分割が到達して、そこで止まります)。

Wedge0202

これで計算した結果の動画です。

最終収束ステップの画像です。

Wedge0203

要素がつぶれてしまっています。収束したとはいえ、このような状態になってしまっては、それ以上先のステップへは進めないことが分かります(進む意味もありません)。

このような方法で、エラーで止まっても、途中までの結果を出力することが出来ました。

この方法の難点は、「try: ・・・ pass」の記述をEficasで編集できないことです。Eficasでcommファイルを出力してから、テキストエディタで追加の編集をする必要があります。

更に面倒なことに、この記述を加えた後のcommファイルは、Eficasで読めません。「RESU=・・・」の前にインデントを入れたままEficasで読もうとしても、「余分なインデントがある」というエラーになって読めません。このインデントを無くして読んだとしても、再編集してcommを書き出すと、try:やexceptなどの記述が消えてしまいます。

そこで、この記述の部分を「コメント」にします。
追加の記述部分の行頭に「#」を入れてコメントアウトし、「RESU=・・・」の前のインデントを無くします。そうすると、Eficasで読めるようになります。

Wedge0204

Eficasで読むと下図の様になります。

Wedge0205


再度JOBを実行するときは、commを書き出した後、テキストエディタで「#」を取り去り、「RESU=・・・」の前にインデントを入れます。

以上です。

今回のメッシュファイル(mmed)とcommファイル(実行できる&Eficasで読めない)は下記です。
wedge02.zip

2012年7月18日 (水)

節点反力を集約する方法

材料引っ張り試験の模擬などを行う際は、拘束した節点の反力を集約することが必要になるのですが、その方法がなかなか分からずにいました。

しかし、公式フォーラムを検索していると、下記のやり取りを参考にすれば出来そうと言うことが分かりましたので、試してみました。

http://www.code-aster.org/forum2/viewtopic.php?id=14243

試してみたのは、以前実施したゴム材料の引っ張り試験です(下記リンク)。
http://salome-meca.cocolog-nifty.com/blog/2012/04/post-51b0.html

このときは、本当は変位制御式で反力を出力したかったのですが、反力の集約方法が分からず、やむをえず荷重制御式にしていました。

このときの条件ですが、下図のとおり、ゴム板の一方の切断面であるfixXのX変位を拘束し、反対側のloadP点に2[N]の荷重をかけていきました(fixY,fixZはそれぞれ対称面なので、これらも面内に拘束です)。
loadP点は、load面とX変位を連結しています。

Reac041

このとき、fixX面に属する節点の反力を集約すれば、loadP点の荷重と一致するはずです。それを試してみました。

ソフトウェアは、DEXCS2011-Salome-A3(64)(SALOME6.3 + Code_Aster STA10.5)の利用です。

基本的には、こちらのモデルを流用しています。
まずMeshなのですが、集約する節点を名前で区別するため、節点グループに「fixX」を追加しました。それでもう一度メッシュファイルを書き出しました。

Reac042

次にcommファイルをEficasで編集しました。
反力を出力するため「POST_RELEVE_T」と「IMPR_TABLE」を下図の通り追加しました。

Reac043

POST_RELEVE_TのRESULTANTEは、今回必要な「DX」のみとしていますが、必要に応じて「DY」や「DZ」も指定できるようです。

また、集約する対象の節点は、GROUP_NOに設定します。SALOMEのMeshで名前を付けた節点グループ名です。

このPOST_RELEVE_Tには「reacF」という名前を付けましたが、これをIMPR_TABLEにて、TABLEで指定します。

なお、この2つの詳細については、下記のマニュアルが参考になると思いますが、私もまだちゃんと読んでいませんので、今回の方法以外にもっと効率の良い設定があるかもしれません。

POST_RELEVE_T:u4.81.21.pdf
http://www.code-aster.org/V2/doc/v10/en/man_u/u4/u4.81.21.pdf

IMPR_TABLE:u4.91.03.pdf
http://www.code-aster.org/V2/doc/v10/en/man_u/u4/u4.91.03.pdf

さて、反力の集約に関する部分については以上なのですが、その他にも下記の部分を変更しています。

●計算時間短縮のため、荷重を与えるステップ分割数を、50から20に減らしました。
(DEFI_LIST_REELのNOMBRE)

●収束性向上のため、「STAT_NON_LINE」の「b_meth_newton」-「NEWTON」に、「REAC_ITER = 1」を追加しました。
これの詳細に関しては、こちらの記事を参照ください。

●得られる反力の精度を上げるため、収束許容値を厳しくしました。
「STAT_NON_LINE」の「CONVERGENCE」で「REGI_GLOB_RELA」の値を0.000001にしました。
この値が0.1や0.01など緩めの値だと、荷重と反力に若干のズレが生じます。
それはそれで正常な動作なので別に構わないのですが、今回は関しては誤差があると気持ち悪いので厳しくしています。まあ、気分の問題です。

以上の設定にて計算を実行しました。

その結果、「reac04.resu」ファイル(resuファイル)に下記のように出力されました。

Reac044

荷重はX方向に2.0 [N]を20分割して与えていますが、反力はX方向マイナス側に同じ値が出力されています。これをエクセルに取り込めば、fixXの反力値として使えそうです。

これで、変位制御式の引っ張り試験も模擬できるようになりそうです。機会があれば、また変制御式での反力についても試してみたいと思います。

今回のメッシュファイル(mmed)と、commファイルは、下記です。

reac04.zip

2012年6月26日 (火)

追従荷重(その2)

結果出力の間隔設定」および「追従荷重」にて、MAEDAさんよりコメントでアドバイスを頂きましたので、それで解析内容を改善してみました。

ソフトウェアは、DEXCS2011-Salome-A3(64)(SALOME6.3 + Code_Aster STA10.5)の利用です。

まず、結果出力の間隔設定ですが、「IMPR_RESU」-「RESU」ではなく、「STAT_NON_LINE」-「ARCHIVAGE」のほうに「LIST_INST」を設定しました。今回は、time=0~1を20分割しています。

また。「IMPR_RESU」-「RESU」についても、 FORMAT、UNITE、RESULTATだけにしています。
下図はEficasによる編集内容ですが、かなりすっきりしました。

Folf02b01

次に、収束性の改善と計算ステップ数の低減です。

まず、「STAT_NON_LINE」にて「NEWTON」-「REAC_ITER = 1」を追加しました。これに関しては、公式のマニュアル「R5.03.01」の「2.4.4.1節」が参考になります。「REAC_INCR」と「REAC_ITER」は共に、ニュートン法の求解過程において、どのタイミングで接線剛性を更新する(作り直す)かという指定です。

「REAC_INCR」は、これで指定した数値のインクリメント(タイムステップ)毎に更新します。これは公式マニュアル「U4.51.03」の「3.9.2節」によると、デフォルトで「1」になっており、インクリメントが「1」進むたびに更新します。

「REAC_ITER」は、これで指定した数値のイテレーション毎に更新します。イテレーションは、各インクリメント(タイムステップ)内における繰り返し計算です。これは、同じく公式マニュアル「U4.51.03」の「3.9.2節」によると、デフォルトで「0」になっており、同じインクリメント(タイムステップ)内では更新されません。

結局、デフォルトの状態では、「R5.03.01」の「図2 -2.4.4.1-c: Use of the tangent matrix revalued with each time step」のような状態になっていると思われます(図を勝手に転載できないので、お手数ですがリンクをたどってご覧ください)。

この状態は、イテレーション毎の接線剛性計算のコストは省けますが、非線形性が強いと収束しにくくなります。今回のようにゴム材料の大変形という、非線形性の強い問題では、イテレーション毎に接線剛性を計算しなおしたほうが、タイムステップ幅も大きく出来るので、結局速くなるということなのでしょう。今回も、REAC_ITERを「1」とすることで、MAEDAさんのコメントどおり、80ステップまでタイムステップを低減できました。Eficasによる編集内容は下図の通りです。

Folf02b02

最後にソルバーですが、これも「Mumps」に変更しました(モデルが小さいので、「OUT_OF_CORE」はデフォルトの「NON」のままです)。

これで、前回2000ステップで20分ぐらいかかっていた計算が、80ステップで80秒程度まで低減できました。ただしソルバーについてはデフォルト(MULT_FRONT)のままでも80秒程度でした。このあたりはCode_Asterのバージョンによるのかもしれません。

MAEDAさん、どうもありがとうございました。

なお、今回のメッシュファイル(mmed)とcommファイルは下記です。

folf02b.zip

2012年6月23日 (土)

結果出力の間隔設定

非線形解析は、数多くのステップ(インクリメント)にて、変位や荷重を小刻みに変化させて解析を進めます。しかし、全てのステップにおける結果が必要になるわけではありません。結果出力を数ステップ毎にすれば、結果ファイルのサイズが小さくなり、ディスクスペースが節約できます。

前回のゴム材料の引っ張り解析では、Time=0~1の間を50分割(50ステップ)で計算しており、デフォルトではその全てのステップにおける結果が出力されていました。
今回、これを5ステップ毎(つまり全体で10分割)にして、結果ファイルのサイズ低減を図ってみました。

ソフトウェアは、DEXCS2011-Salome-A3(64)(SALOME6.3 + Code_Aster STA10.5)の利用です。

その方法ですが、まず計算ステップの設定と同様に「DEFI_LIST_REEL」で出力間隔を定義します。今回は、"output"という名前で新しくDEFI_LIST_REELを作成し、DEBUT(開始時間)を"0"、JUSQU A(終了時間)を"1"、NOMBRE(分割数)を"10"としました。これで、time=0~10の間を10分割して出力してくれます。

Hype01b1

次に、「IMPR_RESU」-「RESU」の設定で、各RESU(ここではRESU_1~RESU_6)の全てについて、「b_extrac」の下に「LIST_INST」を追加します。そのパラメータとして先ほど作成した「output」を指定します。

Hype01b2


これで終了です。
計算してParaVISのObject Inspector(Informationタブ)で見ると、確かに10等分で出力されています。

Hype01b3

ファイルサイズも、前回の50ステップ全部出力では約5.5MBでしたが、今回の5ステップ毎出力では約1.1MBと、5分の1になりました。

なお、DEBUT, JUSQU A, NOMBREの値を変えることで、結果出力する時間の範囲を自由に指定できます。例えばtime=0~1で実行される解析の最後の方だけを細かく見たい場合、DEBUTの値を0.9などの値にすれば良いです。ただし、出力される間隔と計算される間隔は同期していないといけないようです。つまり、出力される時間ポイントにおいては、計算も必ず行われる必要があるということです(これがズレていると、結果がちゃんと書き出されないようです)。

それから、上記の「RESU」-「b_extrac」の下の"LIST_INST"の設定ですが、ここを"INST"にすると、出力する時間ポイントを直接指定できます(出力時間指定は、カンマ区切りで指定します)。
この場合、「DEFI_LIST_REEL」の作成は不要ですので、出力したい時間ポイントが限られている場合、こちらのほうが一見手軽のように思えます。
しかし、出力時間ポイントの変更などを行う場合、RESUの出力項目数が多いと、それらについて全て打ち変えていく必要があるので面倒です。DEFI_LIST_RELLを一箇所変えるだけのほうが楽なように思えます。

これらの設定については、Code_asterのマニュアル「U4.71.00」(4.3.5節)にて説明されています。

今回のcommファイルとmmedファイルは下記です。

hype01b.zip

2012年3月31日 (土)

自由度連結の方法(LIAISON_UNIF)と反力の出力

節点自由度の連結方法を試してみましたので、ご紹介します。
今回、下記の情報を参考にさせていただきました。

●OpenCAE Users Wiki - 「SALOME-Mecaの活用技術」の「AFFE_CHAE_MECAの使い方」資料

http://opencae.gifu-nct.ac.jp/pukiwiki/index.php?plugin=attach&refer=SALOME-Meca%A4%CE%B3%E8%CD%D1%B5%BB%BD%D1&openfile=AFFE_CHAR_MECA_LIAISON_1.pdf

●Code_Aster V10の英文マニュアル u4.44.01
http://www.code-aster.org/V2/doc/v10/en/man_u/u4/u4.44.01.pdf

前回、塑性変形の荷重と変位の関係をプロットしましたが、荷重面を拘束していなかったので、正確な結果とは言えませんでした。そこで今回、荷重面の変位を連結して試行してみました。

Code_Asterでは、この自由度連結に「AFFE_CHAE_MECA」の「LIAISON_UNIF」を用います。

###

まず、前回と同じく上面に荷重を与える条件での試行です。

荷重を与える点(press面とZ軸が交差する点)に「pressP」というグループ名を付けました。このグループはGeometryで作成し、その後Meshにて「Mesh_1」を右クリックし、「Create Group from Geometry」を選択し、「Elements」の方にObject BrowserのGeometryの下の「pressP」を指定して作成しています。

07201

下図は、拘束・荷重に関して前回から変更した部分です。

07202

荷重面の荷重方向(つまりZ方向)の変位を「LIAISON_UNIF」で連結しました。ここで荷重に直角方向(X方向、Y方向)の変位については拘束していません。荷重に直角方向への変形を妨げないようにするためです。

荷重は「pressP」点のみに与えました。前回の荷重値が200MPaだったので、

  200 (MPa) × 5 (mm) × 5 (mm) = 5000 (N)

としました。

計算・結果出力についてですが、下図のとおり追加しました。

07203

節点値の計算(CALC_NO)では、FORC_NODA(節点力)とREAC_NODA(節点反力)を追加しました。

結果出力指定(IMPR_RESU)では、RESU_5にFORC_NODA、RESU_6にREAC_NODAのそれぞれ追加して、ポスト処理で参照できるようにしました。

これで計算した結果が下図です。
変位は10倍表示で、4000N(160MPa)で引っ張った状態です。

07204

荷重はpressP点の1点のみに与えていますが、応力や変形は均等になっていることが分かります。

ここで、pressP点について、DEPL(変位)、FORC(節点力)、REAC(節点反力)を出力してエクセルでグラフにしました。

07205


これを見ると、pressP点の荷重値(FORC・青色の線)は、引っ張り側への荷重最高点で10Nしかありません。この時点では外部から4000Nの荷重が作用していますが、残りの3990Nは連結されている他の節点に振り分けられています。そしてpressP点は、それら他の節点からの反力を受けています(REAC・赤色の線)。

変位連結した面のいずれかの節点に荷重を与えた場合、その荷重値は、荷重節点のFORC値からREAC値を引いた値となります。上の例では、「10 - ( -3990) = 4000」です。

これはデータ書き出し後に計算する必要があり、間違いの元になります。したがって、このような場合の荷重値としては、結果からプロットするのではなく、Timeを元にして、与えた荷重値を計算するべきと思います。

なお、pressP点の荷重値ですが、全体荷重4000Nのうち10Nということで、10/4000=1/400となっています。
press面には均等な応力が作用しており、同じ面積の正方形要素面が10x10=100面あります。したがって、要素1面あたりの荷重は全体の1/100です。要素1面には節点が4つあるので、節点1個あたりの荷重値は、1/400となります。pressP点は要素1個のみに所属しており、他の要素と共有されていませんので、全体の1/400(つまり10N)がそのまま荷重値となります。

press面の辺上の節点は2要素で共有されているので20N、面内の節点は4要素で共有されているので、40Nとなります(下図参照)。

07206

この事例のメッシュファイル(mmed)とAsterファイル(comm)は下記にzip圧縮して置いています。

http://salome-meca.cocolog-nifty.com/samples/07-20d.zip

###

次に、press面を変位で制御する方法での試行を行いました。

本来ならLIAISON_UNIFを用いる必要は無く、press面の全節点をDDL_IMPOでZ方向に一律変位させれば済む話ですが、変位を連結した上で、変位条件を与える点を1点のみとしたとき、その点に反力が集中すれば、反力を読み取る点が一つで済むのではないかという目論見です。

commファイルの変更点は下図の通りです。

07207

時間とともに変化するpressP点へのZ方向変位を1辺の初期長さである「5mm」として、DEFI_FONCTIONによる時間変化を
(0,0)
(1, 0.0015) *弾性範囲
(2,0.016)  *塑性範囲
とひずみ値にすることで、5mm×(ひずみ)=(変位)になるようにしました。

DEFI_LIST_REELは、time=1を超えたところ(塑性域のはじめのところ)を細かい刻みにしています。

なお、この試行については、引っ張りの片道だけ(ひずみ0.16まで)にしています。

これで計算させて反力を表示させたところ、反力はpress面全体に均等分布してしまいました。

07208

pressP点の1点に反力がまとめられるという目論見が外れたことになります。

私が普段仕事で使っているFEM解析ソフト(Marc)では、この方法によって1節点に全反力を集約できるので、よく使う手法なのですが、Code_Asterの場合はどうすればよいのか調べる必要がありそうです。

ビーム要素を剛体扱いにしてpress面と結合させれば出来そうな気もしますが、結果出力のためだけにそのようなモデルを作るのもちょっと面倒です。

この事例のメッシュファイル(mmed)とAsterファイル(comm)は下記にzip圧縮して置いています。

http://salome-meca.cocolog-nifty.com/samples/07-20e.zip

ParaviewによるデータのプロットとCSVファイル出力

前回、日本語例題7.2(塑性-結果の検証)を実施しましたが、データをグラフにプロットして結果を確認するとことが未実施でした。

解説文書「07-20.pdf」では、データをPost-Proで表示させて値を読み取っているようですが、STEP数が多いとそれも大変なので、私は「Paraview(ParaViS)を用いてグラフを表示させる」方法と、「CSVファイルへの出力->WindowsのEXCELででグラフ化する」方法の2通りを試してみました。

Paraviewの使い方、特に今回使用した「選択機能」、「時間に沿ってプロットする」については、下記の日本語解説文書を参考にしました。

サイト
http://paraview.org/Wiki/The_ParaView_Tutorial

ファイル名
http://paraview.org/Wiki/images/0/03/ParaViewTutorial38-jp.pdf

まず、プロット対象となる節点(ポイント)を選択する必要がありますので、その準備として、メニューから「View」-「Windows」-「Selection Inspetor」を表示させました。

それから、変形表示のために表示状態にしていた「Warp by Vector」フィルタを非表示にして、元形状「07-20b.med」の方を表示し、ツールバーのドロップダウンリストで表示の方法を「Points」にしました。

これで節点が選択できる状態になったので、ツールバーのアイコンから「Select Points through」をクリックし、画面に表示されているポイントから、press面とZ軸の交差しているところの節点を囲むようにして選択しました。

0720c1

本来は「Selects Points On Surface」で選択すべきかと思いますが、今回の例では目標の点がどうしても選択できませんでした。これについては、メッシュ作成時点で対象節点にグループ名を付けておけば、Selection Inspectorでその名前にて選択できたかもしれません。

選択した節点番号は、「Selection Inspector」-「Point Label」で「Visible」のチェックを入れ、Label Modeを「MED_POINT_ID」にすると表示されました。

対象の節点が選択できたら、次にグラフの作成です。

メニューから「Filters」-「Data Analysis」-「Plot Selection Over Time」を選択してから、「Object Inspector」-「Properties」タブ上の「Apply」ボタンをクリックすると、Object browserに「PlotSelectionOverTime」が追加され、表示領域にグラフが表示されました。

0720c2

このグラフは初期状態で、横軸が時間、縦軸が選択した節点のパラメータになっていましたので、これをカスタマイズしていきました。

Object Browserで「PlotSelectionOverTime」を選択した状態で、「Object Inspector」-「Display」タブを表示させ、「X Axis Data」に「RESU_DEPL(2)」、Line Seriesの中で「RESU_SIEF_NOEU(2)」のみにチェックが入った状態にしました。

RESU_DEPL(2)は変位(Z成分)、RESU_SIEF_NOEU(2)は応力(Z成分)です。これらのカッコ内の数字は各成分の順になっており、0がX成分、1がY成分、2がZ成分になります。

プロット点は●にしました(Maker Style -> Circle)

0720c3

更に軸見出しやフォントなどを調整しました。これは下図の赤丸の部分(グラフの左上)をクリックすると、「View Settings」画面が立ち上がりますので、そこで行いました。

0720c4

これで「応力-変位グラフ」が完成しました。解説文書(07-20.pdf)の荷重ー変位のグラフと同じになっていると思います(縦軸の符号が逆ですが)。

#######################

完成したグラフをキャプチャするなどしてレポートに貼り付けてもいいのですが、データやグラフを加工するなら慣れたエクセルで行いたいと思いましたので、データのCSV出力も行いました。

これはまず、「PlotSelectionOverTime」のグラフの右上のView分割ボタン(日)を押して、新たに作られたViewで「Spreadsheet View」を選択します。

0720c5

これでスプレッドシートが現れますので、Showingが「PlotSelectionOverTime」になっていることを確認します。また、精度(有効数字)を必要に応じて変えます。それから「File」-「Export」でcsvファイル名を指定して書き出します。

0720c6

これでCSVファイルが出力できました。ファイルの内容はグラフの内容にかかわらず、結果出力されている全ての節点値が書き出されています。これをエクセルに取り込んで編集とグラフ作成を行いました。

下図は、横軸を変位(Z成分)から計算したひずみ(Z成分)として、縦軸を応力(Z成分)としたものです。

0720c7

これはZ方向の応力とひずみですが、解説文書の相当応力-相当ひずみのグラフと近い結果になっています(全く同じではありません)。今回の事例が単軸変形に近い状態であるためと思われます。

今回は以上です。

2012年3月17日 (土)

Paraviewでのアニメーション

Salome-Mecaでの解析結果の可視化には、Paraview(SALOMEのメニューではParaViS)が使えます。

前回実施した日本語例題「07-10」についても後半の結果はParaviewで表示させてみましたが、やはりPost-Proより高機能だし、説明書を見なくても直感的に使える部分が多いです。

とはいえ、細かい部分は説明を見ないと分かりにくい部分もあります。これについては、下記URLに詳しい日本語解説があります。

http://www.paraview.org/Wiki/The_ParaView_Tutorial

(Japanese Translation of the ParaView Tutorial for version 3.8.)というリンクです。

これを読めばほとんど問題なく使えると思いますので、ここでは使用法の説明は省略し、使ってみた例を紹介します。

前回の日本語例題(塑性変形-負荷を変化)の結果について、変形前形状をメッシュのワイヤーフレームで、変形後の形状をサーフェース&Mises応力のコンターで表示し、それらを半分のカット図にしてみました。

0710b1

Paraviewの表示は、元の状態に色々なフィルターを多段接続してカスタマイズしていけますが、これは元の形状に「Warp By Vector」を接続して変形とMises応力のコンター表示をさせ、更に「Clip」フィルターでカットしています。そして表示させるのは「Clip」の方であり、元の形状は非表示にしています。

それから元の形状にもう一つ「Clip」フィルターを繋ぎ、変形前のワイヤーフレームも半分にカットしました。これらフィルターの接続状況は画面左のブラウザーで確認できます。

細かいことを言えば、ワイヤーフレームの色をデフォルトの白から黒に変えたり、カット面を非表示にするなどの変更もしています。

アニメーション表示も簡単で、単に画面上部のアニメーション制御のボタン類を押すだけで表示できます。ここで私が感心したのは、アニメーションの表示モードの「Real Time」です。これは、「View」-「Windows」-「Animation View」で画面下がアニメーション制御画面になり、そこで選択できます。

0710b2

通常パラパラ漫画形式のアニメーションでは、コマ間隔が等間隔になりますので、「07-10」のように、時間軸方向にコマが不均等に並んでいる場合、細かい出力間隔の所ほどスロー再生となってしまします。しかし「Real Time」モードでは、時間の進みを一定にしてアニメーションを進行させてくれます。実時間に対する倍数も設定できます(解析時間の1秒を実時間の10秒で再生するなど。今回はそうしてます)。

さて、このアニメーションは動画出力も出来ます。

「File」-「Save animation」で表示される画面にて、出力する動画の設定を行いました。

0710b3

出力時間は10秒、フレームレート15fps、解像度は表示されているまま(874x481)としました。これで「Save Animation」を押すだけです。

ただし、出力されるのは「ogv」形式の動画となります。DEXCS-Salomeでは、OSに「動画プレーヤー」が入っており、これでogv形式動画を見る事が出来ますが、Windows機やWEBに出したい場合は変換する必要があります。

これも都合良く、DEXCS-SalomeのOS(Ubuntu)に「PiTiVi」という動画編集ソフトが有り、これでAVI形式に変換できました。

0710b4

「クリップの取り込み」で作成したogvファイルをライブラリに追加し、それをタイムライン部分にドラッグ&ドロップするだけで、編集はしていません。そのまま「プロジェクト」-「プロジェクトのレンダリング」で、ファイルの出力を行いました。

ここでサイズを640×480とし、フレームレートを15fps、コンテナをAvi muxerとしました。音声は入っていないのですが「無し」という設定が無かったので、そのままです。
(後で気が付きましたが、アスペクト比を全く考慮せずサイズ変更していました。今後はちゃんと計算が必要です)。

これで出力ファイル名をaviの拡張子で設定して、「レンダリング」を押すと、avi形式の動画ファイルが書き出されました。

これはWindows機でもメディアプレーヤーで見る事が出来ます。後は、WMVやMP4に変換したり、WEBにアップロードしたり好きなように出来ますが、ブログに貼り付けるならYouTubeにアップロードするのが簡単そうです。

このように、Paraviewは色々と凝ったことが出来そうですが、凝り出すとどんどんと時間を食っていきそうなので、程々にしておく必要がありそうです。

今回は以上です。

より以前の記事一覧