OpenSeesPyによる立体フレームモデル解析

今回は，三次元フレームモデルを作って，構造解析をします．
また，より複雑なモデル作成にも使えるよう，部材管理，部材番号管理のコツも公開します！！
全体の流れは次のとおりです．

目次

問題設定
部材管理，材料管理のコツ
プログラムの構成

1. 問題設定

下の図のような構造とします．部材の断面などはきちんと検討しておらず，適当に設定しました．

荷重は，各梁に$10~\mathrm{kN/m}$の等分布荷重を載荷します．

2.部材管理，材料管理のコツ

部材管理のコツ

OpenSeesPyでは，各部材要素(element)，節点(node)に固有の番号(Tag)が割り当てられています．例えば，線形弾性の柱梁部材を定義する際は，

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

といった形で，eleTagとeleNodesに部材番号と節点番号を入力するという仕組みになっています．
今回の建物のような単純な形状であれば頭の中やノートに書いて簡単に管理することができますが，より複雑な建物になるとどの番号がどこを示しているのか解りづらくなります．
例えば，後で2階の梁に荷重をかけたり，節点の変位を出力するといった場合に，要素番号を直接参照していてはプログラムの可読性が下がってしまいます．

そこで，Pythonの辞書型を使って各部材，場合によっては節点も管理しやすい番号で参照することにしつつ，OpenSeesPyの入力上ではプログラム内で作成された順番に番号をつける，という方法で管理することで，プログラムの可読性を上げ，後で参照する際にも便利になります．

材料管理のコツ

OpenSeesPyでは部材を定義する際に，段面積やヤング率，断面二次モーメントなどのパラメータを入力します．今回はH型鋼と角形鋼管を使用しますが，H型鋼と角形鋼管クラスを定義し，実際に材料を使用して部材を定義する際に使用することとします．

class box_beam():
    """
    角型鋼管の断面情報
    - wy 幅 m
    - wz 幅 m
    - t 厚さ m
    - E ヤング係数N/mm^2
    - G せん断弾性係数 N/mm^2
    """
    def __init__(self, wy, wz, t, E=205000, G=79000, ):
        self.wy = wy
        self.wz = wz
        self.t = t
        self.E = E*10**6
        self.G = G*10**6

    def A(self):
        return self.wy*self.wz - (self.wy - 2*self.t)*(self.wz - 2*self.t)
    
    def Iy(self):
        return self.wy*self.wz**3/12 - (self.wy-2*self.t)*(self.wz-2*self.t)**3/12
    
    def Iz(self):
        return self.wz*self.wy**3/12 - (self.wz-2*self.t)*(self.wy-2*self.t)**3/12
    
    def Zy(self):
        return self.Iy()/self.wz*2
    
    def Zz(self):
        return self.Iz()/self.wy*2
    
    def S(self):
        """
        中心線で囲まれる領域の面積
        """
        return (self.wy - self.t) * (self.wz - self.t)
    
    def Jxx(self):
        """
        ねじり定数
        """
        a = self.wy - self.t
        b = self.wz - self.t
        return 2*self.t*self.S()**2/(a + b)
    
    def section(self):
        return [self.A(), self.E, self.G, self.Jxx(), self.Iy(), self.Iz()]
    
    def section_JFE(self):
        # JFE鋼構造便覧断面表に掲載の単位で出力
        # ただし，便覧の数値はrなどが効いているため，数値は微妙に異なる．
        return [self.A()*10**4 #cm^2
                , self.E/10**6 # N/mm^2
                , self.G/10**6 # N/mm^2
                , self.Jxx()*10**12 # mm^4
                , self.Iy()*10**8 # cm^4
                , self.Iz()*10**8 # cm^4
                  ]



class H_beam():
    """
    H型鋼の断面情報
    - h せい(m)
    - b 幅 (m)
    - tw ウェブ厚 (m)
    - tf フランジ厚 (m)
    - E ヤング係数N/m^2
    - G せん断弾性係数 N/m^2
    """

    def __init__(self, h, b, tw, tf, E=205000*10**6, G=79000*10**6, ):
        self.h = h
        self.b = b
        self.tw = tw
        self.tf = tf
        self.E = E
        self.G = G
    
    def A(self):
        """
        断面積
        """
        return (self.h - 2*self.tf)*self.tw + 2*(self.b * self.tf)

    def Iy(self):
        return (self.b*self.h**3 - (self.b - self.tw)*(self.h - 2*self.tf)**3) / 12
    
    def Iz(self):
        return (2 * self.b**3 * self.tf + (self.h - 2*self.tf)*self.tw**3) / 12
    
    def Zy(self):
        return self.Iy()/self.h*2
    
    def Zz(self):
        return self.Iz()/self.b*2
    
    def Jxx(self):
        return 1/3 * (2 * self.b * self.tf**3 + (self.h - 2*self.tf)*self.tw**3)
    
    def section(self):
        # N, m単位型
        return [self.A(), self.E, self.G, self.Jxx(), self.Iy(), self.Iz()]
    
    def section_JFE(self):
        # JFE鋼構造便覧断面表に掲載の単位で出力
        # ただし，便覧の数値はrなどが効いているため，数値は微妙に異なる．
        return [self.A()*10**4 #cm^2
                , self.E/10**6 # N/mm^2
                , self.G/10**6 # N/mm^2
                , self.Jxx()*10**12 # mm^4
                , self.Iy()*10**8 # cm^4
                , self.Iz()*10**8 # cm^4
                  ]

3. プログラムの構成

節点の定義

節点番号管理の辞書node_dictを作り，各節点番号と辞書のキーf-x-yと対応させます．例えば，2階のX1Y1通りの接点なら2-1-1となります．
場合によってはxやyといった文字を入れて2f-x1-y1といったように管理してもわかりやすいでしょう．(windowsとmacでライブラリの名称が違うので適宜コメントアウトなど解除して実行してください．)

import openseespymac.opensees as ops　# mac
# import openseespy.opensees as ops　# windows
import opsvis as opsv
import numpy as np
import matplotlib.pyplot as plt

# 固定・自由，座標系番号
class opensees_constants:
    def __init__(self):
        self.free = 0
        self.fixed = 1

        self.X = 1
        self.Y = 2
        self.Z = 3
        self.RX = 4
        self.RY = 5
        self.RZ = 6
opc = opensees_constants()

ops.wipe()
ndm = 3 # 次元
ops.model('basic', '-ndm', ndm, '-ndf', int(ndm*(ndm+1)/2))

numBayX = 1
numBayY = 1
floor = 2
spanX = 6.4
spanY = 6.4
floor_height = [0, 4.3] + [4.3 + 3.7*_ for _ in range(1, floor + 1)]

# 柱梁サイズの設定 (mm)
C1 = np.array([500, 500, 16]) # mm
C2 = np.array([500, 500, 22]) # mm
C1 = box_beam(*C1/1000)
C2 = box_beam(*C2/1000)

B45 = np.array([440, 300, 11, 18])
B45 = H_beam(*B45/1000)

node_dict = dict()
node_cnt = 0
beam_dict = dict()
column_dict = dict()
ele_cnt = 0


# 節点の作成
for f in range(floor + 1):
    for x in range(numBayX + 1):
        for y in range(numBayY + 1):
            node_dict[f'{f}-{x}-{y}'] = node_cnt
            ops.node(node_cnt, spanX*x, spanY*y, floor_height[f])
            node_cnt += 1

部材の配置

# 拘束条件 (CONSTRAINTS)
# 3Dでの固定支持: tag, DX, DY, DZ, RX, RY, RZ
# すべての基礎節点を全拘束
# 鉛直方向の拘束は Z軸 (DZ) に移動
for key in node_dict.keys():
    if key[0] == '0' and key[1] != '0':
        ops.fix(node_dict[key], opc.fixed, opc.fixed, opc.fixed, opc.fixed, opc.fixed, opc.fixed)

# ここでは X軸方向 [1.0, 0.0, 0.0] を指定。
column_geom = 1
beam_x_geom = 2
beam_y_geom = 3
ops.geomTransf('Linear', column_geom, 0.0, 1.0, 0.0)  # 柱用
ops.geomTransf('Linear', beam_x_geom, 0.0, 0.0, 1.0) # 梁X用
ops.geomTransf('Linear', beam_y_geom, 0.0, 0.0, 1.0) # 梁Y用


# 柱の配置
for f in range(floor):
    for x in range(numBayX+1):
        for y in range(numBayY+1):
            column_dict[f'{f}-{f+1}-{x}-{y}'] = ele_cnt
            i_end = f'{f}-{x}-{y}'
            j_end = f'{f+1}-{x}-{y}'
            
            ops.element('elasticBeamColumn', ele_cnt, node_dict[i_end], node_dict[j_end], *C1.section(), column_geom)
            ele_cnt += 1

# Y梁の配置
for f in range(1, floor + 1):
    for x in range(numBayX + 1):
        for y in range(numBayY):
            beam_dict[f'{f}-{x}-{y}_{x}-{y + 1}'] = ele_cnt
            i_end = f'{f}-{x}-{y}'
            j_end = f'{f}-{x}-{y+1}'
            ops.element('elasticBeamColumn', ele_cnt, node_dict[i_end], node_dict[j_end], *B45.section(), beam_y_geom)
            ele_cnt += 1

# X梁の配置
# 梁の配置
for f in range(1, floor + 1):
    for x in range(numBayX):
        for y in range(numBayY + 1):
            beam_dict[f'{f}-{x}-{y}_{f}-{x + 1}-{y}'] = ele_cnt
            i_end = f'{f}-{x}-{y}'
            j_end = f'{f}-{x + 1}-{y}'
            ops.element('elasticBeamColumn', ele_cnt, node_dict[i_end], node_dict[j_end], *B45.section(), beam_x_geom)
            ele_cnt += 1

荷重の設定と解析実行

ここは2次元の時とさほど変わりません．

# 10 kN/m^2を仮定して荷重をかける
ops.timeSeries('Linear',1,'-factor',1.0)
ops.pattern('Plain', 1, 1)

for key in beam_dict.keys():
    ends = key.split('_')
    # print(ends)
    i_end = ends[0]
    j_end = ends[1]
    beam_tag = beam_dict[key]
    Wx = 0
    Wy = 0
    # Wz = (-10*10**3+1.800)*spanX/2 # N/m
    Wz = -6000*spanX/2 # N/m

    ops.eleLoad('-ele',beam_tag, '-type', '-beamUniform', Wx, Wz, Wy) # 座標系ごちゃごちゃだが，これでglobal z方向負に載荷

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

解析モデル，結果の表示

opsv.plot_model()

解析モデル

変形

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection='3d')

opsv.plot_defo(fig_wi_he=(50,50), ax=ax)
ax.view_init(elev=20, azim=-30)

曲げモーメント図

opsv.section_force_diagram_3d('My',sfac=0.00001)

今回の記事は以上です．