4節点アイソパラメトリック要素のFEM例題 その1

今回は4節点アイソパラメトリック要素のFEMを使って以下の例題を解いてみます

E,νは要素内で一定かつ等方性とします

iso-parametric 4-node sample1.wxm
http://f.st-hatena.com/images/fotolife/r/ryooji_f/20160206/20160206224048.png
n : 節点数
NN : 形状関数
J : Nxのヤコビ行列
J_i : Jの逆行列(invert)
j : Nxのヤコビアン(det(J))
形状関数やヤコビアン等を一通り代入します(画面出力は省略)



B : Bマトリックス
D : Dマトリックス
Ke : 要素剛性マトリックス
Keの導出については4節点アイソパラメトリック要素 その4を参照ください(画面出力は省略)



幾何学境界条件を与えます
全体節点変位ベクトルUを%o22式に,未知節点変位ベクトルUaを%o24式にそれぞれ示します



力学的境界条件を与えます
全体節点荷重ベクトルFを%o29式に示します


http://f.st-hatena.com/images/fotolife/r/ryooji_f/20160207/20160207145349.png
F = K.U の連立方程式を,Uaについて解いた結果を%o33式に示します
(globalsolve:trueの宣言により,解は変数に代入されます)
ここでν= 0とした場合のu2(=u3)は%o34式となり,公称歪P/(E*H)をLで積分した値と一致します



%o33式の結果を使って未知のF(反力)を計算します(%o35〜38)
%o39式より,x方向の反力が外力(+P)と釣合っていることが解ります
また%o40式より,y方向の反力が釣合っていることが解ります



具体的に応力ベクトルを計算し,相当応力をプロットしてみます
寸法や材料定数など適当な値を代入します
応力ベクトルtを%o49式に示します
平面応力状態におけるMises相当応力が自然座標rの関数となっていることが解ります(%o50)



求めたMises相当応力を横軸にr1,縦軸にr2として%t51にプロットします