【コード公開】平面トラスの有限要素解析③：新規関数編

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

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

今回からは平面トラスの有限要素解析を実際にMATLABを使ってコードを作成していきます。

私は当時、プログラミングを勉強していた時にいつもこのようなことを思っていました。

学生トライ
学生トライ

コード例が載っていればそれを読み解くことでMATLABの勉強も有限要素法の勉強も同時にできて効率が良いのになぁ…。

学生トライ
学生トライ

ネットで検索すると部分的にはコードが載っているけど、

一貫して載っていないから結局あまり参考にならないなぁ…。

こういった当時の私の悩みを解決するために、

「ここに載っているコードをそのままコピーすればプログラムが正常に動く」

というものを公開します。

※コードの書き方が最適かどうかは個人の意見によるので議論しないこととします。

量が多いため、今回は「新規で作成する関数」を解説し、次回でメインプログラムを解説します。

ぜひ併せて最後までお読み下さい。

手順①：要素剛性行列を求める

ここまで順番に記事を見てくれていた人には毎度おなじみのこのモデルについてのコードを作成していきます。

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

物性値と部材の長さを求める関数を定義する

まずは下記の要素剛性行列を求めるために各物性値を設定します。

MATLAB　有限要素法　トラス

このうち、ヤング係数 E と部材の断面積 A は具体的な数値をメインプログラム内で入力するので今回は割愛します。

部材の長さ L は三平方の定理で求めます。

MATLAB　有限要素法　トラス

この図の場合、部材の長さ L は以下のようになります。

L = sqrt { (Xj – Xi)2 + (Yj – Yi)2 }

「sqrt」はExcelでもおなじみですが「SQuare RooT」のことで、斜めの長さを求める関数です。

２つの節点の座標を入力し、平面トラスの部材の長さを出力する関数（function）

「PlaneTrussElementLength」を定義します。

MATLABコードは以下のようになります。

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

上記のように新たな関数を定義する場合は、書き出しを「function…」とします。

要素剛性行列を求める関数を定義する

要素剛性行列を求めるのに必要な定数

  • ヤング係数 E
  • 断面積 A
  • 部材の長さ L
  • 要素座標系と全体座標系との成す角度 θ

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

MATLABコードは以下のようになります。

function y = PlaneTrussElementStiffness(E,A,L,theta)
    x = theta*pi/180;
    C = cos(x);
    S = sin(x);
    y = E*A/L*[ C*C  C*S -C*C -C*S;  C*S  S*S -C*S -S*S;
               -C*C -C*S  C*C  C*S; -C*S -S*S  C*S  S*S ];
end

手順②：全体剛性行列を求める

次に要素剛性行列を結合（アセンブル）し、全体剛性行列を求めます。

各部材の節点 i と節点 j 、そして要素剛性行列 k を入力して全体剛性行列を出力する関数「PlaneTrussAssemble」を定義します。

MATLABコードは以下のようになります。

function y = PlaneTrussAssemble(K,k,i,j)
    K(2*i-1,2*i-1) = K(2*i-1,2*i-1) + k(1,1);
    K(2*i-1,2*i)   = K(2*i-1,2*i)   + k(1,2);
    K(2*i-1,2*j-1) = K(2*i-1,2*j-1) + k(1,3);
    K(2*i-1,2*j)   = K(2*i-1,2*j)   + k(1,4);
    K(2*i,2*i-1)   = K(2*i,2*i-1)   + k(2,1);
    K(2*i,2*i)     = K(2*i,2*i)     + k(2,2);
    K(2*i,2*j-1)   = K(2*i,2*j-1)   + k(2,3);
    K(2*i,2*j)     = K(2*i,2*j)     + k(2,4);
    K(2*j-1,2*i-1) = K(2*j-1,2*i-1) + k(3,1);
    K(2*j-1,2*i)   = K(2*j-1,2*i)   + k(3,2);
    K(2*j-1,2*j-1) = K(2*j-1,2*j-1) + k(3,3);
    K(2*j-1,2*j)   = K(2*j-1,2*j)   + k(3,4);
    K(2*j,2*i-1)   = K(2*j,2*i-1)   + k(4,1);
    K(2*j,2*i)     = K(2*j,2*i)     + k(4,2);
    K(2*j,2*j-1)   = K(2*j,2*j-1)   + k(4,3);
    K(2*j,2*j)     = K(2*j,2*j)     + k(4,4);
    y = K;
end

手順③：剛性方程式を解き、変形を求める

次に変形が拘束されていない自由度のみを抽出して剛性方程式を解き、変形を求めます。

拘束されていない自由度の抽出はメインプログラムの中で計算するので新規関数は今回は作成しません。

剛性方程式の解の算出については新規の関数ではなく、MATLABで既に用意されているコマンドである「 \ 」もしくは「inv( )」を使用します。

具体的には次回のメインプログラム編で解説します。

手順④：変形から各部材の応力（軸力）を求める

変形が求まったら最後に部材の応力（軸力）を求めて完了です。

  • ヤング係数 E
  • 断面積 A
  • 部材の長さ L
  • 要素座標系と全体座標系との成す角度 θ
  • 変形ベクトル u

を入力し、平面トラス部材の応力度（軸応力度）を出力する関数（function）「PlaneTrussElementStress」を定義します。

MATLABコードは以下のようになります。

function y = PlaneTrussElementStress(E,L,theta,u)
  x = theta*pi/180;
  C = cos(x);
  S = sin(x);
  y = E/L*[-C S C S]*u;
end

これにより得られた応力度（軸応力度）に断面積 A を掛けると応力（軸力）が求まります。

新規作成した関数を保存するには

これまでに作成した新規関数、そして次回解説するメインプログラムは、すべて同じフォルダに保存しておきましょう。

理由は、ディレクトリが異なるとプログラムが正常に動かなくなるためです。

ここでは新規作成した関数の保存方法について解説します。

手順としては以下の3ステップで完了です。

  1. 画面左上の「新規スクリプト」をクリック
  2. エディタ画面にコードをコピー&ペースト
  3. 「Ctrl+S」を押し、予め決めておいたフォルダ内に保存

具体的にOctaveの画面でも確認しておきましょう。

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

MTLABおよびOctaveでは「.m」という拡張子でデータが保存されます。（通称、mファイル）

ファイル名については「mファイルの名前＝関数の名前」としておくと後々便利です。

次回、今回作成した新規関数を使ってメインプログラムの例を掲載します。

まとめ：新規関数を作成して分かりやすいコードにしよう

今回は平面トラスの有限要素解析にあると便利な新規関数を実際のコードを用いて解説しました。

実際のところ、新規関数を定義せずにメインプログラム内ですべて書いても解析自体は可能です。

ですが、そうするとどうしてもメインプログラムのコードが長くなり、全体の流れが分かりにくくなってしまいます。

今回解説したコードを新規関数として保存しておけば、メインプログラムが短くなり、パッと見で全体の流れが理解しやすくなります。

次回はいよいよメインプログラムの解説となります。

実際のコード＋コードの中で使用しているコマンドについてもひとつずつ解説していきます。

かなりのボリュームになる予感がしますが、次回もぜひお読み下さい。

コメント

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