FEM for St.Venant-Kirchhoff Material

サンブナン体の有限要素法による解法

Last Update:2009年05月15日

$[{K_L}_e]^{p,q}_{i,j}=\int_V_e\frac{\partial S_{gh}}{\partial u^q_j}\[B_{gh}\]^p_idV=\int_V_e\bar{C}_{ghef}\[B_{ef}\]^q_j\[B_{gh}\]^p_idV\sime\sum_{\alpha}\omega_{\alpha}\bar{C}_{ghef}\[B_{ef}\]^q_j\[B_{gh}\]^p_iJ_{\alpha}$

ここで $\bar{C}$ は簡単な構成式テンソルであったので $\bar{C}_{efgh}$ を配列として持たずに直接的に乗算を行うことができる。

$\bar{C}_{ghef}\[B_{ef}\]^q_j\[B_{gh}\]^p_i=(\lambda\delta_{ef}\delta_{gh}+\mu\delta_{ge}\delta_{hf}+\mu\delta_{gf}\delta_{he})\[B_{ef}\]^q_j\[B_{gh}\]^p_i\\ \qquad\qquad\qquad\qquad\qquad=\mu(\[B_{gh}\]^q_j\[B_{gh}\]^p_i+\[B_{hg}\]^q_j\[B_{hh}\]^p_i)+\lambda\[B_{gg}\]^q_j\[B_{hh}\]^p_i$

ここで、次のようにプログラムの中で変数が対応づけられているとする。

ndim:空間の次元の数、nnoel:要素あたりの節点数
lambda= $\lambda$ 、myu= $\mu$
detwei= $\omega_{\alpha}J_{\alpha}$
emat[p][q][i][j]= $\[K_{NL}_e\]^{p,q}_{i,j}$ 、stress[g][h]= $S_{gh}$ 、

$\[B_{gh}\]^p_i$

初期剛性行列を作るルーティンは次のようになる。

// 初期剛性行列を足しこむ
for(unsigned int inoel=0;inoel<nnoel;inoel++){
for(unsigned int jnoel=0;jnoel<nnoel;jnoel++){
   for(unsigned int idim=0;idim<ndim;idim++){
   for(unsigned int jdim=0;jdim<ndim;jdim++){
      {　// muの２つの項を足しこむ
         double dtmp1 = 0.0;
         for(unsigned int gdim=0;gdim<ndim;gdim++){
         for(unsigned int hdim=0;hdim<ndim;hdim++){
            dtmp1 += disp2strain[inoel][idim][gdim][hdim]*disp2strain[jnoel][jdim][gdim][hdim]
               +disp2strain[inoel][idim][gdim][hdim]*disp2strain[jnoel][jdim][hdim][gdim];
         }
         }
         emat[inoel][jnoel][idim][jdim] += detwei*myu*dtmp1;
      }
      {  // lambdaの項を足しこむ
         double dtmp1=0.0, dtmp2=0.0;
         for(unsigned int gdim=0;gdim<ndim;gdim++){
            dtmp1+=disp2strain[inoel][idim][gdim][gdim];
            dtmp2+=disp2strain[jnoel][jdim][gdim][gdim];
         }
         emat[inoel][jnoel][idim][jdim] += detwei*lambda*dtmp1*dtmp2;
       }
   }
   }
}