こんにちは、トライです。
以前、平面トラスの有限要素法（以下、FEM）の基本的な理論とMATLABの実際のコードを公開しました。
今回は第2弾として、「立体トラス」のFEMの基本的な理論を解説します。
平面トラスのような２次元のFEMも様々な面で役に立ちますが、実務レベルにまで落とし込もうとすると3次元のFEMとする必要があります。
今回は、平面トラスの記事で述べたFEMの基本的な理論は理解している前提で解説します。まだ平面トラスの記事を見ていない人はまずそちらからご覧ください。
また立体トラスについても平面トラスと同様に、実際に私が作成したMATLABコードをすべて公開します。
立体トラスのFEMのフローについて
上の図のように、節点 i と節点 j から構成されるトラス要素について考えます。
立体トラスのFEMのフローは平面トラスのフローと全く同じで、以下のようになります。
- 各部材の要素座標系における要素剛性行列を算出する
- 各部材の全体座標系における要素剛性行列を算出する【座標変換】
- 架構の全体剛性行列を算出する【アセンブル】
- 剛性方程式を解き、荷重に対する変形を算出する
- 変形から各部材の応力（軸力）を求める
要素座標系における要素剛性行列
平面トラスの回でも解説しましたが、ピン接合のトラス架構は
部材に曲げ応力やせん断応力は発生せず、軸力のみで外力に抵抗する架構です。
つまり、曲げ剛性やせん断剛性はなく、軸剛性のみを持つ架構ということです。
そのため、一次元トラス部材の剛性行列 kbar は以下のように表されます。
これが立体（三次元）トラスになると材軸方向（x方向）だけでなく、部材の強軸方向（y方向）と弱軸方向（z方向）の成分も出てきます。
ですがトラス部材の変形は材軸方向の伸び縮みのみで、部材の強軸方向と弱軸方向の変形（曲げ変形やせん断変形）は生じません。
そのため、要素座標系における剛性行列は以下の式になります。
uiy, uiz や ujy, ujz の係数となる行列の2,3行目と5,6行目、2,3列目と5,6列目はすべて0になります。
全体座標系における要素剛性行列【座標変換】
次に、先ほど求めた「要素座標系における要素剛性行列」を「全体座標系における要素剛性行列」に変換します。
各部材の座標軸の向きはバラバラなため、そのまま単純に足し合わせることはできません。
そのため、各部材の局所座標系で算出した剛性行列を架構全体の座標系（全体座標系）に変換します。
座標変換をするためには、座標変換行列 T を用います。
三次元座標において、ある軸に対する座標変換行列 T は以下の式で表されます。
トラス部材は i 端と j 端から構成されているので、上記の行列が対角に並びます。
よって、三次元トラスの座標変換行列 T とその転置行列 TT は以下のようになります。
局所座標系での剛性行列 kbar に左からTT、右から T を掛けると、全体座標系における要素剛性行列は以下のようになります。
式中の変数の定義を以下にまとめます。
- Cx = cosθx（全体座標系のX軸に対する方向余弦）
- Cy = cosθy（全体座標系のY軸に対する方向余弦）
- Cz = cosθz（全体座標系のZ軸に対する方向余弦）
式だけを見ると煩雑ですが、
- 1,4行目および1,4列目はcosθxの項
- 2,5行目および2,5列目はcosθyの項
- 3,6行目および3,6列目はcosθzの項
となっています。
つまり、1行1列目はcosθx・cosθx、1行2列目はcosθx・cosθyとなります。
立体トラスの新規関数コード
ここからは立体トラスのFEMのMATLABコードを公開します。
※今回は新規関数のみを公開し、次回メインプログラムを公開・解説します。
部材長さLを求める関数：SpaceTrussElementLength.m
i 端と j 端の節点の座標を入力し、立体トラスの部材の長さを出力する関数（function）
「SpaceTrussElementLength」を定義します。
function y = SpaceTrussElementLength(x1,y1,z1,x2,y2,z2)
y = sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1) );
endただの三次元の三平方の定理なので容易に理解できるかと思います。
要素剛性行列kを求める関数：SpaceTrussElementStiffness.m
次に、各部材の要素剛性行列を求めます。
- ヤング係数 E
- 断面積 A
- 部材の長さ L
- 全体座標系の各軸に対する部材の方向余弦 cosθx, cosθy, cosθz
を入力し、立体トラス部材の要素剛性行列を出力する関数（function）「SpaceTrussElementStiffness」を定義します。
function y = SpaceTrussElementStiffness(E,A,L,Cx,Cy,Cz)
w = [Cx*Cx Cx*Cy Cx*Cz ; Cy*Cx Cy*Cy Cy*Cz ; Cz*Cx Cz*Cy Cz*Cz];
y = E*A/L*[w -w ; -w w];
end全体剛性行列Kを求める関数：SpaceTrussAssemble.m
次に、全体座標系における各部材の要素剛性行列をアセンブル（結合）し、全体剛性行列を出力する関数（function）「SpaceTrussAssemble」を定義します。
function y = SpaceTrussAssemble(K,k,i,j)
K(3*i-2,3*i-2) = K(3*i-2,3*i-2) + k(1,1);
K(3*i-2,3*i-1) = K(3*i-2,3*i-1) + k(1,2);
K(3*i-2,3*i) = K(3*i-2,3*i) + k(1,3);
K(3*i-2,3*j-2) = K(3*i-2,3*j-2) + k(1,4);
K(3*i-2,3*j-1) = K(3*i-2,3*j-1) + k(1,5);
K(3*i-2,3*j) = K(3*i-2,3*j) + k(1,6);
K(3*i-1,3*i-2) = K(3*i-1,3*i-2) + k(2,1);
K(3*i-1,3*i-1) = K(3*i-1,3*i-1) + k(2,2);
K(3*i-1,3*i) = K(3*i-1,3*i) + k(2,3);
K(3*i-1,3*j-2) = K(3*i-1,3*j-2) + k(2,4);
K(3*i-1,3*j-1) = K(3*i-1,3*j-1) + k(2,5);
K(3*i-1,3*j) = K(3*i-1,3*j) + k(2,6);
K(3*i,3*i-2) = K(3*i,3*i-2) + k(3,1);
K(3*i,3*i-1) = K(3*i,3*i-1) + k(3,2);
K(3*i,3*i) = K(3*i,3*i) + k(3,3);
K(3*i,3*j-2) = K(3*i,3*j-2) + k(3,4);
K(3*i,3*j-1) = K(3*i,3*j-1) + k(3,5);
K(3*i,3*j) = K(3*i,3*j) + k(3,6);
K(3*j-2,3*i-2) = K(3*j-2,3*i-2) + k(4,1);
K(3*j-2,3*i-1) = K(3*j-2,3*i-1) + k(4,2);
K(3*j-2,3*i) = K(3*j-2,3*i) + k(4,3);
K(3*j-2,3*j-2) = K(3*j-2,3*j-2) + k(4,4);
K(3*j-2,3*j-1) = K(3*j-2,3*j-1) + k(4,5);
K(3*j-2,3*j) = K(3*j-2,3*j) + k(4,6);
K(3*j-1,3*i-2) = K(3*j-1,3*i-2) + k(5,1);
K(3*j-1,3*i-1) = K(3*j-1,3*i-1) + k(5,2);
K(3*j-1,3*i) = K(3*j-1,3*i) + k(5,3);
K(3*j-1,3*j-2) = K(3*j-1,3*j-2) + k(5,4);
K(3*j-1,3*j-1) = K(3*j-1,3*j-1) + k(5,5);
K(3*j-1,3*j) = K(3*j-1,3*j) + k(5,6);
K(3*j,3*i-2) = K(3*j,3*i-2) + k(6,1);
K(3*j,3*i-1) = K(3*j,3*i-1) + k(6,2);
K(3*j,3*i) = K(3*j,3*i) + k(6,3);
K(3*j,3*j-2) = K(3*j,3*j-2) + k(6,4);
K(3*j,3*j-1) = K(3*j,3*j-1) + k(6,5);
K(3*j,3*j) = K(3*j,3*j) + k(6,6);
y = K;
end部材の軸応力度σを求める関数：SpaceTrussElementStress.m
全体剛性行列 [K] が求まったら外荷重ベクトル { f } を定義し、剛性方程式 { f } = [K]{ u } を解いて変位ベクトル {u} を求めます。
- ヤング係数 E
- 部材の長さ L
- 全体座標系の各軸に対する部材の方向余弦 cosθx, cosθy, cosθz
- 変位ベクトル { u }
を入力し、立体トラス部材の軸応力度 σ を出力する関数（function）「SpaceTrussElementStress」を定義します。
function y = SpaceTrussElementStress(E,L,Cx,Cy,Cz,u)
y = E/L*[-Cx -Cy -Cz Cx Cy Cz]*u;
endまとめ：基本的な理論とフローは平面トラスと同じ
今回は立体トラスのFEMの基本的なフローを解説し、MATLABの新規関数コードを公開しました。
公開した新規関数コードは以下の手順でメインプログラムと同じフォルダに保存しておきましょう。
- 「新規スクリプト」をクリック
- エディタ画面で公開したコードをコピー＆ペースト
- 新規関数を保存 ※mファイル名は関数名と同じにする
次回、立体トラスのメインプログラムを公開し、その中身について解説します。
平面トラスの際はMATLABで用意されている関数を中心に解説しましたが、立体トラスではメインプログラムを作成した際の考え方に重きを置いて解説します。
コメント