OpenSees Developer : TIPs : elasticBeamColumnでの剛性マトリクス

OpenSees のコマンド、elementでelasticBeamColumnを設定した場合（たとえば、下記の様なコマンド)

element elasticBeamColumn 1 1 3 3600 4227 1080000 1

での各要素での剛性マトリクスは、
ElasticBeam2d.cppのgetTangentStiff( )で 求めています。


坂田弘安　島崎和司著　学芸出版社　建築構造力学 II  P78

「軸力、せん断力、曲げを受ける部際が任意の方向を向いている場合の剛性マトリクス」

の、式(3.53) の[T] と、部材座標系[K] から

全体座標系[K} をもとめる式(3.54)

を用いて、[K] のための式を展開し、
ElasticBeam2d.cpp　（およびLinearCrdTransf2d.cpp) 
の中に、プログラムとして記載しています。

OpenSees のソースコード解析に挑戦してみる : 目次

OpenSees Developer : TIPs : Trussでの剛性マトリクス

OpenSees のコマンド、elementでTrussを設定した場合（たとえば、下記の様なコマンド）

element Truss 1 1 4 10.0 1

での各要素での剛性マトリクスは、
Truss.cppのgetTangentStiff( )で 求めています。


Truss::getTangentStiff(void)
{
   if (L == 0.0) { // - problem in setDomain() no further warnings
      theMatrix->Zero();
     return *theMatrix;
   }
    double E = theMaterial->getTangent();
   // come back later and redo this if too slow
   Matrix &stiff = *theMatrix;

   int numDOF2 = numDOF/2;
   double temp;
   double EAoverL = E*A/L;
   for (int i = 0; i < dimension; i++) {
     for (int j = 0; j < dimension; j++) {
              temp = cosX[i]*cosX[j]*EAoverL;
             stiff(i,j) = temp;
             stiff(i+numDOF2,j) = -temp;
             stiff(i,j+numDOF2) = -temp;
             stiff(i+numDOF2,j+numDOF2) = temp;
     }
   }
   return stiff;
}

たとえば、２次元２自由度の場合は
numDOF2 は 2
dimensionは２
が入っています

cosX[0]がcos
cosX[1]がsin
です

L: 長さ
A:断面積
E:ヤング率

stiff(i,j) = temp;
は(0,0)(0,1)(1,0)(1,1)
stiff(i+numDOF2,j) = -temp;
は(2,0)(3,0)(2,1)(3,1)
stiff(i,j+numDOF2) = -temp;
は(0,2)(0,3)(1,2)(1,3)
stiff(i+numDOF2,j+numDOF2) = temp;
は(2,2)(2,3)(3,3)(3,2)
の、剛性マトリクス成分に値が入ります。

そのため、 下記(3.33)の　[K]　が作成されることになります
(坂田弘安　島崎和司著　学芸出版社　建築構造力学 II  P67参照）

OpenSees のソースコード解析に挑戦してみる : 目次

つぶやき （と 謝辞）

OpenSeesファンの皆様

このブログを見ていただきましてありがとうございます。現在6415ビュー。Googleさんの統計によりますと、日本語て書いてあることもあり、「日本」からの閲覧数が　5024　と　圧倒的に多いです。

また、どの部分に閲覧者の方が興味をもたれているか、とうところでは、分野の目次（ページで書いてます）のビュー数から考えると、455 の　「OpenSeesコマンド」　が一番多いようです。
しかし「C言語での 「有限要素法における平面トラスの解析」」の閲覧数もこのごろ急激にあがってますので、この部分にも興味をもたれている方が多いのかもしれません。




「C言語での 「有限要素法における平面トラスの解析」」は、現在まだ、更新中です。
「OpenSeesコマンド」　の更新は、今お休み中です。

というのは、
「OpenSees のソースコード解析に挑戦してみる」　の解析を進めているためです。
目次にも書いてあるとおり、「苦戦中」です。　メモはだいぶ書きたまってきましたが、Blogにアップできるレベルのものまでになっていません。

「OpenSeesコマンド」　からわかるとおり、構造解析が簡単なtcl スクリプトで、試せてしまう、ということを実現するために、バックエンドのプログラム（OpenSeesのソースコード）は、こんなに大変なことをしているんだ。。。と半ば感動をもって読んでいます。
さすが、カリフォルニア大学バークレー校(UCB)、すごいですね。

<!----- さて、
このBlogを書いている環境ですが、週１回アルバイトで行っている某大学　建築学科　Y研究室で週１、２回ペースで解析を進めています。　本来の仕事の範疇からは、思いっきり外れているのですが、「大学って好きなことをする場所だから」とおっしゃってくれる、心の広いY教授に本当に感謝です。

また
「OpenSees のソースコード解析に挑戦してみる」　には、同じ建築学科の S研究室 S教授とY助教の多大なるコメントと、本に助けられています。　まったくまとまっていない内容についての質問にもかかわらず、（怪我したため、治療中の松葉つえをついて、わざわざ出てきてまでして）アドバイスをしてくださるS教授と、解析のためのアプローチ方法などに、思いもかけない切り口を考えてくれるY助教の助けがなかったら、ここまで読み込めなかったと思います。　このお二人にも大感謝として言いようがないほど。それに、S研のＯさん、行くとお菓子で歓迎してくれて、ありがとう！！

あと、
そもそも「OpenSees」というネタ（というか、言葉）を教えてくれたT教授にも、ありがとうを。

---!> こちらのブログは、あくまで趣味として、まったりと追記しておりますが、
今後も、同じ調子で進めていきますので、よろしくお願いします。

Innovativeな建物がたくさんできますように！

2015/4/24  桜が全部ちって緑がまぶしくなってきた丘で

<有限要素法>マトリクスなど、必要となるデータのメモリーのallocation

<有限要素法> C言語でのフローチャート と　使用するヘッダーファイルなどの定義で、定義しました行列に対して、メモリのアロケーションを行います。

point_num, element_numには、csvファイルから読み込んだNode（節点）の数とelement（要素）の数が入っています

メモリのアロケーションが終わったら、0を代入しておきます。

次のページ　→
<有限要素法>要素剛性マトリクス算出

「有限要素法における平面トラスの解析」プログラムをC言語で書いてみる : 目次

OpenSees Developer : TIPs : 逆行列

(TIPsは、断片的に見つけたことの覚書です)


逆行列[A]-1を行っているクラスはここにありました。

matrix.cpp

matrix class のMatrix::Invert( ) 


この中では、やはり高速数理計算用のライブラリ　LAPACKを使用していました。

//m行n列の行列[A]を行交換による部分枢軸選択を利用してLU分解
DGETRF(&n,&n,Aptr,&ldA,iPIV,&info);

//LU分解後の行列を利用して、逆行列[A]-1を計算
DGETRI(&n,Aptr,&ldA,iPIV,Wptr,&workSize,&info);

OpenSees のソースコード解析に挑戦してみる : 目次

<有限要素法>変位の計算 （まとめ）

<有限要素法>変位の計算　のプログラムをまとめたものは、以下の様になります。
（途中、いらなくなったメモリを解放しています

次のページ　→
<有限要素法>応力の計算とバランスチェック と　

「有限要素法における平面トラスの解析」プログラムをC言語で書いてみる : 目次

<有限要素法>変位の計算

変位の算出 は、以下の連立一次方程式を解くことによって行います

連立方程式を解くには、メモリを大量に使用するため、コンピュータの処理時間が長くなります。
そのため、固定点については解くことをしないようにするため
<有限要素法>固定条件の処理
で求めた、剛性マトリックス等を使用します。


剛性マトリックス　all_valid_ESmatrix_EK
変位：　valid_displacement:
荷重力（外力　節点力）valid_loading_force

ここでは、ガウス消去法(LEQファンクション）を使用して解いています。
後程LAPACKライブラリを使用して、さらに高速化した解法を示します。


//計算を高速化するため、固定されていない成分のみの変位を計算する
LEQ(all_valid_ESmatrix_EK, valid_loading_force, valid_displacement,valid_point_num);

valid_displacement[ ] に計算で得られた変位が入りますが、
これは、固定点を間引いた順番での並びになっていますので、これを間引く前の並びに変更します。constraint_sequence[　]に、固定点を間引いた並びと、間引く前の並びを対応させたものが入っています。


for (int i=0; i<valid_point_num; i++) {
displacement[constraint_sequence[i]] = valid_displacement[i];
}



ガウス消去法での連立方程式の解法ファンクション(LEQ)は以下の通りです。

次のページ　→
<有限要素法>変位の計算　（まとめ）

「有限要素法における平面トラスの解析」プログラムをC言語で書いてみる : 目次