SonofSamlawのブログ

うひょひょ

有限要素法(FEM)とは

 FEMの世界を理解しよう!!!

 


構造解析、電磁界解析、流体解析などで使われる有限要素法の解説

 

書きかけです

 

  

 FEMに関しての誤解

 FEMを理解して、プログラムまでいくのは、容易なことではありません。いいかげんな回答が多すぎます。ある回答でいいかげんなものがありました。プログラムを使ったことも作ったこともない者のいうことはばかげている。

 

> 節点というのは、要素の頂点のことを指します。
> 四面体要素であれば、4つの節点があります。
> したがって、節点はトラスの接合部にあたります

 

などと言っている。要素の話になぜ「トラス」が出てくるのか? よくわかっている者にはわけわからない文章ですね。さらには、

 

> ”基礎的な有限要素法では、要素はトラスである”

 

とのことですが、「基礎的有限要素法」とは何か? 要素は何もしません。要素では、ただただ節点の値が補完され、それが積分される区域なのです。このようなことを言っている者がFEMの正しい理解をさまたげています。 要素はトラスでもラーメンでもない、節点の間の空間にすぎません。要素は物体ではありません。当然、モーメントの伝達などナンセンスです。要素という剛体があるわけではありません。

 私は十年以上もFEMの研究をしていたが、要素がトラスだのラーメンだの、モーメントだのという議論は聞いたことがない(名著培風館「有限要素法ハンドブック」やツェンキーヴィッツマトリックス有限要素法」参照)。ただ、「軸方向の力のみを伝達する要素」という特殊な要素(トラス要素)はありますね。

 

■多い誤解について

 わたしが昔、部外者と話しているときにわかったのですが、よくある誤解は2つです。

 各要素の中で微分方程式を解き、それらをつなげていく。

 

 各要素の連結で全体ができているとして、その結合部分のみに解くべき変数    が存在する。

 くべき変数とは、変位・応力・ひずみですね。上の例の誤解もこのへんからきています。「節点はトラスの接合部にあたります」などは、変形しない要素同士がどのようにつなげてあるか? という疑問から出てきたものでしょう。しかし、もしそうだとしたら節点とはなんなんでしょうか?要素と要素の接続部分が問題となるべきであって、節点などどうでもいいはずです。完全に剛体である要素がつながっていると考える。そこでそれらの接続の仕方が問題となるというわけです。

 

■正しくは?

 ここで正しいFEMの姿を説明しておきましょう。 上の間違ったイメージとはまったく違います。たぶん誰もが予想もしなかったものだと思います。

 対象の物体内に有限個の点を設定します。そしてその点、つまり節点のみでの関数の値を変数とし、他の点ではこれらの節点での変数の値から補間するのです。この同じ節点列から補間される区域が要素」なのです。つまり求める関数は連続である座標によるのではなく、有限な節点での変数の値の列によって決められます。関数を有限な点での変数の値によって決定する。その間の点の関数の値は要素の周囲の節点の関数の値より補完するのであります。でありますから、「要素同士の接続の仕方」などという問題はありえません。要素は単なる補間区域にすぎないのです。

 求める関数を有限な節点における値の列として、以下の変分式や残差式を解くのです。つまり、関数の関数であるものが、単に有限の数の変数の関数となってしまうのです。そこで最終的に、この変数で偏微分して、それが0とする変数の数の連立方程式になり、これを解くことにより、そして補間することにより、求める関数が求まるわけであります。



■まず変分法からの説明 

 FEMはなんといっても変分法を抜かしては考えられません。ここから理解していくのが筋でしょうね。

 関数の関数である汎関数Fの変分=0、つまり式の値を最小(停留)にする関数fを離散的な点で求める方法です。

 ある汎関数Fの変分=0、つまりFの値を最小にする(停留させる)関数fはオイラーの方程式と呼ばれる微分方程式の解となります。汎関数とは式の中にある変数が数値ではなく関数なのでして、定積分の形をしています。

 逆に、ある微分方程式があれば、変分=0となる関数fが、はじめの微分方程式の解であるような汎関数Fをつくれます。

 有限要素法では初めのの微分方程式(材料力学、流体、熱分布、電磁界などはすべて微分方程式の解として求められる)にたいおうする汎関数をつくり、変分=0とする関数fを求めようとします。

 このあたりのことは、変分法の本を見てください。

 さて、汎関数Fの変分=0とする関数fを、有限のサンプル点で求める方法が有限要素法です。汎関数の変分=0とするような関数を決定するために、関数fをサンプル点での値のみ未知なものとしてその他の点では周りのサンプル点の値から補間することにします。すると関数fはNをサンプル点数とすると、N個の変数によって表せますね。これは、一般に節点数と呼ばれています。そして、この節点の囲む領域を要素と呼びます。要素内の関数fの値は、周囲の節点のfの値から補間されます。

 すると、汎関数FはN個の変数(f1からfN)の普通の関数になってしまいます。

 つまり、汎関数Fの変分=0の条件は、有限要素法では関数F(f1、・・・、fN)とすると、さらにそれがf1からfNの2次しきであるなら、

   Fをf1で偏微分=0
       ・
       ・
       ・
   FをfNで偏微分=0

というN個の連立方程式になります。

このことについては、

http://note.chiebukuro.yahoo.co.jp/detail/n180985

参照。

 この連立方程式(たいてい数十万、数百万)を解く事ではじめの微分方程式の数値解を得ようとするものなのであります。

 現在ではガラーキン法など、上の説明と違う方法が一般的ですが・・・

■ より、おおざぱな説明です

 汎関数Fは解くべき微分方程式から推測して作ります。Fの変分が0と
なるような関数fが満たすべき微分方程式が、解くべき方程式となるよう
なものとしてです。変分とは、「関数による微分」を言います。これも
どうやるかは、「変分法」の本などみてください。
 
 次に、汎関数Fの変分=0となるような関数fを求めることになります。
まず関数fを空間の有限の点(サンプル点)の値(f1~fN)とその
間は、周囲のfiで補完することにします。
 そうすると、関数fは、f1~FNというN個の数列で表せることに
なります。これを汎関数Fに代入し、変分=0の条件を満たす、
f1~fNを見つけます。それは、上のとおりN個の連立方程式を解く
ことになります。
 この解f1~fNと補完されたものが解関数となります。

 

 以上のように解きたい方程式をEulerの方程式とする変分問題に直し,直接法によって解く方法をRayleigh-Ritz法とよぶ。

 


 例として弾性問題の微分方程式を書いておく。

 

   G*∇^2*u+G*∂e/∂x*(1+2ν/(1-2ν))+X=0      
      G*∇^2*v+G*∂e/∂y*(1+2ν/(1-2ν))+Y=0

      G*∇^2*w+G*∂e/∂z*(1+2ν/(1-2ν))+Z=0

 

   これは、各位置での変位(u,v,w)を未知関数とする方程式である。ここで、

   u,v,w:各点のx,y,z方向の変位

   X,Y,Z:外力のx,y,z成分

       e:体積ひずみ=∂u/∂x+∂v/∂y+∂w/∂z
     

       ν=εd/εl
    εd: 横ひずみ
    εl: 縦ひずみ
       E:: ヤング率

 

 ここで、未知変数はu,v,wである。上の微分方程式の解を求めるのが目標となる。ひずみ、応力などは求められたu,v,wから算出する。

 

垂直ひずみと変位の関係   
   
εx=∂u/∂x
   εy=∂v/∂y
   εz=∂w/∂z


せん断ひずみと変位の関係
   
γxy=∂u/∂y+∂v/∂x
   γyz=∂v/∂z+∂w/∂y
   γzx=∂w/∂x+∂u/∂z

 

により、変位からひずみが計算されます。このひずみより、次の関係により、

 

よ応力が計算されるというわけですね。

 

 

■しかし、

 変分法偏微分方程式の近似解法としては,理論的,数学的,物理的な背景が堅牢で理解しやすいのであるが,等価な変分問題を持つような微分方程式で無いと適用できない。つまり、問題となる微分方程式に対応する変分式がない場合である。元となる微分方程式に少しでも余分な項が入ったりすると、それに対応する変分式が見つけられなかったりする。そこででてくるのがRitz法、ガラーキン法である。

 


 http://nkl.cc.u-tokyo.ac.jp/11s/intro/CS02.pdf#search='Ritz法'

によれば、ガラーキン法が最も精度がよい– 汎関数がある問題については,変分法とガラーキン法は答えが一致する。 多くの商用コードでガラーキン法を使用。

 

■そこで重みつき残差法です。

 

 微分方程式F(f)=0において、

 

  I=∫F(f)dv=0

 

を考えてみる。この場合、I=0だからといって、F(f)=0になるとはかぎらない。ある部分では+、ある部分ではーであるので総和として0になっていることもあるから。そこで、重みつき、ということになる。これは、そのようなことが許されないようにしてしまうのである。F(f)=0でない限り、I=∫F(f)dv=0が成り立たないようにしてしまうのである。

 

 

 

 梁要素の剛性マトリックス

 http://d.hatena.ne.jp/ryooji_f/20110614/1308062579