このプログラムのソース

入力パラメータ Go Back

Nums of M 質量の数(自由度の数)(1〜10の整数)
Nums of S ばねの数(1〜10の整数)
Nums of D 減衰(速度に比例)の数(ばねと並列に接続し1〜10の整数)
Mass  (M0,M1,...M9) 質量(kg)
Force 質量に作用する周期的な外力を Fsin(2πft) とする、振幅 F(N)
Stiff.(K0,K1,...K9) ばね定数(N/m)
Damp. (D0,D1,...D9) 速度に比例する減衰(N*sec/m)
Connect to ばねが接続される質量の番号(0 〜 質量の数−1 の整数)  
固定端の場合は −1とする
Display M# 計算結果を表示する質量の番号(0 〜 質量の数−1 の整数)
      

周波数応答(伝達関数)とインパルス応答


計算方法及びインプリメンテーション

 計算プログラムの主要な部分の抜粋を以下に示します。


 剛性マトリックスの生成

  質量の変位と速度に関する力のつり合いから、連立1次方程式の係数として、

  ばねの接続条件により決定されます。

// Assemble Structure Stiffness Matrix [k0]

      for (i=0; i<ns; i++){
          e1=(int)(stiff[i][1]);/* connected end to mass #e1 */
          e2=(int)(stiff[i][2]);/* connected end to mass #e2 */
          if (e1==-1)
              k0[e2][e2]=k0[e2][e2]+stiff[i][0];
          else if (e2==-1)
              k0[e1][e1]=k0[e1][e1]+stiff[i][0];
          else{
              k0[e1][e2]=k0[e1][e2]-stiff[i][0];
              k0[e2][e1]=k0[e2][e1]-stiff[i][0];
              k0[e1][e1]=k0[e1][e1]+stiff[i][0];
              k0[e2][e2]=k0[e2][e2]+stiff[i][0];
          }
       } 
// Assemble Structure Damping Matrix [c0]

      for (i=0; i<nd; i++){
          e1=(int)(stiff[i][1]);/* connected end to mass #e1 */
          e2=(int)(stiff[i][2]);/* connected end to mass #e2 */
          if (e1==-1)
              c0[e2][e2]=c0[e2][e2]+damp[i][0];
          else if (e2==-1)
              c0[e1][e1]=c0[e1][e1]+damp[i][0];
          else{
              c0[e1][e2]=c0[e1][e2]-damp[i][0];
              c0[e2][e1]=c0[e2][e1]-damp[i][0];
              c0[e1][e1]=c0[e1][e1]+damp[i][0];
              c0[e2][e2]=c0[e2][e2]+damp[i][0];
          }

       }



 剛性マトリックスの係数が複素数となるため以下の方法で、逆マトリックスを計算

  連立1次方程式を [ A + jB ]{ x }={ P }
     ここで、
      [ A ] ばねによる係数マトリックス
      [ B ] 速度に比例する減衰による係数マトリックス
      { x } 質量の変位ベクトル
      { P } 質量に働く外力ベクトル
  とすると変位ベクトルは、

      { x }=[ C + jD ]{ P }=[ E + jF ]

  [ C + jD ]を求めるために、マトリックス[ A ]及び[ B ]を以下のように配置

  し、逆マトリックスを計算

     [A   - B    I    0 ]
     [B    A    0    I ]
  ここで、
     [ I ]は単位マトリックス

  以下のプログラムでは、m1[][]に[ A ]、[ B ]及び[ I ]を配置し、
  逆マトリックスをガウスの消去法により計算

    for (j=0; j<nm; j++){
      for (i=0; i<2*nm; i++){
        i0[ i ][ j ]=0;
        i0[ j ][ j ]=1;
      }
    }
    for (i=0; i<nm; i++){
      for (j=0; j<nm; j++){
       if (i==j)                  // w1=2πf 
         m1[i][i]=k0[i][i]-w1*w1*mass[i];// system matrix [ A ]
       else{
         m1[i][j]=k0[i][j];
         m1[i][i]=k0[i][i]-w1*w1*mass[i];// system matrix [ A ]
       }
     }
     for (k=0; k<nm; k++){
         m1[nm+i][k]=w1*c0[i][k];// system matrix [ B ]
     }
   }

   for (i=0; i<nm; i++){
     for (j=0; j<nm; j++){
       m1[i][nm+j]= -m1[nm+i][j];
       m1[nm+i][nm+j]= m1[i][j];
     }
   }
   for (i=0; i<2*nm; i++){
     nj= 2*nm;
     for (j=nj; j<4*nm; j++){
       m1[i][j]=0;
     }
     m1[i][2*nm+i]=1;
   }
   for (k=0; k<2*nm; k++){
     pivi=m1[k][k];
     k1=k+1;
     for (j=k1; j<4*nm; j++){
       m1[k][j]=(m1[k][j])/pivi;
     }
     for (i=0; i<2*nm; i++){
       if(i !=k){
         w=m1[i][k];
         for (j=k1; j<4*nm; j++){
           m1[i][j]=m1[i][j]-w*m1[k][j];
         }
       }
     }
   }
// end of inversion
以下省略


 参考文献「高速フーリエ変換」E.O Brigham 宮川洋/今井秀樹訳

     「ディジタル信号処理入門」城戸健一 丸善

     「Cによる科学技術計算」小池慎一 CQ出版

     「有限要素法の数値計算」K.J Bathe/E.L.Wilson 菊池文雄訳

     「演習有限要素法」三好俊郎/白鳥正樹 サイエンス社

     「工業振動学」S.Timoshenko 谷下市松/渡辺茂訳

     「有限要素解析の基礎」R.H.Gallagher 川井忠彦 川島/中沢/藤谷訳


Go back to home page