« 2012年8月 | トップページ | 2012年10月 »

2012年9月

2012年9月16日 (日)

荷重を与えて変形させた状態で固定する方法(境界条件の切り替え)

「解析対象物に荷重を与えて変形させ、その変形の状態で固定する。
その後、更に荷重や変位条件を与えていく」

という解析をしたい場合があります。
この場合、解析の途中の段階で、荷重条件から変位条件への切り替えが必要になります。

途中から適用される条件も荷重であれば、荷重対象面(または辺・点)を指定しておき、最初の値はゼロで解析の途中から値が増えるように設定すれば良いのですが、変位条件の場合、解析の途中で「変位=0」等という与え方は(通常の方法では)出来ません。これをすると、「変形前の変位=0の点」に戻ってしまいます(多くの場合、変位増分が大きすぎてエラーになります)。

Code_Asterの非線形解析において、このような境界条件の切り替えをしたい場合はどうするかというと、どうやら「STAT_NON_LINE」を2つ続けて記述し、解析を2段構えにすると良いようです。
まず第1段階で引っ張り荷重を与える計算を行い、続く第2段階にて変位条件を「TYPE_CHARGE='DIDI'」のオプション付きで与えます。

これは公式ドキュメント「U4.51.03(STAT_NON_LINE)」の3.2.3項で説明されていますが、「ETAT_INIT/NUME_DIDI」で指定された時点(デフォルトでは計算の再開時点)において、その時点での位置を基準にして変位条件を適用するものです(正しくは“Dirichlet条件”と言うようですが、構造解析においてはおおよそ“変位条件”という理解で良いようです。)。

変位条件なので、固定(変位0)に限らず、指定した時点からの相対変位で変形させていくことも出来るようです。

-------------------------------------------------------
さて、この「TYPE_CHARGE='DIDI'」について検証をしてみました。

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

解析モデル(メッシュ)です。

Bound0201

長さ100mm × 幅10mm × 厚み4mm
の単純なゴム板です。要素は六面体2次要素で、一辺2mmのサイズで分割しました。

解析手順ですが、今回は先に結果を見て頂くことにします。

【第1段階】X方向が長手方向となるゴム板の一方の端面(fix)を固定し、他端面(dispX)に4[N/mm^2]の面荷重(合計では160[N])を与えて引っ張ります。

Bound0202


【第2段階】dispX面を第1段階が終了した時点の位置で固定し、ゴム板の上面(fz)に0.5[N/mm^2]の面荷重(合計で500[N])をZ方向マイナス向きに与えてたわませます。

Bound0203


解析結果の動画です。

この解析のcommファイルについて、ポイントをEficasの編集画面でご説明します。

●ゴム物性

Bound0204

Mooney-Rivlin形式としました。

●境界条件

Bound0205

・CHAR0:fix面の拘束(第1段階、第2段階共通)

・CHAR1:dispX面の拘束(第2段階のみ適用)

・fx:dispX面への荷重(第1段階のみ適用、fonc_fxで徐々に増加)。

・fz:fz面への荷重(第2段階のみ適用、fonc_fzで徐々に増加)。

●非線形計算部分(STAT_NON_LINE)

・STAT_NON_LINE を2段重ねにしました。

Bound0206

・第1段階は、time = 0~1
・第2段階は、time = 1~2

・第1段階は単純に引っ張るだけです。

・第2段階で、第1段階の計算完了時の状態(応力や変形)を取り込むため、「ETAT_INIT」を指定しました。これも色々とオプションがありますが、今回のように単に第1段階の完了時の状態を指定するだけなら、「EVOL_NOLI」で良いようです。

・境界条件は、fxに代わってCHAR1を、「TYPE_CHARGE='DIDI'」のオプション付きで指定しました。これで、ETAT_INITで指定した時点の状態に対し、変位条件が与えられます。

・更に、ゴム板上面への荷重「fz」を追加しました。

これらの他、「自動時間ステップ調整」、「結果の出力間隔の指定」、「反力の集約値を出力」などの機能を使っています(詳細は各記事とcommファイルをご覧ください)。
時間ステップの指定については、第1段階と第2段階でそれぞれ別に定義しています。

設定については以上です。

結果(変形図)は上の方で示した通りですが、反力についても下図の通り、与えた条件通りの値が出ています。
(X方向の反力は、Z方向の荷重の増加につれて増えています)。

Bound0207


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

さて、今回の主題は以上なのですが、今回のように複数の段階を続けて計算する場合、第2段階の計算で何らかの不具合があって異常終了すると、もう一度第1段階からやり直しになってしまいます。第1段階に多くの計算時間を必要とするケースの場合、これはかなり非効率です。

このようなとき、「POURSUITE」を使うと、複数の段階をそれぞれ分けて、ひとつづつ計算出来ます。

「POURSUITE」の使い方ですが、
OpenCAE Users Wiki」の「オープンCAE初心者勉強会

にて、MK2さんが「SALOME-Mecaの活用技術

を公開されており、この中の「ポスト処理」-「Post_Lin_ver10.pdf 」で詳細に説明されています。

この「SALOME-Mecaの活用技術」は非常に役に立つ情報が多く、私はいつもこちらの情報を参考にさせて頂いております。また、現在はGoogleの「https://sites.google.com/site/codeastersalomemeca/」の方にも活用情報をまとめられつつあります。

【余談】
あと、実戦的な情報源としては「公式フォーラム」の情報が充実していると思います。「こんなことは出来ないのかな?」というような疑問があったとき、ここでキーワード検索してみると、おおむね何か出てきます。Code_Asterを始めたばかりの私が悩むようなことは、既に多くの方が色々と考えてきたことなのでしょう。しかし英語やフランス語なので、細かなニュアンスとかがちょっと分かりにくいですね。

さて、「POURSUITE」を用いて複数の計算を分けて実行する場合、計算データを格納する「.base」フォルダをASTKで指定する必要があります。この「.base」フォルダですが、計算が正常に完了した場合には最新の状態に更新されるようです。そのため、例えば第2段階だけについて設定を変えて何度も実行したいというような場合、第2段階の実行時には、終了後にデータが更新されないよう、「R」のチェックを外しておく必要があります。

Bound0208

ただしこの場合には、第2段階の計算結果は「.base」には保存されません(medやresuなどのファイルには、もちろん出力されます)。別名でバックアップを取っておく等の方法でも良いかもしれません

今回のファイルは下記のzipにまとめています。

bound02.zip

bound02.mmed (メッシュファイル)
bound02.comm (2段階分を通して計算)
bound02-2a.comm(2段階を分けて計算・第1段階)
bound02-2b.comm(2段階を分けて計算・第2段階)

なお、第1段階終了時にも結果を可視化したいため、結果出力の処理を行っています。これについても、結果出力の部分だけを別のcommにするなど、様々な実行パターンが考えられます。

2012年9月 8日 (土)

COMP_INCRとCOMP_ELASの違い(弾塑性編)

前回の検証事例「非線形解析の設定と、出力される応力・ひずみの種類について(弾塑性編)」では、検証の実施に「COMP_INCR」のみを用い、「COMP_ELAS」は省略しました。

その理由ですが、私が参考にした公式サイトの解説ドキュメント

公式サイトの「Support」-「Code_Aster trainings」のページにある「Non linear capabilities

では、「COMP_ELAS」は弾性挙動をするという説明があったためです。つまり、弾塑性解析にCOMP_ELASを用いた場合、負荷過程と除荷過程で同じS-S曲線上を通ってしまい、永久ひずみが残らなくなります。

この辺りの説明は、もちろん公式ドキュメントにも記載されていますが、読解がなかなか難しいです。
上述した「Code_Aster trainings」のシリーズのような、分かりやすい解説資料が出てきたのは有り難いです。

さて、今回はその「COMP_INCR」と「COMP_ELAS」の、弾塑性解析における挙動の違いを確認しました。

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

解析モデルですが、2種類の挙動を一度に並べて比較するため、前回の1要素のみのモデルをコピーして2要素にしました。

下図のように、SalomeのMeshにて、「Modification」-「Translation」で、Y方向に2mmの距離で要素をコピーしました。

Nl0501


このとき、「Generate Groups」にチェックを入れると、元のメッシュに設定されていた各グループも一緒にコピーされました。ただし、コピーされたグループ名は長い名前になるので、手作業で短い名前に直しています(“fixX2”など)。

このように2つの要素とグループが出来ました。

Nl0502


向かって左側の要素を「COMP_INCR」用、右側を「COMP_ELAS」用とし、これらを要素ごとに分けて設定するため、Volumeのグループ(incrとelas)も追加しました。

次にcommファイルのEficasによる編集です。

材料は、前回の「公称応力 - 工学ひずみ」の定義による弾塑性材料です。

Nl0503


境界条件は、各要素の各方向3面を面内拘束します。
荷重は、X方向に2N[N]を、time=0->0.5で与えて、.time=0.5->1.0でゼロに戻します。

Nl0504


そして、STAT_NON_LINEにて、左側の要素に「COMP_INCR」、右側の要素に「COMP_ELAS」を設定しました。

Nl0505

さて、計算結果です。

まず荷重が最大となった時点です(コンターは、荷重方向の応力です)。

Nl0506


荷重が増える方向では、両者は同じです。

次に、荷重をゼロに戻した状態です。

Nl0507


COMP_INCRでは永久ひずみが残り、COMP_ELASでは初期形状に戻ってしまうことが分かります。

動画です。

荷重と変位の関係をグラフにしたものです。

Nl0508


COMP_ELASは、一本の線上を往復しています。

このように、弾塑性解析における「COMP_INCR」と「COMP_ELAS」の挙動の違いを確認することが出来ました。

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

nl05.zip

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月 | トップページ | 2012年10月 »