ここでは、Openseespyの時刻歴応答解析を用いて、地震波形の応答スペクトルの描画を行ってみたいと思います。この応答スペクトルは前回までに述べてきた構造物の固有周期などの諸元などを考慮することで、その最大地震応答の簡易評価などに用いることができます。
地震波形の応答スペクトルは、様々な固有周期を有する1質点系モデルのある地震波に対する応答値をプロットしたグラフになります。したがって、これまでに述べてきたOpenseespyによる1質点系モデルの時刻歴応答解析をある固有周期帯の間で繰り返し行うことで、描くことができます。
地震を受ける構造物は、地震に対してその特性に応じた様々な減衰特性を有することが知られており、応答スペクトルを様々な減衰定数を対象に描くことができれば、様々な構造物の応答評価に使うことができます。そこで、Openseespyで実施する応答解析に用いる1質点系モデルにおいて今回減衰係数を新たにモデル化します。ある減衰定数hにおける1質点系モデルの挙動をモデル化しようとする場合、以下の式で考慮したい減衰定数hを減衰係数Cに置き換えて1質点系モデルに対してダッシュポットとして付与することで、減衰を考慮することが可能になります。
上記を踏まえ、今回は、h=0.02、0.05、0.1で加速度応答スペクトルをかいてみたいと思います。対象地震波は前回検討と同様にBCJ-L2とします。Openseespyの解析を用いて応答スペクトルを描くためのコードは、以下の通りとなります。
import openseespy.opensees as ops
import numpy as np
import pandas as pd
# ============================================================
# 1. 入力地震波の読み込み
# ============================================================
dt = 0.01 # 時刻刻み
csv_file = "BCJ-L2.csv"
df = pd.read_csv(csv_file)
time = df.iloc[:, 0].values
acc = df.iloc[:, 1].values
# ============================================================
# 2. 加速度応答スペクトル計算
# ============================================================
T_list = np.linspace(0.05, 5.0, 200) #T=0.05～5sまで200刻みで
h = 0.05 # 減衰 5%の場合
Sa = [] # 加速度応答スペクトル（加速度）
#解析対象範囲の周期ごとに一質点系モデルの応答解析を実行
for T in T_list:
ops.wipe()
ops.model('basic', '-ndm', 1, '-ndf', 1)
# 1質点系モデルのパラメータ
m = 1
w = 2 * np.pi / T
k = m * w**2
c = 2 * h * w * m
# 1質点モデルを作成
ops.node(1, 0.0)
ops.node(2, 0.0)
ops.fix(1, 1)
ops.mass(2, m)
ops.uniaxialMaterial("Elastic", 1, k)
ops.uniaxialMaterial("Viscous", 2, c, 1.0)
ops.element("zeroLength", 1, 1, 2,
"-mat", 1, 2,
"-dir", 1, 1)
# 地震動の設定
ops.timeSeries("Path", 1, "-dt", dt, "-values", *acc)
ops.pattern("UniformExcitation", 1, 1, "-accel", 1)
# 解析設定
ops.system("BandGeneral")
ops.numberer("Plain")
ops.constraints("Plain")
ops.integrator("Newmark", 0.5, 1/6) #線形加速度法
ops.algorithm("Linear")
ops.analysis("Transient")
# 時刻歴解析
max_acc = 0.0
for i in range(len(acc)):
ops.analyze(1, dt)
rel_acc = ops.nodeAccel(2, 1)
abs_acc = rel_acc +acc[i] #相対加速度から絶対加速度へ変換
max_acc = max(max_acc, abs(abs_acc))
Sa.append(max_acc) #加速度応答最大値をSaに保存
result = pd.DataFrame({"T":T_list,"Sa":Sa})
result.to_csv("spectrum.csv") #CSV形式にて結果出力 上記のコードでは、対象地震波BCJ-L2波を読み込んだ後に、応答スペクトルを描く範囲の固有周期のリストを作成し、作成したリストに基づいて、一質点系モデルに対して、前ブログと同様に時刻歴応答解析をOpenseespyを用いて繰り返し実行しています。（実際のところ一質点系モデルの応答解析程度なら、一からプログラムを書いた方が早いのかもしれませんが、ここではOpenseespyの演習ということであえてこのように記載しています。）繰り返し計算により得られた各固有周期を有する1質点モデルの加速度応答の最大値を加速度応答スペクトルSaとしてリストに格納し、最後にcsv形式でまとめて出力しています。
ここでの重要なポイントは、以下のようにOpenseespyにより得られた応答値である相対加速度に対して、地動加速度を足すことで絶対加速度に変換していることです。
abs_acc = rel_acc +acc[i] #相対加速度から絶対加速度へ変換
応答スペクトルは絶対加速度を書いていますので、この処理を忘れてしまうと正しい結果を得ることができません。
上記のプログラムにより得られた各減衰定数における加速度応答スペクトルを、応答解析計算における離散化誤差のないNigam-jennings法による正解値と比較してみます。
上記のグラフでは、実線が今回のプログラムで得られた結果で、点線がnigam-jennings法による正解値になります。それぞれの結果を見てみますと、固有周期が1s以下の範囲では、各減衰定数hに対してnigam jennings法による正解値を表現できているのに対し、特に2s以上の長周期領域で正解値よりも大きな値となってしまっており、長周期領域において十分な精度で応答スペクトルが得られていないことが分かります。
地震応答解析における離散化誤差のないnigam-jennings法と差が出ているという点から、試しに応答解析における計算の刻み幅dtを0.01→0.001により細かくしたうえで、上に記載したものと全く同じプログラムでふたたび計算をやり直してみました。すると、今度は下記の通りnigam-jennings法による正解値と概ね整合した値が得られていることがわかります。したがって、上記がスペクトルの解析のずれに影響していたことが考えられます。
参考に今回結果の比較に使ったnigam-gennings法による応答スペクトルの計算プログラムを以下に示しておきます。
import math
import numpy as np
import pandas as pd
# Nigam・Jennings法による応答計算
def nigam_jennings(period, acc0, dt, h):
w = 2.0 * math.pi / period
w_dt = w * dt
dw = math.sqrt(1.0-h**2.0) * w
dw_dt = dw * dt
c = math.cos(dw_dt)
s = math.sin(dw_dt)
e = math.e ** (-1.0*h*w_dt)
c0 = h / math.sqrt(1.0-h**2.0)
c1 = c - c0 * s
c2 = (2.0 * h**2.0 - 1.0) / (w**2.0 * dt)
c3 = s / dw
c4 = (2.0 * h) / (w**3.0 * dt)
c5 = dw * s + h * w * c
c6 = 1.0 / (w**2.0 * dt)
c7 = 1.0 / (w**2.0)
a11 = e * (c0 * s + c)
a12 = e * s / dw
a21 = -1.0 * w * e * s / math.sqrt(1.0-h**2.0)
a22 = e * c1
b11 = e * ((c2+h/w)*c3 + (c4+c7)*c) - c4
b12 = -1.0 * e * (c2*c3 + c4*c) - c7 + c4
b21 = e * ((c2+h/w)*c1 - (c4+c7)*c5) + c6
b22 = -1.0 * e * (c2*c1-c4*c5) - c6
dis = 0.0
vel = 0.0
dis_max = 0.0
vel_max = 0.0
acc_max = 0.0
for n in range(len(acc0)-1):
d = a11*dis + a12*vel + b11*acc0[n] + b12*acc0[n+1]
v = a21*dis + a22*vel + b21*acc0[n] + b22*acc0[n+1]
a = -1.0 * (2.0*h*w*v + w**2.0*d)
dis = d
vel = v
acc = a
if abs(dis_max) < abs(d) : dis_max = abs(d)
if abs(vel_max) < abs(v) : vel_max = abs(v)
if abs(acc_max) < abs(a) : acc_max = abs(a)
acc = -1.0*(2.0*h*w*vel+w**2.0*dis)
return period, dis_max, vel_max, acc_max
y0n1_acc = list(pd.read_csv("bcj-L2.csv")["acc"])
swit=1
h=0.05
max_acc = []
max_vel = []
max_dis = []
temp=0
beta=0.25
dt=0.01
T=[]
stepNum =len(y0n1_acc)
for i in range(501):
temp += 0.01
T.append(np.round(temp,2))
for item in T:
period,m_dis,m_vel,m_acc = nigam_jennings(item, y0n1_acc, dt, h)
#save results
max_acc.append(m_acc)
result = pd.DataFrame({"T":T,"Sa":max_acc})
result.to_csv("spectrum_nigam.csv")