OpenSeesPyを用いた1質点系の時刻歴応答解析その1

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

この記事では，OpenSeesPyを用いた1質点系の弾性時刻歴応答解析を実装する方法について解説します．

1. 例題

今回考える問題は次の通りです．地面につながれた節点1，質量$M = 200 \mathrm{kg}$の質点がある節点2がヤング率$E = 205000\mathrm{N/mm^2}$，断面積$A = 200\mathrm{mm^2}$の部材でつながれています．さらに，減衰定数$h=0.05$のダンパーがついています．この時地面に地震動が作用した時の質点の動きを時刻歴応答解析を用いてシミュレーションします．

2. 応答解析のコード

2.0 固定，座標系クラスの作成

openseespyを使っていると，拘束条件で固定が$0$, 自由は$1$であったり，座標系が$x, y, Rz$の順番に$1, 2, 3$であったりするなど，コーディングの際や後で読み返す際にややこしいと感じますよね．

そこで，今回は次のクラスを定義することで，この後のコードの可読性を上げたいと思います．

# import ilbrary
import openseespy.opensees as ops
import numpy as np

class opensees_constants:
    def __init__(self):
        self.free = 0
        self.fixed = 1

        self.X = 1
        self.Y = 2
        self.ROTZ = 3

opc = opensees_constants()

2.1 モデルの設定

まずモデルの設定をします．空間の次元が$x, y$の2次元で，節点の自由度が$x, y, Rz$の三次元なので以下のようになります．

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

続いて，解析モデルの設定を行います．節点1(bottom)が支点，節点2(top)が質点sの位置です．今回はzero_length_elementと呼ばれる，長さ0の仮想的なばねを用いて支点と質点をつなぐので，どちらも座標を0とします．

bottom = 1
top = 2
ops.node(bottom, 0.0, 0.0)
ops.node(top, 0.0, 0.0)

続いて，節点の拘束条件を設定します．ここで，2.0で定義したopensees_constantsクラスを用いてコードの可読性を上げています．

# FIX
# node_Tag, x, y, Rz
ops.fix(top, opc.free, opc.fixed, opc.fixed)
ops.fix(bottom, opc.fixed, opc.fixed, opc.fixed)

さらに，$y$, $Rz$方向の変位はtopとbottomで同じになるはずなので，拘束をかけてあげます．これにはops.equalDOF(node1_Tag, node2_Tag, *Args)を使います．Argsにはnode1とnode2で等しい変位になる変位のタグを入れます．ここでは，$y$と$Rz$なので，コードは次のようになります．

# y, Rz方向の変位はtopとbottomで同じなので拘束する
ops.equalDOF(bottom, top, *[opc.Y, opc.ROTZ])

次に，質点に質量を与えます．$M = 200 \mathrm{kg}$なので，次のようにします．

重量を重力加速度で割って入れていましたが，ops.massコマンドは質量で入力するので訂正しました．

mass = 200 # 質量
ops.mass(top, mass, 0.0, 0.0)

続いて，部材の材料の設定をします．ヤング率が$E = 205000\mathrm{N/mm^2}$
なので，次のようになります．

# Define material
elastic_mat_tag = 1
mat_type = 'Elastic'
E = 205000*(10**3)**2
mat_props = [E]
ops.uniaxialMaterial(mat_type, elastic_mat_tag, *mat_props)

続いて，部材の定義をします．今回はzeroLengthエレメントを用います．これはopenseespy上でよく使われる，長さが0の部材のことです．

# zero length elementを指定する
beam_tag = 1
ops.element('zeroLength', beam_tag, bottom, top, '-mat', elastic_mat_tag, '-dir', 1, '-doRayleigh', 1)

最後に減衰の設定を行います．減衰定数が$h = 0.05$なので，次のように設定します．まず，eigen_1 = ops.eigen('-fullGenLapack', 1)により，固有円振動数$\omega_i$の2乗の値を計算します．(固有値を計算しています.)そして，この値のルートをとって固有円振動数を計算します．

ops.rayleigh(alpha_m, beta_k, beta_k_init, beta_k_comm)によってレイリー減衰を適用できます．レイリー減衰では，減衰力$D$は次のように計算します．
$$D = \alpha_M\times M + \beta_K\times K_{curr} + \beta_{K_{init}}\times K_{init} + \beta_{K_{comm}}\times K_{commit}$$
ここで，$\alpha_{M}$は質量比例型減衰の係数，$\beta_{K}$は瞬間剛性比例型減衰の係数，$\beta_{K_{init}}$は初期剛性比例型減衰の係数です．今回は瞬間剛性比例型の減衰としたいので，(線形弾性なので初期剛性比例型でも結果は同じですが)
$$\beta_{K} = \frac{2h}{\omega_1}$$
とします．そのほかの係数はすべて0とします．


# 減衰の設定
h = 0.05
eigen_1 = ops.eigen('-fullGenLapack', 1)
angular_freq = eigen_1[0] ** 0.5 #固有値の0.5乗, omega
alpha_m = 0.0
beta_k = 2 * h / angular_freq
beta_k_comm = 0.0
beta_k_init = 0.0

ops.rayleigh(alpha_m, beta_k, beta_k_init, beta_k_comm)

3. まとめ

今回はここまでです．次回は地震動の波形を読み込み，応答解析を実行したいと思います．
今回のコードをまとめると次の通りです．

# import ilbrary
import openseespy.opensees as ops
import numpy as np

class opensees_constants:
    def __init__(self):
        self.free = 0
        self.fixed = 1

        self.X = 1
        self.Y = 2
        self.ROTZ = 3

opc = opensees_constants()

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

bottom = 1
top = 2
ops.node(bottom, 0.0, 0.0)
ops.node(top, 0.0, 0.0)

# FIX
# node_Tag, x, y, Rz
ops.fix(top, opc.free, opc.fixed, opc.fixed)
ops.fix(bottom, opc.fixed, opc.fixed, opc.fixed)

# y, Rz方向の変位はtopとbottomで同じなので拘束する
ops.equalDOF(bottom, top, *[opc.Y, opc.ROTZ])

G = 9.8 #重力加速度
mass = 200 / G # 質量
ops.mass(top, mass, 0.0, 0.0)


# Define material
elastic_mat_tag = 1
mat_type = 'Elastic'
E = 205000*(10**3)**2
mat_props = [E]
ops.uniaxialMaterial(mat_type, elastic_mat_tag, *mat_props)

# zero length elementを指定する
beam_tag = 1
ops.element('zeroLength', beam_tag, bottom, top, '-mat', elastic_mat_tag, '-dir', 1, '-doRayleigh', 1)

# 減衰の設定
h = 0.05
eigen_1 = ops.eigen('-fullGenLapack', 1)
angular_freq = eigen_1[0] ** 0.5 #固有値の0.5乗, omega
alpha_m = 0.0
beta_k = 2 * h / angular_freq
beta_k_comm = 0.0
beta_k_init = 0.0

ops.rayleigh(alpha_m, beta_k, beta_k_init, beta_k_comm)