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