構造モデルの時刻歴応答解析においては、その減衰についてある仮定をおいて設定することが一般的です。実際の減衰はどのようなものか、なかなか得体が知れないものではあるのですが、一般的に減衰が全くない構造物はなく何らかの仮定をおいて設定されます。そうした仮定の一つとして、解析時に構造モデル全体に与える減衰モデルとしてレーリー減衰が定義されており、Openseespyにおいても実装されています。そこで、それらを用いて、多質点モデルの時刻歴応答解析において減衰を考慮する検討を行ってみたいと思います。
レーリー減衰についての詳しい説明はここでは省きますが、レーリー減衰は、構造モデルの剛性マトリクスおよび質量マトリクスにそれぞれ比例係数α、βを考慮し、その線形結合として表現される減衰です。
したがって、比例係数α、βの値によって、剛性マトリクスのみに比例する減衰（剛性比例減衰、α=0とした場合）、質量マトリクスのみに比例する減衰（質量比例減衰、β=0とした場合）、剛性、質量マトリクスの両者に比例する減衰（レーリー減衰）の大きく分けて3つの減衰が解析上表現できます。それぞれの減衰は、構造物のそれぞれの周期帯において以下のグラフに示すように表現されます。
図より、剛性比例減衰は固有周期に比例する減衰に、質量比例減衰は固有周期に反比例する減衰となっており、レーリー減衰はその線形結合で表現されていることが分かります。よって例えばその1次周期が長周期となる超高層ビルなどを解析の対象にした場合、仮に剛性比例減衰を用いると、長周期領域で減衰定数を過大に評価してしまうため、主要な振動体での構造物の振動を過小評価してしまう可能性があるといえます。一方で、レーリー減衰は、建物の2つの固有周期帯をターゲットにして減衰定数を指定することができるため、そのような減衰の過大評価を行うことを避けることができます。いずれにせよ、Openseespyではレーリー減衰におけるα、βのパラメータを調節することで、上記の3つのパターンの減衰をいずれも表現することができます。
それでは、実際にOpenseespyにおいて減衰を設定し応答解析を行ってみます。ここでは、以下に示す実際の建物に近い構造特性を有する5層の多質点せん断モデルを解析対象モデルとして設定します。
本構造物の固有値解析により得られた、1次の固有周期（T1=0.70s（ω1=8.97rad））、2次の固有周期（T2=0.28s（ω2=22.56rad））のそれぞれをターゲットとしてレーリー減衰のパラメータを設定します。各次における設定減衰定数としては、鉄筋コンクリート構造物の減衰定数として標準的に用いられるh=0.03（3%）を想定します。すると、レーリー減衰のα、βはそれぞれ以下の式により求めることができます。
#レーリー減衰の係数を算定する h1=0.03 h2=0.03 alpha=2*omega[0]*omega[1]*(h1*omega[1]-h2*omega[0])/(omega[1]**2-omega[0]**2) beta=2*(h2*omega[1]-h1*omega[0])/(omega[1]**2-omega[0]**2) op.rayleigh(alpha,0,beta,0)
omegaは固有値解析により求めた1次、2次の固有円振動数です。また、今回は初期剛性に対して比例係数を設定するため、reyleigh関数の3番目の位置に係数βを設定しています。仮に、2番目の位置に係数βを設定した場合、各ステップの剛性マトリクスに対してβが考慮される瞬間剛性比例型の減衰が設定されます。特に材料の非線形性などを考慮し、時間ステップによって剛性マトリクスの値が変わるような非線形解析などを行う場合は、この機能が有効にすることで応答解析結果が変わります。
上記のコマンドを用いて多質点モデルの地震応答解析に対してレーリー減衰を適用してみます。解析を実行することで得られるモデル頂部の応答は、下記の通りとなります。グラフより、減衰が0（h=0）の時は応答が収束していませんが、レーリー減衰を考慮することで、振動解析途中でモデルの応答が減衰しており、設定した減衰が効いていることが確認できます。
上記の応答において、減衰が正しく適用されているか応答スペクトルの観点から確認してみたいと思います。対象モデルは1次モードの有効質量比が80%を超えており支配的といえます。したがって、モデルの頂部節点では特に1次モードの振動による変位が最も卓越することが想定されます。そこで、まず、対象モデルの1次固有周期（T=0.7s）時の変位応答スペクトルと時刻歴応答解析によるモデルの最大応答相対変位の関係性をみてみます。
今回設定した減衰定数と同じ条件として、減衰定数h=0.03とh=0.10におけるBCJ-L2の変位応答スペクトルをかいてみると以下の通りとなります。
モデルの1次固有周期と変位応答スペクトルの交点をみることで、各減衰定数想定時の最大変位応答SDが算定できます。算定した最大変位応答SDに対して、本対象モデルの1次モードにおける刺激係数β=1.38を乗じることで、h=0.03とした場合の1次モードによる最大変位応答は18.07cm、h=0.10とした場合のモデル頂部の最大変位応答は10.35cmとそれぞれ算定することができます。多質点モデルの時刻歴応答解析におけるモデル頂部節点の変位時刻歴に上記の算定値を点線でプロットしたものを下記に示します。
図より、赤破線で示す算定値は、1次モードの応答のみを考慮したものであるため、高振動数における変位波形のスパイクによる最大変位振幅などは評価できていませんが、それらを除き概ねそれぞれの減衰考慮時の最大時刻歴応答と対応していることが分かります。したがって、本応答解析において、各減衰定数におけるレーリー減衰が概ね正しく設定できていると考えられます。
最後に本検討に用いた多質点せん断モデルの応答解析部分のpythonコードを示します。
import numpy
import math
import openseespy.opensees as op
import pandas as pd
import matplotlib.pyplot as plt
#unit
m_ft=39.3700787#m→inch
kN_pond=0.22480894387096#kN→kips
g = 386.0885827#acceleration of gravity(in/s2)
ratio = 0.393700787 #cm/s2→in/s2
#basic
op.wipe()
op.model('basic', '-ndm', 2) #2D3dof
# 節点座標の定義
m_inch=39.3700787 #単位系m→inchへの変換
op.node(1,0,0*m_inch) #節点1の座標設定(節点番号，x座標，y座標)
op.node(2,0,0*m_inch) #節点2の座標設定
op.node(3,0,0*m_inch) #節点3の座標設定
op.node(4,0,0*m_inch) #節点4の座標設定
op.node(5,0,0*m_inch) #節点5の座標設定
op.node(6,0,0*m_inch) #節点6の座標設定
#boundary 基礎固定
op.fix(1,1,1,1)
op.fix(2,0,1,1)
op.fix(3,0,1,1)
op.fix(4,0,1,1)
op.fix(5,0,1,1)
op.fix(6,0,1,1)
#各層の質量を設定
op.mass(2,12000*kN_pond/g,12000*kN_pond/g,0)
op.mass(3,12000*kN_pond/g,12000*kN_pond/g,0)
op.mass(4,12000*kN_pond/g,12000*kN_pond/g,0)
op.mass(5,12000*kN_pond/g,12000*kN_pond/g,0)
op.mass(6,6000*kN_pond/g,6000*kN_pond/g,0)
#各層のせん断ばねを設定
op.uniaxialMaterial("Elastic", 1, 1200000*kN_pond/(m_ft))
op.uniaxialMaterial("Elastic", 2, 1100000*kN_pond/(m_ft))
op.uniaxialMaterial("Elastic", 3, 900000*kN_pond/(m_ft))
op.uniaxialMaterial("Elastic", 4, 650000*kN_pond/(m_ft))
op.uniaxialMaterial("Elastic", 5, 400000*kN_pond/(m_ft))
op.element("zeroLength", 1, 1, 2,"-mat", 1,"-dir", 1)
op.element("zeroLength", 2, 2, 3,"-mat", 2,"-dir", 1)
op.element("zeroLength", 3, 3, 4,"-mat", 3,"-dir", 1)
op.element("zeroLength", 4, 4, 5,"-mat", 4,"-dir", 1)
op.element("zeroLength", 5, 5, 6,"-mat", 5,"-dir", 1)
dt = 0.01
motion = pd.read_csv('BCJ-L2.csv') #地震波の読み込み
motion = list(motion['acc'])
ratio = 0.393700787 #cm/s2→in/s2
values = [round(float(i)*ratio,4) for i in motion]
op.timeSeries('Path', 1, '-dt', dt, '-values', *values)
#地動加振力の設定
op.pattern('UniformExcitation', 1, 1, '-accel', 1)
#固有値解析の実施
lamda=op.eigen('-fullGenLapack',5)
#固有円振動数を求める。
omega= list(map(lambda x:x**0.5,lamda))
#1次2次周期に対してレーリー減衰を設定する
h1=0.03
h2=0.03
#レーリー減衰の係数を算定する
alpha=2*omega[0]*omega[1]*(h1*omega[1]-h2*omega[0])/(omega[1]**2-omega[0]**2)
beta=2*(h2*omega[1]-h1*omega[0])/(omega[1]**2-omega[0]**2)
op.rayleigh(alpha,0,beta,0)
op.algorithm('Linear')
op.system('BandGen')
op.constraints('Transformation')
op.integrator('Newmark', 0.5, 0.25)
op.analysis('Transient')
op.test('NormUnbalance', 0.001, 10, 0, 2)
analysis_time = (len(values)) * dt #解析時間の算定
outputs = {
"time": [],
"rel_accel": [],
"rel_vel": [],
"rel_disp": [],
}
analysis_time = (len(values)) * dt #解析時間の算定
while op.getTime() < analysis_time:
curr_time = op.getTime()
op.analyze(1, dt)
outputs["time"].append(curr_time)
outputs["rel_disp"].append(op.nodeDisp(6,1)/ratio) #モデル頂部の応答記録
outputs["rel_vel"].append(op.nodeVel(6,1)/ratio)
outputs["rel_accel"].append(op.nodeAccel(6, 1) / ratio)
#応答解析結果の出力
df=pd.DataFrame(outputs)
df.to_csv("result.csv")