STEP３：弾性方程式の離散化（後編）

弾性方程式の離散化

前編で述べたようにある構造物に仮想変位が生じているとしてもその構造物がつりあいの状態にあるとみなすことができるので、 \begin{equation} \delta W + \delta U = 0 \label{10} \end{equation} が成り立ちます。このとき外力\(\,\vec{T_i}\,\)によって仮想変位\(\,\vec{u}\,\)を生じたとすると、この外力の仮想仕事量は、 \begin{equation} \delta \boldsymbol{W}=\int (\delta \vec{u})^{\mathrm{T}}\vec{T} \, d\Gamma \end{equation} ですから、両式を式\eqref{10}に代入すると、 \begin{eqnarray} &\int_{S^{e}} (\delta \vec{u})^{\mathrm{T}}\vec{T} \, d\Gamma+\int_{V^{e}} \{\delta \vec{\varepsilon}^{\mathrm{T}}\vec{\sigma}-(\delta \vec{u})^{\mathrm{T}}\vec{F}\} \, d\Omega=0 \notag \\ &\int_{S^{e}} (\boldsymbol{N}\delta \vec{u^{e}})^{\mathrm{T}}\vec{T} \, d\Gamma+\int_{V^{e}} \{(\boldsymbol{B}\delta \vec{u^{e}})^{\mathrm{T}}\vec{\sigma}-(\boldsymbol{N}\delta \vec{u^{e}})^{\mathrm{T}}\vec{F}\} \, d\Omega=0 \notag \\ &(\delta \vec{u^{e}})^{\mathrm{T}}\int_{S^{e}} \boldsymbol{N}^{\mathrm{T}}\vec{T} \, d\Gamma+(\delta \vec{u^{e}})^{\mathrm{T}}\int_{V^{e}} (\boldsymbol{B}^{\mathrm{T}}\vec{\sigma}-\boldsymbol{N}^{\mathrm{T}}\vec{F}) \, d\Omega=0 \notag \\ &(\delta \vec{u^{e}})^{\mathrm{T}}\left\{ \int_{S^{e}} \boldsymbol{N}^{\mathrm{T}}\vec{T} \, d\Gamma+\int_{V^{e}} (\boldsymbol{B}^{\mathrm{T}}\vec{\sigma}-\boldsymbol{N}^{\mathrm{T}}\vec{F}) \, d\Omega \right\} =0 \end{eqnarray} となります。これは任意の仮想変位についてなり立たなければならないので、 \begin{equation} \int_{S^{e}} \boldsymbol{N}^{\mathrm{T}}\vec{T} \, d\Gamma+\int_{V^{e}} (\boldsymbol{B}^{\mathrm{T}}\vec{\sigma}-\boldsymbol{N}^{\mathrm{T}}\vec{F}) \, d\Omega=0 \label{30} \end{equation} が得られます。応力を節点変位で表した式 \begin{equation} \vec{\varepsilon}= \begin{bmatrix} \boldsymbol{B}_1 & \boldsymbol{B}_2 & \cdots & \boldsymbol{B}_{20} \end{bmatrix}\vec{u^{e}}=\boldsymbol{B}\vec{u^{e}} \label{25} \end{equation} を代入すると式\eqref{30}は、 \begin{eqnarray} \int_{S^{e}} \boldsymbol{N}^{\mathrm{T}}\vec{T} \, d\Gamma+\int_{V^{e}} (\boldsymbol{B}^{\mathrm{T}}\boldsymbol{DB}\vec{u^{e}}-\boldsymbol{N}^{\mathrm{T}}\vec{F}) \, d\Omega=0 \\ \therefore \int_{V^{e}} \boldsymbol{B}^{\mathrm{T}}\boldsymbol{DB} \, d\Omega \cdot \vec{u^{e}}=\int_{V^{e}} \boldsymbol{N}^{\mathrm{T}}\vec{F} \, d\Omega - \int_{S^{e}} \boldsymbol{N}^{\mathrm{T}}\vec{T} \, d\Gamma \label{32} \end{eqnarray} となります。したがって \begin{equation} K^{e}=\int_{V^{e}} \boldsymbol{B}^{\mathrm{T}}\boldsymbol{DB} \, d\Omega \, , \quad \vec{F^{e}}=\int_{V^{e}} \boldsymbol{N}^{\mathrm{T}}\vec{F} \, d\Omega \, , \quad \vec{T^{e}}=-\int_{S^{e}} \boldsymbol{N}^{\mathrm{T}}\vec{T} \, d\Gamma \end{equation} と置くことにより、式\eqref{32}は、 \begin{equation} K^{e}\vec{u^{e}}=\vec{F^{e}}+\vec{T^{e}} \label{34} \end{equation} と表され、これが要素\(\,e\,\)について弾性方程式を離散化した式となっています。

数値積分法

各要素の節点変位は式\eqref{34}で定義されましたが、この積分方程式をコンピュータで算出するために、Gaussの数値積分を活用します。式\eqref{34}は直行座標系の式であるため、まずは前編で定義したBマトリクスの成分を局所座標系で表現する工夫が必要となります。これに伴って、各要素の体積(あるいは面積)積分の被積分領域も局所座標系に変換します。

局所座標への変換

節点\(\,i\,\)の形状関数\(\,N_i\,\)は、 \begin{equation} N_i=N_i(x,y,z) \, , \quad x=x(\xi,\eta,\zeta) \, , \quad y=y(\xi,\eta,\zeta) \, , \quad z=z(\xi,\eta,\zeta) \end{equation} であるため、1次偏導関数は、 \begin{eqnarray} \dd{N_i}{\xi}=\dd{N_i}{x}\dd{x}{\xi}+\dd{N_i}{y}\dd{y}{\xi}+\dd{N_i}{z}\dd{z}{\xi}\notag \\ \dd{N_i}{\eta}=\dd{N_i}{x}\dd{x}{\eta}+\dd{N_i}{y}\dd{y}{\eta}+\dd{N_i}{z}\dd{z}{\eta}\\ \dd{N_i}{\zeta}=\dd{N_i}{x}\dd{x}{\zeta}+\dd{N_i}{y}\dd{y}{\zeta}+\dd{N_i}{z}\dd{z}{\zeta} \notag \end{eqnarray} と表わされます。したがってこれをマトリクス表示すると、 \begin{eqnarray} \begin{bmatrix} \dd{N_i}{\xi}\\ \dd{N_i}{\eta}\\ \dd{N_i}{\zeta} \end{bmatrix}&=& \begin{bmatrix} \dd{x}{\xi} & \dd{y}{\xi} & \dd{z}{\xi} \\ \dd{x}{\eta} & \dd{y}{\eta} & \dd{z}{\eta} \\ \dd{x}{\zeta} & \dd{y}{\zeta} & \dd{z}{\zeta} \end{bmatrix} \begin{bmatrix} \dd{N_i}{x} \\ \dd{N_i}{y} \\ \dd{N_i}{z} \end{bmatrix} \notag \\ &=&\begin{bmatrix} J_{11} & J_{12} & J_{13} \\ J_{21} & J_{22} & J_{23} \\ J_{31} & J_{32} & J_{33} \end{bmatrix} \begin{bmatrix} \dd{N_i}{x} \\ \dd{N_i}{y} \\ \dd{N_i}{z} \end{bmatrix} \notag \\ &=& \boldsymbol{J} \begin{bmatrix} \dd{N_i}{x} \\ \dd{N_i}{y} \\ \dd{N_i}{z} \end{bmatrix} \label{37} \end{eqnarray} と表すことができます。上式で\(\,(x,y,z)\,\)直交座標座標を形状関数を用いて\(\,(\xi,\eta,\zeta)\,\)局所座標で表現すると、\(\,(x,y,z)\,\)は形状関数を用いて式\eqref{22}によって定義されるため、式\eqref{37}のJacobianマトリクスの成分は、 \begin{alignat}{1} J_{11}=\sum_{i=1}^{20}\dd{N_i}{\xi}x_i \quad J_{12}=\sum_{i=1}^{20}\dd{N_i}{\xi}y_i \quad J_{13}=\sum_{i=1}^{20}\dd{N_i}{\xi}z_i \notag \\ J_{21}=\sum_{i=1}^{20}\dd{N_i}{\eta}x_i \quad J_{22}=\sum_{i=1}^{20}\dd{N_i}{\eta}y_i \quad J_{23}=\sum_{i=1}^{20}\dd{N_i}{\eta}z_i \\ J_{31}=\sum_{i=1}^{20}\dd{N_i}{\zeta}x_i \quad J_{32}=\sum_{i=1}^{20}\dd{N_i}{\zeta}y_i \quad J_{33}=\sum_{i=1}^{20}\dd{N_i}{\zeta}z_i \notag \end{alignat} となります。したがってJacobianの逆行列が定まれば式\eqref{37}からすぐさま、\(\boldsymbol{B}\)マトリクスの成分が求まることになります。そこで逆行列と行列式を求めると、 \begin{eqnarray} \boldsymbol{J}^{-1}&=&\frac{1}{|\boldsymbol{J}|} \begin{bmatrix} J_{22}J_{33}-J_{23}J_{32} & -(J_{21}J_{33}-J_{23}J_{31}) & J_{21}J_{32}-J_{22}J_{31} \\ -(J_{12}J_{33}-J_{13}J_{32}) & J_{11}J_{33}-J_{13}J_{31} & -(J_{11}J_{32}-J_{12}J_{31}) \\ J_{12}J_{23}-J_{13}J_{22} & -(J_{11}J_{23}-J_{13}J_{21}) & J_{11}J_{22}-J_{12}J_{21} \end{bmatrix}^{\mathrm{T}} \notag \\ &=& \frac{1}{|\boldsymbol{J}|} \begin{bmatrix} J_{22}J_{33}-J_{23}J_{32} & J_{13}J_{32}-J_{12}J_{33} & J_{12}J_{23}-J_{13}J_{22} \\ J_{23}J_{31}-J_{21}J_{33} & J_{11}J_{33}-J_{13}J_{31} & J_{13}J_{21}-J_{11}J_{23} \\ J_{21}J_{32}-J_{22}J_{31} & J_{12}J_{31}-J_{11}J_{32} & J_{11}J_{22}-J_{12}J_{21} \end{bmatrix} \end{eqnarray} \begin{equation} |\boldsymbol{J}|=J_{11}J_{22}J_{33}-J_{11}J_{23}J_{32}+J_{12}J_{23}J_{31}-J_{12}J_{21}J_{33}+J_{13}J_{21}J_{32}-J_{13}J_{22}J_{31} \end{equation} となるので、これらを次式 \begin{equation} \begin{bmatrix} \dd{N_i}{x} \\ \dd{N_i}{y} \\ \dd{N_i}{z} \end{bmatrix}=\boldsymbol{J}^{-1} \begin{bmatrix} \dd{N_i}{\xi} \\ \dd{N_i}{\eta} \\ \dd{N_i}{\zeta} \end{bmatrix} \end{equation} に代入することで\(\,\boldsymbol{B}\,\)マトリクスの成分は、 \begin{eqnarray} \dd{N_i}{x}=\frac{1}{|\boldsymbol{J}|}\left\{(J_{22}J_{33}-J_{23}J_{32})\dd{N_i}{\xi}+(J_{13}J_{32}-J_{12}J_{33})\dd{N_i}{\eta}+(J_{12}J_{23}-J_{13}J_{22})\dd{N_i}{\zeta}\right\} \notag \\ \dd{N_i}{y}=\frac{1}{|\boldsymbol{J}|}\left\{(J_{23}J_{31}-J_{21}J_{33})\dd{N_i}{\xi}+(J_{11}J_{33}-J_{13}J_{31})\dd{N_i}{\eta}+(J_{13}J_{21}-J_{11}J_{23})\dd{N_i}{\zeta}\right\} \\ \dd{N_i}{z}=\frac{1}{|\boldsymbol{J}|}\left\{(J_{21}J_{32}-J_{22}J_{31})\dd{N_i}{\xi}+(J_{12}J_{31}-J_{11}J_{32})\dd{N_i}{\eta}+(J_{11}J_{22}-J_{12}J_{21})\dd{N_i}{\zeta}\right\} \notag \end{eqnarray} と定まり、これが局所座標系で表わされた式となっています。 次に、被積分領域を局所座標へ変換することについて論じます。まず、被積分領域が体積である場合、直行座標系で任意の座標を、 \begin{equation} \vec{r}=x\vec{i}+y\vec{j}+z\vec{k} \, , \quad x=x(\xi,\eta,\zeta) \, , \quad y=y(\xi,\eta,\zeta) \, , \quad z=z(\xi,\eta,\zeta) \end{equation} とすると局所座標の任意の座標\((d\vec{r_{\xi}},d\vec{r_{\eta}},d\vec{r_{\zeta}})\)は、 \begin{eqnarray} d\vec{r_{\xi}}=\dd{\vec{r}}{\xi}d\xi=\left(\dd{x}{\xi}\vec{i}+\dd{y}{\xi}\vec{j}+\dd{z}{\xi}\vec{k}\right)d\xi \notag \\ d\vec{r_{\eta}}=\dd{\vec{r}}{\eta}d\eta=\left(\dd{x}{\eta}\vec{i}+\dd{y}{\eta}\vec{j}+\dd{z}{\eta}\vec{k}\right)d\eta \\ d\vec{r_{\zeta}}=\dd{\vec{r}}{\zeta}d\zeta=\left(\dd{x}{\zeta}\vec{i}+\dd{y}{\zeta}\vec{j}+\dd{z}{\zeta}\vec{k}\right)d\zeta \notag \end{eqnarray} ですので、直行座標系の体積\(\,d\Omega\,\)を局所座標系に変換したものは、\(\,(x,y,z)\,\)座標において、 \begin{eqnarray} d\Omega &=& d\vec{r_{\xi}}\times d\vec{r_{\eta}} \cdot d\vec{r_{\zeta}} \notag \\ &=& \begin{vmatrix} \vec{i} & \vec{j} & \vec{k} \\ \dd{x}{\xi} & \dd{y}{\xi} & \dd{z}{\xi} \\ \dd{x}{\eta} & \dd{y}{\eta} & \dd{z}{\eta} \end{vmatrix} \cdot \begin{bmatrix} \dd{x}{\zeta} \\ \dd{y}{\zeta} \\ \dd{z}{\zeta} \end{bmatrix} = \begin{vmatrix} \dd{x}{\xi} & \dd{y}{\xi} & \dd{z}{\xi} \\ \dd{x}{\eta} & \dd{y}{\eta} & \dd{z}{\eta} \\ \dd{x}{\zeta} & \dd{y}{\zeta} & \dd{z}{\zeta} \end{vmatrix}d\xi d\eta d\zeta = |\boldsymbol{J}|d\xi d\eta d\zeta \label{45} \end{eqnarray} と表わされます。 次に被積分領域が面積である場合、直行座標系で \begin{equation} \vec{r}=x\vec{i}+y\vec{j}+z\vec{k} \, , \quad x=x(\xi,\eta,-1) \, , \quad y=y(\xi,\eta,-1) \, , \quad z=z(\xi,\eta,-1) \notag \end{equation} と表わされるベクトルは、\((\xi,\eta,\zeta)\)空間における局所座標平面\(\,\zeta=-1\,\)を表します。そこで表面荷重のかかる要素面をこれにとると、平面\(\,\zeta=-1\,\)上の点では、\(N_9～N_{20}=0\)ですから、この面に応じた\(\,(x,y,z)\,\)座標は、 \begin{equation} x=\sum_{i=1}^{20}N_ix_i=\sum_{i=1}^8 N_ix_i \quad , \quad y=\sum_{i=1}^{20}N_iy_i=\sum_{i=1}^8 N_iy_i \quad , \quad z=\sum_{i=1}^{20}N_iz_i=\sum_{i=1}^8 N_iz_i \notag \end{equation} と表すことができます。このとき、 \begin{eqnarray} d\vec{r_1}=\dd{\vec{r}}{\xi}d\xi=\left(\dd{x}{\xi}\vec{i}+\dd{y}{\xi}\vec{j}+\dd{z}{\xi}\vec{k}\right)d\xi \notag \\ d\vec{r_2}=\dd{\vec{r}}{\eta}d\eta=\left(\dd{x}{\eta}\vec{i}+\dd{y}{\eta}\vec{j}+\dd{z}{\eta}\vec{k}\right)d\eta \end{eqnarray} と定義した\(\,d\vec{r_1}\,\)と\(\,d\vec{r_2}\,\)は\(\,(\xi,\eta,\zeta)\,\)空間における局所座標平面\(\,\zeta=-1\,\)上の\(\,\xi\,\)方向のベクトルと、\(\,\eta\,\)方向のベクトルとなりますから、 \(\,d\vec{r_1}\,\)と\(\,d\vec{r_2}\,\)で張られる面の法線方向の成分\(\,d\vec{\Gamma^{p}}\,\)は\(\,(x,y,z)\,\)空間において、 \begin{eqnarray} d\vec{\Gamma^{p}} &=& d\vec{r_1}\times d\vec{r_2} \notag \\ &=& \begin{vmatrix} \vec{i} & \vec{j} & \vec{k} \\ \dd{x}{\xi} & \dd{y}{\xi} & \dd{z}{\xi} \\ \dd{x}{\eta} & \dd{y}{\eta} & \dd{z}{\eta} \end{vmatrix}d\xi d\eta \notag \\ &=& \begin{vmatrix} \vec{i} & \vec{j} & \vec{k} \\ \displaystyle \sum_{i=1}^8 \dd{N_i}{\xi}x_i & \displaystyle \sum_{i=1}^8 \dd{N_i}{\xi}y_i & \displaystyle \sum_{i=1}^8 \dd{N_i}{\xi}z_i \\ \displaystyle \sum_{i=1}^8 \dd{N_i}{\eta}x_i & \displaystyle \sum_{i=1}^8 \dd{N_i}{\eta}y_i & \displaystyle \sum_{i=1}^8 \dd{N_i}{\eta}z_i \end{vmatrix}d\xi d\eta d\zeta \notag \\ &=& \begin{vmatrix} \vec{i} & \vec{j} & \vec{k} \\ J_{11} & J_{12} & J_{13} \\ J_{21} & J_{22} & J_{23} \end{vmatrix}d\xi d\eta \end{eqnarray} で表されます（ただし方向は\(\,d\vec{r_1}\,\)から\( \,d\vec{r_1}\,\)に右ねじを回して進む向きです。そこでこの方向を\(\,p\,\)方向としています）。これをベクトル表示したものを、 \begin{equation} d\vec{\Gamma^p}= \begin{bmatrix} J_{12}J_{23}-J_{13}J_{22} \\ J_{13}J_{21}-J_{11}J_{23} \\ J_{11}J_{22}-J_{12}J_{21} \end{bmatrix}d\xi d\eta= \begin{bmatrix} V^p_x\\ V^p_y\\ V^p_z \end{bmatrix}d\xi d\eta \label{48} \end{equation} と定義しておきます。 このとき\(\,p\,\)方向、\(\,\xi\,\)方向、\(\,\eta\,\)方向それぞれの単位ベクトルは\(\,(x,y,z)\,\)座標において \begin{eqnarray} \vec{n_p}&=&\frac{1}{\sqrt{ (V^p_x)^2 + (V^p_y)^2 +(V^p_z)^2 }} \begin{bmatrix} V^p_x \\ V^p_y\\ V^p_z \end{bmatrix}= \begin{bmatrix} n_{px} \\ n_{py}\\ n_{pz} \end{bmatrix} \\ \vec{n_\xi}&=&\frac{1}{\sqrt{ J^2_{11} + J^2_{12} +J^2_{13} }} \begin{bmatrix} J_{11} \\ J_{12} \\ J_{13} \end{bmatrix}= \begin{bmatrix} n_{\xi x} \\ n_{\xi y}\\ n_{\xi z} \end{bmatrix} \\ \vec{n_\eta}&=&\frac{1}{\sqrt{ J^2_{21} + J^2_{22} +J^2_{23} }} \begin{bmatrix} J_{21} \\ J_{22} \\ J_{23} \end{bmatrix}= \begin{bmatrix} n_{\eta x} \\ n_{\eta y}\\ n_{\eta z} \end{bmatrix} \label{51} \end{eqnarray} です。また、 \begin{equation} d\vec{\Gamma^\xi}=\vec{n_\eta}\times\vec{d\Gamma^p} \quad , \quad d\vec{\Gamma^\eta}=\vec{d\Gamma^p}\times\vec{n_\xi} \end{equation} ですから、式\eqref{48}～\eqref{51}を代入すると、 \begin{eqnarray} d\vec{\Gamma^\xi}&=& \begin{bmatrix} n_{\eta x} \\ n_{\eta y}\\ n_{\eta z} \end{bmatrix}\times \begin{bmatrix} V^p_x\\ V^p_y\\ V^p_z \end{bmatrix}d\xi d\eta = \begin{bmatrix} n_{\eta y}V^p_z-n_{\eta z}V^p_y \\ n_{\eta z}V^p_x-n_{\eta x}V^p_z \\ n_{\eta x}V^p_y-n_{\eta y}V^p_x \end{bmatrix}d\xi d\eta = \begin{bmatrix} V^\xi_x \\ V^\xi_y \\ V^\xi_z \end{bmatrix}d\xi d\eta \label{53}\\ d\vec{\Gamma^\eta}&=& \begin{bmatrix} V^p_x\\ V^p_y\\ V^p_z \end{bmatrix}d\xi d\eta \times \begin{bmatrix} n_{\xi x} \\ n_{\xi y}\\ n_{\xi z} \end{bmatrix}= \begin{bmatrix} n_{\xi z}V^p_y-n_{\xi y}V^p_z \\ n_{\xi x}V^p_z-n_{\xi z}V^p_x \\ n_{\xi y}V^p_x-n_{\xi x}V^p_y \end{bmatrix}d\xi d\eta = \begin{bmatrix} V^\eta_x \\ V^\eta_y \\ V^\eta_z \end{bmatrix}d\xi d\eta \label{54} \\ \end{eqnarray} となり、これを\(\,\xi-\eta\,\)平面に係る表面荷重の演算に用います。

Gaussの数値積分

Gaussの数値積分は多項式\(\,f(x)\,\)に対して、 \begin{equation} \Delta=\int^1_{-1}f(x)dx-\sum_{i=1}^n f(x_i)w_i \end{equation} の\(\Delta\)を最小にするように幾何対称を利用して表2.1のようにGauss point \(\,x_i\,\)とその重み\(\,w_i\,\)を定めるため、多項式の積分はGaussの数値積分によって厳密な解を得ることができます。 3次元問題においては\(\,n=3\,\)で十分であり、次の積分： \begin{equation} I = \int^1_{-1}\int^1_{-1}\int^1_{-1}f(\xi,\eta,\zeta)d\xi d\eta d\zeta \end{equation} を求めるにあたり、まず、\(\,\eta\,\)、\(\,\zeta\,\)を固定した\(\,\xi\,\)のみの関数とみて重み\(\,w_i\,\)を考慮して、次に示すように内側の積分を行います。 \begin{equation} \int^1_{-1}f(\xi,\eta,\zeta)d\xi=\sum^n_i f(\xi_i,\eta,\zeta)w_i=\psi(\eta,\zeta) \end{equation} すると、\(I\)は、 \begin{equation} I = \int^1_{-1}\int^1_{-1}\psi(\eta,\zeta)d\eta d\zeta \label{58} \end{equation} となります。同様に、\(\,\zeta\,\)を固定して\(\,\eta\,\)のみの関数とみると、 \begin{equation} \int^1_{-1}\psi(\eta,\zeta)d\eta=\sum^n_j \psi(\eta_j,\zeta)w_j=\phi(\zeta) \label{59} \end{equation} となります。この式に式\eqref{58}\eqref{59}を代入すると \begin{equation} I=\sum^n_k \sum^n_j \sum^n_i f(\xi_i,\eta_j,\zeta_k)w_i w_j w_k \end{equation} としてGaussの数値積分を表す式が得られます。表１にGauss Pointとその重みを示します。

要素剛性マトリクスと荷重ベクトルのGauss積分表示

前節では離散化された弾性方程式を式\eqref{34}として示しました。この要素剛性マトリクスと荷重ベクトルを局座標変換してGauss積分にて表示し、コンピュータプログラムにセットできる状態まで、より具体的な展開を行います。

要素剛性マトリクス

要素剛性マトリクスについては被積分領域が体積ですから、式\eqref{45}の変換により、 \begin{eqnarray} K^{e} &=& \int_{V^{e}} \boldsymbol{B}^{\mathrm{T}}\boldsymbol{DB} \, d\Omega \notag \\ &=& \int^1_{-1} \int^1_{-1} \int^1_{-1} \boldsymbol{B}^{\mathrm{T}}\boldsymbol{DB}|\boldsymbol{J}|d\xi d\eta d\zeta \notag \\ &=& \sum^3_{k=1} \sum^3_{j=1} \sum^3_{i=1} \left\{\boldsymbol{B}^{\mathrm{T}}\boldsymbol{DB}\right\}_{ijk}|\boldsymbol{J}|w_i w_j w_k \label{62} \end{eqnarray}\label{64} と表わされます。\(\boldsymbol{B}\)マトリクスが\(\,6\times 60\,\)、\(\boldsymbol{D}\)マトリクスが\(\,6\times 6\,\)の行列ですから、上式\eqref{62}の{ }の中身は\(\,60\times 60\,\)の行列となっています。この行列を \begin{equation} K^e = \begin{bmatrix} K^e_{mn} \end{bmatrix} \quad (m=1, \, 60 \quad n=1, \, 60) \end{equation} と表すと、その成分は、 \begin{equation} K^e_{mn}=\sum^3_{k=1} \sum^3_{j=1} \sum^3_{i=1} \left\{ \sum^6_{p=1} \left( B_{mp}^{\mathrm{T}} \cdot \sum^6_{q=1} D_{pq} B_{qn} \right) \right\}_{ijk}|\boldsymbol{J}|w_i w_j w_k \end{equation} で与えられます。このとき材料マトリクス\(\boldsymbol{D}=\left[D_{ij}\right]\)とは、 \begin{equation} \begin{cases} D_{ij}=\frac{E\nu}{(1+\nu)(1-2\nu)} & (i=1～3 & j=1～3 & i\neq j) \\ D_{ij}=\frac{E(1-\nu)}{(1+\nu)(1-2\nu)} & (i=1～3 & i=j) \\ D_{ij}=\frac{E}{2(1+\nu)} & (i=4～6 & i=j) \\ D_{ij}=0 & (others) \end{cases} \end{equation} と表されるものです。

荷重ベクトル

荷重ベクトルも被積分領域が体積ですから、式\eqref{45}の変換により、 \begin{eqnarray} \vec{F^{e}} &=& \int_{V^{e}} \boldsymbol{N}^{\mathrm{T}}\vec{F} \, d\Omega \notag \\ &=& \int^1_{-1} \int^1_{-1} \int^1_{-1} \boldsymbol{N}^{\mathrm{T}}|\boldsymbol{J}|d\xi d\eta d\zeta \cdot \vec{F} \notag \\ &=& \left[ \sum^3_{k=1} \sum^3_{j=1} \sum^3_{i=1} \boldsymbol{N}^{\mathrm{T}}_{ijk}|\boldsymbol{J}|w_i w_j w_k \right] \cdot \vec{F} \end{eqnarray} と表わされます。この行列を \begin{equation} \vec{F^e} = \begin{bmatrix} F^e_{m1} \end{bmatrix} \quad (m=1, \, 60) \end{equation} と表すと、その成分は、 \begin{equation} F^e_{m1}=\left[\sum^3_{n=1} \left\{\sum^3_{k=1} \sum^3_{j=1} \sum^3_{i=1} \left( N_{mn}^{\mathrm{T}} \right)_{ijk}|\boldsymbol{J}|w_i w_j w_k \right\}\right] \cdot F_{n1} \end{equation} で与えられます。なお、この荷重ベクトルは重力などの体積力によるものです。

表面荷重ベクトル

表面荷重ベクトルは、 \begin{eqnarray} \vec{T^{e}} &=& -\int_{S^{e}} \boldsymbol{N}^{\mathrm{T}} \vec{T} \, d\Gamma \notag \\ &=& -\int^1_{-1} \int^1_{-1} \boldsymbol{N}^{\mathrm{T}}\boldsymbol{N}d\vec{\Gamma^\xi} \cdot \vec{T^\xi}-\int^1_{-1} \int^1_{-1} \boldsymbol{N}^{\mathrm{T}}\boldsymbol{N}d\vec{\Gamma^\eta} \cdot \vec{T^\eta}-\int^1_{-1} \int^1_{-1} \boldsymbol{N}^{\mathrm{T}}\boldsymbol{N}d\vec{\Gamma^p} \cdot \vec{T^p} \end{eqnarray} のように局座標で表現できます。この被積分領域は面積ですから、式\eqref{48}、\eqref{53}、\eqref{54}の変換により、\(\xi-\eta\,\)平面に係る表面荷重は、 \begin{eqnarray} \vec{T^{e}} = &-\int_{S^{e}} \begin{bmatrix} N_1I \\ N_2I \\ \vdots \\ N_8I \end{bmatrix} \begin{bmatrix} V^\xi_x \\ V^\xi_y \\ V^\xi_z \end{bmatrix} \begin{bmatrix} N_1 & N_2 & \cdots & N_8 \end{bmatrix} d\xi d\eta \cdot \begin{bmatrix} \vec{T^\xi_1} \\ \vec{T^\xi_2} \\ \vdots \\ \vec{T^\xi_8} \end{bmatrix} \\ &-\int_{S^{e}} \begin{bmatrix} N_1I \\ N_2I \\ \vdots \\ N_8I \end{bmatrix} \begin{bmatrix} V^\eta_x \\ V^\eta_y \\ V^\eta_z \end{bmatrix} \begin{bmatrix} N_1 & N_2 & \cdots & N_8 \end{bmatrix} d\xi d\eta \cdot \begin{bmatrix} \vec{T^\eta_1} \\ \vec{T^\eta_2} \\ \vdots \\ \vec{T^\eta_8} \end{bmatrix} \notag \\ &-\int_{S^{e}} \begin{bmatrix} N_1I \\ N_2I \\ \vdots \\ N_8I \end{bmatrix} \begin{bmatrix} V^p_x \\ V^p_y \\ V^p_z \end{bmatrix} \begin{bmatrix} N_1 & N_2 & \cdots & N_8 \end{bmatrix} d\xi d\eta \cdot \begin{bmatrix} \vec{T^p_1} \\ \vec{T^p_2} \\ \vdots \\ \vec{T^p_8} \end{bmatrix} \notag \end{eqnarray} で与えられます。