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

2012年6月

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日 (土)

追従荷重(Follower force)

追従荷重は、その荷重が作用している要素面の変形(回転)に追従して、荷重の向きも変化するタイプの荷重のことです。従動荷重ともいいます。
Code_Asterの荷重はデフォルトでは追従しないタイプになっていますが、これは、「STAT_NON_LINE」-「EXCIT」で、TYPE_CHARGEに"SUIV"を指定することで追従荷重にすることが出来ます。

以下に、私が試行した例を紹介します。
ソフトウェアは、DEXCS2011-Salome-A3(64)(SALOME6.3 + Code_Aster STA10.5)の利用です。

以前に実施した「超弾性材料の検証」のモデルをそのまま流用し、片持ち梁として曲げることにしました。モデルは幅2.5mm、長さ10mm、厚さ1mmのゴム平板で、1辺0.5mmの立方体で2次メッシュを作成しました。要素数は、5 x 20 x 2 = 200になります。

Folf0201


境界条件を与えるためのグループ名の作成ですが、流用したモデルにはX=0, Y=0, Z=0 の各面にそれぞれfixX, fixY, fixZという名前が既に付いていますので、このうちfixXの面を完全固定面とすることにしました。また、fixXの反対側の端面にはloadという名前が付いています。

また、荷重を与える面ですが、下図のように固定端(fixX)と反対側の端の平板上面に、1要素分の長さで「loadP」という面を作成しました。

Folf0202

最初は、固定端と反対側の端面(「load」面)に、「FORCE_FACE(面荷重)」で下向き(Z方向のマイナス向き)に荷重を与えて曲げようとしていました。これは、通常荷重では上手くいったのですが、追従荷重にすると下記エラーが出て計算できませんでした。

La charge definie dans STAT_NON_LINE en tant que une charge de type suiveuse,
sous le mot-cl? TYPE_CHARGE = 'SUIV' n'est pas une Charge SUIVEUSE.

「FORCE_FACE」は追従荷重に対応していないのかもしれませんので、その代わりに「PRES_REP(圧力)」を用いることにしました。そのため、平板の上面端にこのようなグループが必要になった訳です。

これまでグループはジオメトリに依存したものばかり作成してきましたが、今回は要素面をグループ化して名前をつけることになります。これは下図のように、メニューから「Mesh」-「Create Group」でMeshとNameを指定し、Elements typeをFace、Group typeをStandalone groupにして、該当する5つの要素面を選択します。複数選択は、シフトキーを押しながら各面を順にクリックしていくことで出来ます。選択し終えたら、Contntで「ADD」を押し、「Apply and Close」でグループが作成出来ます。

Folf0203

Eficasによるcommファイルの編集ですが、これも「超弾性材料の検証」のcommファイルを修正して作成していきました。

DEFI_MATERIAUですが、C10=0.66666, C01=0.16666, NU=0.499にしています。「超弾性材料の検証」より硬めの材料に変更しました。

Folf0204


境界条件は、fixX面をFACE_IMPOで完全拘束(DX=0, DY=0, DZ=0)。
loadP面にPRES_REPで0.05の圧力を与えました。
解析時間はtime=0~1とし、計算の時間刻み(DEFI_LIST_REELの"pas")は2000分割としました。

Folf0205


計算の刻み幅が細かすぎると思うのですが、後述する通り収束の許容値をかなり緩くしても、2000分割まで細かくしないと収束してくれませんでした。引張りでは50分割程度で済んだのですが、曲げになると何故このように収束性が悪化するのか不明です。

刻み幅が細かすぎるため、結果出力は100ステップ毎(つまり20分割)にしています。これは、DEFI_LIST_REELの"output"で指定しています。
(出力間隔のコントロールは前回の記事にて説明しています。)

今回のポイントである追従荷重の指定ですが、下図のとおり、、「STAT_NON_LINE」-「EXCIT」で、該当する荷重に「TYPE_CHARGE "SUIV"」を追加します。
これは、Code_Asterのマニュアル「U4.51.03」(3.2.3節)で解説されています。

Folf0206

計算の収束許容値ですが、残差力の相対値での判定ではなかなか収束しなかったので、残差力の絶対値(REGI_GLOB_MAXI)で、0.001という値にしました。これはかなり緩い値なのですが、それでも収束しにくく、試行には苦労しました。

結果出力は、各RESUのb_extracの下に「LIST_INST output」を指定し、出力を20分割としました。それ以外は特に変えていません。

Folf0207

さて、計算結果ですが、まず追従荷重の変形図です。

Folf0208

圧力がゴム板の変形に追従するため、曲がり角度が90度を超えて、荷重面が下向きになっていることが分かります。

次に、「TYPE_CHARGE "SUIV"」 を消去して、同じ大きさの通常荷重にした場合の結果です。

Folf0209

荷重面が上向きのままです。
このように、追従荷重の効果が確認できました。

変形の動画です(追従荷重のみ)

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

folf02.zip

なお、今回は2000ステップもの細かい時間刻みで計算したため、六面体2次要素がわずか200要素だけにもかかわらず、20分ぐらいの計算時間がかかってしまいました。実際は試行錯誤で何度か計算し直したので、かなりの時間を費やしました。商用ソフトではあまり何も考えなくても、デフォルトの設定で同様の計算が20ステップ、約40秒で出来ます。

ここまでCode_Asterで接触、摩擦、塑性、超弾性等などの非線形解析を試してみての感想としては、「収束しにくい!」ということです。まだ検証していない設定項目が色々あるので、それらを問題に応じて詰めていく必要がありそうです。

結果出力の間隔設定

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

前回のゴム材料の引っ張り解析では、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年4月 | トップページ | 2012年7月 »