四角形要素・六面体要素 - 超音波流体屋のプログラム備忘録

形状関数 $N= ^t \{ N_1, N_2, \cdots, N_n \}$の中身、その２。高次要素はやりません。

四角形要素・六面体要素は三角形・四面体のように腕ずくで形状関数を出すというよりかは、
四角形や六面体を予め形状関数が分かっている、(正規)正四角形・立方体にマッピングし直すことで形状関数を導出します。
最終的なK, M行列では体積分が必要であり、数値積分をする必要があるが座標変換のお陰でずいぶん単純になります。

座標変換

四角形や六面体を区間[-1, 1]の正方形・正六面体に座標変換することから始めます。
$\L (x,y,z) \rightarrow (\xi, \eta, \zeta)$

$\L x_i = N_n x_{ni}$
座標変換後の形状関数ははじめから用意されていて、
$\L N_n = \prod_k \frac{1 + S_{kn} \xi_k}{2} = \frac{1}{8} (1 + S_{\xi n} \xi)(1 + S_{\eta n} \eta)(1 + S_{\zeta n} \zeta) $
$\L S_{kn} = \left[ \begin{array}{cccccccc} -1 & 1 & 1 & -1& -1& 1& 1& -1 \\ -1 & -1 & 1 & 1& -1& -1& 1& 1 \\ -1 & -1 & -1 & -1& 1& 1& 1& 1 \end{array} \right]$
正規正六面体の各頂点で自分の物理量の影響が1、それ以外が0になることを確認すると、どうしてこんな形なのか納得できるかなと。

座標変換のヤコビアンが形状関数の微分に必要になります。
$\L J_{ij} = \frac{\partial x_i}{\partial \xi_j} = x_{ni} \frac{\partial N_n}{\partial \xi_j} $
$\L \frac{\partial N_n}{\partial \xi_j} =S_{nj} \prod_k \frac{1 + \bar{\delta}_{jk} S_{nk} \xi_k}{2} $
$\L \bar{\delta}_{ij} = 1-\delta_{ij} = int(i\neq j)$ はクロネッカーのデルタの補集合。

微分

$\L N_{n,i}$の算出には$\xi$の経由が必要です。

$\L N_{n,i} = \frac{\partial N_n}{\partial x_i} = \frac{\partial N_n}{\partial \xi_j} \frac{\partial \xi_j}{\partial x_i} = \frac{\partial N_n}{\partial \xi_j} J_{ij}^{-1} = \left[S_{nj} \prod_k \frac{1 + \bar{\delta}_{jk} S_{nk} \xi_k}{2} \right] \left[x_{ni} \prod_k \frac{1 + \bar{\delta}_{jk} S_{nk} \xi_k}{2} \right]^{-1} $

よって、
$\L [\tilde{K}_{nm}]_{(\xi_k)} = N_{m,i} N_{n,i} $

数値積分に必要なdVを座標変換すると、
$\L dV = |\det J_{(\xi_k)}| d\xi d\eta d\zeta $

数値積分

二次以下の関数$f(\xi)$を[-1 1]の領域で積分するには、ガウスルジャンドルの二点積分公式が有効で、
$\L \int_{-1}^1 f(\xi) d\xi \simeq f(-1/\sqrt{3}) + f(1/\sqrt{3}) $

$\L [K]_{nm} = \int [\tilde{K}]_{nm} dV = \sum_l \left. [\tilde{K}_{nm}] |\det J| \right|_{(\xi_k=S_{kl}/\sqrt{3})} $

質量行列のほうはまじめに数値積分についてはしてもいいけど、メモリ節約の面から、やっぱり各質量を対角要素に割り振って解くことの方が多い。
$\L [M]_{mn} = \delta_{mn} \frac{dV}{\sum_n 1}$
$\L dV = \sum_l |\det J|_{(\xi_k=\omega S_{kl})} $

弾性運動方程式の${\cal K}, {\cal M}$は $N_{n,i}, M_{mn}$から積分前係数行列の導出に従い代入する。