Introduction of FEM

有限要素法の概略

Last Update:2009年11月11日

概要

概要
有限要素法とは？
- メッシュを用いて空間的に離散化された場
  - メッシュの例
弱形式化された偏微分方程式
- ポアソン方程式の弱形式化の例
有限要素法離散化による弱形式偏微分方程式から連立一次方程式の作成
- ポアソン方程式で見る要素単位の積分
有限要素法解析の手続き
参考文献

有限要素法とは？

有限要素法とは，メッシュを用いて空間的に離散化された場の中から弱形式化された偏微分方程式の解を見つける方法です．

メッシュを用いて空間的に離散化された場

物理的な普通の場は無限の次元を持ちますが，コンピュータで無限次元の場は勿論表現できないので，計算機上では有限個の値で場を表現します．これが"有限"要素法と呼ばれる所以です．ではどのようにして，場を有限の変数を使って表現しているのか見てみましょう．

有限要素法による場のデータ表現(有限要素法補間場)はメッシュを用いています．メッシュは解析領域を単純な多角形で（２次元の場合は三角形や４角形，３次元の場合は４面体，６面体）分割したものです．有限要素法補間された場は節点配列と要素配列の２つをセットにして作られます．

節点とは空間上の一点で，その点における場の値を持ちます．節点配列は節点の座標と場の値(温度，変位など)が並んでいる配列です．

節点だけでは節点以外の場がどうなっているのかはわかりません．そこで要素を使います．要素はメッシュの中の１つの多角形で，その要素がどの節点を持っているかという情報を持っています（要は節点Indexの集まりです）．要素の中は補間関数を用いて補間されています．

隣り合う要素間で節点を共有しているので，場が不連続になることはありません( C^0 級)．メッシュは全ての解析領域を埋め尽くすように作られます．ですから，解析領域内の１点は必ずある要素かその境界上にあるはずです．点が含まれる要素内やその辺での補間を求めれば解析領域内の任意の点の値が求まります．要素内やその境界で連続な補間を用いているので，この場は連続であることがわかります．

メッシュの例

下のようなメッシュの要素配列は次のようになる．

要素配列

要素番号	節点Index
0	0 1 3
1	1 4 3
2	4 6 3
3	5 3 6
4	3 5 2

弱形式は偏微分方程式をもっと広い範囲で解が取れるようにした拡張です．強形式では２回が連続に微分可能( C^2 級)という解にかなり強い滑らかさが要求されましたが，弱形式では１階の弱微分が可積分( H^1 級)であればよく，かなり広い範囲で解を取ることができます．有限要素法補間された場は要素の接続部分で微分が不可能であることが多いです．例えば１次元の１時補間は折れ線(区分一次関数)になるので，そこでの微分は定義されません．しかし弱微分ならこのような関数も微分をすることができます．強形式の解が存在すれば，それは弱形式の解(弱解という)にもなっています．弱形式は強形式の自然な拡張になっています．

弱形式でもう一つ重要なことは，それが変分形式であり，積分の形で書かれているということです．積分は領域内での足し算のような演算です，解析領域が要素で充填されているので，各要素毎に積分値を足していけば領域全体の積分値を求めることができます．各要素は単純な形をしているので，要素内での積分は数値的,解析的方法を使って簡単にできます．積分ができるというのは，補間場を作るのと並んでメッシュの重要な機能です．

ポアソン方程式の弱形式化の例

ポアソン方程式を例にとって，弱形式化をやってみます．

$\{\begin{array}{l} -\nabla( \nabla \phi) = f\;\;in\;\Omega\\\ \phi=0\;\;on\;\partial\Omega \end{array}$

上の強形式に $\delta\phi\;\;(\delta\ph=0\;\;on\;\partial\Omega)$ を掛け，領域で積分した後にGreen-Gaussの公式を使うと以下のような弱形式に変換できる．

$\int_{\Omega} \nabla\delta\phi\cdot\nabla\phi d\Omega = \int_{\Omega}\delta\phi f d\Omega$

有限要素法の多くの問題は $\phi$ の属する空間に共役な空間の元である $\b{f}$ (Rieszの表現定理による)と双一次形式の作用素 $\b{a}$ を使って次のように書けます．

$\b{a}(\delta\phi,\phi)=\delta\phi\cdot f\;\;\;\forall \delta\phi$

これはポテンシャル

$\b{W}(\phi)=\frac{1}{2}a(\phi,\phi)+f\cdot\phi$

とした場合のポテンシャルの極値が解であることを意味しています．

$find\;\phi\;st\;\delta\b{W}=0$

はノルムとしての性質を持っていてこれをエネルギーノルムといいます．さらに，作用素の強圧性や有界性などの性質から，弱解は唯一存在することがわかります．よって解はポテンシャル最小値であることがわかります．

$find\;\phi\;st\;minimize\; W$

このポテンシャルは物理的な問題を解くと大抵の場合エネルギーと同じになります．例えば固体力学解析におけるポテンシャルは歪エネルギーを表します．電磁気におけるポテンシャルも電磁場のエネルギーを表します．ポテンシャルが最小であるということは，エネルギー最小化の原理と同じことです．

有限要素法離散化による弱形式偏微分方程式から連立一次方程式の作成

さて，次に弱形式化された偏微分方程式に有限要素法補間場を代入して連立一次方程式を作ってみましょう． $\phi,\delta\phi$ が次のように節点における値 $\phi^b,\delta\phi^a$ と場の基底関数 N^a,N^b との積の和で表されているとしましょう．

$\phi = N^b\phi^b,\;\;\;\delta\phi = N^a\phi^a$

これを前の式に代入すると解くべき連立一次方程式が得られます．

$a(N^a\delta\phi^a,N^b\phi^b)=(N^a\delta\phi^a)\cdot\b{f}\;\;\;\forall \delta\phi^a\\\Leftright \delta\phi^a a(N^a,N^b)\phi^b = \delta\phi^a(N^a\cdot\b{f})\;\;\;\forall \delta\phi^a\\\Leftright a(N^a,N^b)\phi^b = (N^a\cdot\b{f})\\\Leftright A^{ab}\phi^b = f^a\;\;\{A^{ab}=a(N^a,N^b),\;f^a=N^a\cdot\b{f}\}$

得られた連立一次方程式を解くことで解を得ます．実際の係数行列の成分 $A^{ab}=a(N^a,N^b)$ の計算は，要素ごとに行い，それらを足し合わせて作ります．この計算についてポアソン方程式の場合について見てみましょう．

ポアソン方程式で見る要素単位の積分

弱形式化されたポアソン方程式は次のとおりでした．

$\int_{\Omega} \nabla\delta\phi\cdot\nabla\phi d\Omega = \int_{\Omega}\delta\phi f d\Omega$

積分の計算は要素ごとに行った物の足し算です．

$\sum_e^n\int_{\Omega_e} \nabla\delta\phi\cdot\nabla\phi d\Omega = \sum_e^n\int_{\Omega_e}\delta\phi f d\Omega$

要素内で $\phi,\delta\phi$ が次のように節点における値 $\phi^q,\delta\phi^p$ と場の基底関数 N^q,N^p との積の和で表されているとしましょう．

$\phi = N^q\phi^q,\;\;\;\delta\phi = N^p\phi^p$

また，要素内の番号づけ $\{p,q\}$ は節点の番号づけ $\{a,b\}$ に対応しているとします．これを代入すると上の式は次のようになります．

$\sum_e^n\int_{\Omega_e} \nabla (N^p\delta\phi^p)\cdot \nabla (N^q\phi^q ) d\Omega = \sum_e^n\int_{\Omega_e}(N^p\delta\phi^p) f d\Omega\;\;\;\forall\delta\phi^p\\\Leftright \sum_e^n\delta\phi^p\int_{\Omega_e}\nabla N^p\cdot\nabla N^q d\Omega \phi^q = \sum_e^n\delta\phi^p\int_{\Omega_e}N^p f d\Omega\;\;\;\forall \delta\phi^p\\\Leftright\delta\phi^a\sum_e^n $A_e^{pq}$\phi^b=\delta\phi^a\sum_e^n f_e^p\phi^a\;\;\;\forall \delta\phi^p\\\Leftright\sum_e^n $A_e^{pq}$\phi^b=\sum_e^n f_e^p$

ここで

$A_e^{pq}=\int_{\Omega_e}\nabla N^p\cdot\nabla N^q d\Omega$

は要素剛性行列と呼びます．前の式と見比べると，全体剛性行列 $A^{ab}$ はこの足し合わせであることがわかります．足し合わせる際には要素内の節点の番号付けと節点番号の対応づけを使います．今，要素内の節点の番号づけ $\{p,q}$ が全体の節点の番号づけ $\{a,b\}$ に対応しているので， $A_e^{pq}$ は $A^{ab}$ に足し合わされます．このような足し合わせの処理をマージ(Marge)と呼びます．

弱形式化された方程式の解(弱解)は，ある空間(大抵の問題では境界条件を満たすような H^1 空間の部分空間)で，エネルギーノルムを最小にするような解を取ることがわかりました．

H^1 空間の中でも特にその部分空間である有限要素法補間場の中でエネルギーノルムを最小にするような解が有限要素法においての解です．

補間関数とテスト関数が同一の空間である場合は，これをGalerkin法と呼びます．Galerkin法ではエネルギーノルムでもっとも弱解に近いような有限要素法補間場で現される解でもあることがわかります．しかしながら流体のようにエネルギーが存在しない場合は，補間関数とテスト関数を別の空間から選ぶことでより精度の良い解が得られることもあります．

有限要素法解析の手続き

下に示すのは有限要素法解析の全体の流れです．

有限要素法の流れ

有限要素法の擬似コードは以下のようになります．

// 有限要素法の手順の擬似コード

①モデルを作る
②メッシュを生成する
③有限要素法補間場を作る
do i=1，要素数
　　④要素が参照する節点データを取ってくる
　　⑤弱形式化された支配方程式を要素あたりに積分する
　　⑥積分して得られた要素あたりの剛性行列と残差ベクトルを全体へ足し合わせる．
end do
⑦全体剛性行列と残差ベクトルから成る連立一次方程式を解いて，更新ベクトルを得る
⑧更新ベクトルを使って解を更新する．
⑨場を可視化する

参考文献

有限要素法概説 -理工学における基礎と応用-

菊池文雄著