<有限要素法>要素剛性マトリクス算出

各要素剛性マトリクスを算出し、最終的に全体の剛性マトリクス算出するまでのプログラムの説明です。

全体の剛性マトリクスは、各要素（element)から算出された要素剛性マトリクスを全体の剛性マトリクスに加えこんで行くことにより、作成されます。

そのため要素剛性マトリクスは、入力された要素（element)分の回数、作成すること、つまり、

for (element_index = 0; element_index < element_num; element_index++) {
      ：
      ：
      ：
｝
のfor文の中で、各要素剛性マトリクスを算出し、また全体の剛性マトリクスに加えこんで行く形になります。
element_num: 入力される要素の数


最初に、各要素から算出される要素剛性マトリクスは、以下の様に作成されます。

まず、各要素の長さ：element_point_lenを計算します

// Calculate element_point_len
element_point_len_x = (double) (point_x[element_point_1[element_index]] - point_x[element_point_2[element_index]]);
element_point_len_y = (double) (point_y[element_point_1[element_index]] - point_y[element_point_2[element_index]]);

element_point_len = sqrt(pow(element_point_len_x,2.0) + pow(element_point_len_y,2.0));

部品座標系での要素剛性マトリクスのための係数ｋ（ESmatrix_GK_k)は、各要素でのヤング率（element_young）と断面積element_areaの積をelement_point_lenで割ったものです。

このkから部品座標系での要素剛性マトリクス[GK}(ESmatrix_GK)は以下の様になります。

ESmatrix_GK_k = element_young[element_index] * element_area[element_index] / element_point_len;


// Set ESmatrix_GK
for (int i = 0; i < ESM_NUM; i++)
{
for (int j = 0; j < ESM_NUM; j++)
{
ESmatrix_GK[i][j]=0.0;
}
}

ESmatrix_GK[0][0]=ESmatrix_GK_k;
ESmatrix_GK[2][2]=ESmatrix_GK_k;
ESmatrix_GK[0][2]=-1*ESmatrix_GK_k;
ESmatrix_GK[2][0]=-1*ESmatrix_GK_k;


次に、部品座標系から全体座標系への座標変換マトリクス[T] （tranGtoE）を求めます。

cos_points = (double)(element_point_len_x/element_point_len);
sin_points = (double)(element_point_len_y/element_point_len);


// Set tranGtoE
for (int i = 0; i < ESM_NUM; i++)
{
for (int j = 0; j < ESM_NUM; j++)
{
ESmatrix_GK[i][j]=0.0;
}
}


　　　　　　　tranGtoE[element_index][0][0]=cos_points;
tranGtoE[element_index][1][1]=cos_points;
tranGtoE[element_index][2][2]=cos_points;
tranGtoE[element_index][3][3]=cos_points;
tranGtoE[element_index][0][1]=sin_points;
tranGtoE[element_index][1][0]=sin_points * -1;
tranGtoE[element_index][2][3]=sin_points;
tranGtoE[element_index][3][2]=sin_points * -1;




以上より全体座標系での要素剛性マトリクス[K] (ESmatrix_EK)は、以下の式から求められます。

//  Set ESmatrix_EK
double output_matrix[ESM_NUM][ESM_NUM];
bool ret=MxM(ESmatrix_GK, tranGtoE[element_index], output_matrix, ESM_NUM);
ret=WxM(tranGtoE[element_index], output_matrix, ESmatrix_EK[element_index], ESM_NUM);

ここでESM_NUMは、（プログラムの最初に定義していますが）

#define ESM_NUM 4  // 要素剛性マトリックスの成分の数 (2次元のため)


です。

これをすべての要素に対して算出し、最終的に全体の剛性マトリクスに加えこんで行くことにより、全体の剛性マトリクス（all_ESmatrix_EK）が作成されます。


なお
上記で使用している　MxM( ), WxM( )　のファンクションについては
<有限要素法>マトリクス計算のファンクション
を参照してください。

　　　　//--------------------------
// 全体の剛性マトリクス算出
//--------------------------
int all_element_point[ESM_NUM];
all_element_point[0] = element_point_1[element_index]*2;
all_element_point[1] = element_point_1[element_index]*2+1;
all_element_point[2] = element_point_2[element_index]*2;
all_element_point[3] = element_point_2[element_index]*2+1;

for (int i=0; i<ESM_NUM; i++) {
for (int j=0; j<ESM_NUM; j++) {
all_ESmatrix_EK[all_element_point[i]][all_element_point[j]]
=all_ESmatrix_EK[all_element_point[i]][all_element_point[j]]  + 　　　　ESmatrix_EK[element_index][i][j];
}
}

次のページ→

<有限要素法>要素剛性マトリクス算出 (まとめ）

「有限要素法における平面トラスの解析」プログラムをC言語で書いてみる : 目次