こんにちは、トライです。
今回は立体トラスの有限要素法（以下、FEM）のメインプログラムのMATLABコードを公開し、コードの内容について解説します。
【注意】今回の記事で公開したメインプログラムのコードだけをコピーして解析実行ボタンを押してもプログラムは動きません。
以下の記事で公開した新規関数をあらかじめmファイルとしてメインプログラムと同じフォルダに保存しておいてください。
保存時の注意点などもまとめていますので、まだ以下の記事を読んでいない人はまずはこちらを読んでみて下さい。
立体トラスのFEMのMATLABコード公開
今回は以下のような立体トラスのFEMをMATLABで実施します。
例題の条件を整理します。
- 支持条件 節点１、節点２、節点３：ピン支持
- 部材のヤング係数 E = 200×106 [kN/m2]
- 部材断面積 部材①：0.001 [m2]、部材②：0.002 [m2]、部材③：0.001 [m2]
- 外力 節点４にX方向に12kNの集中荷重
この架構のFEM解析をMATLABで実施するためのコード例は以下のようになります。
clear
dim = 3; % 次元数(dimension:3次元)
dof = 3; % 節点の自由度数(degrees of freedom:立体トラスよりDX,DY,DZの3種類)
E = 200e6; % ヤング係数 単位[kN/m^2]
A = [1e-3; 2e-3; 1e-3]; % 断面積(Area) 単位[m^2]
XYZ = [0 -4 0; -3 0 0; 0 4 0; 0 0 5]; % 節点座標[x,y,z] 単位[m]
CON = [1 4; 2 4; 3 4]; % 接続条件(connection) [節点i, 節点j]
BOUN = [1 1 1 1; 2 1 1 1; 3 1 1 1]; % 境界条件 [節点番号, X方向拘束, Y方向拘束, Z方向拘束(0:自由, 1:固定)]
LOAD = [4 12 0 0]; % 外力[節点番号, X方向集中荷重, Y方向集中荷重, Z方向集中荷重] 単位[kN]
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) = SpaceTrussElementLength(XYZ(CON(i,1),1),XYZ(CON(i,1),2),XYZ(CON(i,1),3), ...
XYZ(CON(i,2),1),XYZ(CON(i,2),2),XYZ(CON(i,2),3));
end
L0 = L; % 初期部材長さL0
% 各部材の方向余弦Cxyzを定義 [cosθx cosθy cosθz]
Cxyz = zeros(nm,dim);
for i = 1:nm
for j = 1:dim
Cxyz(i,j) = ( XYZ(CON(i,2),j) - XYZ(CON(i,1),j) ) / L(i);
end
end
% 要素剛性行列ke、全体剛性行列Kを算出
K = zeros(nd);
for i = 1:nm
ke{i} = SpaceTrussElementStiffness(E,A(i),L(i),Cxyz(i,1),Cxyz(i,2),Cxyz(i,3));
K = SpaceTrussAssemble(K,ke{i},CON(i,1),CON(i,2));
end
% 外力ベクトルを定義
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;
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
% 軸応力度σと軸力Nを算出
for i = 1:nm
sigma(i) = SpaceTrussElementStress(E,L(i),Cxyz(i,1),Cxyz(i,2),Cxyz(i,3),u{i});
N(i) = sigma(i) * A(i);
end
% 座標更新
for i = 1:nn
for j = 1:dim
XYZ(i,j) = XYZ(i,j) + U(dof*(i-1) + j);
end
end
% 変形図を描写
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 );
Z0(i,j) = XYZ0( CON(i,j), 3 );
end
plot3(X0(i,:),Y0(i,:),Z0(i,:),'o-', 'color', 'blue','MarkerSize',8, ...
'MarkerFaceColor','white','MarkerEdgeColor','blue','LineWidth',4);
hold on
end
scale = 100; % 変形をscale倍して描写
plotXYZ = XYZ0 + scale * [XYZ - XYZ0];
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 );
Z(i,j) = plotXYZ( CON(i,j), 3 );
end
plot3(X(i,:),Y(i,:),Z(i,:),'o-', 'color', 'red','MarkerSize',8, ...
'MarkerFaceColor','white','MarkerEdgeColor','red','LineWidth',4);
hold on
end
% 節点番号を挿入
for i = 1:nn
text( XYZ0(i,1), XYZ0(i,2), XYZ0(i,3), ['n=',num2str(i)],'HorizontalAlignment','left',...
'VerticalAlignment','top','FontSize',30 );
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) ), ...
1/2*( XYZ0(CON(i,1),3) + XYZ0(CON(i,2),3) ),
['m=',num2str(i)],'HorizontalAlignment','left',...
'VerticalAlignment','middle','FontSize',30,'color','black' );
hold on
end
axis equal
xlabel('x(m)','FontSize',30); ylabel('y(m)','FontSize',30); zlabel('z(m)','FontSize',30);
grid on
set(gca,'FontSize',30);
hold off公開コードの解説
今回公開したコードについて、平面トラスの記事で解説していない部分や平面トラスと記述が異なる部分について解説します。
各部材の方向余弦cosθの算出方法
平面トラスの時は三角関数「atan」を用いて全体座標系に対する部材の角度θを求めました。
ですが、今回は
- 全体座標系に対する部材の角度θをX, Y, Z軸それぞれに対して求める必要があること
- 最終的に必要なものは角度θではなくcosθであること
- コードの書き方の正解はひとつではないこと
こういったことを考慮し、平面トラスとは違うコード内容としました。
【44-50行目】
% 各部材の方向余弦Cxyzを定義 [cosθx cosθy cosθz]
Cxyz = zeros(nm,dim);
for i = 1:nm
for j = 1:dim
Cxyz(i,j) = ( XYZ(CON(i,2),j) - XYZ(CON(i,1),j) ) / L(i);
end
end各部材の方向余弦を返す「Cxyz」を定義するときに頭の中で考えたことを説明します。
まずは「行列の行には何が入るか、列には何が入るか」といった完成形をイメージします。
今回の例題の場合、定義したい方向余弦行列 Cxyz の完成形はこのような感じになります。
このように完成形を書き出してみると以下のようなことが見えてきます。
- 行には部材番号が来る。すると行の数は部材の本数、つまりnm※1になるな。
- 列には各部材のcosθx, cosθy, cosθzが入る。ということは列の数は全体座標系の軸の本数、つまりdim※2になるな。
- だから行列Cxyzのサイズは(nm)×(dim)だ。
※1 nmは部材数（number of materials）今回は部材が3本なので nm = 3）
※2 dimは次元数（dimension）今回は立体つまり三次元なので dim = 3）
このように考え、最初にCxyzを「Cxyz = zeros(nm, dim);」と定義しました。
変数のサイズは最初に必ず定義しておかないといけないの？
定義しなくても解析自体は回ります。
ただし、forループが進むにつれてベクトルや行列のサイズが変わるようなコードだと計算速度が遅くなってしまいます。
変数の完成形をイメージする訓練という意味でも、最初に変数のサイズは定義しておきましょう。
plot3：三次元の図を描写
平面トラスではMATLABに標準装備されているコマンド「plot」を用いました。
それに対し今回のような三次元の架構に対しては「plot」ではなく「plot3」を用います。
【94-104行目】
% 変形図を描写
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 );
Z0(i,j) = XYZ0( CON(i,j), 3 );
end
plot3(X0(i,:),Y0(i,:),Z0(i,:),'o-', 'color', 'blue','MarkerSize',8, ...
'MarkerFaceColor','white','MarkerEdgeColor','blue','LineWidth',4);
hold on
endplot3…の行で最低限必要なコマンドは「plot3( X0(i, :), Y0(i, :), Z0(i, 🙂 )」のみです。
これ以降は図の見た目を整えるためのコマンドです。
今回公開したコードに記述したものを簡単に解説すると
- ‘o-‘：節点位置を〇、部材を線で描写する
- ‘color’, ‘blue’：部材を青色にする
- ‘MarkerSize’, 8：節点の〇のサイズを8ptにする
- ‘MarkerFaceColor’, ‘white’：節点の〇を白色で塗りつぶす
- ‘MarkerEdgeColor’, ‘blue’：節点の〇の線を青色にする
- ‘LineWidth’, 4：部材の線の太さを4ptにする
これらを設定すると、以下のような結果図が描写されます。
axis, label, grid, set：軸の設定
この他にも公開したコードには変形図を描写するのに多くのコマンドを使用しています。
ひとつずつ解説していきます。
【134-138行目】
axis equal
xlabel('x(m)','FontSize',30); ylabel('y(m)','FontSize',30); zlabel('z(m)','FontSize',30);
grid on
set(gca,'FontSize',30);
hold off- axis equal：各軸の幅（縦横比）を統一する。
- xlabel(‘x(m)’, ‘FontSize’,30)：X軸のタイトルを「x(m)」とし、フォントサイズを30ptにする。
- grid on：座標軸にグリッドラインを表示する。
- set(gca, ‘FontSize’, 30)：座標軸のフォントサイズを30ptにする。
- hold off：既存の図
私が変形図を描写する際は、これらのコマンドはほぼすべて毎回定義しています。
例えば「axis equal」を定義せずに変形図を描写するとこのような図が出てきます。
同じ1mのグリッドでもY軸のグリッドのほうが長く描かれていますね。
また、節点番号の表記「nn=3」の下半分が消えてしまっています。
実際の架構形状を正しく把握するためにも、また図中の表記がちゃんと表示されるようにするためにも、「axis equal」は常に入力しておくことを勧めます。
今回紹介したもの以外にも、MATLABの図の描写に関しては本当に様々なメニューが用意されています。
MATLABのヘルプセンターで「plot」や「plot3」と検索してみると
へぇー。こんなにたくさん表現の種類が用意されてるんだ。
と驚くと思います。ぜひ以下のリンクから一度確認してみて下さい。
text：節点番号や部材番号の挿入
結果図を描写後、図中に説明テキストを挿入します。
【118-133行目】
% 節点番号を挿入
for i = 1:nn
text( XYZ0(i,1), XYZ0(i,2), XYZ0(i,3), ['n=',num2str(i)],'HorizontalAlignment','left',...
'VerticalAlignment','top','FontSize',30 );
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) ), ...
1/2*( XYZ0(CON(i,1),3) + XYZ0(CON(i,2),3) ),
['m=',num2str(i)],'HorizontalAlignment','left',...
'VerticalAlignment','middle','FontSize',30,'color','black' );
hold on
end- text( XYZ0(i,1), XYZ0(i,2), XYZ0(i,3), …)：説明テキストを座標[ XYZ0(i,1), XYZ0(i,2), XYZ0(i,3) ]に配置する。
- [ ‘n=’, num2str(i) ]：「n=」＋「数値 i 」を表記する。（「n=1」「m=2」など）
- ‘HorizontalAlignment’, ‘left’：配置点がテキストの左端となるような位置に表示する。
- ‘VerticalAlignment’, ‘top’：配置点がテキストの上端となるような位置に表示する。
- ‘FontSize’,30：テキストのフォントサイズを30ptにする。
「HorizontalAlignment」は配置点に対するテキストの水平方向の配置を指定します。
配置の指定が ‘left’ なら配置点がテキストの左端になるようにテキストを配置します。
「VerticalAlignment」は配置点に対するテキストの垂直方向の配置を指定します。
配置の指定が ‘middle’ なら配置点がテキストの中央になるようにテキストを配置します。
配置位置のレパートリーはMATLABのヘルプセンターにまとめられているので抜粋しておきます。
例として平面トラスの記事で取り扱った架構でテキスト位置を変更した場合の例を示します。
左の図のコード
% 節点番号を挿入
for i = 1:nn
text( XYZ0(i,1), XYZ0(i,2), ['node',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) ), ...
['element',num2str(i)],'HorizontalAlignment','center',...
'VerticalAlignment','middle','FontSize',18,'color','black' );
hold on
end右の図のコード
% 節点番号を挿入
for i = 1:nn
text( XYZ0(i,1), XYZ0(i,2), ['node',num2str(i)],'HorizontalAlignment','center',...
'VerticalAlignment','middle','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) ), ...
['element',num2str(i)],'HorizontalAlignment','right',...
'VerticalAlignment','bottom','FontSize',18,'color','black' );
hold on
end
左の図では、部材の ‘VerticalAlignment’ を ‘middle’ としたために、部材の線と「element1」という文字が重なってしまっています。
逆に右の図では、節点の ‘HorizontalAlignment’ を ‘center’ としたために、節点の〇と「node1」という文字が重なってしまっています。
このように、MATLABでは図の表現方法はかなり細かく設定することができます。
本記事では紹介しきれなかった設定方法についてはMATLABヘルプセンターで「text」のプロパティを確認してみて下さい。
まとめ：コード作成時は変数の完成形をイメージしよう
今回は立体トラスのメインプログラムのコードを公開し、記述内容について解説しました。
繰り返しになりますが、コードを書く上で大切なことは
「各変数の完成形をイメージすること」
これに尽きます。
完成形を正しく理解できていないと正しいコードは書けないからです。
本ブログで私が作成したひとつの正解のコードを公開することで
- まだ自分でコードが書けない
- コードを書いてみたものの、解析が回らない。
- 自分で書いたコードのどこが間違っているのか分からない
といった人にとっての手掛かり、きっかけとなれば幸いです。
また今回でトラス架構の理論解説およびコードの公開は終わりです。
次回は軸力に加え、せん断力や曲げ応力も負担する「骨組架構」について紹介します。
コメント