条件設定

2013年1月 3日 (木)

ゴムバンドの摩擦力による動力伝達(問題点)

「ゴムバンドの摩擦力による動力伝達」について、前々回(前編)前回(後編)の2回に分けてご紹介しましたが、また問題点が残っていました。

これについて改良を試みましたが、まだ解決出来ていません。現在の状況をご報告します。

利用ソフトウェアは、Salome-Meca2012.2(SALOME 6.5 + Code_Astr 11.2.0)で、OSはdebian6です。

【問題点1】

ゴムバンドの応力分布にムラがあり、ゴムバンドの動きに合わせて動いてしまう。

ゴムバンドが二つのローラーに掛けられて定常状態で回転しているとき、応力や歪はメッシュ(つまりゴムバンド)について動いていくのではなく、空間上の「位置」によって定まるようになるはずです。これは、パイプを流れる液体の流体解析において、液体の状態がその「位置」によって決まるのと同じです。

ただし、流体解析は実際に空間に対してメッシュを作成するのですが、今回は構造解析ですので、ゴムバンドについてメッシュを作成しています。それでも、定常状態における応力は位置によって決まるはずなのです。つまり、ゴムバンドを基準として見た場合、ある決まった応力分布が、ゴムバンド上を一定速度で移動していくことになります。

しかし、前回の解析では、応力の大きい部分がゴムバンドの移動につれて動いています。

Rband0301


【問題点2】

従動ローラーの回転が、駆動ローラから遅れている。

駆動ローラーはtime=2~4の間に2回転するので、従動ローラーも2回転するはずですが、下図の通り、やや遅れています。すべりが大きいのかもしれません。

Rband0302


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

これらの問題を解決するため、色々試行錯誤してみました。

まず【問題点1】ですが、一番疑わしいのは、メッシュが不均質なために応力にムラが生じているのではないかという事です。そこで、ゴムバンドは出来るだけ形の揃った四辺形(2次要素)でメッシュを作成しなおしました。
また、ローラーについても、表面の要素辺サイズが揃うようにしました。

Rband0303

それから怪しいのが、剛体移動防止のために設定している「弱いバネ」です。これにゴムバンドが引っ張られているのかもしれません。
前回のモデルでは、バネ定数は全て10のマイナス3乗でしたが、これを10のマイナス6乗まで小さくしました(それ以上小さくすると、不安定になり収束しませんでした)。

Rband0304


次に【問題点2】です。
前回は、摩擦力の計算アルゴリズム(ALGO_CONT)を「STANDARD」にしていました。これでは、ローラーを回転させる段階において、「COEF_FROT」(デフォルト値:100)を、「10」まで下げないと計算が収束しませんでしたので、その値にしています。これが小さすぎて、摩擦力が小さめになっている可能性があります。

そこで今回は、ALGO_FROTを「PENALISATION」に変えました。これまでの経験から、すべりの大きな摩擦計算では、「STANDARD」よりも「PENALISATION」のほうが収束しやすかったためです。その上で、「COEF_PENA_FROT」を「500」にしました。ゴムの物性値がヤング率換算で5[MPa]程度ですので、二桁大きな値になります。

Rband0305

また、初回インクリメントでなかなか収束しにくい事が多かったのですが、それに関しては「SEUIL_INIT」に「0.1」を設定する事でかなり改善しました。これは初期の接触面圧を与えるもので、剛体移動の抑制に効くようです。この値は大きいほど効果があるようですが、反面どのような悪影響(副作用)があるか分かりませんので、出来るだけ小さな値で済むようにしました。今回は0.1程度より大きくすると収束性改善効果が感じられましたので、0.1にしました。

これについては、公式ドキュメント「u2.04.04」の4.2節に解説があります。
http://www.code-aster.org/V2/doc/default/en/man_u/u2/u2.04.04.pdf

これによると、継続して実行するSTAT_NON_LINEにおいては「CONTACT_INIT='INTERPENETRE'」も指定したほうが良さそうなので、これも入れていますが、これの効果はいまひとつ分かりません。

なお、計算のステップは下記の通りで、

(1)time 0 -> 1
 駆動ローラーを固定(drp1, drp2を固定)し、従動ローラー中心(center)に右向きの荷重を与え、ゴムバンドを所定の軸荷重で張っていく。

(2)time 1 -> 2
 ゴムバンドを張った状態で、従動ローラーをcenter位置でピン固定し、軸荷重は除去。

(3)time 2 -> 4
 駆動ローラーを反時計回りに2回転させる(time=1あたり1回転)。
 このとき、回転始め~1/4回転迄の間で、従動ローラーのトルクを増加させて所定の値にする。

前回は(2)と(3)の間でJOBを前編と後編に区切っていたのですが、(2)の最初で収束しにくい状況だったため、試行の効率を考え、今回は(1)と(2)の間で区切りました。

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

以上のような変更を行って計算したのですが、最初に書いた通り、残念ながら解決には至りませんでした。

まず【問題点1】ですが、やはりゴムバンドの応力にはムラが生じており、ゴムバンドと一緒に移動してしまいます。

Rband0306


このゴムバンドの「応力ムラ」についてよく見ると、どうやら初期状態からバンドを引っ張った状態(Time=1~2)のときの応力状態を「引きずって」いるようです。

Rband0307


このときに応力の大きい部分が、回転しても応力が大きいままです。
しかしこれはおかしい事です。
塑性変形しているなら分かりますが、これは超弾性材料なので、変形途中の応力状態に影響されるはずがありません。あくまでも初期状態(真円)のときが応力ゼロとして計算されるべきです。

なにか、塑性変形的な設定をしてしまっているのかと思い、STAT_NON_LINEの中の「COMP_INCR」を「COMP_ELAS」に変えても結果は同じでした(塑性解析では違いが出ます)。

また、摩擦の影響かと思って、Time=0~2までを「摩擦なし」で計算させてみましたが、それでも同じでした。
低減積分要素を用いても変わらず。

ということで行き詰りました。

どうも、STAT_NON_LINEを連続して実行するときの、前段階からの状態の引継ぎ設定【ETAT_INIT=_F(EVOL_NOLI=RESU,),】あたりが怪しそうなのですが・・・。

それから【問題点2】についても、前回の結果と同じで、従動ローラーの回転が遅れています。

Rband0308

しかしこれは、もしかしたら正しいのかもしれません。せん断変形などによって、実際にずれていくということも考えられます。もう少し薄いバンドで、屈曲変形の影響を排除して確かめる必要が有りそうです。

以上のように、まだまだ検討の余地が有ります。

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

rband03.zip

(ZIPファイルの内容は下記の通りです)

メッシュファイル: rband03.mmed
commファイル(前編): rband03a.comm
commファイル(後編): rband03b.comm

2012年12月29日 (土)

Eficasのマニュアルと、便利な自動翻訳

CAE懇話会の「Salome-Meca活用研究会」で、Sさんから教えて頂いた情報です。

メニューなどが英訳されている「Salome」に比べ、「Eficas」はメニューやボタンがフランス語で、使い方もいまひとつ分かりにくいです。

しかし、SalomeのWizardで出来ない設定はどうしてもEficasを使う必要があるので、私は今まで手探りでなんとなく使ってきました。

しかし、どうやら英文のマニュアルがちゃんとあるようです(私が今まで気が付かなかっただけですが・・・)。

マニュアルの場所は、Salome-Meca2012.2においては下記です。WWWブラウザから下記のローカルファイルを開くと閲覧できます。

/opt/SALOME-MECA-2012.2-LGPL/SALOME/Salome-V6_5_0p1-OT-LGPL-x86_64/modules/EFICAS_V6_5_0p1/share/doc/salome/gui/EFICAS/index.html

このマニュアルの場所は、Salome-Mecaのバージョンによって異なるようですが、概ね似たような場所に保管されているようです。

内容はHTMLで、ブラウザから見る形式になっています。

英語で書かれていて日本人には分かりにくいので、そこはブラウザの自動翻訳を活用すると、分かりやすくなります。

Sさんによると「Google Chrome」がお勧めとの事です。ただし、自動翻訳のためにはインターネットに繋がっている状態である必要があります。

2012年12月 9日 (日)

ゴムバンドの摩擦力による動力伝達(後編)

前回(前編)の続きで、2本のローラー間にゴムのバンドを掛けて、その摩擦力でローラーの回転と動力を伝達させるという事例です。

(前編の記事掲載以後、様々な改良を行っていますが、今回はとりあえず改良前のモデルの後編の説明です)

利用ソフトウェアは、Salome-Meca2012.2(SALOME 6.5 + Code_Astr 11.2.0)で、OSはdebian6です。

前編では、ローラーに所定の軸荷重を与えてゴムバンドを張り(time=0~1)、その位置で従動ローラの中心を固定し、軸荷重および安定化の為の拘束を除去したところまで(time=1~2)でした。

後編では、time=2~4で、駆動ローラーを反時計回りに2回転させます。つまり、time=1あたり1回転になります。また、回転はじめから1/4回転までの間で、従動ローラーの内側突起部の側面(trq)に負荷トルクとしての荷重を与えていき、1/4回転からは一定のトルク値を保持します。

以下、後編のcommファイル(rband01b.comm)の説明です。

■計算の継続

Rband0117


前回のデータを読み込み、継続して計算をする場合には「POURSUITE」を用います。
Eficasで新規にcommファイルを作成する際に「POURSUITE」を追加すると、前回のcommファイルを指定するようにメッセージが表示され、ファイルを選択するBOXが開きますので、そこで前回の「rband01a.comm」を選択します。
また、「POURSUITE」から始まるcommファイルを開いたときも同様で、前回のcommファイルを指定します。
この作業によって、前回のcommファイルで定義済みの変数や関数、パラメータなどがそのまま使えるようになります。

■駆動ローラーの回転変位設定(drrot)

Rband0118



time=2~4で、駆動ローラー上の2ヶ所の点(drp1, drp2)に変位を与え、駆動ローラーを反時計回りに回転させる設定です。time=1あたり1回転とします。
2ヶ所の点のx,yの変位を、INST(time)の関数として、それぞれ関数(FORMULE)として定義し、「AFFE_CHAR_MECA」で境界条件「drrot」として設定しました。

このように変形体を回転させる方法についての詳細は「回転体とゴム板の接触」をご参照ください(参照先の制御節点はX軸上に並んでいますが、今回はY軸上に並んでいますので、式の内容は異なります)。

2ヶ所の点は駆動プーリの内径ライン上に位置するので、振幅は内径の半径である「44mm」になります。

■負荷トルクの設定(trq)

Rband0119_2


従動ローラーの内径側に出ている突起の側面(trq面)に面圧荷重を与え、トルクを発生させます(突起の詳細など、モデルについては前編をご覧下さい)。

従動ローラー中心からtrq面中央までの距離が43[mm]、圧力値が1.16[N/mm^2]、trq面の面積が2[mm^2]なので、トルクは、
    43 x 1.16 x 2 = 100 [N・mm] となります。

厳密に言えば、trq面への荷重は半径方向に分布しているので、上記のような計算式ではなく積分計算をする必要がありますが、ここはtrq面が微小であるとして、単純な計算で済ませています。

■接触の設定

Rband0120



接触の設定「DEFI_CONTACT」は既に前編のcommファイル(rband01a.comm)で設定済みなので、それと同じ設定で良ければ、このcommファイルで再設定する必要はありません。しかし今回の試行では、前編の設定のままでは後編の回転部分においてどうしても計算が収束しませんでした。そのため、前編(cont)と異なる名前(cont2)にて再設定しました。

前編と違うところは、「COEF_FROT」の値を「10」にした所です。前編はデフォルト値の「100」でしたので、その1/10になっています。「COEF_FROT」は公式ドキュメント「U4.44.11」でもあまり詳しい記述が無いのですが、おそらく、接触する両面間に相対変位が生じたときに、摩擦力が変位と共に立ち上がるときの傾きのようなものではないでしょうか(これについては自信がありません)。もしそうだとすると、この値が小さくなると、相対変位に対する摩擦力の増加が緩やかになり、収束しやすくなる反面、現実と比べて十分な摩擦力が得られない可能性もあります。

■タイムステップの設定(time=2~4)

Rband0121



pas3aは、回転開始のtime=2から、1/4回転してトルクが所定値に達するtime=2.25までを180分割と細かくし、それからトルクが一定値を保持するtime=4までは315分割と粗くしています。

output3は結果の出力間隔ですが、これは2回転の間を均等に360分割しています。つまり、駆動ローラーが2°回転する毎に結果が出力されます。

pas3bは時間ステップ幅の自動調整です。これは前編とあまり変わりませんが、ステップ幅を増やす判断に用いるNewton法の収束回数基準値を、9回と多めにしました。

■非線形解析の設定(time=2~4)

Rband0122



境界条件は、今回定義した「drrot(駆動ローラー回転変位)」と「trq(従動ローラートルク荷重)」を追加しています。従動ローラーセンターのピン固定は継続です。

ここで、「trq」については追従荷重(TYPE_CHARGE='SUIV')を指定し、ローラーが回転しても、常にtrqは突起側面に垂直になるようにしています。この設定で、常に一定のトルクが作用するようにしています。下の動画をご参照ください。

動画を見られない方は、下の静止画をご覧ください。

Rband0124

また、追従荷重の詳細については、「追従荷重」をご参照ください。

それから、収束条件を残差の相対値「REGI_GLOB_RELA」で、値を「1e-4」にしました。
デフォルトは「1e-6」なのですが、それでは収束しなかったため、緩めの判定値にしました。

■結果の計算と出力

Rband0123



これは特に変わった事はしていません。

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

Rband0125



baseについてですが、前編でデータを「rband01a.base」フォルダに書き出していました。
したがって後編ではそこから読み出しても良いのですが、それでは後編の計算が終わったときに上書きされてしまいます。
後編の計算を何度も試行錯誤したかったので、「rband01b.base」フォルダを別に作り、試行のたびに「rband01a.base」の内容を上書きコピーしてから計算をスタートさせるという方法にしました。

以下は結果の図です(Mises応力)。ローラーの応力は不要なので非表示にしています。

●time=2

Rband0126


引っ張っただけなので、上下対称の変形・応力状態です。

●time=2.25

Rband0127


駆動ローラーが反時計回りに1/4回転した状態で、トルクが所定値に達しています。
両ローラー間でゴムバンドが張られている所(スパン)では、トルクによって、上側が強く張られていて応力も大きい事が分かります。

●time=3

Rband0128


駆動ローラーが反時計回りに1回転した状態です。
しかし従動ローラーは1回転に若干足りていない事が分かります。
これは、ローラーとゴムバンドの「すべり」が原因と思われます。

●time=4

Rband0129


駆動ローラーが反時計回りに1回転した状態です。
従動ローラーは更に遅れています。トルクが一定状態でも「すべり」があるようです。

それに、応力の分布状態に「ムラ」があります。これについては改善が必要です。

改めて結果の動画を載せます。

このように、一応この種の解析が可能である事が分かりました。
しかし、まだ下記のような問題が残っており、実用にはそれらの解決が必要です。

  • 応力や変形が若干いびつになっている。これはバネの影響が残っているか、メッシュが粗すぎるなどの原因が考えられる。
  • 計算が不安定で、タイプステップが制限値まで細かく分割されても収束せず、JOBが止まる事が多い。
  • 計算時間が長く掛かる。これは、ローラーを剛体面扱いできず、変形体として解いている事が大きく影響していると思われる。剛体サーフェスの導入か、せめて剛体結合時にも大回転が出来るようにしてほしい。
  • 接触・摩擦計算については更なる高速化と、安定性が必要。
  • これらの問題については、開発元の対応がない限りどうしようもない事もありますが、自分で改善できそうなところは、出来るだけ詰めていきたいと思います。

今回のデータは下記です(前編と同じです)。

rband01.zip

メッシュファイル: rband01.mmed
commファイル(前編): rband01a.comm
commファイル(後編): rband01b.comm

2012年12月 3日 (月)

ゴムバンドの摩擦力による動力伝達(前編)

今回は、2本のローラー間にゴムのバンドを掛けて、その摩擦力でローラーの回転と動力を伝達させるという事例です。設定内容が複雑なため、前編と後編の2回に分けることにします。

利用ソフトウェアは、Salome-Meca2012.2(SALOME 6.5 + Code_Astr 11.2.0)で、OSはdebian6です。接触・摩擦計算が高速化されたCode_aster11.2を使いたかったため、いつものDEXCS-Salomeではなく、この環境になりました。

また、モデル化は2次元(平面ひずみ)としました。はじめから3次元で行うと、試行錯誤に大変な時間がかかるためです。

まず、ローラーとゴムバンドのモデルです。

Rband0101


ゴムバンド(band)は真円の状態が初期形状です。その中に、中空のローラー2個を配置します。向かって左側が駆動ローラー(dr)、右側が従動ローラー(dn)です。drはdriving、dnはdrivenの略です。

寸法は、
・ゴムバンド(外径:260mm、内径:240mm)
・ローラー (外径:100mm、内径:88mm)

従動ローラーの拡大図です。

Rband0102


センター位置をピン固定するためのアームが1本出ているのと、トルクを与えるための突起(高さ2mm)が、内側へ出ています。

上図のように各グループを設定しましたが、それぞれの役割は下記の通りです。

【面グループ】
・band:ゴムバンドの領域設定
・dr, dn:各ローラーの領域設定

【線グループ】
・bandc:ゴムバンドの内面(接触面の定義)
・drc, dnc:各ローラーの外面(接触面の定義)
・trq:従動ローラーに負荷トルクを与えるための圧力荷重面

【点グループ】
・bandp1, bandp2, dnp1:不安定性対策の弱いバネ、および拘束の設定箇所
・drp1, drp2:駆動ローラーの回転変位制御点
・center:従動ローラーの軸荷重点、および回転中心

このモデルについて、下記の手順で計算します。

(1)time 0 -> 1
 駆動ローラーを固定(drp1, drp2を固定)し、従動ローラー中心(center)に右向きの荷重を与え、ゴムバンドを所定の軸荷重で張っていく。

(2)time 1 -> 2
 ゴムバンドを張った状態で、従動ローラーをcenter位置でピン固定し、軸荷重は除去。

(3)time 2 -> 4
 駆動ローラーを反時計回りに2回転させる(time=1あたり1回転)。
 このとき、回転始め~1/4回転迄の間で、従動ローラーのトルクを増加させて所定の値にする。

文章で説明しても分かりにくいので、まずは最終計算結果の動画をご覧下さい。

実は、このような解析が、私がSALOME-Meca導入当初からやりたかった事です(実際はもっと複雑ですが)。のんびりやっていたため、導入・試行を始めてから1年以上かかってしまいました。
今年の春頃から、SALOMEや Code_Asterの様々な機能を検証してきたのですが、その多くがこの解析のための前準備です。

次に、メッシュ分割です。

Rband01mesh_2

全て三角形2次要素で分割しました。
band, dr, dn はそれぞれ別にMesh作成を行い、Compoundでまとめています。
ただし、座標が重なる節点についてのマージは行っていません(回転によってズレていくため)。
もちろん、各グループ名をジオメトリから設定しています。

さて、ここからが今回(前編)の設定内容なのですが、これまでのようにEficasの画面スナップショットを用いると画像が膨大な量になってしまうため、commファイルの内容を用いる事にします(それでもかなり多いですが)。
前編のcommファイル名は「rband01a.comm」です。

前編では、上記の解析手順の(2)まで、つまり、バンドを張った状態で従動ローラーを固定し、それまで与えていた軸荷重を抜くところまでです。残りの回転部分については、後編にて継続計算(POURSUITE)で行います。

■メッシュの修正(エッジの方向)

Rband0103

接触や面圧荷重を受けるエッジ(線)の方向を整える設定です。

■弱いバネ要素の作成

Rband0104


初期の不安定な状態において計算の収束性を向上させるため、ゴムバンド2ヶ所、従動ローラー1ヶ所に弱いバネを作成します。駆動ローラー近くで変位が小さい所(bandp1)をspring1、従動ローラー近くで大きく変位する所(bandp1,dnp1)をspring2と二つに分け、バネ定数などの設定を分けられるようにしています。

■要素の種類の設定

Rband0105

ゴムバンド(band)、ローラー(dr, dn,trq)は平面ひずみ「D_PLAN」、バネ(spring1, spring2)は「2D_DIS_T」としました。ここで、負荷トルクのための圧力荷重面(trq)も「D_PLAN」に含める必要がありました(その理由は分かりません)。

■バネ定数の設定

Rband0106

弱いバネのバネ定数の設定です。種類は「K_T_D_N」です。
ここの値の設定は、何度も値を変えての試行錯誤が必要で、かなり苦労しました。
バネが弱すぎると初期の不安定性を抑えきれないし、強すぎるとゴムバンドが大きく伸びたときにバネの影響が出てしまいます。最終的に、spring1, spring2共に“1e-3”で落ち着きました。

私が使っている商用ソフトでは、バネ定数をTimeの関数にして、バネを徐々に消滅させていく事が出来るのですが、Code_Asterにはそのような機能が無いようです(“VALE_F”はあったのですが、Timeの関数には出来ませんでした)。この機能はぜひ入れて欲しいところです。

■材料物性の設定

Rband0107

ゴムバンドはMooney-Rivlin材料としました(「超弾性材料の検証」を参照)。

また、ローラーは剛体にしたかったのですが、「LIAISON_SOLIDEによる剛体は大回転出来ない」の結果から、やむなく変形体として、大きなヤング率(1000[GPa])を設定しました。

■接触の設定

Rband0108

Code_aster 11.0.28以降では、接触・摩擦解析が高速化」の結果から、Code_Aster11.2.0を用い、
 ・ALGO_RESO_FROT=’NEWTON’
 ・ALGO_RESO_CONT=’NEWTON’
を設定しました。

また、ゴムバンドを「マスター」、2つのローラーを「スレーブ」としています。
硬い材料であるローラーをマスターにしたかったのですが、複数のZONEに共通している面がある場合、それをマスターにしないといけないようです(今回の場合はbandです)。

摩擦係数は0.4で、他にはあまり調整したところはありません。

■タイムステップの設定(time=0~1)

Rband0109


pas1aは、従動ローラーを右(X方向)に引っ張り、ゴムバンドを張っていく段階(time=0~1)のタイムステップ定義です。解析初期(time=0.5まで)の不安定なところの刻みを多くしています。

output1は、出力の時間刻み設定です。
これに関しては「結果出力の間隔設定」をご参照ください。
今回は、time=0~1を均等に50分割しています。

pas1bは、pas1aを基本として、タイムステップの自動調整を追加したものです。これについては「時間ステップ幅の自動調整機能」をご参照ください。
今回のような解析は非常に不安定で収束性が悪く、この機能が無ければ、現実的な時間で計算する事が難しかったと思います。

最小時間幅(PAS_MINI)の値は小さめにしています。

また、時間幅を増やす場合の増加割合(PCENT_AUGM)のデフォルトは100(%)なのですが、時間幅が元の2倍になると再び収束しなくなって、もう一度分割される、という非効率な現象が良く起こっていたので、これを50(%)にしました。

■タイムステップの設定(time=1~2)

Rband0110

time=1~2では境界条件の切り替えを行うだけなので、時間刻み幅を「1」としてます。出力の刻み幅も「1」です。
一応、収束しなかったときのために自動調整機能も入れていますが、実際は動作せず、時間刻み幅は「1」のままでした。

■境界条件の設定

Rband0111


・fix_dr:駆動ローラーの変位固定
・fix_y_c:従動ローラーのセンター位置を、上下(y方向)に固定
・force_dn:従動ローラーを右に引っ張る軸荷重

       10[N]としています。
       平面ひずみ要素を用いているので、厚み(Z方向)が1mmあたりの値になります。

・fix_x_c:従動ローラーのセンター位置を、左右(x方向)に固定
・fix_y_b:弱いバネを設定した節点を、y方向に固定
      バネがあるので本来この条件は不要なはずですが、これを入れないと収束しませんでした。

・force:force_dnを増減させるための関数です。time=0~1では増加します。特に解析初期であるtime=0.5までの増加量を少なくしています。
time=1~2でゼロになります。

■非線形解析の設定(time=0~1)

Rband0112


time=0~1では、駆動ローラーを固定し、従動ローラーに軸荷重(force_dn)を与えてゴムバンドを引っ張っていきます。そのため、「fix_x_c」は使いません。

COMP_INCRでは、ローラー(ELAS)とゴムバンド(ELAS_HYPER)を異なるRELATIONに設定しました。また「異なる「DEFORMATION」設定の混在」も可能なのですが、今回は両方とも「GROT_GDEP」にしています。

NEWTON=_F(REAC_ITER=1,), は、ニュートン法の求解過程において、イテレーション毎に接線剛性を更新するものです。デフォルトではインクリメント(タイムステップ)毎になっているので、非線形性の強い計算では収束性が悪くなります。詳しくは「追従荷重(その2)」を参照ください。

■非線形解析の設定(time=1~2)

Rband0113

time=1~2では、ゴムバンドを引っ張った状態で従動ローラーのセンター位置を固定(fix_x_cの条件を追加)し、軸荷重(force_dn)を除去します。また、後編でのローラー回転に備え、安定化のためのy方向固定(fix_y_b)も除去します。

ここで、従動ローラーのセンター一固定(fix_x_cの条件)については、初期位置から変位した場所で固定するため、「TYPE_CHARGE='DIDI'」の設定が必要になります。これに関しては「荷重を与えて変形させた状態で固定する方法(境界条件の切り替え)」を参照ください。

また、time=1からの継続計算なので、「ETAT_INIT=_F(EVOL_NORI=RESU,)」を設定しています。

■結果の計算と出力設定です。

Rband0114


これは特に問題ないと思います。

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

以上の設定で、ASTKを用いて計算を実行しました。

Rband0115


この後に後編で続きの計算(ローラー回転)を行うので、データを「rband01a.base」フォルダに保存する設定にしておきます。

time=2の時点での結果です(Mises応力)。

Rband0116


ゴムバンドが伸ばされて応力が生じている事が分かります。
現段階では単に引っ張っただけですので、変形や応力状態は上下対称になっています(メッシュによる誤差を除いて)。

なお、ローラーも変形体なので応力が出力されているのですが、評価対象ではないため非表示にしています。そのほうが、ゴムバンドの応力が見やすくなります。

前編はここまでです。

今回のデータは下記です
rband01.zip

メッシュファイル: rband01.mmed
commファイル(前編): rband01a.comm
commファイル(後編): rband01b.comm

2012年11月 3日 (土)

回転体とゴム板の接触

前回の検証で、LIAISON_SOLIDEによる剛体結合は大回転できないことが分かりました。

そこで今回、変形体を用いて回転体の解析を行うことにしました。しかし硬くて変形しない部品を単に回転させるだけでは面白くないので、ゴムとの接触も含めた例題を作成してみました。

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

まず、解析モデルを下図に示します。

Rigid0301


上側の穴の開いた長細い部品が回転体です。これは穴を中心に回転し、下方のゴム板に接触します。本来はこの部品を剛体面または剛体結合として扱いたかったのですが、それが出来なかったため、変形体として扱うことにしました。

ゴム板は荷重が作用していない状態で両端(fix1,fix2)を拘束されており、回転体の接触によって下側へ押されます。

回転体の回転のさせ方ですが、まず、Geometryにおいて下図のように回転体を平面で分割し、分割面と穴の内面の交わる線(2本)を、それぞれ「line1」、「line2」というグループ名にしました。

Rigid0302


この2本の線に変位を与え、回転体を動かします。本来ならば穴の内面の点全てに変位を与えたほうが良いのでしょうが、今回は柔らかいゴムに接触するだけなので、この2本の変位拘束だけで荷重を十分支持できると判断しました。荷重が大きい場合は、回転体の弾性率も大きくする必要があります。

寸法は特に重要ではないのですが、一応下図に示します。
穴の内径は12mm(半径6mm)です。

Rigid0303


次にメッシュですが、回転体は四面体2次要素、ゴム板は六面体2次要素で作成しました(それぞれ別にメッシュ作成し、Compoundで一つにまとめています)。
ゴム板は曲げ変形するため、厚み方向に2層以上は欲しいところですが、今回は例題ですので計算時間を優先し、1層にしています。
もちろん、Geometryから各グループを作成しています。

Rigid0304

続いて、Eficasによるcommファイルの編集です。

複数の面を含む面グループに接触の設定をする場合、「ORIE_PEAU_3D」で面を揃える必要があります(これを結構忘れます)。

Rigid0305

材料物性ですが、回転体は構造用鋼としました。剛体として考えるなら、もっと大きな弾性率のほうが良いかもしれません。ゴムはMooney-Rivlinです。

Rigid0306

接触の設定ですが、摩擦は「なし」にしました。
FORMULATIONはcontinueです。
面の設定ですが、本来は弾性率の大きい回転体の方をマスターにするべきなのでしょうが、メッシュが回転体のほうが細いため、ゴム板のほうをマスターにしました。このあたりの設定によっても、収束性が異なってくるかもしれません。

Rigid0307

時間ステップですが、0秒~1秒の間で回転体が1回転するとして、その1秒間を計算します。
初期時間分割は40分割、結果出力も40分割、ただし、自動時間ステップ調整機能をを用います。

Rigid0308

さて、ここからが今回の本題である、回転させるための境界条件設定ですが、まずは下準備として、line1,line2に対する変位を時間の関数で定義することが必要です。

X-Y平面で見ると、line1, line2の初期配置は下図の通りです。

Rigid0309

この図では穴の中心座標が(0,0)となっていますが、実際には異なります。しかし、与える変位量としては同じことです。
穴の半径 r=6(mm)、周期 T=1秒なので、変位は下式の通りになります。

【line1】
X変位 = 6 * ( cos(2π(t-0.5)) + 1)
Y変位 = 6 * ( sin(2πt) )

【line2】 
X変位 = 6 * ( cos(2πt) - 1)
Y変位 = 6 * ( sin(2π(t-0.5)) )

これらの式を、Eficasで「FORMULE」に関数として、それぞれ設定していきます。
INSTは時間(=t)、piはπ(円周率)です。

Rigid0310

そして、定義した関数を「AFFE_CHAR_MECA_F」にて、それぞれの変位として設定します。
なお、line1,line2共に、Z方向の変位は0としました。これは、AFFE_CHAR_MECA_Fの中に入れても良いし、私のようにAFFE_CHAR_MECAを別に作っても良いようです。

STAT_NON_LINEの設定は下図です。
特に変わったことはしていません。回転体(solid)とゴム板(rubber)は別々に「COMP_INCR」を設定しました。

Rigid0311

一応、fix1,fix2のY方向(上下方向)の反力値も出力するようにしています。

Rigid0312

このような設定で計算を行いました。

結果の動画です。

静止画です

Rigid0313


回転体が回転して、ゴム板を押す様子が計算できました。

ゴム板の固定点(fix1, fix2)の反力も妥当な形で出ています。

Rigid0314


今回、X-Y平面上で4つの自由度を拘束しているため、ちゃんと変位を計算して回転させているつもりでも。計算誤差などによってわずかな変形と応力が発生しているようです。今回は実質的に影響が無いくらい小さいのでOKですが、この点にも注意する必要がありそうです。
本来は自由度3つ(1点のX,Yと、もう一点のXかY)の方が良いのですが、それだと回転の向きが定まらない時点があり、結局4つの自由度拘束になっています。何か、うまい方法は無いものでしょうか。

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

rigid03.zip

2012年10月21日 (日)

面荷重(FORCE_FACE)を荷重の合計値で与えるには

3次元要素の面への荷重(「AFFE_CHAR_MECA」-「FORCE_FACE」)は、単位面積当たりの荷重値で入力します。SI単位なら、[N/m^2]になります。

様々なケースの構造解析をしていると、荷重を対象の面全体の合計値で与えたい場合もあるのですが、残念ながらCode_Asterではそのようなコマンドは無いようです。

荷重を与える面が円や四角など単純な場合は、あらかじめ面積を求めておき、合計荷重値をその面積で割った値を入力すれば良いのですが、計算できないような形状の場合はどうすればいいのでしょうか。

そのような疑問が生じ、公式フォーラムなどを少し検索してみたのですが、簡単な方法は無いようです。
たとえば、

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

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

では、「POST_ELEM」の機能を流用して、面グループの面積を求める方法が紹介されていますが、これらは本解析とは別にJOBを実行する必要があり、面倒です。Pyrhonスクリプトを駆使すれば、1JOBにまとめることも出来るようですが、それは更に面倒です。

結局、Code_Asterだけで解決するのはあきらめ、CAD(Geometry)の段階で対象面の面積を求めておくことにしました。形状を3次元CADで作成しているのであれば、そのCADの機能を使えば良いし、Salomeなら、「Geometry」に面積を求める機能があります。

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

解析モデルは下図の通りです。
断面が1[mm] × 1[mm]、長さ10[mm]の角柱を作り、一方の端面の一辺にR=1[mm]のフィレットをつけます。

Forcef0101

フィレットを付けた方の面を荷重面としてグループ名を「force」にし、他端を固定面としてグループ名を「fix」にします。
この荷重面の長手方向(X方向)に合計1[N]の荷重を与えて引っ張るものとします。

このとき、荷重面がフィレット無しの四角のままなら、断面積は1[mm^2]なので、単に1[N/mm^2]を与えれば済みます。しかし、フィレットを作成したことにより断面の面積が変わっていますので、その断面積を求めておく必要があります。
(今回の「force」面の場合は明らかに、π/2[mm^2]なのですが、一応「わからない」ということにしておきます)。

面の面積の求め方ですが、Geometryのメニューから「Measure」-「Basic Properties」で、下図のようなウィンドウが出ます。

Forcef0102


ここで、オブジェクトブラウザーで該当する面グループを選択すると、その面の面積が、「Surface is:」の欄に表示されます。「force」面の面積は、1.570795[mm^2]と表示されています(これはπ/2であり、計算値と一致します)。
また、「Length」は面を囲む線の合計長さのようです。

このモデルで単純な線形解析を行いました。
メッシュ図です。

Forcef0103


四辺形2次要素を用いました。

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

Forcef0104


FORCE_FACEのFXの値を
1÷1.570795 = 0.63662031
としました。

また、fix面のX方向反力の合計値を出力するようにしています。

Forcef0105


この方法の詳細については、「節点反力を集約する方法」をご覧ください。

結果ですが、まずfix面のX方向反力の合計値は「resu」ファイルに出力され、下記のようになりました。

#TABLE_SDASTER
INTITULE         RESU     NOM_CHAM         NUME_ORDRE   INST       
DX
Reaction_Force   RESU     REAC_NODA                   1  0.00000E+00
-9.99999E-01 

反力はほぼ1[N]となり、設定どおりの力が作用していることが分かりました。

下図は、X方向の応力図です。

Forcef0106_2


角柱の上面の応力が大きく、下面が小さいことが分かります。
変形も、よく見ると「お辞儀」しています。

これは、下図のように荷重が角柱の上面側に偏るためです。

Forcef0107


「FORCE_FACE」は対象の面に(デフォルトで)均等分布しますが、荷重は法線方向ではなくグローバル座標の成分になることに注意する必要があります。
法線方向の荷重(つまり圧力)にしたい場合は「PRES_REP」になります。

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

forcef01.zip

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月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のみです。応力はコーシー応力(真応力)になります。

今回は以上です。