こんにちは、トライです。
今回は平面トラスの最終回として、平面トラスの有限要素解析のメインプログラムコードを載せ、コードの内容について解説します。
私の経験から、プログラミングの勉強やコードの作成でつまずく理由のひとつとして
「完成形がイメージできていない」ことが挙げられます。
ここでいう完成形とは、例えば
この行列は〇〇が行に来て●●が列に来るから、〇行●列のサイズの行列になるな
といった感じです。
計算自体はコンピュータがやってくれますが、コンピュータに指示を出すのは人間です。
そしてコンピュータに指示を出す手段がプログラミングです。
指示を出す人間がプログラムの完成形、行列で言えばサイズがいくつで、どのような種類の値がどの位置になければ来なければならないのかを理解していないと正しいコードは書けません。
そんなこと言ったって勉強し始めたばかりで、完成形のイメージがつかないよ…
という人のためにも、本記事で私が作成したMATLABのコードを無料公開します。
前回の記事と合わせれば
「ここに書かれているコードをコピー&ペーストすればプログラムが動く」
というものになっていますので、ぜひ最後までお読み下さい。
平面トラスの有限要素解析のMATLABコード公開
以下のモデルに対してMATLABを用いて有限要素解析を実施します。
早速ですが、私が作成した平面トラスの有限要素解析のMATLABメインプログラムコードを公開します。
※2021.6.19 13:55 部材角θの算出方法を修正、変形図を出力する変数 plotXYZ を追加
clear
dim = 2; % 次元数
dof = 2; % 自由度数（今回はDX,DYの2種類）
E = 210e6; % ヤング係数 単位[kN/m^2]
A = [1e-4; 1e-4; 1e-4]; % 断面積 単位[m^2]
XYZ = [0 0; 4 0; 2 3]; % 節点座標[x,y] 単位[m]
CON = [1 2; 1 3; 2 3]; % 接続条件 [節点i, 節点j]
BOUN = [1 1 1; 2 0 1]; % 境界条件 [節点番号, X方向の拘束, Y方向の拘束(0:自由, 1:固定)]
LOAD = [3 5 -10]; % 集中荷重 [節点番号, X方向集中荷重, Y方向集中荷重] 単位[kN]
nn = size(XYZ, 1); % 総節点数
nm = size(CON, 1); % 総部材数
nd = dof * nn; % 総自由度数
nl = size(LOAD, 1); % 総外力数
XYZ0 = XYZ; % 変形前の節点座標
% 変形を拘束する自由度の算出
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
% 変形を許容する自由度の算出
icount = 0;
for i = 1:nd
if fix(i) == 0
icount = icount + 1;
free(icount) = i;
end
end
% 部材長さLの算出
for i = 1:nm
L(i) = PlaneTrussElementLength(XYZ(CON(i,1),1),XYZ(CON(i,1),2), ...
XYZ(CON(i,2),1),XYZ(CON(i,2),2));
end
L0 = L; % 初期部材長さL0
% 全体座標系に対する部材角θの算出
for i = 1:nm
theta(i) = atan( (XYZ(CON(i,2),2)-XYZ(CON(i,1),2)) /...
(XYZ(CON(i,2),1)-XYZ(CON(i,1),1)) )*180/pi;
end
% 要素剛性行列ke, 全体剛性行列Kの算出
K = zeros(nd);
for i = 1:nm
ke{i} = PlaneTrussElementStiffness(E,A(i),L(i),theta(i));
K = PlaneTrussAssemble(K,ke{i},CON(i,1),CON(i,2));
end
% 外力ベクトルfの算出
f = zeros(nd,1);
for i = 1:nl
f( dof*(LOAD(i,1))-1, 1 ) = LOAD(i,2);
f( dof*(LOAD(i,1)) , 1 ) = LOAD(i,3);
end
% 変形を許容する自由度を抽出し、剛性方程式を解く
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)-1);
U(dof*CON(i,1)) ;
U(dof*CON(i,2)-1);
U(dof*CON(i,2)) ];
end
% 各部材の応力度σ, 軸力Nの算出
for i = 1:nm
sigma(i) = PlaneTrussElementStress(E,L(i),theta(i),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
% 変形前と変形後の架構形状を描画
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,:),'o-', 'color', 'blue','MarkerSize',8, ...
'MarkerFaceColor','white','MarkerEdgeColor','blue','LineWidth',2);
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);
end
plot(X(i,:),Y(i,:),'o-', 'color', 'red','MarkerSize',8, ...
'MarkerFaceColor','white','MarkerEdgeColor','red','LineWidth',2);
hold on
end
% 節点番号を挿入
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
axis equal;
xlabel('x(m)','FontSize',18); ylabel('y(m)','FontSize',18);
xlim([-0.5 4.5]); ylim([-0.5 3.5]);
title('Plane Truss Analysis','FontSize',18)
set(gca,'FontSize',18);
hold off公開コードの解説
ここからは先ほどのメインプログラムについて解説していきます。
すべてのコマンドを解説するのは難しいので、今回解説していないものの中で使い方が分からないコマンドについては
「MATLAB＋スペース＋コマンド名」
でGoogle検索すればMathWorksのヘルプセンターに解説が載っているので確認してみて下さい。
変形を許容する自由度を算出するコード
まずはじめに、新規関数編で触れなかった
「変形が許容されている自由度を抽出し、剛性方程式を解く」
ためのコードについて解説します。
【19-36行目】
% 変形を拘束する自由度の算出
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
% 変形を許容する自由度の算出
icount = 0;
for i = 1:nd
if fix(i) == 0
icount = icount + 1;
free(icount) = i;
end
endこのコードは以下のような考え方で作成しました。
- 変形可能な行を0、変形が拘束されている行を1とするベクトルfixを作成
- ベクトルfixで、0が存在する行番号（何行目か）をベクトルfreeに記録する
1. により、以下のようなベクトルfixが得られます。
fix =
1
1
0
1
0
0境界条件BOUNで「節点1はX方向とY方向を拘束、節点2はY方向を拘束」と指定したので、ベクトルfixの1, 2, 4行目が「1」、それ以外が「0」となっています。
2. 「0」が存在する行は3, 5, 6行目なので、[ 3 5 6 ] というベクトルfreeを返すコードを記述します。
実際に得られたベクトルfreeを確認すると以下のようになっています。
free =
3 5 6clear：ワークスペースから変数データを削除する
何回か続けて解析を実行すると、前に実行したプログラムの変数データが残っていたために、後で実行した解析がうまくいかなくなることがあります。
そういったエラーを回避するために、各解析の実行前にすべての変数データをクリアしておく必要があります。
そのためのコマンドが「clear」です。
コマンド「clear」はメインプログラムの最初につける癖をつけておきましょう。
実際にコマンド「clear」を実行すると以下のようにワークスペースに記録されていた変数がすべて削除されます。
全部じゃなくて特定の変数だけ削除することはできないの？
「clear＋半角スペース＋削除したい変数」と入力すれば特定の変数だけ削除することができます。
例えば、断面積 A のみをワークスペースから削除したい場合は、コマンドウインドウで
「clear A」と入力し、Enterキーを押します。
すると先ほどまでワークスペースの一番上にあった変数Aのみ削除されています。
%：コメントの挿入
公開したコード内では、以下のように各変数の前の行もしくは右にコード内にコメントを挿入しています。
【3行目】
dim = 2; % 次元数MATLABではプログラムの解説やコメント、メモ等を挿入する場合、先頭に「%」を付けます。
「%」以降は文字が緑色になり、プログラムには影響しない文字列となります。
コード作成時に自分の頭を整理するため、また時間を空けて見た際にすぐに書いてある内容が理解できるように、定義した変数の説明や、何をするためのコードなのかは必ず加えておきましょう。
…：コードを複数行に分けて記載
1行あたりのコードが長すぎる場合は「…」（ピリオド3個）というコマンドを使って複数行にコードを分割することができます。
【38-42行目】
% 部材長さLの算出
for i = 1:nm
L(i) = PlaneTrussElementLength(XYZ(CON(i,1),1),XYZ(CON(i,1),2), ...
XYZ(CON(i,2),1),XYZ(CON(i,2),2));
endコードが長すぎると全体を見るためにいちいち横にスクロールしなければなりません。
画面内にすべてのコードが納まるように「…」を活用しましょう。
{ }：セル配列
MATLABはMatrix Laboratory というだけあり、基本的には行列計算です。例えば、
あるベクトル A の i 行目を返したいときは「A( i ) 」
ある行列 B の i 行 j 列目を返したいときは「 B( i, j ) 」
と記述します。
上記のほかにもMATLABでは「ある行のある列の値」だけでなく、「ある行列そのもの」を返すこともできます。
例えば今回の場合、部材が3本なので、要素剛性行列が3つ存在します。
このうち、どれかひとつの要素剛性行列の値を確認したい場合、参考にしているテキストでは
- 要素①の要素剛性行列：k1
- 要素②の要素剛性行列：k2
- 要素③の要素剛性行列：k3
と変数名を分けていますが、これでは正直、汎用性がありません。
なので今回はMATLABのコマンドのひとつ「変数 { i } 」を用いました。
今回のコードに対し、コマンドウィンドウで 「ke」 と入力し、Enterキーを押してみると
ke =
{
[1,1] =
5250 0 -5250 0
0 0 0 0
-5250 0 5250 0
0 0 0 0
[1,2] =
1792.1 2688.2 -1792.1 -2688.2
2688.2 4032.2 -2688.2 -4032.2
-1792.1 -2688.2 1792.1 2688.2
-2688.2 -4032.2 2688.2 4032.2
[1,3] =
1792.1 -2688.2 -1792.1 2688.2
-2688.2 4032.2 2688.2 -4032.2
-1792.1 2688.2 1792.1 -2688.2
2688.2 -4032.2 -2688.2 4032.2
}と返ってきます。これはどういうことかというと、
- 1番目（部材①）の要素剛性行列が ke{1}（コマンドウィンドウ上では [1,1] ）
- 2番目（部材②）の要素剛性行列が ke{2}（コマンドウィンドウ上では [1,2] ）
- 3番目（部材③）の要素剛性行列が ke{3}（コマンドウィンドウ上では [1,3] ）
ということを表しています。
また、コマンドウィンドウ内で「ke{1}」、「ke{2}」、「ke{3}」と入力してもそれぞれ以下のように返ってきます。
>> ke{1}
ans =
5250 0 -5250 0
0 0 0 0
-5250 0 5250 0
0 0 0 0
>> ke{2}
ans =
1792.1 2688.2 -1792.1 -2688.2
2688.2 4032.2 -2688.2 -4032.2
-1792.1 -2688.2 1792.1 2688.2
-2688.2 -4032.2 2688.2 4032.2
>> ke{3}
ans =
1792.1 -2688.2 -1792.1 2688.2
-2688.2 4032.2 2688.2 -4032.2
-1792.1 2688.2 1792.1 -2688.2
2688.2 -4032.2 -2688.2 4032.2plot：解析結果の2次元プロット
次に設定した解析条件や解析結果がどうなっているかを確認するために、結果を図示（プロット）します。
今回公開したコードでは、以下のような結果図が出力されます。
ちゃんと節点１はX, Y方向、節点2はY方向の変位が拘束されていることが視覚的に確認できましたね。
今回はコマンドの紹介も兼ねて出来るだけ多くのコマンドを盛り込みました。
実際、公開したコードのうち、91-137行目までは結果図を出力するためのコードです。
図の出力についてはMATLABは本当にたくさんのコマンドが用意されているため、今後も例を挙げながら少しずつ紹介していきます。
解析結果の検証
最後に作成したプログラムから得られた結果が正しいかを検証しておきます。
手計算で解くと支点の反力や部材の軸力は以下のように求められます。
以上で手計算による反力や軸力の算出が終わりました。
作成したプログラムでもコマンドウィンドウで反力ベクトル F と軸力ベクトル N の値を確認してみましょう。
※解説を緑文字で入れました。
F =
-5.0000e+00 ←節点1のX方向反力
1.2500e+00 ←節点1のX方向反力
8.8818e-16 ←節点2のX方向反力（≒0）
8.7500e+00 ←節点2のX方向反力
5.0000e+00 ←節点3のX方向集中荷重
-1.0000e+01 ←節点3のX方向集中荷重
N =
5.8333 -1.5023 -10.5162 ←左から部材1, 部材2, 部材3の軸力手計算で算出した値とバッチリ合っていますね。
あれ？節点2のX方向反力が0じゃないけど。
MATLABで連立方程式を解くと、どうしてもピッタリ0とはならず、非常に小さい値（ゴミ値）が残ってしまうことがあります。
ですが、例えば10のマイナス10乗といったような非常に小さい値なので、実質的には0と判断しても全く問題ありません。
まとめ：具体的なコードを参考に効率的な学習を
今回公開したコードはあくまで私が執筆時点で
- 汎用性のあるコードか（定数だけ変えれば他の架構形式にも適用可能か）
- 無駄なコードがないか
といった点を意識して書いたものです。
当然、これが最適な記述方法だとは思いませんが、極力分かりやすく記述したつもりです。
MATLABの機能も含めて、学習ポイントを以下にまとめます。
- ベクトルや行列の完成形をイメージしてコードを書く
- メインプログラムの1行目には「clear」をつけておく
- 先頭に「%」を付けるとそれ以降はコメント化（文字は緑色）される
- 「…」を付けるとコードを複数行に分けて記述できる
- 解析結果を「plot」で視覚的に確認する
- 簡単な架構なら手計算でプログラムの結果が正しいか検証する
もし掲載したMATLABコードに関して質問があればコメント欄に記載いただければ回答させていただきます。
私はプロのプログラマーではないので、あまりに専門的な質問は答えられないかもしれません。その点はご了承下さい。
ちなみにMATLABの公式サイトでもユーザーからの質問や公式サイトからの回答が閲覧できます。
こちらも非常に参考になるのでお気に入り登録をしておくと良いでしょう。
最後までお読みいただき、ありがとうございました。
コメント