OpenSeesPyによる建築構造解析の基礎2

このページに掲載されている内容によって生じた損失や損害に対して一切責任を負いません．
もしコードの誤りなどありましたらご連絡ください．

OpenSeesPyの解析結果の可視化

前回作成した構造解析の解析モデルや結果を可視化する手法について説明します．

1. 可視化に用いるライブラリ

OpenSeesPyの解析結果の可視化に用意されているライブラリは次の2種類があります．

1. opsvis (opensees visualization) https://opsvis.readthedocs.io/en/latest/
2. vfo (visualization for openseespy) https://vfo.readthedocs.io/en/latest/

どちらもOpenSeesPyの解析モデルや解析結果の可視化のために開発されたライブラリです．

一つ目のopsvisはmatplotlibをもとに開発されているもので，質点系やトラス，柱梁モデル等幅広い解析モデルの可視化に特化しているのが特徴です．

それに対して二つ目のvfoは二次元，三次元の柱梁モデル等の解析結果の可視化に特化しており，表示した解析結果をマウスを使ってCAD上のモデルのように視点を移動したりズームアウト，ズームインが簡単にできるのが特徴です．

今回の記事では，前回(https://qiita.com/srk1-bot/items/6a1a636905382f6518e4)の記事で使用した構造解析のプログラムをもとにopsvisを用いて解析モデルと解析結果の可視化を行います．

2. opsvisの使い方

インストールはコマンドプロンプトに次の命令を打てばできます．

py -m pip install opsvis

一部，Pythonのバージョンによってはインストールしても以下のプログラムでエラーが出る場合があります．

(自分の場合は，anacondaのPython3.8.10にインストールして実行したらエラーが出ました．pythonのアップデートをすれば解決すると思われますが，Pythonのアップデートがうまくいかず，いまだに解決していません．どなたか教えてください．．．)

ライブラリのインポートは次のように行います．

import openseespy.opensees as ops # openseespyのインポート
import opsvis as opsv # opsvisのインポート

import matplotlib.pyplot asplt

2.1 解析モデルの描画

まず，前回の記事(https://qiita.com/srk1-bot/items/6a1a636905382f6518e4)のコードのうち前半部分の解析モデルの設定のコードは次の通りです．それぞれの文法は前回の記事をご覧ください．

# 初期化
ops.wipe()
# 解析モデルの設定
ops.model('basic', '-ndm', 2, '-ndf', 3)

# 節点の定義
p0 = [0.0, 0.0]
p1 = [0.0, 2.0]
p2 = [1.0, 2.0]
ops.node(0, *p0)
ops.node(1, *p1)
ops.node(2, *p2)

# 境界条件
const0 = [1, 1, 1]
ops.fix(0, *const0)

# 材料の定義
ops.uniaxialMaterial('Elastic', 0, 3000.0)

# geotransf
ops.geomTransf('Linear', 0) # よくわからん

# 部材の定義
A = 5.0 # 断面積
E = 3000 # ヤング率
Iz = 2 # 断面二次モーメント
ops.element('elasticBeamColumn', 0, *[0, 1], A, E, Iz, 0)
ops.element('elasticBeamColumn', 1, *[1, 2], A, E, Iz, 0)

このコードで定義したフレームの解析モデルをopsvisを用いて描画します．モデルの描画には，
opsv.plot_model()メソッドを用います．節点をを表示したりしなかったり，部材を表示したりしなかったりと，細かい設定をする引数が用意されていますが，基本的には引数を入れずに実行しても動くはずです．上のコードを実行したのち，次のコードで解析モデルを描画します．

opsv.plot_model()
plt.show()

すると，次の画像が表示されるはずです．

.ipynbファイルであれば一行目のopsv.plot_modelのみで動くはずですが，.pyファイルではplt.show()と命令してあげると，別ウィンドウで画像が表示されるはずです．

2.2 荷重の描画

ここではモデルに加えた荷重の描画方法を説明します．荷重の描画にはopsv.plot_loads_2d()メソッドを使用します．このメソッドも引数がたくさん設定されており，描画に関する細かい設定ができるようになっていますが，ここでは特に気にせずに実行すればいいと思います．

荷重を設定したのち，構造解析を実行し，設定した荷重を描画するコードは次の通りです．

# 荷重の定義
ops.timeSeries('Constant', 0) # 時間依存性タイプの定義
ops.pattern('Plain', 0, 0) # 荷重タイプの定義
ops.load(2, *[0, -10, 0]) # 荷重の定義

# define and run analysis 

# create SOE
ops.system("BandSPD")

# create DOF number
ops.numberer("Plain")

# create constraint handler
ops.constraints("Plain")

# create integrator
ops.integrator("LoadControl", 1.0)

# create algorithm
ops.algorithm("Linear")

# create analysis object
ops.analysis("Static")

# perform the analysis
ops.analyze(1)


# result of the analysis

N1=ops.basicForce(0)
N2=ops.basicForce(1)
print("N0=",N1,"N1=",N2)

# 荷重の描画
opsv.plot_loads_2d()
plt.show()

# 荷重の描画
opsv.plot_loads_2d()
plt.show()

これにより次の画像が出力されると思います．

2.3 変形の描画

変形の描画には，opsv.plot_defo(sfac)メソッドを使用します．このメソッドにも様々な引数が用意されていますが，ここではsfacだけを理解しておけば大丈夫です．sfac = (float)はscale factorの略で，変形を実際の何倍して描画するかを指定する引数です．ここではsfac = 100として描画してみましょう．

opsv.plot_defo(sfac = 100)
plt.plot()

すると次の画像が出力されます．

2.4 曲げモーメント図，軸力図，せん断力図の描画

曲げモーメント図，軸力図，せんだん力図の描画には，opsv.section_force_diagram_2d(sf_type, sfac, nep=17)メソッドを使用します．主要な引数は次の三つです．

1. sf_type | (str) 図の種類を指定するメソッド．曲げモーメント図，せんだん力図，軸力図の時がそれぞれsf_type='M', sf_type='V', sf_type='N'に対応する．
2. sfac | (float)何倍の大きさで曲げモーメント図を描画するかここではsfac=0.01とする．
3. nep=17 | (int)部材の両端を含んで何等分点で部材力を評価するかを指定する整数で，デフォルトでは17とせってされている．

# Bending Moment
opsv.section_force_diagram_2d('M', sfac = 0.01)
# shear Force
opsv.section_force_diagram_2d('V', sfac = 0.01)
# Axial Force
opsv.section_force_diagram_2d('N', sfac = 0.01)

これを実行すると次の3つの画像が出力される．

曲げモーメント図

せんだん力図

軸力図

3 まとめ

今回使用したコードを以下にまとめておく．

import openseespy.opensees as ops
import opsvis as opsv
import matplotlib.pyplot as plt
ops.wipe()
ops.model('basic', '-ndm', 2, '-ndf', 3)

p0 = [0.0, 0.0]
p1 = [0.0, 2.0]
p2 = [1.0, 2.0]
ops.node(0, *p0)
ops.node(1, *p1)
ops.node(2, *p2)

# 境界条件
const0 = [1, 1, 1]
ops.fix(0, *const0)

# 材料の定義
ops.uniaxialMaterial('Elastic', 0, 3000.0)

# geotrans
ops.geomTransf('Linear', 0) # よくわからん

# 部材の定義
A = 5.0 # 断面積
E = 3000 # ヤング率
Iz = 2 # 断面二次モーメント
ops.element('elasticBeamColumn', 0, *[0, 1], A, E, Iz, 0)
ops.element('elasticBeamColumn', 1, *[1, 2], A, E, Iz, 0)

ops.timeSeries('Constant', 0)
ops.pattern('Plain', 0, 0)
ops.load(2, *[0, -10, 0])

# define and run analysis 

# create SOE
ops.system("BandSPD")

# create DOF number
ops.numberer("Plain")

# create constraint handler
ops.constraints("Plain")

# create integrator
ops.integrator("LoadControl", 1.0)

# create algorithm
ops.algorithm("Linear")

# create analysis object
ops.analysis("Static")

# perform the analysis
ops.analyze(1)


# result of the analysis

N1=ops.basicForce(0)
N2=ops.basicForce(1)
print("N0=",N1,"N1=",N2)

# -----------^-^-------------
# 解析モデルの描画
opsv.plot_model()
plt.show()

# -----------^-^-------------
# 荷重の描画
opsv.plot_loads_2d()
plt.show()

# -----------^-^-------------
# 変形の描画
opsv.plot_defo(sfac = 100)
plt.show()

# -----------^-^-------------
# Bending Moment
opsv.section_force_diagram_2d('M', sfac = 0.01)
plt.show()
# shear Force
opsv.section_force_diagram_2d('V', sfac = 0.01)
plt.show()
# Axial Force
opsv.section_force_diagram_2d('N', sfac = 0.01)
plt.show()