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