積分前係数行列の導出 - 超音波流体屋のプログラム備忘録

前 : 有限要素法・導出 にて省略した質量行列の作り方とか。

積分前の係数行列

$ [M] =\left[ \iiint \{N_n\} {}^t\{N_n\} dV \right] = \iiint [\tilde{M}] dV $
$ [K] =\left[ \iiint \{\nabla N_n\} {}^t\{\nabla N_n\} dV \right] = \iiint [\tilde{K}] dV $
$ [{\cal M}] = \iiint [\tilde{\cal M}] dV $
$ [{\cal K}] = \iiint[\tilde{\cal K}] dV $
における、$[\tilde{M}], [\tilde{K}], [\tilde{\cal M}], [\tilde{\cal K}]$を扱います。
数値積分については各要素の解説の所にて。

添字

$(i, j, k, l, s, t) \in \{x, y, z\}$ : 空間次元に使います。
$1 \leq (n, m) \leq P$ : 形状関数の番号付けに使います。 $P$は要素のもつ接点数。

M行列

なにも問題なし…なんだけど数値積分がめんどくさいので、
普通は$dV$を均等に各接点対角要素に配るだけだけだったりします。

$ [\tilde{M}]_{mn} = N_m N_n $
$ [\tilde{M}] = \left[ \begin{array}{c} N_1 \\ \vdots \\ N_m \end{array} \right] [N_1 \cdots N_n] $

K行列

勾配同士の積が噛むので注意。ただし、$N_{n,x} = \partial N_n / \partial x$
三角形・四面体要素は微分すると形状関数が定数になるので、そのままdVを掛けるだけで$[K]$になります。

$ [\tilde{K}]_{mn} =　N_{m,i} N_{n,i} $

$ [\tilde{K}] = \left[ \begin{array}{ccc} N_{1,x} & N_{1,y} & N_{1,z} \\ & \vdots & \\ N_{m,x} & N_{m,y} & N_{m,z} \end{array} \right] \left[ \begin{array}{ccc} N_{1,x} & & N_{n,x} \\ N_{1,y} & \cdots & N_{n,y} \\ N_{1,z} & & N_{n,z} \end{array} \right] $

高階テンソルについて

弾性運動方程式の話しに行く前に、このページのテンソル表記について。

本導出では当面フォークト表記を使いません。
つまり応力テンソルを${}^t{\bf \sigma} = {}^t[\sigma_{xx} \cdots \tau_{xy}] $とベクトルで扱いません。
当然メモリを食いますが、応力を行列のまま扱えるのは回転やら大変形やらで旨みがあるので。

弾性定数行列 ${\cal C}$ は4階のテンソルで、等方性材料の場合ラメ定数$\lambda, \mu$を用いて以下のようになります。
$ \sigma = {\cal C} : \epsilon $
$ {\cal C}_{ijkl} = \lambda \delta_{ij}\delta_{kl} + \mu (\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) $
$\delta_{ij}$はクロネッカーのデルタ。

高階テンソルのイメージに…

あんまりよろしくないけど、縮退次元的な表現としてこんな表記をしてみます。
部分行列ではないので注意してください。ただし、$\kappa = \lambda + 2\mu $とします。

テンソル積$\sigma_{ij} = {\cal C}_{ijkl} \epsilon_{kl}$ のイメージが少し湧くのではないかなと。

剛性テンソル

テンソル演算と変位のインデックス$i$を保持した状態で、ひずみ・応力を導出する。
$ \epsilon_{kl} = {\cal D}_{kli} u_{i} $
$ {\cal D}_{kli} = \frac{1}{2} \left[ \delta_{ik} \frac{\partial }{\partial x_l} + \delta_{il} \frac{\partial }{\partial x_k} \right] $

$x^d = \frac{\partial }{\partial x}$

これを元に形状関数$ u_{i} = N_{n} u_{in} $を使って、書き直すと
$ \epsilon_{kl} = {\cal B}_{klin} u_{in} $
$ {\cal B}_{klin} = \frac{1}{2} \left[ \delta_{ik} \frac{\partial N_n}{\partial x_l} + \delta_{il} \frac{\partial N_n}{\partial x_k} \right] = \frac{1}{2} \left[ \delta_{ik} N_{n,l} + \delta_{il} N_{n,k} \right] $

先の${\cal C}$テンソルを掛けて、応力は
$ \sigma_{uv} = {\cal C}_{uvkl} \epsilon_{kl} = {\cal C}_{uvkl} {\cal B}_{klin} u_{in} $

仮想仕事式は$ \delta \epsilon : \sigma $よりテンソル表記は、
$ \delta　\epsilon_{uv} : \sigma_{uv} = \delta u_{jm} {\cal B}_{uvjm} {\cal C}_{uvkl} {\cal B}_{klin} u_{in} $

以上より、要素剛性テンソルは
$ [\tilde{\cal K}]_{ijmn} = {\cal B}_{uvjm} {\cal C}_{uvkl} {\cal B}_{klin} $

質量テンソル

形状関数も三階のテンソルを使って書き直すと、
$ u_k = {\cal N}_{kin} u_{in} $
$ {\cal N}_{kin} = \delta_{ik} N_n $

慣性項の仮想仕事式は以下のとおり。
$\delta u_j \rho \ddot{u}_i = \delta u_{jm} \tilde{\cal M}_{ijnm} \ddot{u}_{in} =$

$ \tilde{\cal M}_{ijnm} = \rho {\cal N}_{kjm} {\cal N}_{kin} = \rho \delta_{jk} N_m \delta_{ik} N_n = \rho \delta_{ij} M_{nm}$

行列に戻したいときは

$ I = i + 3n, u_I = u_{in} $など普通に添字の次元を落とせば良い。
全体行列に編みこむ時で十分だと思うが。