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])