トピックス プログラム、シミュレーション 科学・技術・数学

有限要素法の中身

投稿日:

有限要素法(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)は、右辺の行列の青枠の行と列を掛けた和なので、k×u1+(ーk)×u2+0×u3+0×u4-k2u2となります。残りの行列は赤枠の部分になります。

p.8 行列の式を「ガウスの前進消去法」で手順(1)~(5)で計算していきます。青枠内の(0 0 P )を行列の中に入れてた後、対角線を「1」、対角線の下半分を「0」にしようとしています。

p.9 元に戻して、u4u3に、u3u2の式に代入していきます。この操作のことを「後退代入」と呼びます。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がどんなことをしているか、流れは掴んでいただけましたでしょうか? 転置行列の処理辺りは、私もまだ十分にできていませんので、わかり次第またアップします。

-トピックス, プログラム、シミュレーション, 科学・技術・数学

Copyright© 進化するガラクタ , 2019 All Rights Reserved Powered by STINGER.