こんにちは、トライです。
今回は平面骨組要素の有限要素法（以下、FEM）のメインプログラムのMATLABコードを公開し、コードの内容について解説します。
【注意】今回の記事で公開したメインプログラムのコードだけをコピーして解析実行ボタンを押してもプログラムは動きません。
前回の記事「新規関数編」で公開した新規関数をmファイルとしてメインプログラムと同じフォルダに保存しておいてください。
まだ前回の記事を読んでいない人はまずはこちらから
平面骨組要素のFEMのMATLABコード公開
前回の「新規関数編」でも登場した以下のラーメン架構に対してFEMを実施します。
解析モデルの条件は以下のように設定します。
- 部材のヤング係数 E：210×106 [kN/m2]
- 部材の断面積 A：部材①②③ともに0.02 [m2]
- 部材の断面二次モーメント I：部材①②③ともに0.00005 [m2]
- 支持条件：節点１および節点４が固定支持
- 外力 節点２にーX方向に20kN、節点３に時計回りに12kNの集中荷重
MATLABによるFEMのコード例は以下のようになります。
clear
dim = 2; % 次元数(dimension:2次元)
dof = 3; % 節点の自由度数(degrees of freedom:平面骨組よりDX,DY,θの3種類)
E = 210e6; % ヤング係数 単位[kN/m^2]
A = [ 2e-2; 2e-2; 2e-2 ]; % 断面積 単位[m^2]
I = [ 5e-5; 5e-5; 5e-5 ]; % 断面二次モーメント 単位[m^2]
XYZ = [ 0 0; 0 3; 4 3; 4 0 ]; % 節点座標[x,y] 単位[m]
CON = [ 1 2; 2 3; 3 4 ]; % 接続条件(connection) [節点i, 節点j]
BOUN = [ 1 1 1 1; 4 1 1 1 ]; % 境界条件(boundary condition) [節点番号, X方向変位, Y方向変位, 節点回転(0:自由, 1:固定)]
LOAD = [ 2 -20 0 0;
3 0 0 12 ]; % 外力[節点番号, X方向集中荷重, Y方向集中荷重, 節点モーメント荷重] 単位[kN,kNm]
nn = size(XYZ, 1); % 節点数(number of nodes)
nm = size(CON, 1); % 部材数(number of materials)
nd = dof * nn; % 総自由度数 number of degrees of freedom
nl = size(LOAD, 1); % 外力数 number of loads
XYZ0 = XYZ; % 初期座標
% 変形を許容する自由度を0、変形を拘束する自由度を1とするベクトルfixを定義
fix = zeros(nd,1);
for i = 1:size(BOUN,1)
for j = 1:dof
if BOUN(i,j+1) == 1
fix(dof*(BOUN(i,1)-1)+j) = 1;
end
end
end
% 変形を許容する自由度の番号を返すベクトルfreeを定義
icount = 0;
for i = 1:nd
if fix(i) == 0
icount = icount + 1;
free(icount) = i;
end
end
% 各部材の長さL[m]を算出
for i = 1:nm
L(i) = PlaneFrameElementLength(XYZ(CON(i,1),1),XYZ(CON(i,1),2), ...
XYZ(CON(i,2),1),XYZ(CON(i,2),2));
end
L0 = L; % 初期部材長さL0
% 各部材の要素座標軸と全体座標軸との角度θに基づくcosθ, sinθの算出
% C = [cosθ1 cosθ2 ... cosθnm], S = [sinθ1 sinθ2 ... sinθnm]
C = zeros(nm,1); S = zeros(nm,1);
for i = 1:nm
C(i) = ( XYZ(CON(i,2),1) - XYZ(CON(i,1),1) ) / L(i);
S(i) = ( XYZ(CON(i,2),2) - XYZ(CON(i,1),2) ) / L(i);
end
% 要素剛性行列ke, 全体剛性行列Kの算出
K = zeros(nd);
for i = 1:nm
ke{i} = PlaneFrameElementStiffness(E,A(i),I(i),L(i),C(i),S(i));
K = PlaneFrameAssemble(K,ke{i},CON(i,1),CON(i,2));
end
% 外力ベクトルfの算出
f = zeros(nd,1);
for i = 1:nl
f( dof*(LOAD(i,1))-2, 1 ) = LOAD(i,2);
f( dof*(LOAD(i,1))-1, 1 ) = LOAD(i,3);
f( dof*(LOAD(i,1)) , 1 ) = LOAD(i,4);
end
% 変位を許容した自由度に対して剛性方程式を解き、変位ベクトルUと節点力ベクトルFを求める
U = zeros(nd,1);
U(free,1) = K(free,free)\f(free);
F = K * U;
% 各部材の変位ベクトルuの算出
% u = [節点iのX方向変位 節点iのY方向変位 節点iの回転角 節点jのX方向変位 節点jのY方向変位 節点jの回転角]'
for i = 1:nm
u{i} = [ U(dof*CON(i,1)-2);
U(dof*CON(i,1)-1);
U(dof*CON(i,1)) ;
U(dof*CON(i,2)-2);
U(dof*CON(i,2)-1);
U(dof*CON(i,2)) ];
end
% 各部材の応力ベクトルForceの算出
% Force = [ 節点i位置の軸力 節点i位置のせん断力 節点i位置の曲げモーメント 節点j位置の軸力 節点j位置のせん断力 節点j位置の曲げモーメント]'
for i = 1:nm
Force{i} = PlaneFrameElementForces(E,A(i),I(i),L(i),C(i),S(i),u{i});
end
% 節点座標の更新
for i = 1:nn
for j = 1:dim
XYZ(i,j) = XYZ(i,j) + U(dof*(i-1) + j);
end
end
% 変形前と変形前の架構形状を描画
figure(1);
cla
for i=1:nm
for j = 1:size(CON,2)
X0(i,j) = XYZ0((CON(i,j)),1);
Y0(i,j) = XYZ0((CON(i,j)),2);
end
plot(X0(i,:),Y0(i,:),'b-','LineWidth',2);
hold on
end
scale = 100; % 変形をscale倍して描画
plotXYZ = XYZ0 + scale * [XYZ - XYZ0]; % 変形図の描写用にplotXYZを定義
for i=1:nm
for j = 1:size(CON,2)
X(i,j) = plotXYZ((CON(i,j)),1);
Y(i,j) = plotXYZ((CON(i,j)),2);
end
plot(X(i,:),Y(i,:),'r-','LineWidth',2);
hold on
end
% 節点番号を挿入
for i = 1:nn
text( XYZ0(i,1), XYZ0(i,2), ['nn=',num2str(i)],'HorizontalAlignment','left',...
'VerticalAlignment','top','FontSize',18 );
hold on
end
% 要素番号を挿入
for i = 1:nm
text( 1/2*( XYZ0(CON(i,1),1) + XYZ0(CON(i,2),1) ), ...
1/2*( XYZ0(CON(i,1),2) + XYZ0(CON(i,2),2) ), ...
['nm=',num2str(i)],'HorizontalAlignment','center',...
'VerticalAlignment','middle','FontSize',18,'color','black' );
hold on
end
axis equal;
xlabel('x(m)','FontSize',18); ylabel('y(m)','FontSize',18);
title('Plane Frame Analysis','FontSize',18)
set(gca,'FontSize',18);
hold off
公開コードの解説
FEMフローの確認
前回の記事でも述べた平面骨組のFEMのフローを再確認します。
- 各部材の要素座標系における要素剛性行列を算出する
- 各部材の全体座標系における要素剛性行列を算出する【座標変換】
- 架構の全体剛性行列を算出する【アセンブル】
- 剛性方程式を解き、荷重に対する変形を算出する
- 変形から各部材の応力（軸応力・せん断応力・曲げ応力）を求める
FEMフローと公開コードの対応
ここからはFEMのフローと公開コードの対応関係を確認します。
公開したコード内では
「1. 各部材の要素座標系における要素剛性行列を算出する」と
「2. 各部材の全体座標系における要素剛性行列の算出する【座標変換】」
はまとめて実施しています。
対応コードは1行目の「clear」から58行目の
「ke{i} = PlaneFrameElementStiffness(E,A(i),I(i),L(i),C(i),S(i));」までです。
はじめに要素剛性行列を構成する各係数（E, A, I, L）を設定し、各部材のsinθとcosθを求め、要素剛性行列を算出しています。
「3. 架構の全体剛性行列を算出する【アセンブル】」に対応したコードは
59行目の「K = PlaneFrameAssemble(K,ke{i},CON(i,1),CON(i,2));」です。
新規関数編で全体剛性行列を作成する新規関数を定義したおかげで、メインプログラムではこの1行を for ループの中に入れるだけで全体剛性行列が求まります。
「4. 剛性方程式を解き、荷重に対する変形を算出する」に対応したコードは
変形を許容する自由度を返すベクトル free を定義した21行目から38行目までの範囲と、
62行目から72行目の「U(free,1) = K(free,free)\f(free);」です。
FEMの剛性方程式は、変形を許容する自由度に対してのみを抽出して解きます。
「5. 変形から各部材の応力（軸応力・せん断応力・曲げ応力）を求める」 に対応したコードは
各部材の変位ベクトル u を算出した75行目から84行目までの範囲と、
各部材の応力ベクトル Force を求めた86行目から90行目までです。
公開コードの検証
ここからは、本記事で公開したコードに誤りがないかを検証していきます。
先ほどの例は3次不静定架構（反力数6＋部材数3＋剛接合数2ー2×節点数4＝3）のため、手計算で応力を求めるのは骨が折れます。
そのため、以下のような静定架構（反力数3＋部材数3＋剛接合数2ー2×節点数4＝0）でコードの検証を実施します。
手計算による応力の算出
一級建築士の学科Ⅳ（建築構造）で出題される構造力学の問題もこれくらいのレベルです。
一級建築士試験を受験する予定の人は一度ご自身でも解いてみましょう。
はじめに支点の反力を求めます。
X方向の力のつり合い式、Y方向の力のつり合い式、節点1回りのモーメントのつり合い式は以下のようになります。
ΣX = 0：H1 – 10kN = 0 ∴H1 = 10kN
ΣY = 0：V1 + V5 – 20kN = 0
ΣM1 = 0：-10kN・3m + 20kN・2m – V5・4m = 0 ∴V5 = 10kN, V1 = 17.5kN
次に各部材の応力を求めます。
断面切断位置における力のつり合い式
ΣX = 0：Q + 10kN = 0 ∴Q = -10 (kN)
ΣY = 0：N + 17.5kN = 0 ∴N = -17.5 (kN)
ΣMx = 0：-M – 10kN・xm = 0 ∴M = -10x (kNm)
断面切断位置における力のつり合い式
ΣX = 0：N – 10kN + 10kN = 0 ∴N = 0 (kN)
ΣY = 0：-Q + 17.5kN = 0 ∴Q = -17.5 (kN)
ΣMx = 0：-M – 10kN・3m + 17.5kN・xm = 0 ∴M = 17.5x – 30 (kNm)
断面切断位置における力のつり合い式
ΣX = 0：Q = 0 ∴Q = 0 (kN)
ΣY = 0：N + 2.5kN = 0 ∴N = -2.5 (kN)
ΣMx = 0：M = 0 ∴M = 0 (kNm)
力のつり合い式
ΣX = 0：-N = 0 ∴N = 0 (kN)
ΣY = 0：Q + 2.5kN = 0 ∴Q = -2.5 (kN)
ΣMx = 0：M – 2.5kN・xm = 0 ∴M = 2.5x (kNm)
以上により求められた結果を応力図としてまとめます。
FEM解析結果の確認
手計算の結果を踏まえて、次にFEM解析によって得られた結果を確認します。
はじめに公開したコードから、以下の部分を変更しました。
【はじめに公開したコードの7-13行目】
A = [ 2e-2; 2e-2; 2e-2 ]; % 断面積 単位[m^2]
I = [ 5e-5; 5e-5; 5e-5 ]; % 断面二次モーメント 単位[m^2]
XYZ = [ 0 0; 0 3; 4 3; 4 0 ]; % 節点座標[x,y] 単位[m]
CON = [ 1 2; 2 3; 3 4 ]; % 接続条件(connection) [節点i, 節点j]
BOUN = [ 1 1 1 1; 4 1 1 1 ]; % 境界条件(boundary condition) [節点番号, X方向変位, Y方向変位, 節点回転(0:自由, 1:固定)]
LOAD = [ 2 -20 0 0;
3 0 0 12 ]; % 外力[節点番号, X方向集中荷重, Y方向集中荷重, 節点モーメント荷重] 単位[kN,kNm]【検証用コードの7-13行目】
A = [ 2e-2; 2e-2; 2e-2; 2e-2 ]; % 断面積 単位[m^2]
I = [ 5e-5; 5e-5; 5e-5; 5e-5 ]; % 断面二次モーメント 単位[m^2]
XYZ = [ 0 0; 0 3; 2 3; 4 3; 4 0 ]; % 節点座標[x,y] 単位[m]
CON = [ 1 2; 2 3; 3 4; 4 5 ]; % 接続条件(connection) [節点i, 節点j]
BOUN = [ 1 1 1 0; 5 0 1 0 ]; % 境界条件(boundary condition) [節点番号, X方向変位, Y方向変位, 節点回転(0:自由, 1:固定)]
LOAD = [ 2 -10 0 0;
3 0 -20 0 ]; % 外力[節点番号, X方向集中荷重, Y方向集中荷重, 節点モーメント荷重] 単位[kN,kNm]プログラムを回し、得られた結果について確認していきましょう。
まず節点力ベクトル F を確認してみると
支点の反力については
- 節点1のX方向の節点力が10kN
- 節点1のY方向の節点力が17.5kN
- 節点5のY方向の節点力が2.5kN
となっているので、手計算の結果と合っていますね。
10-14みたいな非常に小さな値は≒0と考えてOKです。
外力については
- 節点2のX方向の節点力が-10kN
- 節点3のY方向の節点力が-20kN
になっていて、入力した外力と一致しているね。
次に部材の応力ベクトル Force を確認すると
例として曲げモーメントMの項を見てみると
- 部材①の節点1でM=0kNm、節点2でM=30kNm
- 部材②の節点2でM=30kNm、節点3でM=5kNm
- 部材③の節点3でM=5kNm、節点4でM=0kNm
- 部材④は節点4, 5ともにM=0kNm
となっており、手計算の結果と一致していますね。
軸力N、せん断力Qについても手計算の結果と一致しています。
最後に変形図も確認しておきます。（節点番号と部材番号の表記も少し変更しました。）
節点1はピン支持なので変位なし、節点5はローラー支持なのでX方向に変位しています。
変形図を確認してもおかしなところはありませんね。
まとめ
今回は平面骨組要素のFEMのメインプログラムのMATLABコードを公開し、FEMフローとの対応確認、手計算による検証、FEMの結果の見方について解説しました。
繰り返しになりますが、骨組要素（梁要素）はFEMで最も重要な要素のひとつです。
基本的な理論とFEMのフローはトラスでも梁でも同じです。
また、一度理解してしまえばそれほど難しい話でもないので、なんとなく理解しておくと良いと思います。
次回は最も汎用性の高い立体骨組要素（三次元梁要素）の解説記事を作成予定です。
そちらもぜひご覧下さい。
コメント