« 2012年7月 | トップページ | 2012年9月 »

2012年8月

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月13日 (月)

ゴム(超弾性)の物性値を同定する方法(Code_Asterでの検証)

昨日、ゴム材料の物性同定をエクセルにて行いましたが、その検証をCode_asterにて行いました。

検証のための解析モデルには、先日実施した「節点反力を集約する方法」のゴム平板の単軸引っ張りモデルを流用しました。

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

流用したモデルの変更点は下図の通り、ゴム物性と、引っ張り条件のみです。

Hype051

ゴム物性(DEFI_MATERIAU)は、同定したC10,C01,C20に変更し、引っ張り条件を変位制御式(DDL_IMPO)に変更して、引っ張り変位を15にしました。元の長さが10なので、工学ひずみが1.5になります。

上図にはありませんが、出力分割数も30にして、工学ひずみ0.05毎に出力されるようにしています。

下図が結果です。

Hype052


同定計算の真応力は、公称応力に伸長比を掛けて換算した値です。
Code_Asterの結果ですが、節点反力(REAC_NODA)の集約値を、元の断面積で割った値を公称応力としました。また、SIEF_NOEUは真応力になるようです。

これによって、エクセルによる同定計算と、Code_Asterの計算が一致することが確認できました。

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

hype05.zip

2012年8月12日 (日)

ゴム(超弾性)の物性値を同定する方法

Code_Asterでゴム(超弾性)を扱う場合の材料物性として、公式ドキュメントU4.43.01」では、「DEFI_MATERIAU」にて「ELAS_HYPER」を指定すると、Signorini形式の材料が使えると書かれています。この超弾性材料の設定に関する理論は「R5.03.19」で解説されていますが、ひずみエネルギー密度関数の定数として、C10, C01, C20 の3つを使う形式のようです。このうち、C20を省くと「Mooney-Rivlin」形式、C01とC20の両方を省くと「Neo-Hooke」形式になります。
「Ogden」形式は、ドキュメントには記載がないのでサポートされていないのかもしれません。

いずれの形式を使うにせよ、実験データからこれらの定数を同定する必要があります。商用ソフトでは、実験データから材料定数を同定する機能(カーブフィッティング)が用意されています(プリプロセッサに組み込まれていることが多いようです)。しかし、SalomeやCode_Asterにはこの機能は無いようです。公式フォーラムでは「MACR_RECAL」(U4.73.02)で出来るというような情報もありますが、なんだか難しそうです。

そこで今回は、多くの人が持っているであろう「EXCEL」で、超弾性の物性を同定してみます。
なお、ゴム物性をきちんと把握するためには、単軸の引っ張りだけではなく、1軸拘束2軸引っ張り、2軸均等引っ張りなどの2軸引っ張り試験データも必要ですが、ここでは単軸引っ張りだけで簡便に扱います。

まずは予備知識として、ひずみエネルギー密度関数(W)や、伸張比(λ)、ひずみの不変量(I1, I2)などに関する基礎理論を理解することが必要ですが、これについては、ここでは書ききれません。
様々な書籍や、WEBでも検索すれば多くの情報が得られると思います。
私がよく参考にしているのは下記の書籍です。

「設計のための高分子の力学」 著:深堀美英、技報堂出版(2000年) (ISBN 4-7655-0389-5)

ここで、今回の同定に用いる式は下記の通りです(詳細は上記の書籍などをご参照ください)。

(1)ひずみエネルギー密度関数 W

Hype041
ゴム物性を定義する式です。変形によるエネルギーは、ひずみの不変量I1,I2の関数になるという考えです。ひずみの不変量は、材料のひずみの度合いを表す指標のようなものです。これはI1,I2,I3の3種類ありますが、上式ではI1とI2の二つを使っています。

この式は自分の目的に応じて好きなものを使って良いのですが、今回用いる式は、Code_Asterで用いるC10,C01,C20の各項に加え、汎用性を考えてC11の項も追加しました。

(2)ひずみの不変量I1, I2を、伸張比λから計算する。

Hype042_2
単軸引っ張りの場合において、伸張比λから、ひずみの不変量I1とI2を計算する式です。これは、変形モード(単軸引っ張り、2軸均等引っ張りなど)によって異なってきます。
伸張比λは、工学ひずみεに1を足したものです(λ=ε+1)。

(3)ひずみエネルギー密度関数Wを、I1,I2で偏微分した式

Hype043

Hype044

これは、(1)の式をI1, I2で偏微分すると得られます。

(4)公称応力σを求める式

Hype045
この式と、(2),(3)式によって、公称応力σを、伸張比λから計算できます。工学ひずみをεとすると、λ=ε+1なので、公称応力と工学ひずみの関係が得られます。

これらの式をエクセルのシートにまとめたものが下記のファイルです。

hype04-mooney4.xls

画面のスナップショットです。

Hype046


まず、実験データ(工学ひずみと公称応力)を左側に入れます。
各行とも上述の式に従い、伸張比、I1,I2, dW/dI1, dW/dI2, 計算した公称応力が表示されます。
(行が足りない場合は数式を含めてコピーして増やし、余る場合は削除してください)。

このとき、左上の定数(C10, C01, C11, C20)を適当に打ちかえると、計算された公称応力が様々に変化することが分かります。しかし、計算結果を実験に合わせるのは、手動での打ち変えで調整するやり方では難しいです。そこで、エクセルの「ソルバー」機能を使います。

ソルバー機能はエクセルのデフォルトのインストールでは組み込まれていないかもしれませんので、無い場合は組み込みます(「エクセル ソルバー アドイン」で検索してください)。私の使っているエクセル2010では、「データ」タブの「分析」にボタンが追加されました。

ソルバーによる同定の方針ですが、公称応力の「実験値」と「計算値」の「差の2乗」を最小化します。エクセルシート上で、各行の残差の2乗を計算し、それらの合計をまとめます(K3セル)。
ソルバーの目的セルをその(K3)にして、目標値を「最小値」にします。
変数セル(変化させるセル)は、C10(C1), C01(C2), C20(C4)を複数選択します。Code_AsterではC11は使わないので、C3セルの値はゼロにしておき、変数セルには含めません

C10, C01, C20の初期値を適当に入れ、「解決」を押すと、適切な数値が入ります。
ここで、「制約の無い変数を非負数にする」にチェックを入れると、定数は全て正の値になります。定数に負の値が含まれると、カーブフィットの精度は上がりますが、変な極値がでて挙動が不安定になる場合があるので、注意が必要です。

同定が完了しました。

Hype047

このようにして、単軸試験におけるゴム物性を同定できました。
これを発展させると、2軸試験を含めた複数カーブの同時フィットも出来ます。
ただし、エクセルのソルバーは指数関数は苦手のようで、Ogdenのカーブフィットは上手くいきませんでした。

なお、今回のエクセルファイルですが、説明用に手早く作成したものなので、数式の間違いなどがあるかもしれません。
業務には使用しないよう、お願いいたします。

2012/8/13 追記

*同定値のCode_Asterによる検証はこちら

*今回の同定対象データは、手違いにより、ひずみ0.5~1.0の間が0.1刻みになっていました(他の区間は0.05刻み)。計算内容や同定結果には影響はありません。

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年8月 4日 (土)

異なる「DEFORMATION」設定の混在

前回の「ゴムと金属の接触・摩擦解析(平面ひずみ)」において、「STAT_NON_LINE」-「COMP_INCR」の「DEFORMATION」の設定を、材料毎に分けて設定しました。ゴムを「GROT_GDEP」、鉄を「PETIT」にしたのですが、同じ解析計算の中で異なる「DEFORMATION」が混在することになります。

これまで使ってきた何種類かの商用ソフトではこのようなことは出来なかったので、それで良いのかどうか、今回、簡単な剛体回転の例題で確かめてみました。

ソフトウェアは、DEXCS2011-Salome-A3(64)(SALOME6.3 + Code_Aster STA10.5)の利用です。
解析モデルは下図のとおりです。

Nl0101


100mm角の正方形Faceを二つ並べ、平面ひずみ解析を行うために四角形2次要素でメッシュを切りました。左側をfaceA、右側をfaceBとしています。
それぞれ、左上の頂点をX,Y共に拘束し、右下の頂点にX方向「-100mm」の強制変位を与えることにしました。これで、時計回りに45度の剛体回転が起こるはずです。

Eficasでの境界条件設定は下図の通りです。

Nl0102


fixA, fixB点をX,Y共に拘束。
dispA,dispB点に、X方向 -100mmの変位を、0~1秒で、20分割で与えます。

STAT_NON_LINEの設定です。

Nl0103


COMP_INCR_1は、faceA(左側)で、DEFORMATIONを「GROT_GDEP」に設定。
COMP_INCR_2は、faceB(右側)で、DEFORMATIONを「PETIT」に設定。

また、この例題は剛体回転で反力が出ないため、収束条件を反力の絶対値としています。

計算結果は下図です。

Nl0104

右側(PETIT)は、回転に伴って正方形が膨らんでいることが分かります。これは微小変位設定で大きな回転をさせたときに起こります。それに対し、右側(GROT_GDEP)はそのようなことはなく、元の形状のままで剛体回転しています(コンターは変位量にしているので色が付いていますが、大きさは変わっていないことは分かります)。

これは、その動画です。右側の正方形が回転に伴って膨らんでいます。

このように、同じ計算の中で違うDEFORMATIONを混在させても良さそうです。ちゃんと確認するには、それぞれの応力やひずみの値を検証する必要がありますが、それはかなり手間がかかりますので、またの機会にしたいと思います。

それ以前に、計算オプションや出力オプションによって、どのような種類の応力やひずみが出力されるのか、把握する必要がありそうです。

PETITはもちろん工学ひずみと公称応力だと思いますが、あとはオプションや材料の組み合わせで、「グリーン・ラグランジュひずみ&第2ピオラ・キルヒホッフ応力」になるか、「対数ひずみ&コーシー応力」になるかが違ってくるのではと思います。

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

2012年8月 3日 (金)

ゴムと金属の接触・摩擦解析(平面ひずみ)

今回は、超弾性材料(ゴム)と金属の接触・摩擦問題について試行しました。

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

Code_Asterの接触・摩擦(Contact)に関するドキュメントとしては、

[U2.04.04] Instructions for use of contact

[U4.44.11] Operator DEFI_CONTACT

等があります。しかしこれらは理論も難しく、私の外国語の読解力不足もあり、なかなか全てを把握するのは難しいです。これから徐々に理解を進めていきたいと思っています。

さて、今回の解析モデルは下図のとおりです。

Wedge0101

クサビ状のゴムを、鉄のV字状の溝に押し込み、引き抜きます。
クサビの角度(片側)は、ゴムが24°で、鉄の溝が20°であり、ゴムのほうが若干広角度になっています。
モデルは「平面ひずみ」で、左右方向は対称なので、半分のみモデル化しました。

このモデルですが、SalomeのGeometryにて外形線を作成してから、「New Entry」-「Build」-「Surface」にて外形線を複数選択(シフトキーを押しながらクリック)して、Faceを作成しました。また角のフィレットは、「Operations」-「Fillet 2D」で、Planar Faceに該当Faceを、Vertexesにフィレットを作成する角部の点を指定して作成しました。

境界条件等を与えるためのグループは下図の通り作成しました。

Wedge0102

ゴムの「dispR」は変位制御でゴムを下方に押し込むための変位拘束エッジ、「fixXR」は対称条件による左右方向(X方向)の拘束エッジ、「contR」は接触を定義するエッジです。鉄の「fixS」は固定点として全自由度拘束するエッジ、「contS」は接触を定義するエッジです。
また、各Faceには材料定義のために、「Rubber」、「Steel」のグループ名を付けました。

メッシュは、三角形2次要素としました。柔らかくて変形が大きいゴムの方のメッシュを、スレーブ側想定として細かくしています。

Wedge0103

メッシュについても、「Mesh」-「Create Groups from Geometry」によって、ジオメトリのグループ名を引き継いでいます。後工程での汎用性を考え、Nodes, Edges, Facesの全てに、出来るだけのグループを作成しました。

その後、「Mesh」-「Build Compound」にてMeshesにMesh_1とMesh_2の両方のメッシュを複数選択(シフトキーを押しながらクリック)で指定し、メッシュをまとめました(Compound_Mesh_1)。このCompound_Mesh_1をオブジェクトブラウザで右クリックし、「Export to MED file」にて、メッシュファイルに出力しました(wedge01.mmed)。

次に、Eficasによるcommファイルの編集です。今回は項目がかなり多いですが、順に上から見ていきます。

●材料と要素の定義です。

Wedge0104

材料定義(DEFI_MATERIAU)は、ゴムはMooney-Rivlin材料として、鉄は線形弾性材料としました。これらの材料を「AFFE_MATERIAU」にてそれぞれの要素グループに割り当てました。AFFE_MODELEにおける要素の種類は全て「D_PLAN」(平面ひずみ要素)としました。

Wedge0105

「MODI_MAILLAGE」-「ORIE_PEAU_2D」は各接触面におけるエッジの方向を揃えるために入れています。これはSolomeのAsterの線形解析Wizardを用いて、複数エッジにわたるグループに対して荷重を設定したときにも自動生成されるのですが、これまでの私の非線形解析試行では抜けていたようです。これ無しで計算すると接触処理でエラーになります。

●今回の主題である接触定義の部分(DEFI_CONTACT)の部分です。

Wedge0106_2


「FORMULATION」は4種類あって、デフォルトは「DISCRETE」ですが、今回の試行では「CONTINUE」のほうが大きな押し込み変位まで計算できましたので、こちらにしています。これらの理論は下記ドキュメントで解説されていますが、難しいです(MarcのDiscreteおよびAnalyticalと同じようなものと考えて良いのでしょうか。それならCONTINUEのほうが滑らかな気がします)。

[R5.03.50] Discrete formulation of the contact-friction

[R5.03.52] Contact elements derived from a hybrid continuous formulation

「FROTTEMENT」(摩擦)は、クーロン摩擦(COULOMB)です。

収束計算においての繰り返し計算回数をデフォルトから増やしました。これらは、Code_Asterの計算エラー時のメッセージに「この値を増やせ」という意味の記述があったので、それに従いました。

・ITER_GEOM_MAXI(デフォルト10)-> 50
   幾何的なリサイクル
・ITER_FROT_MAXI(デフォルト10)-> 50
   摩擦力計算のリサイクル?

「ZONE」は接触面の定義です。鉄の面(contS)をマスター(GROUP_MA_MAIT)に、ゴムの面(contR)をスレーブ(GROUP_MA_ESCL)にしました。摩擦係数は0.3です。

●境界条件の定義です。

Wedge0107

fixXRはゴムの対称面の拘束、fixSは鉄の底面と側面の完全拘束です。
ゴムの押し込みは、dispR面を変位制御で、下方(Y方向マイナス向き)に1.6mmとしました。本当はもっと多く押し込みたかったのですが、これ以上の押し込み量では計算が出来ませんでした。
この1.6mmの変位を、0秒~0.5秒間で与え、0.5秒~1.0秒の間で0に戻します(つまり
引き抜きます)。DEFI_FONCTION(depl_imp)で定義しています。

●時間ステップの定義です。

Wedge0108


自動時間ステップを使用しています。自動時間ステップの詳細については「時間ステップ幅の自動調整機能」をご参照ください。

初期ステップ幅は、0.01(1秒間を100分割)としました。押し込み50ステップ、引き抜き50ステップになります。
摩擦計算はイテレーションが多くなりそうですので、時間幅を増やす条件 (ADAPTATION)において、NEWTON法のイテレーション回数の判定基準を「15」と大きめにしています。

●結果出力間隔を指定する設定です。

Wedge0109

結果出力間隔の設定については「結果出力の間隔設定」をご覧ください。
今回は、初期時間ステップと同じ、0.01(1秒を100分割)としています。

●非線形解析の設定です。

Wedge0110


接触(CONTACT)には、上で定義した「contact」を設定します。

それから、「COMP_INCR」について「COMP_INCR_1」と「COMP_INCR_2」の2つを定義しています。
はじめは「COMP_INCR_1(RELATION:ELAS_HYPER)」の内容を全要素に対して設定していましたが、それで計算すると、Steelの材料がELASであるためか、エラーが発生しました(“ELAS_HYPERが無い”というような意味のメッセージが出ました)。
そこで、「COMP_INCR_1」はRubberグループに限定し、「COMP_INCR_2(RELATION:ELAS)」を新たに追加してSteelグループに設定しました。

これで計算は一応出来たのですが、この設定の仕方で良いのかどうかは、よく分かりません。「COMP_INCR_2」ではうっかり「DEFORMATION」の指定を忘れたのですが、この場合、SteelのDEFORMATIONはデフォルトの「PETIT」になったようです(messファイルより)。Steelは微小変形で変位も小さいのでそれで構わないのですが、そもそも、GROT_GDEP(大回転・大変位・微小ひずみ)とPETIT(微小回転・微小変形・微小ひずみ)を同じ計算の中で混ぜて良いのでしょうか。その場合、出力される応力やひずみは何になるのかも気になります。
しかし、今回はひとまず保留です。

「ARCHIVAGE」-「LIST_INST」には、上で定義した結果出力の間隔リストである「output」を指定しています。

●計算する応力等の指定です。

Wedge0111_2


●反力の集約値を出力するための設定です。

Wedge0112


節点反力の集約値を出力するための方法については、「節点反力を集約する方法」をご参照ください。

変位制御エッジである「dispR」のY方向の節点反力(REAC_NODA)を集約しています。
dispRを指定するには、要素グループ(GROUP_MA)ではなく節点グループ(GROUP_NO)で指定する必要があるようです。

結果はresuファイルに出力されます。

●ParaVIEWで表示するmedファイルへの、結果書き出しの設定です。

Wedge0113


RESU_1~5の項目が無ければ、計算したデータが全て書き出されますので、普通はそのほうが良いです(MAEDAさんからの情報)。

しかし、今回は全部書き出されると具合が悪いことがありました。

「DEFI_CONTACT」において、「FORMULATION」を「CONTINUE」に設定した場合、節点変位(DEPL)に「LAGS_C」という成分(自由度)が追加されます。これはスレーブ側の接触圧力らしいのですが、2次元解析では「DX,DY,LAGS_C」の順となり変位の第3成分になります。これをポストプロセッサ(ParaVIEW)で読み込んだときに、このLAGS_CがZ方向変位として扱われるので、その状態で変形図を出すと、無いはずのZ方向変位が生じてしまいます。

Wedge0116

そこで、RESU_1でNOM_CHAM、NOM_CMPを用いて、DEPLの出力をDXとDYのみに限定しています。
そのため、他の応力や反力などについても、RESU2~RESU5で出力指定しています。

私は最初はこの情報を知らず、結果表示でZ方向に変な変位が出て驚いたのですが、公式フォーラムを検索して、下記のやりとりから解決法を得ました。
http://www.code-aster.org/forum2/viewtopic.php?id=13405

以上でcommファイルの編集は終わりです。

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

さて、上記の設定で計算させた結果が、下の動画です。

ゴムが押し込まれ、引き抜かれる様子が分かります。
押し込み時に鉄のR部分にゴムが引っ掛かるようで、その部分での変形が大きく、あまり大きく押し込めないようです。形状の設定がまずかったかもしれません。

摩擦なし(DEFI_CONTACTで、FROTTEMENTを削除)の条件でも計算し、resuファイルからエクセルに反力値を読み込んで、変位と反力の関係を比較しました。

Wedge0114

摩擦係数0.3の結果は、摩擦なしよりも反力がかなり大きくなっています。また、押し込み時と引き抜き時で違う経路を通ります。押し込み時のほうが大きな反力となっていますが、この行き・帰りの線で囲まれた部分が摩擦で失われた仕事になります。
クサビということで、引き抜き時に反力がマイナスになる(引き抜き力が必要になる)結果を期待したのですが、そのようにはなりませんでした。その様な結果を出すには、もう少し形状を調整し、また、荷重制御にて解析する必要があるかもしれません。

摩擦なしの線は一本に見えますが、行き・帰りで同じ経路をたどるため、2本の線が重なって見えています。

最大押し込み時の変形図でも、摩擦の有無で変形の違いがあるようです。

Wedge0115

何回か試行しましたが、摩擦なしの場合は「CONTACT」の「FORMULATION」を「DISCRETE」にしても、かなり大きな変位で押し込むことが出来ています。

このように、とりあえずはゴムの接触・摩擦解析を行うことが出来ました。しかし、設定内容に分からないところや、あやふやな所が数多くあります。まだこれから確認すべき点は多いようです。

今回のメッシュファイル(mmed)とcommファイルは下記です(摩擦係数0.3のケースのみ)。
wedge01.zip

« 2012年7月 | トップページ | 2012年9月 »