2020/05/24

クラウドベースの3次元構造解析ツール SkyCiv

Openseesの学習が少し停滞しておりますので、更新滞っております。

Openseesとはあまり関係ありませんが、使えそうなクラウドベースの構造解析ツールをたまたま見つけたので紹介します。その名も、「Sky Civ」です。以下、ページにて詳細を確認できます。(私は、本ソフトの製作者の回し者ではありませんので、悪しからず。。)

https://skyciv.com/

通常構造解析ツールは自分のパソコンにインストールしてそのパソコンでライセンスを取得してしようする形式が基本ですが、これですと、基本的に解析結果を確認できるのでは会社のパソコンのみとなります。
一方、設計者が外出先でちょっと条件を変えた解析結果をすぐ見たいといったニーズはあるかと思います。今回紹介するSky Civはクラウドベースの3次元構造解析ツールで、いつでも、どこでも、ただパソコンさえあれば、構造解析を行うことが可能です。構造設計者にとっては、いつでも仕事に打ち込める環境が整ってしまうということで、ますますワーカホリックが加速されるツールといえます。
本ツールは線形解析、非線形解析、固有値解析に対応しています。
クラウドベースのツールということで、professionalバージョンで109ドル(現在)のサブスクリプションモデルを採用しています。なかなかいいお値段しますが、通常の汎用構造解析ツールと比べると安いのか。。今回は、Freeversion(3次元応力解析の場合、解析モデルのノード10節点以下の制約あり)を少し触ってみました。



解析モデルはウェブブラウザ上で動作するクラウドツールの操作画面にて、節点をクラウド上の表シートにコピペで入力し、表示された節点同士をマウスでつなぐことで簡単に作成できます。その操作性は直感的でかなり良好です。



部材断面も様々なタイプの形状が既に組み込まれており、変更したい寸法部分をクリックして、寸法を入力することで直感的に入力できます。サクサク動きます。



応力解析結果も簡単に確認できました。
計算書も簡易なものですが、ボタン一つでpdfのファイルとして自動出力されます。


また、解析ファイルはクラウド上に保存されますので、どんなPC環境でも常に呼び出すことが可能です。リモートワークにも最適です。今回対象としたモデルはシンプルなフレームモデルだったからか、一般的な構造解析ツールと比較しても、全体的に動きが軽快で操作が快適な印象を受けました。

ソフトウェアを自分のパソコンにインストールしなくてもクラウド上でどこでも解析が行え、結果をどこでも自由に確認できることから、チームで設計を行う上でのコミュニケーションツールとして役立ちそうです。MicrosoftTeamsとの解析結果の連携し、解析結果やコメントを設計チームでスムーズに共有することも可能となっているようです。また、Pythonなどの自作プログラムと組み合わせた高度な活用を想定し、JSON形式によるAPIも装備しています。

さらに、比較的解析時間の大きい複雑なモデルを想定した場合、こうしたクラウド上のツールがスパコン等と連携することで手早く解析できるなどの可能性が考えられ、こうしたクラウドベースのツールの今後の進化にも注視する必要がありそうです。

2020/05/06

トラス要素から梁要素への拡張(Openseesによる3次元応力解析)

前回構築したOpenseesによる3次元フレーム要素の弾性応力解析プログラムについて、補足的に説明します。プログラムのソースはあまり3次元トラス要素のプログラム構成と変わりません。そこで、前回までに構築した3次元トラス要素プログラムの一部を修正したポイントを以下にかいていきます。扱う部材がトラス要素→梁要素となったことで以下のように部材の曲げモーメント、せん断力を考慮する必要がでてきます。修正したポイントはその点に関連するところがほとんどです。



・節点固定度、荷重項を3自由度→6自由度に

トラス要素では部材端がピンのみでしたが、梁要素では部材端が剛接の場合も考えられます。そこで、節点については回転方向の固定、作用力も考慮できるようにする必要がありました。部材の回転を考えない場合、部材固定度、節点荷重も3自由度にて表現されますが、回転の概念が入ってくることで、6自由度となります。そこで、プログラムのインプットにおいて6自由度を考慮できるようにリストの構成等を修正しました。

・要素をトラス要素→梁要素に

トラス要素は軸力のみを伝達することを考えるため、断面積のみを設定する仕様ですが、フレーム要素にした場合曲げも伝達することになるため、部材各方向の断面二次モーメント、ねじれモーメントを設定する必要があります。Openseesでは弾性梁要素として、elasticBeamColumnコマンドにて以下のように部材特性を設定します。

element('elasticBeamColumn', eleTag, *eleNodes, Area, E_mod, G_mod, Jxx, Iy, Iz, transfTag)

eleTag:部材番号
eleNodes:材端の座標
Area:断面積
E_mod:弾性係数
G_mod:せん断弾性係数
Jxx:ねじり断面二次モーメント
Iy:部材座標系y軸方向に対する断面二次モーメント
Iz:部材座標系z軸方向に対する断面二次モーメント
transfTag:部材要素の座標変換に関する設定(後述)

・部材要素の座標変換に関する設定の追加

全体座標系→部材座標系の変換に以下のコマンドを新たに入力する必要がありました。(通常トラス要素の応力解析でもこのような座標変換は行われていることから、あくまでもOpenseesならではの仕様なのかもしれません。)

transfTag = 1
vecxz =0,0,1
geomTransf('Linear', transfTag, vecxz)
’Linear’:線形での座標変換 transfTag:座標変換を定義するタグ(それぞれの部材要素に設定) vecxz:ユーザーが指定する基本ベクトルであり、部材座標系のx軸と並行であってはならない。 (ここでは解析の対象モデルに鉛直部材がないことを想定して、全体座標系のZ軸方向の単位ベクトルを設定)

上記の設定を各部材要素に対してタグを設定することにより当てはめていきます。

・解析結果の出力

トラス要素の応力解析では部材に軸力のみが生じましたが、梁要素では曲げ及びせん断力が各部材の材端に生じるため、各要素の作用力を出力するbasicforceコマンドを適用すると、1部材につき、6個の数値がリストとして出力されます。そこで、それらのデータを整理し、軸力、曲げモーメント、せん断力のそれぞれについて、各部材のアウトプットをグラデーションによる色分布で図化出力できるように修正しました。前回投稿で記載したものです。

plotlyによる図化出力においては、各部材の線グラフのデータにhovertextとして作用力の数値データを埋め込み、マウスをhoverさせることにより、各部材の作用力が直感的に確認できるようにしました。
ただ、このplotly、表示機能を盛り込みすぎると、動作速度が大変遅くなってしまうのが難点ではあります。(本来はグラフ描画ツールなので。。)みたい結果をその都度絞りながら、結果を表示させていくことが現実的な使い方かもしれません。


以上を考慮することにより、Openseesを使いながら、3次元の任意形状の構造モデルを対象とした梁要素の応力解析プログラムが構成できました。トラス要素の構造解析プログラムとあまり構成が変わらないので、ソースコードの記載は割愛します。上記のポイントを押さえながら比較的簡単に拡張できると思いますので、試してみてください。

2020/05/05

Grasshopperによる3次元ジオメトリのOpenseesによる応力解析

今回は、Openseesにより構成した3次元の応力解析プログラムの実務的な利用について考えてみます。特に複雑な形状モデルを対象とする場合、3次元の形状データから、目視や手入力による作業で応力解析プログラムのインプットを作成することを考えると、膨大で複雑な数値データの入力が必要となることから、大変面倒な作業が必要となってしまいます。そこで、近年建築分野で普及が進んでいる3DcadであるRhinoceros+Grasshopperにより作成した形状データをもとにプログラムによるOpenseesの解析入力ファイルの作成にトライしてみます。また、Rhinoceros+Grasshopperによる入力データの作成および解析結果のKaramba3Dによる解析結果の簡単な検証を試みます。

・Rhinoceros+Grasshopperとは
RhinocerosはRobert McNeel & Associates社が開発した3次元CADであり、NURBS局面などの複雑なジオメトリを柔軟に扱えることから、近年、建築分野でも多く活用されています。またGrasshopperはRhinocerosをスクリプトベースで自由に操ることができるツールで、視覚的なアイコンをつなぐことで容易にプログラムが構成できるため、コンピュテーショナルデザインなどに活用されています。さらに近年は世界中のエンジニア・プログラマーにより開発された最適化や応力解析などの多くの外部モデュールが公開されており、その多くがフリーで利用可能となっています。そうしたツールは以下のサイトfood4Rhinoで利用可能となっています。

https://www.food4rhino.com/

・Karamba3Dとは
Grasshopperには、多くの外部機能がPythonのモジュールのように提供されており、Karamba3DはRhinoceros上のジオメトリを用いた応力解析を容易におこなうことができる応力解析モジュールです。近年は有料化されていますが、(unlimitedlicenceで1150ユーロ)多くの構造エンジニアにも使用され始めています。また、最適化ツールなどの機能拡張が日々行われています。

https://www.karamba3d.com/

まず、応力解析の対象とするトラス屋根モデルを以下のように懸垂円筒状のトラスモデルとして作成してみました。作成に当たってはトラス要素のジオメトリを簡便に作成できる外部モデュール(Lunchbox)を使用しています。Lunchboxについても先述したfood4Rhinoにてフリーで弾雨ロードすることが可能です。


Grasshopper上では、上記のようなアイコンをつなげるだけの簡単なコーディングで以下のような複雑なトラスの構造モデルが作成できます。立体的な幾何学を扱うにはとっても便利なツールです。



このような複雑な形状の構造解析モデルの座標や要素の情報を一つ一つinputファイルに記述するのは例えExcelファイルによる記述であってもとても大変です。

しかし、幸いなことに、Grasshopper上ではPythonで書いたプログラムをコンポーネントとして扱うことができるため、作成したジオメトリデータ(line、座標データなど)とPythonを連携させて使用することできます。そこで、Grasshopper上で作成したジオメトリデータをOpenseesの構造解析用のインプットデータとしての節点、要素の数値データに置換するプログラムをPythonで作成し、そのプログラムをつかってOpenseesによる3次元解析プログラムのインプットデータを作成しました。プログラムといっても、Rhinoceros上のジオメトリの節点座標と、各梁要素の材端座標を拾って、input用のデータ書式に再構成するだけのシンプルなものとなります。

例題として、作成した解析モデルの頂部の中心節点に鉛直荷重が加わった場合の応力解析をOpenseesにより行います。
前回可視化に利用したPlotlyにより、Openseesによる応力解析結果を視覚化するとこんな感じになります。


軸力図

軸力図については圧縮材が赤色、引張材が青色となるように設定しました。


最大曲げモーメントの分布図


そして変位図


Karambaによる変位図

上記のOpenseesによる解析結果とKarambaによる解析結果とを比較すると、変形性状は概ね対応していることがわかります。上記では、Grasshopper→Openseesによる応力解析により複雑な形状の構造物の応力解析ができました。Pythonによるコーディングにおいて幾何学形状を扱うことは労力がかかります。そうしたことから、このように幾何学の扱いに長けており、デザイナーによるコーディングも可能なRhinoceros+Grasshopperをデザインなどで使用し、それらの分析にOpenseesを活用するといったスキームは、その間をつなぐ環境がしっかりと構築されれば、建築設計実務でも十分活用できるのではと考えます。

2020/04/29

Plotlyによる応力解析結果の出力

前回まではmatplotlibを用いた2次元画像による応力解析結果の可視化を試してみました。
しかし、3次元の構造解析モデルの応力解析結果を確認する場合、汎用解析ツールのUIのように
同じく3次元空間でグリグリとオブジェクトを回しながら結果確認ができると実務的に便利です。
そこで今回は、Pythonにも対応しているインタラクティブなデータ可視化ライブラリPlotlyを用いてOpenseesPyによる3次元トラスの応力解析結果の出力を行います。

Plotlyはオープンソースのデータ可視化ライブラリです。
Pythonだけでなく、JAVAscriptやMATLAB等にも対応しています。
同じデータ可視化ライブラリであるmatplotlibは主にグラフデータの画像出力に対応していましたが、
Plotlyでは出力したグラフデータをブラウザ上でマウスなどを動かしながらインタラクティブな状態で確認できます。
Pythonでの使用にはOpenseesPyの場合と同じく、ライブラリのインストールおよびプログラムへのインポートが必要となります。
インストール法や使用方法の詳細につきましては、以下サイトをご覧ください。

https://plotly.com/python/

今回もこれまでに構築した3次元トラス解析プログラムをベースに検討を進めます。
Plotlyを使ってデータをグラフ化するプロセスを以下に簡単に示します。

・ライブラリのインポート

import plotly.offline as po
import plotly.graph_objs as go
上記によりPlotlyをプログラムにインポートします。

・初期化

po.init_notebook_mode()
上記によりPlotlyを初期化します。

・グラフの定義

graph=go.Scatter3d(x=node_x, y=node_y, z=node_z, name='beam%s'%(int(ele_no[i])), mode='lines+text', line=dict(width=int(area[i]),color=color_tag),
hovertext='AxialForce %s'%(round(temp[0],2)), showlegend=True)
上記スクリプトの形式で描画するグラフを定義します。
なお、上記スクリプトはモデル図兼応力図を描くためのスクリプトになります。

モデルの各梁要素を線グラフとして描画します。線の太さは梁要素の断面積のデータから、線の色は応力解析により得られた梁要素の軸力から決定するように定義しています。なお、軸力についてはグラフ上でカーソルを合わせると確認できる「hovertext」として定義しています。

・レイアウトの定義とグラフの描画

data.append(graph)
layout = go.Layout(xaxis=dict(showgrid=False))
fig2 = dict(data = data, layout=layout)
po.plot(fig2, filename="model")
定義したグラフデータを保存したdataリストとグラフのレイアウト等を定義したlayoutリストをfig2に呼び出して、po.plotにて描画しています。


上記のPlotlyによるデータ図化機能を応用して、結果としてグリグリ動かせるモデル軸力図・変形図が出力できました。
ブラウザ上の画面で動かせる形式になっており、
マウスを梁上に移動させると軸力も表示させることができます。


変位図を出力したり、画面上で選択することで、モデル図応力図と変位図の2つのグラフを同時に出現させたりすることもできます。
マウスで図を回転させたり、拡大したりすることが非常に直感的に行えます。



以上、汎用の構造解析ツールのような操作感覚で、解析結果を確認することができました。
今後さらに複雑な立体トラスモデルの応力解析結果を確認する場合には、非常に役立ちそうなツールだと思います。
最後に今回使用したソースコードを以下に示します。
#plotlyによる描画
po.init_notebook_mode()
trace1 = go.Scatter3d(x=cor_x, y=cor_y, z=cor_z, mode='markers+text', marker_color='blue', marker_size=2, name="Node", text=no, showlegend=True)
data = [trace1]

#モデル図・軸力図の描画
for i in range(len(ele_no)):
    temp2=[]
    temp3=[]
    node_x=[]
    node_y=[]
    node_z=[]
    temp2=nodeCoord(int(n1[i]))
    temp3=nodeCoord(int(n2[i]))
    node_x.append(temp2[0])
    node_x.append(temp3[0])
    node_y.append(temp2[1])
    node_y.append(temp3[1])
    node_z.append(temp2[2])
    node_z.append(temp3[2])
    
    temp=basicForce(int(ele_no[i]))#各部材の応力値をRGB値に変換
    if temp[0] >=0:
        if max_tens == 0:
            color_tag='rgb(0,0,'+str(255*temp[0])+')'
        else:
            color_tag='rgb(0,0,'+str(255*temp[0]/max_tens)+')'     
    else:
        if max_comp == 0:
            color_tag='rgb('+str(-255*temp[0])+',0,0)'
        else:  
            color_tag='rgb('+str(-255*temp[0]/max_comp)+',0,0)'
            
    graph=go.Scatter3d(x=node_x, y=node_y, z=node_z, 
                               name='beam%s'%(int(ele_no[i])), mode='lines+text', line=dict(width=int(area[i]),color=color_tag),
                               hovertext='AxialForce %s'%(round(temp[0],2)), showlegend=True)

    data.append(graph)

#変位図の描画
node_x=[]
node_y=[]
node_z=[]
temp2=[]
temp3=[]
for i in range(len(ele_no)):
    temp2=[x + y * mag_factor for (x, y) in zip(nodeCoord(int(n1[i])),nodeDisp(int(n1[i])))]
    temp3=[x + y * mag_factor for (x, y) in zip(nodeCoord(int(n2[i])),nodeDisp(int(n2[i])))]
    node_x.append(temp2[0])
    node_x.append(temp3[0])
    node_y.append(temp2[1])
    node_y.append(temp3[1])
    node_z.append(temp2[2])
    node_z.append(temp3[2])
graph=go.Scatter3d(x=node_x, y=node_y, z=node_z, 
                               name="Disp", mode='lines', line_color='blue',line=dict(width=2, dash='dot'),
                               showlegend=True)
data.append(graph)

layout = go.Layout(xaxis=dict(showgrid=False))
fig2 = dict(data = data, layout=layout)
po.plot(fig2, filename="model")

2020/04/12

Matplotlibによる応力・変形図の図化出力(3D対応)

前回紹介したMatplotlibによる図化出力ですが、以下の手順で3次元トラスの応力解析結果の図化出力にも拡張できました。蛇足ですが説明します。

図化出力のプログラムの頭に以下のコードを記載します。
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

また、グラフのスクリプトを記述する際には、以下を記載します。
ax1 = fig.add_subplot(2,1,1, projection='3d')

以上を考慮するとMatplotlibによるグラフにおいてz軸の概念が導入され、2次元のグラフ描画と同様のルールでグラフの表示カスタマイズ(軸最大値の設定、表示する角度の設定など)も可能となります。Matplotlibによる3次元グラフ描画の詳細については以下のサイトに記載があります。


上記を踏まえて組んだプログラムにより、以下のように3次元のグラフが描画できました。意外といろいろなことができるんですね。Matplotlib。


参考として、今回の図化出力に用いた3次元トラスの応力解析結果の描画部分のソースコードを以下に示します。

#変位図、応力図を描く    
fig = plt.figure(figsize=(6, 9))
max_x=5
max_y=5#グラフの表示範囲を設定
max_z=5

ax1 = fig.add_subplot(2,1,1, projection='3d')
ax2 = fig.add_subplot(2,1,2, projection='3d')
for e in getEleTags():
    color_tag=0
    nodes = eleNodes(e)
    Nnodes = len(nodes)

    xyz = np.zeros((Nnodes,3), dtype=np.double)
    uu   = np.zeros((Nnodes,3), dtype=np.double)
    for i in range(Nnodes):
        xyz[i,:] = nodeCoord(nodes[i])#return coordinate
        uu[i,:] = nodeDisp(nodes[i])#return current displacement

    x = np.zeros(12,dtype=np.double)
    y = np.zeros(12,dtype=np.double)
    z = np.zeros(12,dtype=np.double)    
    u = np.zeros(12,dtype=np.double)
    v = np.zeros(12,dtype=np.double)
    t = np.zeros(12,dtype=np.double) 

    conec = [0, 1]

    x = xyz[conec,0]
    y = xyz[conec,1]
    z = xyz[conec,2]
    x[2::3] = np.nan
    y[2::3] = np.nan
    z[2::3] = np.nan

    mag_factor = 20.#変形図の拡大率
    u = xyz[conec,0] + mag_factor*uu[conec,0]
    v = xyz[conec,1] + mag_factor*uu[conec,1]
    t = xyz[conec,2] + mag_factor*uu[conec,2]
    u[2::3] = np.nan
    v[2::3] = np.nan
    t[2::3] = np.nan
    
    temp=basicForce(e)#各部材の応力をプロット
    if temp[0] >=0:
        if max_tens == 0:
            color_tag=temp[0]
        else:
            color_tag=temp[0]/max_tens*1       
        ax1.text(np.mean(x),np.mean(y),np.mean(z),round(temp[0],2))
        ax1.plot(x,y,z,color=[0,float(1-color_tag),1],lw=area[e-1])
    else:
        if max_comp == 0:
            color_tag=temp[0]
        else:        
            color_tag=-temp[0]/max_comp*1
        ax1.text(np.mean(x),np.mean(y),np.mean(z),round(temp[0],2))
        ax1.plot(x,y,z,color=[1,0,float(1-color_tag)],lw=area[e-1])

    ax2.plot(x,y,z,marker="o",color="0.5")       
    ax2.plot(u,v,t,color="0.",linestyle="dashed")#変形図をプロット

#各グラフのフォーマットの設定
ax1.set_title("Stress Diagram")
ax1.view_init(elev=40, azim=180)
ax1.set_xlim(-1,max_x)
ax1.set_ylim(-1,max_y)
ax1.set_zlim(-1,max_z)
ax1.grid(False)
ax1.set_xticks([]) 
ax1.set_yticks([]) 
ax1.set_zticks([]) 
ax2.set_title("Displacement Diagram")
ax2.view_init(elev=40, azim=180)
ax2.set_xlim(-1,max_x)
ax2.set_ylim(-1,max_y)
ax2.set_zlim(-1,max_z)
ax2.grid(False)
ax2.set_xticks([]) 
ax2.set_yticks([]) 
ax2.set_zticks([]) 
plt.show()

2020/04/11

Matplotlibによる応力・変形図の図化出力

今回はOpenSeesPyによる解析結果の図化出力して遊んでみます。
図化出力の対象は前回までに構築した2次元のトラス解析プログラムの応力解析による結果とします。OpenSeesPyにより出力される節点変位、部材軸力などの数値データを図上に整理します。

図化出力にはpythonのグラフ描画ライブラリMatplotlibを用いることにします。Matplotlibはグラフのみならず様々な図の出力に応用できるようなのです。
私はまだそんなに玄人ではありませんので、少しずつ勉強しながら使っていくことにします。

今回はとりあえず応力解析結果として、部材の軸力図と変形図を出力することにしました。

軸力図は以下のような定義とします。
軸力図における部材の軸力分布はグラフの線色で表現します。
圧縮軸力は赤線、引張軸力は青線で表現、また各部材の軸力の度合いは色のグラデーションで表現します。
さらに、設定した構造部材の断面積を線の太さで表現することにしました。

これらを踏まえて、圧縮軸力となる部材と引張軸力となる部材でそれぞれ場合分けし、軸力図の出力部分は以下のようなコードとしました。

if temp[0] >=0: if max_tens == 0: color_tag=temp[0] else: color_tag=temp[0]/max_tens*1 ax1.text(np.mean(x),np.mean(y),round(temp[0],2)) ax1.plot(x,y,color=[0,float(1-color_tag),1],lw=area[e-1]) else: if max_comp == 0: color_tag=temp[0] else: color_tag=-temp[0]/max_comp*1 ax1.text(np.mean(x),np.mean(y),round(temp[0],2)) ax1.plot(x,y,color=[1,0,float(1-color_tag)],lw=area[e-1])

変形図については、各節点座標とOpenSeesPyのアウトプットコマンドnodeDispで得られた各節点変位の和をグラフ上にプロットします。変位図のプロットの際には、設定した拡大率(mag_factor)により実際の変位を拡大表示します。

以上、組んだプログラムを用いて以下のような応力図が出力できました。


以下のような2次元の任意形状のトラスの応力解析結果の出力に対応します。


このように、色々な形状を設定しながら遊んでみるとトラス構造物の応力の流れが直感的に把握できます。しかし、煩雑なモデルになればなるほどExcelでちまちまと節点や部材の情報を入力するのが面倒になってきます。

そこで今後、CADなどのプログラムと連携させ、より複雑な形状を簡易にインプットできるようなシステムの構築にも挑戦できたらと思います。

今回紹介した図化出力プログラムのソースコードは以下の通りです。
正直、Matplotlibはまだまだ全然使いこなせている感じがしません。
それでも。。
本来はグラフを作成するツールで、枠線、軸線を消したり、2段にまとめたレイアウトとすることによりシンプルな図化出力ができました。
今後も引き続き勉強していきます。

#変位図、応力図を描く    
fig = plt.figure()
max_x=9
max_y=5#グラフの表示範囲を設定

#グラフレイアウトの設定 
ax1 = fig.add_subplot(2,1,1)
ax2 = fig.add_subplot(2,1,2)

#各要素ごとに出力
for e in getEleTags():
    color_tag=0
    nodes = eleNodes(e)
    Nnodes = len(nodes)

    xyz = np.zeros((Nnodes,2), dtype=np.double)
    uu   = np.zeros((Nnodes,2), dtype=np.double)
    for i in range(Nnodes):
        xyz[i,:] = nodeCoord(nodes[i])#return coordinate
        uu[i,:] = nodeDisp(nodes[i])#return current displacement

    x = np.zeros(12,dtype=np.double)
    y = np.zeros(12,dtype=np.double)
    u = np.zeros(12,dtype=np.double)
    v = np.zeros(12,dtype=np.double)

    conec = [0, 1]

    x = xyz[conec,0]
    y = xyz[conec,1]
    x[2::3] = np.nan
    y[2::3] = np.nan

    mag_factor = 10.#変形図の拡大率
    u = xyz[conec,0] + mag_factor*uu[conec,0]
    v = xyz[conec,1] + mag_factor*uu[conec,1]
    u[2::3] = np.nan
    v[2::3] = np.nan
    
    temp=basicForce(e)#各部材の応力をプロット
    if temp[0] >=0:
        if max_tens == 0:
            color_tag=temp[0]
        else:
            color_tag=temp[0]/max_tens*1       
        ax1.text(np.mean(x),np.mean(y),round(temp[0],2))
        ax1.plot(x,y,color=[0,float(1-color_tag),1],lw=area[e-1])
    else:
        if max_comp == 0:
            color_tag=temp[0]
        else:        
            color_tag=-temp[0]/max_comp*1
        ax1.text(np.mean(x),np.mean(y),round(temp[0],2))
        ax1.plot(x,y,color=[1,0,float(1-color_tag)],lw=area[e-1])

    ax2.plot(x,y,marker="o",color="0.5")       
    ax2.plot(u,v,color="0.",linestyle="dashed")#変形図をプロット

#軸ラベルを消す
ax1.tick_params(labelbottom=False,
                labelleft=False,
                labelright=False,
                labeltop=False,
                bottom=False,
                left=False,
                right=False,
                top=False)
ax2.tick_params(labelbottom=False,
                labelleft=False,
                labelright=False,
                labeltop=False,
                bottom=False,
                left=False,
                right=False,
                top=False)

#各グラフのフォーマットの設定
ax1.set_title("Stress Diagram")
ax1.set_xlim(-1,max_x)
ax1.set_ylim(-1,max_y)
ax1.spines["right"].set_color("none")
ax1.spines["left"].set_color("none")
ax1.spines["top"].set_color("none")
ax1.spines["bottom"].set_color("none")
ax2.set_title("Displacement Diagram")
ax2.set_xlim(-1,max_x)
ax2.set_ylim(-1,max_y)
ax2.spines["right"].set_color("none")
ax2.spines["left"].set_color("none")
ax2.spines["top"].set_color("none")
ax2.spines["bottom"].set_color("none")
plt.show()

2020/04/05

2次元トラス解析プログラムの構築(OpenSeesPy)

OpenSeesなどを用いて建築物などの構造解析を行う場合、様々な形状の構造物の応力解析を行う必要があり、そうした構造物を想定する場合どうしても入力する必要のある情報(構造物の節点、部材構成など)の量が増えることになります。そのため、構造物の情報を整理するプログラムをOpenSeesPyモジュールと併用しながら構成すると便利です。そこで、今回は以前の例題で用いたコードを抽象化し、構造物の情報を外部ファイルからインプットしながら整理しOpenSeesPyモジュールに受け渡すソースを組み立て、任意形状の入力に対応した2次元のトラス解析プログラムを構築します。

想定するプログラムのinput、outputは以下の通りです。
input:構造物の節点座標、部材構成、部材断面、節点荷重の位置、場所
output:各節点の座標、部材応力

前回までに用いたトラス構造物を想定した場合、inputファイルはcsvファイルで以下のように構成します。



<ポイント>

・input.csvの読み込み

作成したcsvファイルを扱うため、Pandasモジュールを使用します。Pandasはpythonにて、Excel、csvファイルのデータフレームを扱えるオープンソースのライブラリです。このモジュールを用いることで、csvファイルのデータを読み込むことが可能になります。m使用の際には、別途PCにPandasモジュールをダウンロードし、プログラムの中でインポートする必要があります。
データの読み込み部分は以下のように構成します。

#csvファイルから構造物のデータを読み込み
ip = pd.read_csv('input.csv', header=1)
no=ip['no']
・・・・

上記1行目では、csvファイルのヘッダーを一行読み飛ばし、2行目から各列の名目とデータを取得しています。
2行目はpythonのリストに'no'列(上図A列)のデータを読み込んでいます。
リストに読み込んでしまえば、リストのデータをOpenSeesPyのコマンドに受け渡すことにより、解析ファイルを生成することができます。

・NaN行の非カウント

OpenSeesPyのインプットを作成する際、データフレーム上のNaN行を読み飛ばす必要があります。そこで各インプットデータ作成ループの前に、NaNがない場合に限り、必要ループ数をカウントする構文を追記しました。NaNの判定は、以下のように同じ値を比較することで行えます。

count1=0
for i in range(len(fix_node)):#NaNの場合をカウントしない
    if fix_node[i] == fix_node[i]:
        count1+=1

・解析結果の出力

以下のようにfor文を用いて、input.csvの節点数、部材数に応じて動的に結果出力を行えるようにしました。

# output_results
print("<Node Displacement>")
for i in range(len(no)):
    print("Node%s x:%s y:%s" %(i,nodeDisp(int(no[i]),1),nodeDisp(int(no[i]),2)))
print("<Axial Force of each element>")
for i in range(count2):
    print("N%s=%s" %(i,basicForce(int(ele_no[i]))))

前回対象としたトラス構造物をインプットとした場合、出力される結果は以下の通りになります。しかし、数値だけだと解析結果のイメージが直感的につかめません。そこで、次回はpythonのグラフ描画モジュールmatplotlibを用いて、解析結果を図化出力することにトライしてみようと思います。


今回の検討に用いたソースは以下の通りです。赤字が今回追記した部分になります。


from openseespy.opensees import *

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# ------------------------------
# Start of model generation
# -----------------------------

# remove existing model
wipe()

# set modelbuilder
model('basic', '-ndm', 2, '-ndf', 2)

#input information
ip = pd.read_csv('input.csv', header=1)
no=ip['no']
cor_x=ip['cor_x']
cor_y=ip['cor_y']
fix_node=ip['fix_node']
fix_x=ip['fix_x']
fix_y=ip['fix_y']
ele_no=ip['ele_no']
n1=ip['n1']
n2=ip['n2']
area=ip['area']
load_point=ip['point']
dir_x=ip['dir_x']
dir_y=ip['dir_y']
dat_no=len(ip)#データフレーム数を取得

# create nodes
count=0
for i in range(dat_no):#NaNの場合をカウントしない
    if no[i] == no[i]:
        count+=1
for i in range(count):
    node(int(no[i]), float(cor_x[i]), float(cor_y[i]))

# set boundary condition
count1=0
for i in range(dat_no):#NaNの場合をカウントしない
    if fix_node[i] == fix_node[i]:
        count1+=1
for i in range(count1):
    fix(int(fix_node[i]), int(fix_x[i]), float(fix_y[i]))

# define materials
uniaxialMaterial("Elastic", 1, 3000.0)

# define elements
count2=0
for i in range(dat_no):#NaNの場合をカウントしない
    if ele_no[i] == ele_no[i]:
        count2+=1
for i in range(count2):
    element("Truss",int(ele_no[i]),int(n1[i]),int(n2[i]),float(area[i]),1)

# create TimeSeries
timeSeries("Linear", 1)

# create a plain load pattern
pattern("Plain", 1, 1)

# Create the nodal load - command: load nodeID xForce yForce
count3=0
for i in range(dat_no):#NaNの場合をカウントしない
    if load_point[i] == load_point[i]:
        count3+=1
for i in range(count3):
    load(int(load_point[i]),float(dir_x[i]),float(dir_y[i]))

# ------------------------------
# Start of analysis generation
# ------------------------------

# create SOE
system("BandSPD")

# create DOF number
numberer("RCM")

# create constraint handler
constraints("Plain")

# create integrator
integrator("LoadControl", 1.0)

# create algorithm
algorithm("Linear")

# create analysis object
analysis("Static")

# perform the analysis
analyze(1)

# output_results
max_tens=0
max_comp=0

print("<Node Displacement>")
for i in range(count):
    print("Node%s x:%s y:%s" %(i,nodeDisp(int(no[i]),1),nodeDisp(int(no[i]),2)))
print("<Axial Force of each element>")
for i in range(count2):
    print("N%s=%s" %(i,basicForce(int(ele_no[i]))))
    temp=basicForce(int(ele_no[i]))
    if temp[0] >= 0:
        if temp[0] > max_tens:
            max_tens = temp[0] 
    else:
        if abs(temp[0]) > max_comp:
            max_comp = abs(temp[0])

2020/03/29

OpenSeesPyによる静定トラスの解析(材料非線形性の考慮)

通常、構造物の材料は材料非線形性を有しており、部材荷重がある降伏荷重を超えると、部材が降伏することで部材の荷重増加率が低減し、部材の変形が大きく増加する状態になります(降伏現象)。Openseesではこうした部材の材料非線形性を考慮した応力解析が可能です。今回は前回と同形状の静定トラスを対象として、その部材が材料非線形性を考慮したモデルを題材とし、各部材の降伏を追跡した応力解析による検討を行います。

今回考慮する部材の材料非線形性は下図に示すような1つの降伏点を有するバイリニアの材料非線形性とします。部材の降伏荷重sYや降伏後の部材剛性低下率αを定義することによって、部材の材料非線形性を解析において考慮できます。


今回は以下に示すパラメータを有する架空の材料を想定します。
ヤング係数: E = 3000
降伏荷重: sY=10
剛性低下率: α=0.05

OpenSeesPy上で上記の材料特性を記載すると、以下のようになります。

E=3000
sY = 10
alpha = 0.05
uniaxialMaterial("Steel01", 1, sY, E, alpha)

()内は部材タイプ、部材特性の番号、降伏荷重、ヤング係数、剛性低下率の順に記載します。上記のコマンドを入力することで、タグ1が定義されている部材全てにおいて、設定した非線形の材料特性が適用されます。
なお、「uniaxialMaterial」はトラスなどの軸方向のみ剛性を有する部材の材料特性を定義するコマンドです。今回導入するバイリニアの材料特性以外にも様々な部材タイプの材料特性が用意されています。それぞれの材料特性に対し設定すべきパラメータが定義されており、それらを設定することにより様々な材料特性を考慮した解析を行うことが可能です。詳しくは以下のリンク先をご覧ください。また、そのうち整理してみたいと思います。(コンクリートなどの材料特性を想定したモデルも用意されています。世界中の研究者が提案している最新の材料特性モデルなども多く揃えられています。)

https://openseespydoc.readthedocs.io/en/latest/src/uniaxialMaterial.html

基本的には上述したコマンドを加えるだけで解析モデルにおける材料非線形性が定義され、部材の材料非線形性を考慮した応力解析が可能です。さらに、今回はトラス構造物における各部材の降伏に伴う節点変位を追跡するため、設定した節点荷重に至るまでの各荷重係数(λ=0~1)を想定した場合の結果を出力する荷重増分解析を行います。そのため、解析実施部分のスクリプトを以下のように定義します。

integrator("LoadControl", 1.0/Nsteps)

integratorコマンドを入力することで静的荷重増分解析を行うことができます。増分解析を荷重制御にて行う場合,上記のように記述し解析の分割数Nsteps=100を想定します。また各stepにおける頂部節点変位と外力を出力するため、以下のようにfor文を記述することにより、解析の実施及び解析結果の出力を行います。

Nsteps=100
data = np.zeros((Nsteps+1,2))
for j in range(Nsteps):
    analyze(1)
    data[j+1,0] = -nodeDisp(2,2)
    data[j+1,1] = getLoadFactor(1)*150

nodedispコマンドはnodeDisp(対象節点番号、変位方向)を入力することで各節点の変位を取得できます。また、getloadfactorコマンドでは各解析時の荷重係数をそれぞれ取得しています。結果リストをグラフ描画モデュールmatplotlibにより処理することで、以下のような荷重変位関係が取得できます。

グラフより節点の鉛直荷重(Vertical Load)が増加することにより、部材の降伏荷重に達し、その後節点の鉛直変形(Vertical Displacement)が大きく増加していることが確認できます。また、対象トラス構造物における部材1,3の降伏やその後の部材2の降伏を荷重増分解析により適切に追跡できていることが確認できます。


以下に、今回の検討に用いたソースコードを示します。基本は前回のソースコードと同様ですが、赤字が主に変更した部分になります。


from openseespy.opensees import *

# ------------------------------
# モデルの作成
# -----------------------------

# モデルの初期化
wipe()

# モデルの定義
model('basic', '-ndm', 2, '-ndf', 2)

# 節点の定義
node(1, 0.0, 0.0)
node(2, 4.0,  4.0)
node(3, 8.0,  0.0)

# 境界条件の定義
fix(1, 1, 1)
fix(3, 0, 1)

# 部材材料の定義
E=3000
sY = 10
alpha = 0.05
uniaxialMaterial("Steel01", 1, sY, E, alpha)

# 部材の定義
element("Truss",1,1,2,5.0,1)
element("Truss",2,1,3,5.0,1)
element("Truss",3,2,3,5.0,1)

# 時間に依存する荷重係数の定義
timeSeries("Linear", 1)

# 荷重パターンの定義
pattern("Plain", 1, 1)

# 節点荷重の定義
load(2, 0, -150)

# ------------------------------
# 解析条件の定義、解析実行
# ------------------------------

# create SOE
system("BandSPD")

# create DOF number
numberer("RCM")

# create constraint handler
constraints("Plain")

# create integrator
Nsteps = 100
integrator("LoadControl", 1.0/Nsteps)

# create algorithm
algorithm("Newton")

# create analysis object
analysis("Static")

# perform the analysis
data = np.zeros((Nsteps+1,2))
for j in range(Nsteps):
    analyze(1)
    data[j+1,0] = -nodeDisp(2,2)
    data[j+1,1] = getLoadFactor(1)*150

# ------------------------------
# 解析結果のグラフ出力(Matplotlibによる)
# ------------------------------
plt.plot(data[:,0], data[:,1])
plt.xlabel('Vertical Displacement')
plt.ylabel('Vertical Load')
plt.show()

2020/03/28

OpenSeesPyによる静定トラスの解析

今回は図に示すように、静定トラスの頂部に静的な荷重が加わった場合の応力解析をOpenSeesPyを用いて行っていきます。


・解析モデルの設定
modelコマンドにて全体の解析モデルを定義します。
ndmは解析モデルの次元数、ndfは自由度数です。
今回は2次元の静定トラスを想定するので、いずれも2に設定します。

model('basic', '-ndm', 2, '-ndf', 2)

・トラスの節点座標の設定
nodeコマンドにて解析モデル各節点の座標を定義します。
2次元モデルの場合、node(節点番号、x座標、y座標)のように定義します。

・固定されている節点の境界条件の設定
fixコマンドにて解析モデルの各節点の境界条件を定義します。
2次元モデルの場合、fix(節点番号、x方向の境界条件、y方向の境界条件(いずれも0:フリー,1:固定)と定義します。

・トラス部材の材料の条件を設定
今回は弾性のトラス部材を想定するため、UniaxialMaterialコマンドにて部材のヤング係数のみを定義します。

・トラス部材を設定
elementコマンドにて各部材の特性を定義します。
トラス部材の場合、以下のように定義できます。
element(elementtype、部材番号、部材端部の節点番号、断面積、定義した材料番号)

・外力を設定
節点に生じる外力をloadコマンドにて以下のように定義します。
load(設定する節点番号、x方向の外力、y方向の外力)

上記が対象とする問題を想定した場合の主な解析モデル、荷重条件の設定項目になり、解析に関する種々の条件を定義したうえでanalyzeコマンドを実行すると、応力解析が実行されます。解析結果は「Output Commands」に用意されている各種コマンドを用いて確認することができます。例えば、確認したい項目が、節点変位であればnodeDispコマンド、部材応力であればbasicForceコマンドの結果をprintすることで確認できます。今回の場合、basicForceコマンドにより出力された値が手計算の値と一致することから、正しく解析が行われていることが確認できました。



以下、今回の検討に用いたソースコードを示します。


from openseespy.opensees import *


# ------------------------------
# モデルの作成
# ----------------------------

# モデルの初期化
wipe()

# モデルの定義
model('basic', '-ndm', 2, '-ndf', 2)


# 節点の定義
node(1, 0.0, 0.0)
node(2, 4.0,  4.0)
node(3, 8.0,  0.0)

# 境界条件の定義
fix(1, 1, 1)
fix(3, 0, 1)

# 部材材料の定義
uniaxialMaterial("Elastic", 1, 3000.0)

# 部材の定義
element("Truss",1,1,2,5.0,1)
element("Truss",2,1,3,5.0,1)
element("Truss",3,2,3,5.0,1)

# 時間に依存する荷重係数の定義
timeSeries("Linear", 1)

# 荷重パターンの定義
pattern("Plain", 1, 1)

# 節点荷重の定義
load(2, 0, -50)

# ------------------------------
# 解析条件の定義、解析実行
# ------------------------------

# create SOE
system("BandSPD")

# create DOF number
numberer("RCM")

# create constraint handler
constraints("Plain")

# create integrator
integrator("LoadControl", 1.0)

# create algorithm
algorithm("Linear")

# create analysis object
analysis("Static")

# perform the analysis
analyze(1)

# ------------------------------
# 解析結果の確認
# ------------------------------

N1=basicForce(1)
N2=basicForce(2)
N3=basicForce(3)

print("N1=",N1,"N2=",N2,"N3=",N3)

OpenSeesPyの導入②

今回は、実際にpythonにおいてOpenSeesPyを扱う上で必要なポイントについて説明したいと思います。

・Pythonプログラムの作成環境
前回インストールしたopensees anacondaにはプログラムエディタとして、JupyterNotebookとSpyderが付属しています。

・JupiterNotebook
Webブラウザ上で動作し実行結果がすぐにみられる対話型実行環境。

・Spyder
科学者、エンジニア向けにデザインされたPythonの統合開発環境。テキストエディタ形式でプログラムの編集および実行が可能。


Pythonはコンパイルを不要とするインタプリタ言語であるため、いずれもプログラム記述後すぐに実行結果がみられます。本ブログの検討では基本的にSpyderをベースとした環境にてプログラムを構築していきます。

・Pythonプログラム上でのOpenSeesPyのインポート
Pythonプログラムには様々な機能を提供する標準モジュールが「ライブラリ」として用意されており、それらを自身の作成プログラムにimportすることにより、シンプルなコーディングで様々な機能を有するプログラムを作成することが可能になります。「OpenSeesPy」もそうしたライブラリの一つであり、プログラムの冒頭で以下のようにimportすることにより、使用が可能になります。


・OpenSeesの単位系について
OpenSeesはアメリカ合衆国のPacific Earthquake Engineering Research Center (PEER) によって開発されたため、基本的にヤードポンド法をベースとしており、使用の際にはSI単位系に換算したパラメータを用いる必要があります。 

2020/03/22

OpenSeesPyの導入①

Openseesはアメリカ合衆国のPacific Earthquake Engineering Research Center (PEER) が開発したオープンソース型の構造解析ツールです。オープンソースであることから、全てのコードは公開されています。またその内容は世界中のエンジニアの手で日々更新されており、更新履歴は全てweb上で公開されています。
静的な応力を対象とした構造解析に加え、地震波を入力した時刻歴応答解析などが可能で、全てフリーで使えることから世界中の研究者が解析ツールとして利用しています。詳細は以下をご覧ください。

https://opensees.berkeley.edu/

近年、pythonなどが気軽に使えるプログラムとして台頭しており、openseesもpythonプログラム上で導入可能なpythonモジュール「OpenSeesPy」として公開されています。
各コマンドの機能を解説したマニュアルについても、英語版は下記リンク先にて公開されています。

今回はまずその導入方法を紹介したいと思います。

①pythonを導入
pythonの主要なモデュールを含んだanacondaを以下ページより導入します。

https://www.anaconda.com/distribution/#download-section

開いたページの「Download」をクリックします。

②Jupyter Notebookを起動
Jupyter Notebookはpythonコードを書くと、その実行結果がすぐにみられるツールです。同ツールを開き、ツール上でそれぞれの場合において以下のコードを記述します。


・新規インストール
python -m pip install openseespy
python -m pip install --user openseespy
・ソフトウェアの更新
python -m pip install --upgrade openseespy
python -m pip install --user --upgrade openseespy
上記のように入力すると自動的にモデュールのインストールが始まります。

③pythonコード上にモデュールをインポート
以下のコードをPythonコードの冒頭に記載することで、python上でOpenSeesPyモデュールの使用が可能となります。

import openseespy.opensees as ops

上記の操作により、Pythonコーディング上でOpenseesの構造解析、時刻歴応答解析ツールを動かすための環境構築が終了しました。次は、Python上でOpenseesを動かす方法について紹介していきたいと思います。