ここではOpenseesPyでP-Δ効果を考慮した解析について検討してみます。
P-Δ効果とは、構造物の幾何学的非線形性によって生じる付加変形および付加応力のことです。例えば、図に示すように軸力を受ける長柱が水平変位を生じると、その変位に応じて付加的な曲げモーメントが発生します。特に大変形領域ではこの影響が無視できなくなるため、P-Δ効果の考慮が必要となります。OpenSeesPyでは、このような効果を幾何剛性としてモデル化し、動的解析に反映させることができるため、その挙動を確認してみます。
OpenSeesPyでは、各要素（element）に対してgeomTransfコマンドでPDelta TransformationまたはCorotational Transformationを指定することで、部材の幾何学的非線形性を考慮することができます。
ここからは、具体的に1質点モデルを用いて検討してみます。対象として、以下に示す簡単な2次元の1質点モデルを想定します。
まず、1質点モデルの梁要素に対してPDeltaTransforamationを設定するために、以下の通り記述します。通常は"Linear"と定義しているgeomTransfを"PDelta"とします。Corotational Transformationでは部材内の回転角も考慮できる？ようなのですが、こちらを設定する場合はここを同じく"Corotational"と設定します。
m_inch=39.3700787 #単位系m→inchへの変換
kN_pond=0.22480894387096 #kN→kips
E_mod=205000000*kN_pond/((m_inch)**2);G_mod=78846153.8461538*kN_pond/((m_inch)**2)
axial_area=0.463385
I=0.504221
#P-deltaの幾何非線形性を設定
ops.geomTransf('PDelta', 1,1,0,0)
#要素はnonlinearBeamColumn要素にて設定
ops.section('Elastic', 1, E_mod, axial_area*(m_inch)**2, I*(m_inch)**4, G_mod, 0.5)
ops.element('nonlinearBeamColumn', 1, 1,2, 10, 1, 1)また、P-Δ効果を考慮するためには、上記のように部材に幾何学的非線形性を設定した上で、部材に作用する軸力をモデルに対して加えた重力解析を事前に実施する必要があります。以下ではモデルy方向に軸力を加えた10ステップの増分解析を行い、最終ステップの結果を保持したまま解析時間をリセットすることで、解析の最終応力状態を次の解析に引き継ぐことを試みています。
m = 40000 #1質点モデルの重量（kN）
#固定荷重としての軸力設定
ops.timeSeries('Constant', 2)
ops.pattern('Plain', 2, 2)
ops.load(2, 0, -m*kN_pond, 0) #部材に加える荷重は重量で入力
# 重力解析の実行
ops.system('BandGeneral')
ops.numberer('RCM')
ops.constraints('Plain')
ops.integrator('LoadControl', 0.1)
ops.algorithm("Newton")
ops.analysis('Static')
ops.analyze(10)
# その状態を保持（荷重状態を保持したうえで、解析時間をリセット）
ops.loadConst('-time', 0.0) 上記により重力解析の状態を保持した上で、水平方向の荷重をかけた解析を行うことで、その解析上で軸力によるP-Δ効果を再現することができるようになります。
それでは、P-Δ効果を検証するため、まず静的な荷重の載荷を前提としたPushover解析を行います。水平荷重を加えたPushover解析により、幾何学的非線形性を考慮する場合としない場合の比較します。Pushover解析では、各ステップで静的な水平荷重を増分し、最終的にモデル頂部に10000 kNの水平荷重を作用させます。
ここで、少しp-Δ効果について考えてみます。p-Δ効果は軸力による見かけの剛性低下を表す現象なので、今回の単純片持ち梁についてその剛性低下分を簡単に評価してみます。
仮に、高さLのモデル頂部に水平変位Δが生じた場合の、軸力Pによる基部の付加曲げモーメントMaddを考えると以下の通り表せます。
Madd=PΔ基部の回転角Θとすると部材内での回転角を無視できればLΘ≒Δとなるので、
Madd≒PLΘ
また付加水平力は以下の通りとなります。
Fadd=Madd/L≒PΘ≒PΔ/L
したがって、片持ち梁においてP-Δ効果の変形増加分の水平方向の幾何剛性は以下の通り表せることになります。
Kadd≒P/L
上記を踏まえて、p-Δ効果を考慮した解析についてみていきます。解析により得られた各ステップの水平変位および基部曲げモーメントの結果は以下のようになります。
図より、通常の解析ケースに比べて、p-Δ効果を考慮したp-delta、corotationalとして部材を設定した場合、水平変形、基部モーメントともやや増加していることが分かります。また、通常の解析ケースの曲げモーメントをみると、高さ50mの点に水平力をそれぞれ1000kN、10000kN加えているので、それぞれ50000kN、500000kNとなっています。さらに、p-deltaとcorotationalの結果を比較すると、大変形領域では結果にやや差がみられますが、ほぼ同様の結果となっています。
また、当初の単純梁の曲げ剛性3EI/L^3と上記により求めた幾何剛性Kaddの差分をP-Δ効果を考慮した見かけの剛性とみなして、基部の曲げモーメントを求めた結果を併せて示してみます。この手計算による評価と、PDeltaおよびCorotationalコマンドによる解析結果は概ね一致しており、Pushover解析においてP-Δ効果が適切に再現されていることが確認できます。
最後にOpenSeesPyのp-delta、corotationalコマンドを適用した際の時刻歴応答解析結果の比較結果を示します。入力地震波はこれまでの検討と同様にBCJ-L2波とし、1次の剛性比例型減衰3%を設定します。
上に示した結果をみると、p-Δ効果を考慮した場合、特にモデル頂部の応答変形について大きめに評価する傾向が確認できます。p-deltaとcorotationalコマンド設定時の応答をそれぞれ比較すると、若干p-deltaコマンドを設定した方が応答を大きめに評価する傾向がありますが、概ね同程度の応答評価となっていることも分かります。
なお、P-Δ効果を考慮した場合に特に変位応答が大きく増加しているのは、本モデルにおいて40000 kNという比較的大きな軸力を与えているためと考えられます。参考として、軸力を4000 kNとした場合の結果を以下に示します。この場合、P-Δ効果の有無による差は大きく縮小することが確認できます。
最後に本検討で実施した1質点モデルを対象としたp-Δ効果を考慮した応答解析の実行コードを示します。
import openseespy.opensees as ops
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
# ============================================================
# 1. 入力地震波の読み込み
# ============================================================
#unit
m_inch=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
#入力地震動の読み込み
dt = 0.01 # 時刻刻み
csv_file = "BCJ-L2.csv"
df = pd.read_csv(csv_file)
time = df.iloc[:, 0].values
acc = df.iloc[:, 1].values
ops.wipe()
ops.model('basic', '-ndm', 2, '-ndf', 3)
# SDOF のパラメータ
m = 40000
ops.node(1, 0.0,0.0)
ops.node(2, 0.0,50*m_inch)
ops.fix(1, 1,1,1)
ops.mass(2, m*kN_pond/g)
#1質点モデルのばねを設定
E_mod=205000000*kN_pond/((m_inch)**2);G_mod=78846153.8461538*kN_pond/((m_inch)**2)
ops.geomTransf('PDelta', 1,1,0,0) #p-Δを設定する場合
axial_area=0.463385
Jxx=100
I=0.504221
ops.section('Elastic', 1, E_mod, axial_area*(m_inch)**2, I*(m_inch)**4, G_mod, 0.5)
ops.element('nonlinearBeamColumn', 1, 1,2, 10, 1, 1)
#固有値解析
omega=ops.eigen('-fullGenLapack',1)
natural_frequency= list(map(lambda x:x**0.5,omega))
natural_frequency.append(0) #2次固有円振動数は0とする
print(natural_frequency)
#1次2次周期に対してレーリー減衰を設定する
h1=0.03
h2=0.03
#レーリー減衰の係数を算定する
alpha=2*natural_frequency[0]*natural_frequency[1]*(h1*natural_frequency[1]-h2*natural_frequency[0])/(natural_frequency[1]**2-natural_frequency[0]**2)
beta=2*(h2*natural_frequency[1]-h1*natural_frequency[0])/(natural_frequency[1]**2-natural_frequency[0]**2)
ops.rayleigh(alpha,0,beta,0)
#軸力を考慮した増分解析
ops.timeSeries('Constant', 2)
ops.pattern('Plain', 2, 2)
ops.load(2, 0, -m*kN_pond, 0)
ops.system('BandGeneral')
ops.numberer('RCM')
ops.constraints('Plain')
ops.integrator('LoadControl', 0.1)
ops.algorithm("Newton")
ops.analysis('Static')
ops.analyze(10)
ops.loadConst('-time', 0.0) #軸力載荷状態を次の解析に引き継いでリセットする
#時刻歴応答解析
# 地震動の入力
ops.timeSeries("Path", 1, "-dt", dt, "-values", *acc*-ratio)
ops.pattern("UniformExcitation", 1, 1, "-accel", 1)
ops.system("BandGeneral")
ops.numberer("Plain")
ops.constraints("Plain")
ops.integrator("Newmark", 0.5, 0.25)
ops.algorithm("Newton")
ops.test("EnergyIncr", 1.0e-8, 50)
ops.analysis("Transient")
# 応答時刻歴の結果格納
outputs = {
"time": [],
"rel_accel": [],
"rel_vel": [],
"rel_disp": [],
}
for i in range(len(acc)):
ops.analyze(1, dt)
curr_time = ops.getTime()
outputs["time"].append(curr_time)
outputs["rel_disp"].append(ops.nodeDisp(2,1)/ratio/100)
outputs["rel_vel"].append(ops.nodeVel(2,1)/ratio/100)
outputs["rel_accel"].append(ops.nodeAccel(2, 1)/ratio/100)
0 件のコメント:
コメントを投稿