FEM2D 部品

最低限必要な部品。
言語は疑似コードです。

2次節点

class Node2D
{
 //■FEM
 //自身を含む節点を格納する全体リストへの参照
 //FEMクラスでも良い
 List<Node2D> Nodes;
 
 public int GetIndex()
 {
   return this.Nodes.IndexOf(this);
 }
 
 //■要素
 //節点に接続する（節点に構築される）要素
 List<Truss2D> ConnectedElement = new List<Truss2D>();
 
 //■節点
 //初期座標
 Vector2D Initial_pos ;
 //現在座標
 Vector2D Pos ;
   
 Vector2D DX()  
 {
   return this.Initial_pos - this.Pos;    
 } 
 
 //コンストラクタ
 Node2D(List<Node2D> global_list, Vector2D initial_pos)
 {
   this.Nodes = global_list;
   //newして参照をブチ切ること
   this.Initial_pos = new Vector2D(initial_pos);
   this.Pos = new Vector2D(initial_pos);
 } 
}

2次トラス要素

class Truss2D
{
 //■FEM
 //自身が所属するFEMに対する参照があっても良い
 List<Node2D> _nodes;

 //1Dは節点のインデックスを持つのみだったが
 //2DはK行列構築時にトラス要素の傾きが必要になる。
 //なのでクラス内に節点への参照があったほうが楽。
 Node2D StartNode()
 {
   return _nodes.get(StartNodeIndex);
 } 
 Node2D EndNode()
 {
   return _nodes.get(EndNodeIndex);
 }
 
 //EV-SV
 Vector2D GetLV()
 {
   return this.EndNode().Pos - this.StartNode().Pos;
 } 
 float GetLength()
 {
   return GetLV().Length;
 }
 
 //thisトラス要素の傾きを求める(ラジアン)
 double CalcuAngleRadian()
 {
   Vector2D lv = this.GetLV();
   //x軸とLVの成す角
   //戻り値は-pi<Θ<+pi, 時計の3時12時9時が0～pi, 3時6時9時が0～-pi
   double ret = Math.atan2(lv.y,lv.x);
   return ret;
 }
 //thisトラス要素の傾きを求める(度数)
 double CalcuAngleDegree()
 {
   double ret = CalcuAngleRadian()*180/Math.PI;
   return ret;
 }
 
 //端1のインデックス
 public int StartNodeIndex;
 //端2のインデックス
 public int EndNodeIndex;
 //ばね定数
 public double k = 1;
 
 //コンストラクタ
 Truss2D(List<Node2D> global_list)
 {
   _nodes = global_list;
 }
}//Truss2D

マトリックス変位法

 void static_displacement(FEM2D fem, List<int> fix_node_indexes, List<int> displace_indexes)
 {
   //■固定点のコピー
   List<int> copy_fix = copy(fix_node_indexes);
   
   //■全体剛性行列Kのコピー
   double[][] copy_K = make_copy(fem.Global_K);
   //Kに0行0列がある場合、それを固定点として追加
   List<int> zero_row = search_zero_row(copy_K);
   for(int i = 0; i<zero_row.size(); i++)
   {
     if(contain(copy_fix, zero_row[i])){continue;}     
     copy_fix.add(zero_row[i]);
   }

   //■変位ベクトルU更新
   fem._make_u();
   double[] copy_U = make_copy(fem.Global_U);  
   //■力ベクトルF作成
   double[] copy_F = new double[fem.Global_F.length];
   
   //↑↑ここまで行列の大きさは変えてない↑↑

   //変位境界処理
   //削除すべき行と列の対角を1、削除すべき行と列の対角以外の要素を0に
   //それに伴い力行列Fを修正
   boundary_process(copy_K,copy_U,copy_F,copy_fix,displace_indexes);  
   
   //連立方程式を解く
   double[] u = solve(copy_K, copy_F);
   if(u==null)//失敗した場合
   {
     is_singlar = true;    
     return;
   }
   copy_vector(u,fem.Global_U);
   
   //変位適用
   for(int i = 0; i<fem.Node2Ds.size(); i++)
   {
     Node2D tnode = fem.Node2Ds[i];
     
     //固定点でないなら変位適用
     if(!contain(fix_node_indexes, i*fem.ComponentCount))
     {        
       tnode.Pos.x = (float)(tnode.Initial_pos.x + fem.Global_U[i*fem.ComponentCount]);
     }
     if(!contain(fix_node_indexes, i*fem.ComponentCount+1))
     {
       tnode.Pos.y = (float)(tnode.Initial_pos.y + fem.Global_U[i*fem.ComponentCount+1]);
     }         
   }//for
   copy_vector(copy_F,fem.Global_F);    
   
   is_singlar = false;
 }   

変位境界処理

//K行列の削除すべき行と列（固定点、あるいは係数Ｋが飛んだ時
//（例えば2次トラスが水平垂直にならぶと1次に縮退する、その時係数Ｋが０になる項がでる））
//の対角を1、削除すべき行と列の対角以外の要素を0に

//K行列の変更に伴って力ベクトルFを同時に修正

void boundary_process(
 double[][] k, double[] u, double[] f, List<int> fix_indexes, List<int> displacement_indexes)
{
   //変位境界処理
   for(int i = 0; i<k.length; i++)
   {      
     if(contain(fix_indexes, i)||contain(displacement_indexes, i))
     {
       for(int j = 0; j<k.length; j++)
       {
         f[j]-=k[j][i]*u[i];
         k[j][i]=0;
         k[i][j]=0;
       }        
       f[i]=u[i];
       k[i][i]=1;
     }
   }    
}  

全体剛性行列Kの組み立て

要素と節点を組み立てた最初の段階（初期状態）で1回だけ計算すれば良い
逆に要素や節点に変更を加えるのであればその都度再計算が必要

void Assembly_K(FEM2D fem)
{
 init_matrix(fem.Global_K);
 
 for(int i = 0; i<fem.Truss2Ds.size(); i++)
 {
   Truss2D truss = fem.Truss2Ds[i];
   int start_node_index = truss.StartNodeIndex;
   int end_node_index = truss.EndNodeIndex;     
   double k = truss.k;
   double sita = truss.CalcuAngleRadian();
   
   double cos = Math.cos(sita);
   double sin = Math.sin(sita);
   
   if(cos<0.0000001){cos = 0;}
   if(sin<0.0000001){sin = 0;}
   
   double A = cos*cos*k;
   double B = cos*sin*k;
   double C = sin*sin*k;
   
   int start_index = start_node_index*fem.ComponentCount;
   int end_index = end_node_index*fem.ComponentCount;
   
   fem.Global_K[start_index+0][start_index+0] += A;
   fem.Global_K[start_index+0][start_index+1] += B;
   fem.Global_K[start_index+1][start_index+0] += B;
   fem.Global_K[start_index+1][start_index+1] += C;
   
   fem.Global_K[end_index+0][end_index+0] += A;
   fem.Global_K[end_index+0][end_index+1] += B;
   fem.Global_K[end_index+1][end_index+0] += B;
   fem.Global_K[end_index+1][end_index+1] += C;
   
   fem.Global_K[start_index+0][end_index+0] += -A;    
   fem.Global_K[start_index+0][end_index+1] += -B; 
   fem.Global_K[start_index+1][end_index+0] += -B; 
   fem.Global_K[start_index+1][end_index+1] += -C; 
   
   fem.Global_K[end_index+0][start_index+0] += -A;
   fem.Global_K[end_index+0][start_index+1] += -B;
   fem.Global_K[end_index+1][start_index+0] += -B;
   fem.Global_K[end_index+1][start_index+1] += -C;       
 }
}