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行列の削除すべき行と列(固定点、あるいは係数Kが飛んだ時
//(例えば2次トラスが水平垂直にならぶと1次に縮退する、その時係数Kが0になる項がでる))
//の対角を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;
}
}
この記事が気に入ったらサポートをしてみませんか?