【MATLABコード公開】平面骨組要素の有限要素法①～新規関数編～

こんにちは、トライです。

以前、平面トラスおよび立体トラスの有限要素法（以下、FEM）の基本的な理論を解説し、MATLABコードを公開しました。

本記事ではFEMの第3弾として、平面骨組要素、いわゆる梁要素のFEMの基本的な理論と新規作成の関数について解説します。

一般的に建物の構造体は軸力のみを伝達するトラスではなく、軸力・せん断力・曲げ応力を伝達する骨組で構成されています。

骨組要素（梁要素）はFEMで最も重要な要素のひとつです。

トライ

構造設計者が行う建物の構造解析も、この骨組要素を中心に建物をモデル化します。

FEMでよく使用される要素TOP3

第1位：骨組要素（＝梁要素）
第2位：板要素
第3位：トラス要素（以前の記事で解説＆コード公開済）

今回は２次元の平面骨組要素のFEMの基本理論を解説し、FEMの際に使用する新規関数について解説します。

FEMの基本理論については過去2回解説していますが、その際に

「もう少し詳しく述べるべきだったかなぁ」

と反省した点については詳しい解説を追加しました。

また今回も平面トラスや立体トラスの回と同様に、実際に私が作成したMATLABコードをすべて公開します。

次回公開予定のメインプログラム編とあわせれば平面骨組要素のFEMプログラムを動かせるようになります。

ぜひ最後までお読みください。

平面骨組要素のFEMのフローについて
要素座標系における要素剛性行列
全体座標系における要素剛性行列【座標変換】
剛性方程式を解き、荷重に対する変形を算出する
変形から各部材の応力（軸応力・せん断応力・曲げ応力）を求める
平面骨組要素の新規関数コード
まとめ：骨組要素（梁要素）はFEMで最も重要な要素

平面骨組要素のFEMのフローについて

上の図のように、節点 i と節点 j から構成される平面骨組要素について考えます。

図中の各記号の意味を以下にまとめます。

x, y軸：要素座標軸
X, Y軸：全体座標軸
E：ヤング係数
A：断面積
I：断面二次モーメント
L：部材長さ
θ：要素座標軸と全体座標軸のなす角度

平面骨組要素のFEMのフローは平面トラスや立体トラスのフローと全く同じです。

平面骨組のFEMのフロー

各部材の要素座標系における要素剛性行列を算出する
各部材の全体座標系における要素剛性行列を算出する【座標変換】
架構の全体剛性行列を算出する【アセンブル】
剛性方程式を解き、荷重に対する変形を算出する
変形から各部材の応力（軸応力・せん断応力・曲げ応力）を求める

トライ

FEMのフローについては今回で3回目なのでだいぶ定着してきたと思います。

平面トラスや立体トラスの記事をまだ見ていない人は是非そちらもご覧ください。

要素座標系における要素剛性行列

骨組要素（梁要素）はトラス要素とは違い、以下の3種類の力を伝達できる要素です。

軸力
せん断力
曲げ応力

そのため、要素座標系における二次元骨組要素の剛性行列 k_bar は以下のように表されます。

トライ

k_bar の導出過程については、「たわみ角法」というやや複雑な理論を用いる必要があるので本記事では割愛します。

この剛性行列 k_bar を用いて、要素座標系における剛性方程式は以下のように表されます。

N：軸応力、Q：せん断応力、M：曲げ応力

トライ

トラスのときよりもかなり式が複雑ですが、とりあえず

材軸方向の剛性に関する係数はEA/Lでトラスと同じ
材軸直交方向の剛性と節点の回転剛性に関する係数が追加された

くらいの理解でOKです。

12EI/L³や6EI/L²の根拠についても「たわみ角法」を用いなければならないので本記事では割愛します。

全体座標系における要素剛性行列【座標変換】

次に、先ほど求めた「要素座標系における要素剛性行列」を「全体座標系における要素剛性行列」に変換します。

最終的には各部材の要素剛性行列を足し合わせて全体剛性行列を求めたいのですが、

各部材の要素座標軸の向きはバラバラなため、単純に足し合わせることはできません。

そのため、各部材の要素座標系で算出した剛性行列を架構全体の座標系（全体座標系）に変換します。

この作業を座標変換と言います。

座標変換をするためには、座標変換行列 T を用います。

平面骨組要素における座標変換行列 T は以下の式で表されます。

部材は i 端と j 端から構成されているので、高校数学Cで習う回転行列の「θ」に「ーθ」を代入した行列

[ cosθ　sinθ　0 ;　-sinθ　cosθ　0 ;　0　0　1 ] が対角に並びます。

局所座標系での剛性行列 k_bar に左からT^T、右から T を掛けると、全体座標系における要素剛性行列 k が求められます。

トライ

１つにまとめると行列が非常に複雑になるので、今回は２つに分けて表記してあります。

この後に公開する新規関数のコードも見やすくするためにこのような表記にしてあります。

T^T は座標変換行列 T の転置（てんち）行列を意味します。

転置行列とは、ある行列の行と列を入れ換えた行列のことを言います。

例として、2行3列の行列 A の転置行列 A^T は以下のように3行2列の行列となります。

トライ

この転置行列は、MATLABでは変数の後ろに「’」（カンマ）を付けるだけで簡単に計算してくれます。

トライ

またFEMのテキストでは、行数が多いベクトルをそのまま表記すると本のスペースを多く取ってしまうため、この転置行列の表現が多く用いられています。

例えば、n行1列の行ベクトル h について、FEMのテキストでは以下のように表記されることが多いです。

FEMのテキストでの転置表現の使用例

剛性方程式を解き、荷重に対する変形を算出する

全体剛性行列 K が求まったら、次に全体座標系における荷重ベクトル f（以下、荷重ベクトル f）を定義します。

具体的な荷重ベクトル f については、以下のモデルを例に解説します。

まず荷重ベクトル f のサイズについてですが

平面骨組要素は1つの節点の自由度が3（材軸方向、材軸直交方向、回転方向）
このモデルは節点数が4つ

以上より、荷重ベクトル f は3×4=12行の行ベクトルになります。

荷重ベクトル f のサイズは理解したので、各行にどのような値が入るかを解説します。

荷重ベクトル f の構成は以下のようになっています。

f = [ 節点1のX方向集中荷重;　節点1のY方向集中荷重;　節点1のモーメント荷重;　節点2のX方向集中荷重;　・・・　;　節点4のモーメント荷重 ]

今回の例の場合、

節点2のーX方向に20kNの集中荷重
節点3に時計回りに12kNmのモーメント荷重

が掛かっているため、荷重ベクトル f は

f = [ 0;　0;　0;　-20;　0;　0;　0;　0;　12;　0;　0;　0 ]

となります。

トライ

値がプラスかマイナスかについては、以下のルールを理解しておきましょう。

集中荷重は全体座標軸の方向と同じならプラス、反対ならマイナス
モーメント荷重は時計回りがプラス、反時計回りがマイナス

荷重ベクトル f が定義できたら、剛性方程式 f = Ku を解いて、全体座標系における変位ベクトル u（以下、変位ベクトル u）を求めます。

MATLABでは組み込み関数「\」（もしくは「￥」）を用いて

u = K \ f　もしくは　 u = K ￥ f

と入力することで、変位ベクトル u を求めることができます。

トライ

実際には節点2や節点3のように、変形や回転が拘束されていない自由度のみを抽出して剛性方程式を解きます。

詳しくは次回のメインプログラム編もしくは過去に公開した平面トラスや立体トラスのメインプログラム編をご覧ください。

変位ベクトル u の構成は以下のようになっています。

u = [ 節点1のX方向変位;　節点1のY方向変位;　節点1の回転角;　節点2のX方向変位;　・・・　;　節点4の回転角 ]

変形から各部材の応力（軸応力・せん断応力・曲げ応力）を求める

最後に、変位ベクトル u から各部材の応力（軸応力・せん断応力・曲げ応力）を求めて終了となります。

ここで注意点としては、部材の応力（軸応力・せん断応力・曲げ応力）は部材断面に対する方向ごとに定義されているため、

全体座標系ではなく、要素座標系で算出します。

トライ

部材の応力について補足します。

軸応力：部材軸方向に伸びたり縮んだりする変形に抵抗する力
せん断応力：部材軸と垂直方向にずれる変形に抵抗する力
曲げ応力：部材の折れ曲がりに抵抗する力

応力ベクトルは要素座標系における剛性行列 k_bar に右から変位ベクトル u を掛ければいいのですが、先ほど求めた変位ベクトル u は全体座標系における変位ベクトル u です。

そのため、変位ベクトル u の左から座標変換行列 T を掛けて要素座標系における変位ベクトル u_bar に変換する必要があります。

ここで応力ベクトルを F とすると、応力ベクトル F は以下の式で求められます。

平面骨組要素の新規関数コード

ここからは平面骨組要素のFEMのMATLABコード（新規関数）を公開します。

トライ

今回は新規関数のみを公開し、次回メインプログラムを公開・解説します。

部材長さLを求める関数：PlaneFrameElementLength.m

i 端と j 端の節点の座標を入力し、平面骨組要素の長さを出力する関数（function）

「PlaneFrameElementLength」を定義します。

function y = PlaneFrameElementLength(x1,y1,x2,y2)
    y = sqrt( (x2-x1)*(x2-x1)+(y2-y1)*(y2-y1) );
end

トライ

関数名は異なりますが、コードの内容は以前公開した平面トラスの要素長さを求める関数とまったく同じです。

要素剛性行列kを求める関数：PlaneFrameElementStiffness.m

次に、各部材の要素剛性行列を求めます。

ヤング係数 E
断面積 A
断面二次モーメント I
部材の長さ L
全体座標系に対して角度θをもつ部材のcosθ, sinθ（関数内の表記はそれぞれC, S）

を入力し、平面骨組部材の要素剛性行列を出力する関数（function）「PlaneFrameElementStiffness」を定義します。

function y = PlaneFrameElementStiffness(E,A,I,L,C,S)

    w1 = A*C*C + 12*I*S*S/(L*L);
    w2 = A*S*S + 12*I*C*C/(L*L);
    w3 = (A-12*I/(L*L) )*C*S;
    w4 = 6*I*S/L;
    w5 = 6*I*C/L;
    y = E/L*[  w1   w3  -w4  -w1  -w3  -w4;
               w3   w2   w5  -w3  -w2   w5;
              -w4   w5  4*I   w4  -w5  2*I;
              -w1  -w3   w4   w1   w3   w4;
              -w3  -w2  -w5   w3   w2  -w5;
              -w4   w5  2*I   w4  -w5  4*I ];
end

トライ

自分でコードを書く際のポイントとして、

「y = …」以降のように、「;」の後に改行したり、行の中でもスペースを利用して各項の位置を合わせておくことで、後で間違いに気付きやすくなるのでぜひ真似してみて下さい。

全体剛性行列Kを求める関数：PlaneTrussAssemble.m

次に、全体座標系における各部材の要素剛性行列をアセンブル（結合）し、全体剛性行列を出力する関数（function）「PlaneFrameAssemble」を定義します。

function y = PlaneFrameAssemble(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

部材の応力を求める関数：PlaneFrameElementForces.m

全体剛性行列 K が求まったら外荷重ベクトル f を定義し、剛性方程式 f = Ku を解いて変位ベクトル u を求めます。

剛性方程式より求まった変位ベクトル u を要素座標系に変換し、各部材の応力を求めたらFEMはほぼ完了です。

ヤング係数 E
断面積 A
断面二次モーメント I
部材の長さ L
全体座標系に対して角度θをもつ部材のcosθ, sinθ（関数内の表記はそれぞれC, S）
座標変換行列 T
変位ベクトル u

を入力し、平面骨組部材の応力ベクトルを出力する関数（function）「PlaneFrameElementForces」を定義します。

function y = PlaneFrameElementForces(E,A,I,L,C,S,u)
    w1 = E*A/L;
    w2 = 12*E*I/(L*L*L);
    w3 = 6*E*I/(L*L);
    w4 = 4*E*I/L;
    w5 = 2*E*I/L;
    kbar = [ w1   0   0 -w1   0   0;
              0  w2  w3   0 -w2  w3;
              0  w3  w4   0 -w3  w5;
            -w1   0   0  w1   0   0;
              0 -w2 -w3   0  w2 -w3;
              0  w3  w5   0 -w3  w4 ];
    T = [ C  S  0  0  0  0;
         -S  C  0  0  0  0;
          0  0  1  0  0  0;
          0  0  0  C  S  0;
          0  0  0 -S  C  0;
          0  0  0  0  0  1 ];
    y = kbar * T * u;
end

まとめ：骨組要素（梁要素）はFEMで最も重要な要素

今回は平面骨組要素のFEMの基本的なフローを解説し、MATLABの新規関数コードを公開しました。

公開した新規関数コードは以下の手順で、今後作成するメインプログラムと同じフォルダに保存しておきましょう。

新規関数の保存方法

画面左上の「新規のスクリプト」をクリック、もしくは「Ctrl + N」
エディタ画面に本記事で公開したコードをコピー＆ペースト
「ファイルの保存」をクリック、もしくは「Ctrl + S」で新規関数を保存
mファイル名は関数名と同じにしておく
新規関数および今後作成するメインプログラムは同じフォルダに保存する

次回、平面骨組要素のメインプログラムを公開し、その中身について解説します。