有限要素解析の初歩的な質問です

質問

私は学校で有限要素解析を用いた弾塑性特性を示す材料の応力解析を行いたいと思っています。現在、ABAQUS CAEを用いて材力計算の信頼性を図るために、十分に大きい円孔板を左右に強制変位を与え（つまり引っ張ったとき）、その時の円孔周りの応力分布が理論通りになるかを試しています。先日、自分で条件を与え、実際に計算させてみたところ、応力分布形状は、大体あっているようですが、引張り軸と垂直方向の円孔周りの応力値が引張り応力の3.00倍になるという理論値には遠く及ばない2倍くらいになってしまいました。
ちなみに、境界条件の与え方などは間違っていないと思いますが、メッシュの形状が左右対称ではありません。歪んでいます。

さらに、板を1000mmに円孔を1mmにしたりと色々試してみましたが、上手くいきませんでした。
どなたか、この分野に詳しい方、アドバイスをいただけると助かります。
よろしくお願いします。

回答No.6

noname#221368

＞もっと弾性係数が大きい材料を仮定して計算を行った方がいい？
　そんな事はありません。私が言いたかったのは、弾性係数によって応力分布や応力値が変わってくる可能性もあるので、目的の物に合わせた方が安全だという事でした。

＞URLの4.1を見ると何をやりたいのかわかると・・・
　とりあえず読んでみます。

＞平面応力/ひずみ要素厚(Plane stress/strain thickness）を入力するoption
　想像するに、剛性マトリックを計算する際に使用する、要素寸法の一部と思われますが（平面応力にしても平面ひずみにしても、あっておかしくないoptionです）、この前シェル厚というのが出てきました。これとの関係は調べたほうがいいと思います。シェル厚も上記の意味に取れるからです。ソフトの用語は、しょせんそのソフトだけの方言です。

＞エラーについて
　こればっかりは・・・。よく調べてがんばって下さいとしか言えません。がんばって下さい。

回答No.5

noname#221368

補足情報をもとに、知人に訊いてみました。ABAQUSクラスの解析ソフトになると、ベテランでもマニュアルを引きながら使用するのが普通です。知人も、手元にマニュアルがないのでわからない、と言っていた部分もあります。大きなお世話ですが以下、補足情報の各項目についてコメントします。予想の部分には(？)，気になった部分には(！)を付しています。

(1)メッシュ
【Element Shape】Quad：四角形要素
【Tequnique】Free　or　Atructured：わかりませんでした
【Algorithm】Medial axis　or　Advancing frontで、Medial axisを選択(！)
　メッシュが左右非対称になるのは、Medial axis機能がうまく働いていないんだと思います（これにはコンピューターの桁落ち誤差が関係します）。自分もCosmos／Ｍの領域自動メッシュで、似たような目にあいました。最初に1/4モデルを作成し、1/2モデルへ左右折り返しコピーすれば、これは防げると思います。Advancing frontもメッシュ生成ロジックが明快なので、個人的には嫌いではありません。
　Minimize the mesh transition：要素寸法変化を最小化して滑らかにする(？)。付けて正解と思います。
【Element Type】Explicit　or　Standardで、Explicitを選択(？)
　Explicitを選べば、次の【Geometri order】を有効になり、Standardなら内部で勝手にやる、という事かな？と思いました。
【Geometri order】Linear：形状折れ線近似
【Element control】Membrane strain：Finite，Second-order accuracy：No，Hourglass control：Stiffness，Displacement hourglass factor：１，Rotational hourglass scaling factor：１，Out-of-plane displacement hourglass scaling factor：１
　デフォルト値で良いと思います。
【Global seeds】Sizing controls：Approzimate global size：５，Curvecontrol：Deviation factor(h/l)=0.1
　純粋にメッシュ寸法と形状と配置に関する事なので、細かければ応力オーダーに大きな影響は出ないと思います。
【Seed Edge：Biased】円孔付近でメッシュが細かくなるように下底辺において、Biased ratio：10.0，Number of element along the edge：50
　同上。
【Seed Edge：By Number】円孔に沿って、Number of element along the edge：20
　同上。

(2)Shell planar
【断面積分】を解析中と解析前にできる条件を選択：わかりませんでした
【基本Basic】
シェル厚を5(！)：
　知人に確認した結果、シェル厚は板の厚さの寸法だとわかりました。解析対象は厚2mmだったはずです。これでは応力値が1/2.5になります。もっとも結果を応力比で見ているので、問題ないのかもしれませんが、念のため。

弾性係数1.9GPa(！)：
　ポアソン比0.3なので最初鋼材かな？と思いましたが、鋼材は約Ｅ＝200GPaです。2GPa程度となると、概ね鋼材より100倍，コンクリート（約20GPa）より10倍変形しやすい材料です。私は土木出身なので、これら以外の材料をあまり知りませんが（木材ですか？）、弾性係数が違っているという事はないでしょうか？。
　というのはモデルが完全な矩形でないために、このFEMは高次内部不静定問題を解いている事になって、弾性係数によって応力値が変わってきます。これも結果を応力比で見ているので、問題ないのかもしれませんが、念のため。

ポアソン比0.3(！)：
川田技報には弾性係数は載っていませんが、ポアソン比は0です。もしこの技報と比較しているのであれば、ポアソン比によって応力比は確実に異なります。影響度合いは計算しないと不明ですが・・・。

厚さ方向の積分方法：Gauss　or　Simpsonで、Gaussを選択、積分点の数は3：
　標準的と思います。

【高度Basic】
Section Poisson's ratio：use analysis default：わかりませんでした
Thickness modulus：無し：わかりませんでした
Transverse shear：無し：わかりませんでした

(3)その他(！)
　planarには、平面歪／平面応力の選択Optionはないでしょうか？。平面歪と平面応力は全くの別物です。Shellなので当然平面応力と思いがちですが、特殊な使用をするUserのために、このOptionが存在するかもしれません。またShellだから当然平面応力がDefaultだと思いがちですが、たまにソフトのバグで、そうでない場合もあります。

＞引張り軸と垂直方向の円孔周りの応力値が引張り応力の3.00倍になるという理論値には遠く及ばない2倍くらいに・・・
　とありますが、「円孔周りの応力」とは正確には何でしょうか？。例えば、円孔辺上の応力を取ったのであれば、理論的には円孔辺上の応力は0なので（材料表面ですから）、比が大きくなるのは当然です。FEMはあくまで近似計算なので、円孔辺上の応力は0になりません。また引張り応力は、モデル全体で一定ではないはずです。完全な矩形板の場合の引張り応力ですか？。

　シェル厚を5mmとされていると思いますが、これと1.9GPaの単位の整合性はありますか？。そうでないのであれば、少なくとも変形量に関して100万倍の違いがあります。これも結果を応力比で見ているので、問題ないのかもしれませんが、コンピューターの有効桁数は有限です。というわけで念のため。

　Upされた応力図は、本当に引張り軸と垂直方向の応力でしょうか？（怒らないで下さいね）。図の解像度で文字が見えなかったのです。

補足

各項目へのご回答ありがとうございました。
とても助かります。

メッシュについては、ddtddtddtさんのアドバイスをヒントに節点の数を増やし、コピーすることでダイブ左右対称に近づきました。
また、ご指摘いただいた弾性係数ですが、私が最終的に計算したい対象物は建築用の壁材です。擬弾塑性特性を持つセラミックス材料なので鉄鋼材料より、はるかに弾性係数は低くなっています。ご指摘前はあまり気にしませんでしたが、もっと弾性係数が大きい材料を仮定して計算を行った方がいいのでしょうか？

また、平面応力、平面歪についての選択optionですが、一通り探してみた結果、よく分かりませんでした。ただ、平面応力/ひずみ要素厚(Plane stress/strain thickness）を入力するoptionはありました。これのことでしょうか？

また、円孔周りの応力分布についてですが、円孔辺上ではなく、正しくは円孔周りの応力分布です。↓URLの4.1を見ると何をやりたいのかわかると思います。
http://www.aero.osakafu-u.ac.jp/as/lab3/design/2001/03.pdf

川田技報の例は、この問題を有限要素で計算していたので、ご紹介いたしました。

ただ、シェル要素で設定を先日書いた条件の通りに計算してみましたら、エラーが出てしまい、Visualization Moduleで、There is no valid step data available on the database. If the analysis is running, the database must be closed and reopened once the results have been initialized. The requested operation has been canceled.
となってしまいました。

どうしてよいのか、悩みどころです・・・

回答No.4

noname#221368

私事ですが、一昨日、とりあえず納品が終わったので、またここに戻って来ました。
　ABAQUSといえども何か癖を持つかもしれないので、ABAQUSを使用していた知人に訊いてみようと思います。そこで、もしよろしければですが、追加情報をお願いします。

(1)メッシュ
　図を拝見していて気になったのですが、変形が左右対称でない事です。これはメッシュの影響でしょうか？。円孔付近の変形状態を見ると、左の辺を固定して、右に引っ張ったようにも見えます。

(2)支持条件，荷重条件
　以下、左辺を固定せずに左右に引っ張ったという前提で書きます。文面では「左右に強制変位を与え・・・」とありますが、これは正確には、
　1)左右の辺に力を与えて、その結果0.1mmの変位が生じた．
　2)左右の辺に、直接0.1mmの変位を与えた．
のどちらでしょうか？。1)と2)では状況が違います。もし1)のケースなら、支持条件不足という事になります。そのために左右対称な変形にならず、応力精度も悪化したという事も考えられます。

(3)Shell planar
　私の知っている面内応力機能付のShell要素は、以下のものです。同じものでしょうか？。面外曲げ専用の要素と、面内応力を平面応力または平面歪状態で扱う要素を合体させたもので、両者の間に連成はない。パラメータとして、以下の設定が考えられる。

(1)面外曲げのパラメータ
1)弾性係数，ポアソン比，板厚．または、それらを集約し、
2)板剛度で与える．
3)ただし今回のケースでは、面外曲げは生じない．

(2)面内応力のパラメータ
1)弾性係数，せん断弾性係数，ポアソン比，板厚．
2)平面応力と平面歪の切り替えフラグ．ただし最初からどっちかに固定かもしれない．
3)また、せん断弾性係数は、弾性係数（ヤング率）とポアソン比からの自動計算かもしれない．

いちおうこんな処ですが、問題が問題だけに、身近な教官に聞くのが一番早いかもしれないとは思います。

お礼

ddtddtddtさん、毎度毎度のご回答ありがとうございます。
さて、早速ですが追加の情報を付け足します。

(1)メッシュ
メッシュ形状に関してですが、左右対称でないのは引っ張った影響ではなく、メッシュ作製のときからです。参考の為に、解析前におけるメッシュの図を以前示したページに添付しておきました。http://mu4.mp7.jp/ib/?uid=satoyu0117&dir=82​
またこのメッシュを作ったときの条件を下記に示しました。おそらく不必要な情報だとは思いますが・・・

【Element Shape】はQuadであり、【Tequnique】はFreeとAtructuredが選べ、Freeを選びました。【Algorithm】はMedial axisとAdvancing frontが選べ、Medial axisを選びました。更にMinimize the mesh transitionにチェックを付けました。
ただし、これらは全てデフォルト通りです。
【Element Type】のElement LibraryはExplicitとStandardが選べ、Explicitを選択し、【Geometri order】はLinearを選択しました。
【Element control】はMembrane strainがFiniteで、Second-order accuracyがNo、Hourglass controlがStiffness、Displacement hourglass factorが１、Rotational hourglass scaling factorが１、Out-of-plane displacement hourglass scaling factorも１というデフォルトの値を入れました。

【Global seeds】において、Sizing controlsのApprozimate global sizeを５として、CurvecontrolはDeviation factor(h/l)=0.1としました。
また、円孔付近でメッシュが細かくなるようにと、図の下底辺において、【Seed Edge：Biased】をBiased ratioが10.0、Number of element along the edgeが50としました。
また、円孔に沿って【Seed Edge：By Number】のNumber of element along the edge を20としました。
以上がこの図におけるメッシュ作製のときに入れた条件全てです。

長くなってしまったので、分けて書き込ませていただきます。

補足

(1)の続きです。上下ぎゃくになってしまいましたが。。。

(2)支持条件、荷重条件
曖昧な表現を使ってしまい申し訳ございませんでした。
私が使った条件は、(2)の『左右の辺に直接0.1mmの変位を与えた』です。
解析結果の図が左右対称でないのは、上で述べたとおり解析前のメッシュが歪んでいたからです。

(3)Shell planar
すみません。FEM初心者なので分かりませんという言い訳がましくなって、大変恐縮なのですが、ABAQUSに面内応力、面外曲げ応力のパラメータについて詳細に設定できることは初めて知りました。
とりあえず、調べてみまして、要素特性のことかなと思いました。均質の連続体シェルの場合における断面積分を解析中と解析前にできる条件で、基本Basicとしてシェル厚を5、材料特性の弾性係数1.9GPa、ポアソン比0.3、厚さ方向の積分方法はGaussとSimpsonが選べ、Gaussを選び、厚さ方向の積分点の数は3としました。
続いて、高度Basicとして、Section Poisson's ratioはuse analysis defaultを用いました。また、Thickness modulusとTransverse shearの有りor無しも選べましたが、defaultどおり無しにしました。

以上が、補足の追加情報です。
お時間があるときでよろしいので、よろしくお願いします。

回答No.3

noname#221368

低減積分は、まともなモデルを使用していても起こる、不運な異常事態に対処するためにあります。今回のように、見た瞬間問題ありとわかるモデル化に対して（言葉は悪いですが、すいません）、用いる事ができたとしても、その結果はやはり信頼に値しないと思います。
　本当にShear Locking が原因かどうかを確かめるなら、要素を深さ方向に2層以上にする事です。それが根本的な確認方法です。自分の問題だったらそうしたと思いますが、しかし多層にするのは手がかかるので、今回は安易に薦めたくありません。
　じつは、もう一つ疑問があるのです。縦横100mmで厚2mの板ですから、十分に薄板の部類に入ると思います。なのに何故、六面体要素（Solid要素）なのですか？。前半で述べたように薄板に対してSolid要素を用いると、要素の細長比を正常に保つために、平面的には非常に細かいメッシュが必要となり、また深さ方向も2層以上の分割が必要と思えます。こういった面倒さを避けるために、膜応力を考慮できるSheLL要素や板曲げ要素が用意されていると思うのです（塑性範囲もカバーしてるはずです）。そういった要素を用いれば、深さ方向には1層で、平面的には20×20程度の要素数で、十分な結果を出せる気がします。特別な事情がない限り、薄板の扱いは、そういうものだと思っていました。よって、問題はこうなります。

(1) 薄板に対して六面体要素（Solid要素）を使わなければならない、特殊事情があるのか？
(2) たんにメッシュの詳しさと六面体要素の性能評価をしたいだけなら、例題が不適切ではないのか？

です。どうでしょうか？。

お礼

ご回答ありがとうございました。

ddtddtddtさんのご指摘どおり、この例題においてSolid要素を用いる必要性はないと思いました。なぜSolid要素を使ったかと言われれば、あまり深く考えていなかったとしか言えないのですが、私が最終的に解析を行いたいものが釘の刺さっている壁材なので、あまり考えずにSolid要素を用いて解析を行っていました。釘周りにおける応力解析が主ですので、例題も円孔板が適切であると考えました。また、メッシュの詳しさや六面体要素の性能評価をするというよりは、ちゃんと有限要素の計算ができる(扱える)こと。そして、有限要素の計算が理論値とどこまであっているかで、未知の材料の力学特性を有限要素で解析する際の下準備として知りたいというのが目的でした。（すみません。説明不足で・・・）

さて、ddtddtddtのご指摘どおり、薄板（Shell planar 200×200）における解析をQuad Elementを用いて行ってみました。その結果とプログラムの一部を下記URLにアップロードしておきました。

http://mu4.mp7.jp/ib/?uid=satoyu0117&dir=82

解析図の境界条件は、200×200mmの薄板1/2モデルにおける弾性率1.9GPaの円孔の開いた弾性体に対し、横方向(x方向)に0.1mmずつ引っ張った図です。さらに円孔側の下面はx方向のみの変位のみを許す、としました。

しかし、結果はSolid要素のときと変わりなく、精度は20パーセントくらい悪い結果となりました。
やはり、私の解析手法がそこまで悪いのでしょうか・・・
それともShearlockingの影響以外の影響が考えられるのでしょうか？
疑問は深まるばかりです。

何かお気づきの点があればご回答お願いします。

回答No.2

noname#221368

参考URLを読みました。参考はハイアラーキ要素なのですね。FEMから離れて約2年になります。少なくとも2年前の某大学のABAQUSには、ハイアラーキ要素は搭載されてなかったように記憶してます。
　次に、要素数10527とすれば、おおよそ1×1×2mmの寸法になるので、メッシュは厚さ方向に一層なのかな？、と想像しました。つまり、以下は次が前提です。

　　(1)使用要素は、節点変位のみで内部の歪み場が決まる、通常の六面体要素。
　　（六面体要素なのに6自由度？、というのは一寸気になりましたが）
　　(2)メッシュは厚さ方向に一層である。

　これらが違ってたら、御免なさい。

　前提が正しいとして結論を言うと、Shear Lockingでは？、と思いました。荷重条件は等方的でないので、円孔付近ではせん断歪みが発生すると思えます。ところが(2)だとすると、厚さ方向には境界条件も荷重条件も対称なので、板の表裏で同じ節点変位になる気がします。そうすると(1)より、せん断歪みは生じず、Shear Lockingが起こり、せん断変形に対して要素は剛体的になって、結果として実際よりも硬い板になってしまった。そのつけが引張応力に現れたのではないか、と思いました。
　ちなみに参考URLでは、せん断変形モードも独立な一般化変位にできるハイアラーキ要素だから、うまくいったと思えます。

お礼

前回に続き、ありがとうございました。

ddtddtddtさんの前提というか仮定は言われたとおりです。メッシュは厚さ方向に一層で通常の六面体要素を使っています。
shearlockingが原因で解析が上手くいかないのでは。。。とのご指摘ですが、確かに納得することはできたのですが、では、どうすればこのshearlockingの影響をクリアできるのでしょうか？拘束節点を減らした低減積分法というものがあるみたいですが、これはこれでアワーグラス問題などもあり、上手くいかなそうです。先に進まずここにやってきましたが、やはり、有限要素法は経験工学的といわれるだけあって素人には難しいですね・・・

機械系の学科ではないので、周りに聞ける人がいないので、かなりしんどいです。また、ちょっと頑張ってみます。

ご指摘などありましたら、またお願いします。
よろしくお願いします。

回答No.1

noname#221368

ABAQUS は有名な解析ソフトなので、そんなに間違った結果を出すとも思えませんが、状況がわからないため、以下は怒らないで聞いて下さい。とりあえず言える事は、

(1) 理論値は正しいか？。境界条件など解析と違っていないか？
(2) 弾塑性解析の増分量は適当か？
(3) 1.5倍が3倍になってもおかしくないくらい、メッシュは粗くないか？
(4) 何故メッシュを左右対称にしないのか？。歪んでていいのか？

などです。弾塑性解析に代表される非線形計算は、けっこう微妙なものです。

補足

ddtddtddtさん、ご回答ありがとうございます。
質問の内容が十分でなく、すみませんでした。

質問について補足させていただきますと、理論解は弾性論であるので、解析は弾性領域、すなわち降伏応力以下での応力値で行いました。(弾完全塑性体なので、荷重量が小さければ、弾性解が出せます)

ですので、理論値は弾性解ですが有名な無限円孔板の問題ですので、(2)のひずみ増分理論は適当かというご指摘も適応外なので、今の時点では考えなくて良さそうです。（１）のご指摘も近似解となってしまうものの、かなりの精度で理論と等しくなることが分かっており、問題ないと思われます。(参照URL http://www.kawada.co.jp/technology/gihou/pdf/vol21/21_ronbun_01.pdf#search='有限要素%20無限円孔板')
なお、弾塑性解の証明については無限円孔板の問題ができたら、内圧円筒問題で弾塑性解析例があるので、そちらで試すつもりです。

また、メッシュの切り方ですが、もちろんご指摘どおり、左右対称でなくてはならないと思っています。孔の開いていないモノにならば、ちゃんとしたメッシュを与えられるのですが、孔を開けた図形になるとどうも上手くいきません。シードもところどころにおいてるのですが・・・
ちなみに解析に用いたモデルの寸法は、100×100×2mmであり、接点数14664、要素数10527として、6自由度六面体エレメントを用いた解析図を用いましたので、メッシュは十分細かいと思います。（円孔の周りに近づくにつれて、特に細かくなるように設定している）

長々となってすみませんが、今の状況はこんなものです。すなわち、最終的には弾塑性解析にまで到達したいのですが、現在困っているのは弾性解析の部分です。弾塑性解析になると、さらに高い壁があると思われますが、まずは弾性解析をクリアしたいので、よろしくお願いします。