有限要素法・導出 - 超音波流体屋のプログラム備忘録

テンソル演算の嵐なのであまり入門用でないです。
あまり参考にしないように。

微分方程式

有限要素法使って解きたい支配方程式ってだいたいこの辺だろうか。
このページでは3式並列で離散化する。
物理的意味ではなく数理的意味が把握できたら他の微分方程式も応用が効くんじゃないかな。

弾性運動方程式: $\rho \ddot{\bf u} = {\rm div}{\bf \sigma} + \rho {\bf f} $
ポアソン方程式: $\nabla^2 \varphi = - \varrho/\varepsilon $
波動方程式: $\nabla^2 \phi - (1/c^2)\ddot{\phi} = q$

重み付け(仮想仕事の式)

上の式で厄介なのは、空間微分$\nabla$が関わった項。
これを処理する上での有限要素法のトリックは、両辺によくわからない重みを掛けることから始まります。

重みは本来なんでもいいのですが、二階微分の根っこにいる物理量に似たものを掛けるのがガレルキン法というやりかたで、
式三本の弾性運動方程式は内積を取る形で重み$\delta {\bf u}$を、その他には$\delta \varphi, \delta \phi$を掛けます。

弾: $\rho \delta{\bf u} \cdot \ddot{\bf u} = \delta{\bf u} \cdot {\rm div}{\bf \sigma} + \rho \delta{\bf u} \cdot {\bf f} $
ポ: $\delta \varphi \nabla^2 \varphi = - \delta \varphi \varrho/\varepsilon $
波: $\delta \phi \nabla^2 \phi - \delta \phi (\ddot{\phi}/c^2) = \delta \phi q$

重みへの要請は下記の境界条件をみたすことだけ。

1. 不連続でないこと (ディリクレ条件, Dirichlet)
2. 微分可能であること (ノイマン条件, Neumann)

重みに物理的意味を求めることはあまり得策ではない模様ですが、
弾性運動方程式で変位$\delta {\bf u}$を掛けると、全体がエネルギーの次元になるので、仮想仕事の式とか呼ばれるようです。

積分形に戻す

全体に体積分をかけて積分形に戻します。
ここで、divのある項は発散定理が使えるので面積分になります。
重みとdivが混じっている項の処理はGreenの定理(要は部分積分)てのがあってこんな感じになります。
$\iiint w {\rm div} {\bf A} dV = \iint w {\bf A} d{\bf S} - \iiint \nabla w \cdot {\bf A} dV $

弾: $\iiint \left[ \rho \delta{\bf u} \cdot \ddot{\bf u} + \nabla \delta{\bf u} : {\bf \sigma} \right] dV =\iint \delta{\bf u} \cdot {\bf \sigma} d{\bf S} + \iiint \rho \delta{\bf u} \cdot {\bf f} dV $
ポ: $\iiint \nabla \delta \varphi \cdot \nabla \varphi dV = \iint \delta \varphi \nabla \varphi d{\bf S} + \iiint \delta \varphi \varrho/\varepsilon dV $
波: $\iiint \left[\nabla \delta \phi \cdot \nabla \phi + \delta \phi (\ddot{\phi}/c^2) \right] dV = \iint \delta \phi \nabla \phi d{\bf S} - \iiint \delta \phi q dV$

式から二階微分が消えたことに注目されたい。

基礎積分方程式

物理的に意味を持っているものを別のパラメータに書き換える。
仮想変位の勾配→仮想ひずみ $\nabla \delta{\bf u} : {\bf \sigma} = \delta{\bf \epsilon} : {\bf \sigma} = \epsilon_{ij} \sigma_{ij}$
界面への外力と界面に垂直な電界・速度→${\bf t}, E, v$

弾: $\iiint \left[ \rho \delta{\bf u} \cdot \ddot{\bf u} + \delta{\bf \epsilon} : {\bf \sigma} \right] dV =\iint \delta{\bf u} \cdot {\bf t} dS + \iiint \rho \delta{\bf u} \cdot {\bf f} dV $
ポ: $\iiint \nabla \delta \varphi \cdot \nabla \varphi dV = \iint \delta \varphi E dS + \iiint \delta \varphi \varrho/\varepsilon dV $
波: $\iiint \left[\nabla \delta \phi \cdot \nabla \phi + \delta \phi (\ddot{\phi}/c^2) \right] dV = -\iint \delta \phi v dS - \iiint \delta \phi q dV$

形状関数

導入した重みである$\delta$のついた物理量を展開して、数値的に計算できるようにします。
ある要素内における物理量$f$は、周囲節点の物理量と形状関数$N_n({\bf x})$を用いて以下のように表せる。
$ f ={}^t\{N_n\} \{f_n\} $

微分演算では$f$じゃなくて、$N$を微分すればよく要素の形次第で初めから既知。
$ \nabla f = {}^t\{\nabla N_n\} \{f_n\} $

それで、上の基礎積分方程式を構成する大半の積分項はこうなる。$\Theta$は要素内で定数または時間微分演算子など。
$\iiint \nabla \delta f \cdot \nabla f dV = {}^t \lfloor \delta f_n \rfloor \left[ \iiint \{\nabla N_n\} {}^t\{\nabla N_n\} dV \right] \{f_n\} = {}^t \lfloor \delta f_n \rfloor [K] \{f_n\} dV $
$\iiint\delta f \cdot \Theta f dV ={}^t \lfloor \delta f_n \rfloor \left[ \iiint \{N_i\} {}^t\{N_n\} dV \right] \Theta \{f_n\} = {}^t \lfloor \delta f_n \rfloor [M] \Theta \{f_n\} dV $

外力項は形状関数に従って分配する。
$\iint \delta \varphi g dS = {}^t \lfloor \delta f_n \rfloor \left[ \iint \{N_n\} dS \right] g = \{g_n\} dS $
$\iiint \delta f h dV = {}^t \lfloor \delta f_n \rfloor \left[ \iiint \{N_n\} dV \right] h = \{h_n\} dV $

$\delta$のついた物理量、本来はあまり意味のない重み、にかかわらず上の基礎積分(恒等)方程式が成り立つためには、
$\delta$以外のところで方程式が成り立つ必要がある。
つまり、頭についている重みは外して良い。これらを基礎積分方程式に反映すると以下の式となる。

行列$[M], [K]$、弾性運動方程式用の$[{\cal M}], [{\cal K}]$については、
積分前係数行列の導出および各要素の形状関数についての部分にて。

有限要素方程式

一つの要素内においては、以下の話にまとめることができる。

弾: $ \rho [{\cal M}] \{\ddot{\bf u}_n\} + [{\cal K}] \{{\bf u}_n\} = \{t_n\} \frac{dS}{dV} + \{f_n\} $

ポ: $ [K] \{\varphi_n\} = \{E_n\} \frac{dS}{dV} + \{\varrho_n/\varepsilon\} $

波: $ [K] \{\phi_n\} + [M] \{\ddot{\phi}_n\}/c^2 = - \{v_n\} \frac{dS}{dV} - \{q_n\}$

次のステップ

• 積分前係数行列の導出
• 三角形要素・四面体要素
• 四角形要素・六面体要素
• 有限要素法・アルゴリズム (個々の要素行列を全体行列に統合するなど)