【MATLABコード完全公開】立体骨組要素の有限要素法を解説

プログラミング　MATLAB　有限要素法 プログラミング
2021.08.28

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

本ブログではこれまでに

  • 平面トラス要素
  • 立体トラス要素
  • 平面骨組要素（梁要素）

の有限要素法（以下、FEM）の基本的な理論を解説し、MATLABの新規関数コードとメインプログラムのコードを公開してきました。

これらはいわゆる「線材要素」と呼ばれる要素です。

今回は線材要素の締めくくりとして、立体骨組要素のFEMの新規関数およびメインプログラムについて解説します。

前回の平面骨組要素は骨組要素（梁要素）の理論を解説するのに良い例でした。

ですが、実務で最も使うのは間違いなく今回の「立体骨組要素」です。

今回は、私が作成した新規関数およびメインプログラムのMATLABコードをすべて公開します。

このあと公開するコードをコピー＆ペーストすれば立体骨組要素のFEMプログラムを動かせるようになります。

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

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

以下のように、節点 i と節点 j から構成される立体骨組要素について考えていきます。

プログラミング　MATLAB 有限要素法
トライ
トライ

FEMのフローについてはこれまでに3回解説してきているため、今回は割愛します。

FEMのフローを理解していない人はこちらの記事をご覧下さい。

【MATLABコード公開】平面骨組要素の有限要素法①～新規関数編～
骨組要素（梁要素）は板要素、トラス要素に並び、FEMで最も重要な要素のひとつです。 構造設計者が行う建物の構造解析も、この骨組要素を中心に建物をモデル化します。 本記事ではFEMの第3弾として、平面骨組要素、いわゆる梁要素のFEMの基本的な理論と新規作成の関数について解説します。 今回も平面トラス、立体トラスの回と同様に、実際に私が作成したMATLABコードをすべて公開します。 ぜひお読みください。
1st-archiengineer.com
2021.08.13

要素座標系における立体骨組要素の要素剛性行列 kbar は以下の式で表されます。

プログラミング　MATLAB 有限要素法

式中の w1~w10 はそれぞれ以下のようになります。

プログラミング　MATLAB 有限要素法
トライ
トライ

立体骨組要素の場合、1節点あたりの自由度は

  • x方向変位
  • y方向変位
  • z方向変位
  • x軸まわりの回転
  • y軸まわりの回転
  • z軸まわりの回転

の計6つです。

そして、部材は i 端と j 端で構成されるので剛性行列のサイズは12×12となります。

ここで、式中に登場する各変数についてまとめておきます。

  • E：ヤング係数
  • A：断面積
  • L：部材長さ
  • Iz：要素座標系の z軸（断面の縦軸）まわりの断面二次モーメント
  • Iy：要素座標系の y軸（断面の横軸）まわりの断面二次モーメント
  • G：せん断剛性係数
  • J：サンブナンのねじり定数

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

先ほどの要素剛性行列を全体座標系における要素剛性行列に変換します。

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

CXxは部材と全体座標系X軸との方向余弦、CYxは部材と全体座標系Y軸との方向余弦、CZxは部材と全体座標系Z軸との方向余弦を示しています。

この座標変換行列 R を用いて、全体座標系における要素剛性行列 k は

「k = RT kbar R」 となります。

ここで、MATLABコードもあわせて公開します。

部材長さ L を求める新規関数「SpaceFrameElementLength」と

全体座標系における要素剛性行列を求める新規関数「SpaceFrameElementStiffness」を定義します。

【 SpaceFrameElementLength.m】

function y = SpaceFrameElementLength(x1,y1,z1,x2,y2,z2,x3,y3,z3)
    y = sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1) );
end

【SpaceFrameElementStiffness.m】

function y = SpaceFrameElementStiffness(E,G,A,Iy,Iz,J,x1,y1,z1,x2,y2,z2)
    L = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1));
    w1 = E*A/L;
    w2 = 12*E*Iz/(L*L*L);
    w3 = 6*E*Iz/(L*L);
    w4 = 4*E*Iz/L;
    w5 = 2*E*Iz/L;
    w6 = 12*E*Iy/(L*L*L);
    w7 = 6*E*Iy/(L*L);
    w8 = 4*E*Iy/L;
    w9 = 2*E*Iy/L;
    w10 = G*J/L;
    kbar = [ w1   0   0    0   0   0 -w1   0   0    0   0   0;
              0  w2   0    0   0  w3   0 -w2   0    0   0  w3;
              0   0  w6    0 -w7   0   0   0 -w6    0 -w7   0;
              0   0   0  w10   0   0   0   0   0 -w10   0   0;
              0   0 -w7    0  w8   0   0   0  w7    0  w9   0;
              0  w3   0    0   0  w4   0 -w3   0    0   0  w5;
            -w1   0   0    0   0   0  w1   0   0    0   0   0;
              0 -w2   0    0   0 -w3   0  w2   0    0   0 -w3;
              0   0 -w6    0  w7   0   0   0  w6    0  w7   0;
              0   0   0 -w10   0   0   0   0   0  w10   0   0;
              0   0 -w7    0  w9   0   0   0  w7    0  w8   0;
              0  w3   0    0   0  w5   0 -w3   0    0   0  w4 ];
    if x1 == x2 && y1 == y2
        if z2 > z1
            Lambda = [0 0 1 ; 0 1 0 ; -1 0 0];
        else
            Lambda = [0 0 -1 ; 0 1 0 ; 1 0 0];
        end
    else
        CXx = (x2-x1)/L;
        CYx = (y2-y1)/L;
        CZx = (z2-z1)/L;
        D   = sqrt(CXx*CXx + CYx*CYx);
        CXy = -CYx/D;
        CYy = CXx/D;
        CZy = 0;
        CXz = -CXx*CZx/D;
        CYz = -CYx*CZx/D;
        CZz = D;
        Lambda = [CXx CYx CZx ; CXy CYy CZy ; CXz CYz CZz];
   end
   R = [  Lambda  zeros(3) zeros(3) zeros(3) ;
         zeros(3)  Lambda  zeros(3) zeros(3) ;
         zeros(3) zeros(3)  Lambda  zeros(3) ;
         zeros(3) zeros(3) zeros(3)  Lambda   ];
   y = R' * kbar * R;
end

要素剛性行列を結合し、全体剛性行列を求める

次に、全体座標系における各部材の要素剛性行列をアセンブル（結合）し、架構全体の剛性行列 K を算出します。

新規関数「SpaceFrameAssemble」を以下のように定義します。

【 SpaceFrameAssemble.m】

function y = SpaceFrameAssemble(K,k,i,j)
    K(6*i-5,6*i-5) = K(6*i-5,6*i-5) + k(1,1);
    K(6*i-5,6*i-4) = K(6*i-5,6*i-4) + k(1,2);
    K(6*i-5,6*i-3) = K(6*i-5,6*i-3) + k(1,3);
    K(6*i-5,6*i-2) = K(6*i-5,6*i-2) + k(1,4);
    K(6*i-5,6*i-1) = K(6*i-5,6*i-1) + k(1,5);
    K(6*i-5,6*i)   = K(6*i-5,6*i)   + k(1,6);
    K(6*i-5,6*j-5) = K(6*i-5,6*j-5) + k(1,7);
    K(6*i-5,6*j-4) = K(6*i-5,6*j-4) + k(1,8);
    K(6*i-5,6*j-3) = K(6*i-5,6*j-3) + k(1,9);
    K(6*i-5,6*j-2) = K(6*i-5,6*j-2) + k(1,10);
    K(6*i-5,6*j-1) = K(6*i-5,6*j-1) + k(1,11);
    K(6*i-5,6*j)   = K(6*i-5,6*j)   + k(1,12);
    K(6*i-4,6*i-5) = K(6*i-4,6*i-5) + k(2,1);
    K(6*i-4,6*i-4) = K(6*i-4,6*i-4) + k(2,2);
    K(6*i-4,6*i-3) = K(6*i-4,6*i-3) + k(2,3);
    K(6*i-4,6*i-2) = K(6*i-4,6*i-2) + k(2,4);
    K(6*i-4,6*i-1) = K(6*i-4,6*i-1) + k(2,5);
    K(6*i-4,6*i)   = K(6*i-4,6*i)   + k(2,6);
    K(6*i-4,6*j-5) = K(6*i-4,6*j-5) + k(2,7);
    K(6*i-4,6*j-4) = K(6*i-4,6*j-4) + k(2,8);
    K(6*i-4,6*j-3) = K(6*i-4,6*j-3) + k(2,9);
    K(6*i-4,6*j-2) = K(6*i-4,6*j-2) + k(2,10);
    K(6*i-4,6*j-1) = K(6*i-4,6*j-1) + k(2,11);
    K(6*i-4,6*j)   = K(6*i-4,6*j)   + k(2,12);
    K(6*i-3,6*i-5) = K(6*i-3,6*i-5) + k(3,1);
    K(6*i-3,6*i-4) = K(6*i-3,6*i-4) + k(3,2);
    K(6*i-3,6*i-3) = K(6*i-3,6*i-3) + k(3,3);
    K(6*i-3,6*i-2) = K(6*i-3,6*i-2) + k(3,4);
    K(6*i-3,6*i-1) = K(6*i-3,6*i-1) + k(3,5);
    K(6*i-3,6*i)   = K(6*i-3,6*i)   + k(3,6);
    K(6*i-3,6*j-5) = K(6*i-3,6*j-5) + k(3,7);
    K(6*i-3,6*j-4) = K(6*i-3,6*j-4) + k(3,8);
    K(6*i-3,6*j-3) = K(6*i-3,6*j-3) + k(3,9);
    K(6*i-3,6*j-2) = K(6*i-3,6*j-2) + k(3,10);
    K(6*i-3,6*j-1) = K(6*i-3,6*j-1) + k(3,11);
    K(6*i-3,6*j)   = K(6*i-3,6*j)   + k(3,12);
    K(6*i-2,6*i-5) = K(6*i-2,6*i-5) + k(4,1);
    K(6*i-2,6*i-4) = K(6*i-2,6*i-4) + k(4,2);
    K(6*i-2,6*i-3) = K(6*i-2,6*i-3) + k(4,3);
    K(6*i-2,6*i-2) = K(6*i-2,6*i-2) + k(4,4);
    K(6*i-2,6*i-1) = K(6*i-2,6*i-1) + k(4,5);
    K(6*i-2,6*i)   = K(6*i-2,6*i)   + k(4,6);
    K(6*i-2,6*j-5) = K(6*i-2,6*j-5) + k(4,7);
    K(6*i-2,6*j-4) = K(6*i-2,6*j-4) + k(4,8);
    K(6*i-2,6*j-3) = K(6*i-2,6*j-3) + k(4,9);
    K(6*i-2,6*j-2) = K(6*i-2,6*j-2) + k(4,10);
    K(6*i-2,6*j-1) = K(6*i-2,6*j-1) + k(4,11);
    K(6*i-2,6*j)   = K(6*i-2,6*j)   + k(4,12);
    K(6*i-1,6*i-5) = K(6*i-1,6*i-5) + k(5,1);
    K(6*i-1,6*i-4) = K(6*i-1,6*i-4) + k(5,2);
    K(6*i-1,6*i-3) = K(6*i-1,6*i-3) + k(5,3);
    K(6*i-1,6*i-2) = K(6*i-1,6*i-2) + k(5,4);
    K(6*i-1,6*i-1) = K(6*i-1,6*i-1) + k(5,5);
    K(6*i-1,6*i)   = K(6*i-1,6*i)   + k(5,6);
    K(6*i-1,6*j-5) = K(6*i-1,6*j-5) + k(5,7);
    K(6*i-1,6*j-4) = K(6*i-1,6*j-4) + k(5,8);
    K(6*i-1,6*j-3) = K(6*i-1,6*j-3) + k(5,9);
    K(6*i-1,6*j-2) = K(6*i-1,6*j-2) + k(5,10);
    K(6*i-1,6*j-1) = K(6*i-1,6*j-1) + k(5,11);
    K(6*i-1,6*j)   = K(6*i-1,6*j)   + k(5,12);
    K(6*i,6*i-5)   = K(6*i,6*i-5)   + k(6,1);
    K(6*i,6*i-4)   = K(6*i,6*i-4)   + k(6,2);
    K(6*i,6*i-3)   = K(6*i,6*i-3)   + k(6,3);
    K(6*i,6*i-2)   = K(6*i,6*i-2)   + k(6,4);
    K(6*i,6*i-1)   = K(6*i,6*i-1)   + k(6,5);
    K(6*i,6*i)     = K(6*i,6*i)     + k(6,6);
    K(6*i,6*j-5)   = K(6*i,6*j-5)   + k(6,7);
    K(6*i,6*j-4)   = K(6*i,6*j-4)   + k(6,8);
    K(6*i,6*j-3)   = K(6*i,6*j-3)   + k(6,9);
    K(6*i,6*j-2)   = K(6*i,6*j-2)   + k(6,10);
    K(6*i,6*j-1)   = K(6*i,6*j-1)   + k(6,11);
    K(6*i,6*j)     = K(6*i,6*j)     + k(6,12);
    K(6*j-5,6*i-5) = K(6*j-5,6*i-5) + k(7,1);
    K(6*j-5,6*i-4) = K(6*j-5,6*i-4) + k(7,2);
    K(6*j-5,6*i-3) = K(6*j-5,6*i-3) + k(7,3);
    K(6*j-5,6*i-2) = K(6*j-5,6*i-2) + k(7,4);
    K(6*j-5,6*i-1) = K(6*j-5,6*i-1) + k(7,5);
    K(6*j-5,6*i)   = K(6*j-5,6*i)   + k(7,6);
    K(6*j-5,6*j-5) = K(6*j-5,6*j-5) + k(7,7);
    K(6*j-5,6*j-4) = K(6*j-5,6*j-4) + k(7,8);
    K(6*j-5,6*j-3) = K(6*j-5,6*j-3) + k(7,9);
    K(6*j-5,6*j-2) = K(6*j-5,6*j-2) + k(7,10);
    K(6*j-5,6*j-1) = K(6*j-5,6*j-1) + k(7,11);
    K(6*j-5,6*j)   = K(6*j-5,6*j)   + k(7,12);
    K(6*j-4,6*i-5) = K(6*j-4,6*i-5) + k(8,1);
    K(6*j-4,6*i-4) = K(6*j-4,6*i-4) + k(8,2);
    K(6*j-4,6*i-3) = K(6*j-4,6*i-3) + k(8,3);
    K(6*j-4,6*i-2) = K(6*j-4,6*i-2) + k(8,4);
    K(6*j-4,6*i-1) = K(6*j-4,6*i-1) + k(8,5);
    K(6*j-4,6*i)   = K(6*j-4,6*i)   + k(8,6);
    K(6*j-4,6*j-5) = K(6*j-4,6*j-5) + k(8,7);
    K(6*j-4,6*j-4) = K(6*j-4,6*j-4) + k(8,8);
    K(6*j-4,6*j-3) = K(6*j-4,6*j-3) + k(8,9);
    K(6*j-4,6*j-2) = K(6*j-4,6*j-2) + k(8,10);
    K(6*j-4,6*j-1) = K(6*j-4,6*j-1) + k(8,11);
    K(6*j-4,6*j)   = K(6*j-4,6*j)   + k(8,12);
    K(6*j-3,6*i-5) = K(6*j-3,6*i-5) + k(9,1);
    K(6*j-3,6*i-4) = K(6*j-3,6*i-4) + k(9,2);
    K(6*j-3,6*i-3) = K(6*j-3,6*i-3) + k(9,3);
    K(6*j-3,6*i-2) = K(6*j-3,6*i-2) + k(9,4);
    K(6*j-3,6*i-1) = K(6*j-3,6*i-1) + k(9,5);
    K(6*j-3,6*i)   = K(6*j-3,6*i)   + k(9,6);
    K(6*j-3,6*j-5) = K(6*j-3,6*j-5) + k(9,7);
    K(6*j-3,6*j-4) = K(6*j-3,6*j-4) + k(9,8);
    K(6*j-3,6*j-3) = K(6*j-3,6*j-3) + k(9,9);
    K(6*j-3,6*j-2) = K(6*j-3,6*j-2) + k(9,10);
    K(6*j-3,6*j-1) = K(6*j-3,6*j-1) + k(9,11);
    K(6*j-3,6*j)   = K(6*j-3,6*j)   + k(9,12);
    K(6*j-2,6*i-5) = K(6*j-2,6*i-5) + k(10,1);
    K(6*j-2,6*i-4) = K(6*j-2,6*i-4) + k(10,2);
    K(6*j-2,6*i-3) = K(6*j-2,6*i-3) + k(10,3);
    K(6*j-2,6*i-2) = K(6*j-2,6*i-2) + k(10,4);
    K(6*j-2,6*i-1) = K(6*j-2,6*i-1) + k(10,5);
    K(6*j-2,6*i)   = K(6*j-2,6*i)   + k(10,6);
    K(6*j-2,6*j-5) = K(6*j-2,6*j-5) + k(10,7);
    K(6*j-2,6*j-4) = K(6*j-2,6*j-4) + k(10,8);
    K(6*j-2,6*j-3) = K(6*j-2,6*j-3) + k(10,9);
    K(6*j-2,6*j-2) = K(6*j-2,6*j-2) + k(10,10);
    K(6*j-2,6*j-1) = K(6*j-2,6*j-1) + k(10,11);
    K(6*j-2,6*j)   = K(6*j-2,6*j)   + k(10,12);
    K(6*j-1,6*i-5) = K(6*j-1,6*i-5) + k(11,1);
    K(6*j-1,6*i-4) = K(6*j-1,6*i-4) + k(11,2);
    K(6*j-1,6*i-3) = K(6*j-1,6*i-3) + k(11,3);
    K(6*j-1,6*i-2) = K(6*j-1,6*i-2) + k(11,4);
    K(6*j-1,6*i-1) = K(6*j-1,6*i-1) + k(11,5);
    K(6*j-1,6*i)   = K(6*j-1,6*i)   + k(11,6);
    K(6*j-1,6*j-5) = K(6*j-1,6*j-5) + k(11,7);
    K(6*j-1,6*j-4) = K(6*j-1,6*j-4) + k(11,8);
    K(6*j-1,6*j-3) = K(6*j-1,6*j-3) + k(11,9);
    K(6*j-1,6*j-2) = K(6*j-1,6*j-2) + k(11,10);
    K(6*j-1,6*j-1) = K(6*j-1,6*j-1) + k(11,11);
    K(6*j-1,6*j)   = K(6*j-1,6*j)   + k(11,12);
    K(6*j,6*i-5)   = K(6*j,6*i-5)   + k(12,1);
    K(6*j,6*i-4)   = K(6*j,6*i-4)   + k(12,2);
    K(6*j,6*i-3)   = K(6*j,6*i-3)   + k(12,3);
    K(6*j,6*i-2)   = K(6*j,6*i-2)   + k(12,4);
    K(6*j,6*i-1)   = K(6*j,6*i-1)   + k(12,5);
    K(6*j,6*i)     = K(6*j,6*i)     + k(12,6);
    K(6*j,6*j-5)   = K(6*j,6*j-5)   + k(12,7);
    K(6*j,6*j-4)   = K(6*j,6*j-4)   + k(12,8);
    K(6*j,6*j-3)   = K(6*j,6*j-3)   + k(12,9);
    K(6*j,6*j-2)   = K(6*j,6*j-2)   + k(12,10);
    K(6*j,6*j-1)   = K(6*j,6*j-1)   + k(12,11);
    K(6*j,6*j)     = K(6*j,6*j)     + k(12,12);
    y = K;
end

剛性方程式を解き、変位を算出する

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

荷重ベクトル f と変位ベクトル U は具体的に以下のような構成になっています。

プログラミング　MATLAB 有限要素法

式中の各変数についてまとめておきます。

  • X：X軸方向変位
  • Y：Y軸方向変位
  • Z：Z軸方向変位
  • θX：X軸まわりの回転角
  • θY：Y軸まわりの回転角
  • θZ：Z軸まわりの回転角

  • N：軸力
  • QY：Y軸方向のせん断力
  • QZ：Z軸方向のせん断力
  • MX：X軸まわりのモーメント（ねじりモーメント）
  • MY：Y軸まわりのモーメント
  • MZ：Z軸まわりのモーメント

剛性方程式 f = KU の両辺に左から K の逆行列 K-1 をかければ変位ベクトルは

U = K-1 f となります。

トライ
トライ

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

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

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

変形から部材の応力を求める

変位ベクトル U が求まったら応力ベクトル F を求めます。

F = kbarRU

応力ベクトル F の算出式とMATLABコードは以下のようになります。

【 SpaceFrameElementForces.m】

function y = SpaceFrameElementForces(E,G,A,Iy,Iz,J,x1,y1,z1,x2,y2,z2,u)
    L = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1));
    w1 = E*A/L;
    w2 = 12*E*Iz/(L*L*L);
    w3 = 6*E*Iz/(L*L);
    w4 = 4*E*Iz/L;
    w5 = 2*E*Iz/L;
    w6 = 12*E*Iy/(L*L*L);
    w7 = 6*E*Iy/(L*L);
    w8 = 4*E*Iy/L;
    w9 = 2*E*Iy/L;
    w10 = G*J/L;
    kbar = [ w1   0   0    0   0   0 -w1   0   0    0   0   0;
              0  w2   0    0   0  w3   0 -w2   0    0   0  w3;
              0   0  w6    0 -w7   0   0   0 -w6    0 -w7   0;
              0   0   0  w10   0   0   0   0   0 -w10   0   0;
              0   0 -w7    0  w8   0   0   0  w7    0  w9   0;
              0  w3   0    0   0  w4   0 -w3   0    0   0  w5;
            -w1   0   0    0   0   0  w1   0    0   0   0   0;
              0 -w2   0    0   0 -w3   0  w2   0    0   0 -w3;
              0   0 -w6    0  w7   0   0   0  w6    0  w7   0;
              0   0   0 -w10   0   0   0   0   0  w10   0   0;
              0   0 -w7    0  w9   0   0   0  w7    0  w8   0;
              0  w3   0    0   0  w5   0 -w3   0    0   0  w4 ];
    if x1 == x2 & y1 == y2
        if z2 > z1
            Lambda = [0 0 1 ; 0 1 0 ; -1 0 0];
        else
            Lambda = [0 0 -1 ; 0 1 0 ; 1 0 0];
        end
    else
        CXx = (x2-x1)/L;
        CYx = (y2-y1)/L;
        CZx = (z2-z1)/L;
        D = sqrt(CXx*CXx + CYx*CYx);
        CXy = -CYx/D;
        CYy = CXx/D;
        CZy = 0;
        CXz = -CXx*CZx/D;
        CYz = -CYx*CZx/D;
        CZz = D;
        Lambda = [CXx CYx CZx ; CXy CYy CZy ; CXz CYz CZz];
    end
    R = [ Lambda  zeros(3) zeros(3) zeros(3);
         zeros(3)  Lambda  zeros(3) zeros(3);
         zeros(3) zeros(3)  Lambda  zeros(3);
         zeros(3) zeros(3) zeros(3)  Lambda  ];
    y = kbar*R*u;
end
トライ
トライ

新規関数については以上になります。

メインプログラムの公開

ここからはメインプログラムのMATLABコードを公開します。

今回は以下のようなモデルについて考えます。

プログラミング　MATLAB 有限要素法

トライ
トライ

断面積や断面二次モーメントなどは以下のように設定しました。

  • 部材のヤング係数 E：210×106 [kN/m2]
  • せん断剛性係数 G：84×106 [kN/m2]
  • 部材の断面積 A：すべて0.02 [m2]
  • 部材の強軸方向の断面二次モーメント Iy：すべて0.00002 [m2]
  • 部材の弱軸方向の断面二次モーメント Iz：すべて0.00001 [m2]
  • サンブナンのねじり定数 J：すべて0.000005 [m2]
  • 支持条件：節点１~節点４が固定支持
  • 外力　節点７にーX方向に15kN

【メインプログラム】

clear

dim = 3;  % 次元数(dimension:3次元)
dof = 6;  % 節点の自由度数(degrees of freedom:立体骨組よりDX,DY,DZ,θx,θy,θzの6種類)

E  = 210e6;    % ヤング係数     単位[kN/m^2]
G  =  84e6;    % せん断弾性係数 単位[kN/m^2]
A  = [  2e-2;  2e-2;  2e-2;  2e-2;  2e-2;  2e-2;  2e-2;  2e-2 ]; % 断面積 単位[m^2]
Iy = [ 20e-5; 20e-5; 20e-5; 20e-5; 20e-5; 20e-5; 20e-5; 20e-5 ]; % 断面二次モーメント   単位[m^2]
Iz = [ 10e-5; 10e-5; 10e-5; 10e-5; 10e-5; 10e-5; 10e-5; 10e-5 ]; % 断面二次モーメント   単位[m^2]
J  = [  5e-5;  5e-5;  5e-5;  5e-5;  5e-5;  5e-5;  5e-5;  5e-5 ]; % サンブナンのねじり定数 単位[m^2]
XYZ = [ 0  0  0;    0  4  0;    4  4  0;    4  0  0;
        0  0  5;    0  4  5;    4  4  5;    4  0  5  ]; % 節点座標[x,y,z]　単位[m]
CON = [ 1  5; 2  6; 3  7;  4  8;  5  6; 6  7; 7  8; 5  8 ]; % 接続条件(connection) [節点i, 節点j]  
% 境界条件(boundary condition) [節点番号, DX, DY, DZ, θx, θy, θz (0:自由, 1:固定)]  
BOUN = [ 1  1  1  1  1  1  1;
         2  1  1  1  1  1  1;
         3  1  1  1  1  1  1;
         4  1  1  1  1  1  1  ];
% 外力[節点番号, PX, PY, PZ, Mx, My, Mz] 単位[kN,kNm]
LOAD = [ 7  -15 0  0  0  0  0 ];

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) = SpaceFrameElementLength(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} = SpaceFrameElementStiffness( E,G,A(i),Iy(i),Iz(i),J(i), ...
                                        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) );
    K = SpaceFrameAssemble(K,ke{i},CON(i,1),CON(i,2));
end

% 外力ベクトルfの算出
f = zeros(nd,1);
for i = 1:nl
    f( dof*(LOAD(i,1))-5, 1 ) = LOAD(i,2);
    f( dof*(LOAD(i,1))-4, 1 ) = LOAD(i,3);
    f( dof*(LOAD(i,1))-3, 1 ) = LOAD(i,4);
    f( dof*(LOAD(i,1))-2, 1 ) = LOAD(i,5);
    f( dof*(LOAD(i,1))-1, 1 ) = LOAD(i,6);
    f( dof*(LOAD(i,1))  , 1 ) = LOAD(i,7);
end

% 変位を許容した自由度に対して剛性方程式を解き、変位ベクトルUと節点力ベクトルFを求める
U = zeros(nd,1);
U(free,1) = K(free,free)\f(free);
F = K * U;

% 各部材の変位ベクトルuの算出
% u = [DXi DYi DZi θix θiy θiz DXj DYj DZj θjx θjy θjz]'
for i = 1:nm
    u{i} = [ U(dof*CON(i,1)-5);
             U(dof*CON(i,1)-4);
             U(dof*CON(i,1)-3);
             U(dof*CON(i,1)-2);
             U(dof*CON(i,1)-1);
             U(dof*CON(i,1)  );
             U(dof*CON(i,2)-5);
             U(dof*CON(i,2)-4);
             U(dof*CON(i,2)-3);
             U(dof*CON(i,2)-2);
             U(dof*CON(i,2)-1);
             U(dof*CON(i,2)  )  ];
end

% 各部材の応力ベクトルForceの算出
% Force = [ Ni QiY QiZ Mix Miy Miz Nj QjY QjZ Mjx Mjy Mjz ]'
for i = 1:nm
    Force{i} = SpaceFrameElementForces(E,G,A(i),Iy(i),Iz(i),J(i), ...
                                       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),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);
        Z0(i,j) = XYZ0((CON(i,j)),3);
    end
    plot3(X0(i,:),Y0(i,:),Z0(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);
        Z(i,j) = plotXYZ((CON(i,j)),3);
    end
    plot3(X(i,:),Y(i,:),Z(i,:),'r-','LineWidth',2);
    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',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) ), ...
          1/2*( XYZ0(CON(i,1),3) + XYZ0(CON(i,2),3) ), ...
         ['(',num2str(i),')'],'HorizontalAlignment','right',...
          'VerticalAlignment','top','FontSize',18,'color','black' );
    hold on    
end

axis equal; grid on
xlabel('x(m)','FontSize',18); ylabel('y(m)','FontSize',18); zlabel('z(m)','FontSize',18);
xlim([0 4.5]); ylim([0 4.5]); zlim([-0.5 5]);
title('Space Frame Analysis','FontSize',18)
set(gca,'FontSize',18);
hold off

メインプログラムを実行すると、以下のような変形図が作成されます。

節点７がちゃんとーX方向に変形しています。

プログラミング　MATLAB 有限要素法

まとめ

今回はFEMで最も使用される立体骨組要素について具体的なMATLABコードを紹介しました。

実務では汎用有限要素解析ソフトを用いるため、自分でコードを書くことはありません。

ですが、これまで合計9回に渡って解説してきた知識を知った上で汎用ソフトを扱うのか、そうでないかでエンジニアとしての成熟度は大きく変わります。

9回もあると結局どれを見ればいいのか分からない。

全部読む時間もないし…。

という方は、「平面骨組要素」についての２記事だけ読んでもらえればだいたい理解できると思います。

【MATLABコード公開】平面骨組要素の有限要素法①～新規関数編～
骨組要素（梁要素）は板要素、トラス要素に並び、FEMで最も重要な要素のひとつです。 構造設計者が行う建物の構造解析も、この骨組要素を中心に建物をモデル化します。 本記事ではFEMの第3弾として、平面骨組要素、いわゆる梁要素のFEMの基本的な理論と新規作成の関数について解説します。 今回も平面トラス、立体トラスの回と同様に、実際に私が作成したMATLABコードをすべて公開します。 ぜひお読みください。
1st-archiengineer.com
2021.08.13
【MATLABコード公開】平面骨組要素の有限要素法②~メインプログラム編~
今回は平面骨組要素（梁要素）のFEMのメインプログラムのコードを公開します。 1. FEMのフローと公開コードの対応確認、2. 手計算による公開コードの検証、3. FEMの結果の見方の解説を行います。 MATLAB、プログラミング、平面骨組要素（梁要素）のFEMに興味のある人はぜひお読み下さい。
1st-archiengineer.com
2021.08.17

最後までお読みいただきありがとうございました。

コメント

タイトルとURLをコピーしました