「有限要素法(FEM) その2」です。 本日の目標は、次の資料のp.2にあるような片持ち梁(片方が固定されている梁)の右端に力が掛かった場合の変形をFEMでシミュレーションします。 プールの飛び込み台の先端に選手が乗った状態をイメージください。 専用ソフトを使う前に、FEM処理の流れをざっと掴んでください。
資料はこちら → 有限要素法 その2
p.4 昨日の復習です。1つのばねの両端を引張ったときのつり合い式です。
p.5 3つのばねを直列につなげた場合を想定し、一番上が左のばね、2番目が真中、3番目が右のばねになります。各々行列で表したものが右の式です。Fは力、uは伸びた変位量、kがばね定数です。
p.6 3つの行列を足し合わせます。ばねとばねのつなぎ目を「節点」と呼びます。ばね1~3が要素となります。u1は固定されているので変位量0です。ばね3の右端を力Pで引張ります。
p.7 真中に書いてあるように、ばね1と2及びばね2と3の節点では左右の力が相殺されて、外力は0となり、節点4だけがPとなります。これらを左上の式に代入すると左下の式になります。F1(1)は、右辺の行列の青枠の行と列を掛けた和なので、k1×u1+(ーk2)×u2+0×u3+0×u4=-k2u2となります。残りの行列は赤枠の部分になります。
p.8 行列の式を「ガウスの前進消去法」で手順(1)~(5)で計算していきます。青枠内の(0 0 P )を行列の中に入れてた後、対角線を「1」、対角線の下半分を「0」にしようとしています。
p.9 元に戻して、u4をu3に、u3をu2の式に代入していきます。この操作のことを「後退代入」と呼びます。F1(1)=ーPとなりました。
p.11 2次元(平面)に拡張します。長さLで幅Dに板を右方向に引っ張ります。 長さ方向にΔL伸び、幅はΔD縮みます。応力をσ、ひずみをε、ヤング率をE、ポアソン比をνとすると右上の式で表されます。長さ方向(x方向)のひずみεと応力σ、幅方向(y方向)のひずみεと応力σの関係式は下の式になります。ヤング率Eと横弾性係数Gの関係式が右下にあります。横方向と幅方向の弾性係数の関係式になります。
p.12 x方向とy方向のひずみと応力の関係式、せん断ひずみとせん断応力の関係式を行列でまとめたものが真中の式です。弾性定数マトリックス(Dマトリックス)の中にヤング率やポアソン比が収まっています。
p.13 三角形の頂点はi、j及びkで各々が位置座標x、yと変位uが(xi、yi、ui)、(xj、yj、uj)及び(xk、yk、uk)と表し、行列で表すことができます。
p.14~15 複雑な変形をしていますが、最終的には、三角形要素eの変位は形状マトリックス(位置情報)と各節点の変位マトリックスの積で表されます。
p.16 ひずみεは形状マトリックス(Bマトリックス)と変位マトリックスの掛け算になります。
p.17 ひずみε=Bマトリックス(節点の形状情報)×変位、応力σ=Dマトリックス(要素の材料特性)×ひずみεの2つの式をまとめて、応力σ=Bマトリックス×Dマトリックス×変位 つまり、ひずみε=(節点の形状情報)×(要素の材料特性)×変位で表されます。
p.18 Kマトリックスという新しい係数です。応力に面積を掛けて力にするための係数になります。
p.19 B、D及びKマトリックスの関係式です。
p.20 右上のように三角要素が1、2及び3が繋がっていて、節4のところを力Fで引張ります。各々の三角要素について行列が書けます。1つの節においてxとy方向に変位し、また隣の節の影響を受けます。3つの行列を加算した行列が下の式です。
p.21 上の式を簡単に書くと、[K][u]=[F]。 この式から各節の変位[u]を算出し、次にひずみ[ε]を算出そして応力[σ]を求めます。
p.22 以上の流れを示してあります。
p.24 以下壁に左端が固定された「片持ち梁」の事例で計算してみます。片持ち梁の設計条件が中程に記載しました。
p.25 Dマトリックスです。
p.26 Bマトリックスとその転置行列です。
p.27 Kマトリックスの算出方法です。
p.28 Kマトリックスです。
p.29 上述のK1マトリックスは三角要素(1)の節点①、②及び⑧の状態を示しています。水色が節点①と⑧の関係、黄色の部分が②と⑧の関係を示しています。
p.30 三角要素(2)のK2マトリックスです。
p.31 三角要素(3)、(4)、(5)及び(6)のKマトリックスです。
p.32 上述のKマトリックスを加算していきます。
p.33 グレーの部分は節が繋がっていないので0です。
p.34 [K][u]=[F]の[F]の行列です。
p.35 節点①と⑧は拘束されているので、変位uはゼロで、働く力は未知です。
p.36 節点①と⑧の部分を除外したマトリックスでまず変位を計算します。
p.37 FEMの計算結果が下の図です。変位は少し誇張して拡大していますが、片持ち梁の先端ほど変形量は大きく、内部のストレスは固定端(左端)に行くほど大きくなります。(色を濃く着色しています)。
これだけの計算は手計算では無理で、プログラムを組んで計算する必要があります。自分でマクロを組んでも作れそうですが大変です。FEMがどんなことをしているか、流れは掴んでいただけましたでしょうか? 転置行列の処理辺りは、私もまだ十分にできていませんので、わかり次第またアップします。