1. No category

構造振動特論－第 11 回資料
11 第 11 回
１質点１自由度系の非線形地震応答(3)
最終印刷日時：2014/06/26 12:57:00
１質点１自由度系の非線形地震応答(3)
応答計算例
11.1
次いで，応答計算例を示す．構造物が線形弾性応答している時の最大応答復元力 Qe と構造物の降
伏復元力 Qy の比を R とし，式(11.1)で定義する．
R=
Qe
Qy
(11.1)
ここで，R < 1 の場合には構造物が線形弾性範囲で応答する場合と対応し，R > 1 の場合には構造
物が降伏点を超えて非線形領域で応答する場合と対応する．そこで，入力地震動として El Centro
1940NS を用い，構造物の弾性時の固有周期 Te を 0.3 秒と 1.0 秒，R の値を 1.0～5.0 までの範囲で変
化させて，構造物の降伏耐力が応答に及ぼす影響を検討する．ここで，復元力特性としてはノーマ
ルバイリニアーモデルと剛性低下型バイリニアーモデルを用いることとし，骨格曲線における降伏
後の剛性は弾性時の 0.01 倍とする．また，弾性時における減衰定数は 0.05 とする．
図 11-1 に，構造物の弾性時の固有周期 Te を 0.3 秒とした場合における，線形弾性範囲（R = 1）
での応答変位と応答復元力の時刻歴，ならびに復元力特性としてノーマルバイリニアーモデルを用
いた場合（R = 5）での応答変位と応答復元力の時刻歴を示す．
図 11-1 より以下の傾向が確認できる．
(1)
応答変位に関しては，非線形領域で応答する場合（R = 5）の方が線形弾性範囲で応答
する場合（R = 1）と比べて大きくなる．加えて，R = 5 の場合では 5 秒以降で応答が正方向
に偏る傾向がみられる．
(2)
応答復元力に関しては，非線形領域で応答する場合（R = 5）の方が線形弾性範囲で応
答する場合（R = 1）と比べて小さくなり，応答復元力が降伏復元力 Qy の値でほぼ頭打ちと
なる．
同様に，図 11-2 に，復元力特性として剛性低下型バイリニアーモデルを用いた場合（R = 5）で
の応答変位と応答復元力の時刻歴を示す．図 11-2 より，図 11-1 と同様の傾向が確認できるが，剛
性低下型バイリニアーモデルの場合には，R = 1 の場合ならびに R = 5 として復元力特性としてノー
マルバイリニアーモデルを用いた場合と比べて，応答の周期が長くなっている傾向が明瞭に見られ
る．このように，骨格曲線は同一であっても，応答波形が大きく異なる事が確認できる．
なお，図 11-3 は，R = 5 の場合での復元力－変位関係（履歴応答と呼ぶ事もある）を，ノーマル
バイリニアーモデルと剛性低下型バイリニアーモデルで比較して示したものである．図 11-3 より，
どちらのモデルにおいても復元力－変位関係が仮定した復元力特性に応じてループを描いている事
が確認できる．ここで，履歴ループにより囲まれた面積が構造物により吸収されたひずみエネルギ
ーに対応する．
11-1
構造振動特論－第 11 回資料
図 11-1
１質点１自由度系の非線形地震応答(3)
最終印刷日時：2014/06/26 12:57:00
応答時刻歴波形（ノーマルバイリニアーモデル，Te = 0.3s）
11-2
構造振動特論－第 11 回資料
図 11-2
１質点１自由度系の非線形地震応答(3)
最終印刷日時：2014/06/26 12:57:00
応答時刻歴波形（剛性低下型バイリニアーモデル，Te = 0.3s）
図 11-3
11.2
履歴応答の比較（R = 5.0，Te = 0.3s）
古典的な経験則：エネルギー一定則と変位一定則
構造物の耐震設計において，昔からの最大の関心事の１つは，構造物に生じる最大塑性率 μ をあ
る許容値以内に収めるためには降伏耐力 Qy をどのくらいにすればよいか，ということである．これ
に関し，以下の２つの経験則が知られている。
経験則１（エネルギー一定則）
弾性時固有周期 Te が比較的短周期である構造物において，構造物が弾性挙動すると仮定したとき
の最大せん断力 Qe，最大応答変位を δe とする．これに対し，弾塑性領域で応答する場合の降伏せん
11-3
構造振動特論－第 11 回資料
１質点１自由度系の非線形地震応答(3)
最終印刷日時：2014/06/26 12:57:00
断耐力を Qy，弾塑性応答するときの最大応答変位を δmax（= μδy）とすると，式(11.2)の関係が概ね成
立する（図 11-4）．
1
1
Q y ( 2δ max − δ y ) ≈ Qeδ e
2
2
図 11-4
(11.2)
エネルギー一定則
なお，このエネルギー一定則は，前述した「エネルギーのつり合い」とは全く別のものである点
に留意されたい．
経験則２（変位一定則）
弾性時固有周期 Te が比較的長周期である構造物において，構造物が弾性挙動すると仮定したとき
の最大応答変位 δe と弾塑性応答するときの最大応答変位を δmax は，概ね等しい（式(11.3)）．
δ max ≈ δ e
(11.3)
上の経験則を確認するため，固有周期の違いによる最大応答の傾向の違いを調べる．図 11-5 に，
弾性時の固有周期 Te が 0.3 秒と 1.0 秒のケースについて，縦軸に最大応答加速度を取り横軸に最大
応答変位を取ったものを示す．
図 11-5(a)より，Te が 0.3 秒の場合には，R が大きくなる（降伏耐力が小さくなる）につれて最大
応答加速度が減少する一方で最大応答変位が増加する傾向にある事がわかる．特に，R が大きくな
るにつれて最大応答変位の増加が顕著となる．特に，ノーマルバイリニア―モデルの場合，エネル
ギー一定則と概ね対応していることがわかる．
一方，図 11-5(b)より，Te が 1.0 秒の場合には，R が大きくなる（降伏耐力が小さくなる）につれ
て最大応答加速度が減少するものの，最大応答変位もむしろ減少し，R = 5.0 とした場合でも最大応
答変位が弾性時（R = 1.0）と比べて小さくなっている．この傾向は他の地震動においても一般的に
見られる傾向である．
なお，Te が 0.3 秒と 1.0 秒のケースにおいて，復元力特性を剛性低下型バイリニアーモデルとした
場合では，ノーマルバイリニアーモデルの場合と比べて最大応答変位が大きくなる傾向が見られる．
11-4
構造振動特論－第 11 回資料
１質点１自由度系の非線形地震応答(3)
最終印刷日時：2014/06/26 12:57:00
図 11-5 最大応答の比較
11-5
構造振動特論－第 11 回資料
１質点１自由度系の非線形地震応答(3)
最終印刷日時：2014/06/26 12:57:00
演習問題
11.3
与えられた２つの Excel のマクロ（ノーマルバイリニアーモデル，剛性低下型ノーマルバイリニ
アーモデル）を用いて，以下の１自由度系の非線形時刻歴応答解析を実施し，応答性状を比較せよ．
構造物１：質量 100ton，弾性固有周期 Te = 0.4s
構造物２：質量 100ton，弾性固有周期 Te = 0.8s
構造物３：質量 100ton，弾性固有周期 Te = 2.0s
初期剛性の値は，与えられた質量と弾性固有周期から求める事．
減衰定数 h は 0.05 とする．降伏耐力 Qy と降伏後剛性と初期剛性の比は以下の通りとする．
z
Qy = Qe / R: R = 1, 1.5, 2, 3, 4, 5, 10
z
降伏後剛性と初期剛性の比（降伏後剛性低下率）＝0.001 の１通り
Qe の求め方：初めに，Qy = 10000kN（十分に大きい値）としてマクロを実施し，最大応答せん断
力を求める．この値を Qe とする．
解析結果は，今回資料の図 11-5 と以下の表 11-1 のようにしてまとめること．
表 11-1
解析結果の取りまとめ
最大塑性率
R
ノーマルバイリニ
剛性低下型バイリ
アーモデル
ニアーモデル
1.000
1.000
1.000
*.***
*.***
*.***
エネルギー一定則
変位一定則
1.0
1.000
1.5
*.***
2.0
・・・・
なお，エネルギー一定則に基づく最大塑性率μの推定値は式(11.4)より，変位一定則によるμの推定
値は式(11.5)より与えられる．
μ=
1 2
( R + 1)
2
(11.4)
μ=R
(11.5)
注：グラフと表は，構造物１～３で分けて作成すること．なお，図には，ノーマルバイリニアー
モデルと剛性低下型バイリニアーモデルの解析結果を，同じ図に重ねて示すこと．
11-6
Sheet3 - 1
Sheet3 - 2
Sub SDF_Normal_Bi_Linear()
'　１自由度系モデルの非線形地震応答解析（Normal-Bi-Linear Model）（Excel用マクロ）Ver 1.1
'
'
Programmed by K. Fujii at 2012/05/31
'
Last Modified at 2014/06/20
'
'
本プログラムは，復元力特性がNormal-Bi-Linearモデルである１自由度系モデルの非線形地震応答解析プログラムである．
'
'
'
'
Dim mass As Double '系の質量
Dim stiff1 As Double '系の初期剛性
Dim qy As Double '系の降伏耐力
Dim alpha As Double '系の降伏後剛性と初期剛性の比
Dim h1 As Double '系の減衰定数（瞬間剛性比例型）
Dim
Dim
Dim
Dim
amp As Double '地動の入力倍率
nrecord As Integer '加速度の原データ数
istep As Integer '積分刻み
dt As Double '積分時間刻み（dt0＝dt/istep）
Dim beta As Double 'Newmark-β法での定数
Dim
Dim
Dim
Dim
y0 As Double '前ステップでの相対変位
v0 As Double '前ステップでの相対速度
a0 As Double '前ステップでの相対加速度
ag0 As Double '前ステップでの地動加速度
Dim
Dim
Dim
Dim
y1 As Double '現ステップでの相対変位
v1 As Double '現ステップでの相対速度
a1 As Double '現ステップでの相対加速度
ag1 As Double '現ステップでの地動加速度
Dim
Dim
Dim
Dim
dy As Double '現ステップでの相対変位増分
dv As Double '現ステップでの相対速度増分
da As Double '現ステップでの相対加速度増分
dag As Double '現ステップでの地動加速度増分
Dim vmax As Double '最大相対速度
Dim amax As Double '最大絶対加速度
Dim frmas As Double '最大復元力
Dim duc As Double '最大塑性率
Dim eta As Double '累積塑性変形倍率
Dim
Dim
Dim
Dim
Wkmax
Wdmax
Wsmax
Eimax
As
As
As
As
Double
Double
Double
Double
'最大運動エネルギー
'最大減衰吸収エネルギー
'最大ひずみエネルギー
'最大入力エネルギー
'***** 初期設定 *****
Set ws1 = Worksheets("構造物データと解析結果の要約")
Set ws2 = Worksheets("加速度の原記録")
Set ws3 = Worksheets("運動方程式")
Worksheets("運動方程式").Cells.ClearContents
'***** データ読み込み *****
mass = ws1.Cells(2, 3)
stiff1 = ws1.Cells(4, 3)
qy = ws1.Cells(4, 5)
alpha = ws1.Cells(4, 7)
h1 = ws1.Cells(6, 3)
dt0 = ws1.Cells(8, 3)
amp = ws1.Cells(8, 5)
nrecord = ws1.Cells(8, 7)
istep = ws1.Cells(8, 9)
If istep <= 1 Then
istep = 1
dt = dt0
Else
Dim ktan0 As Double '前ステップでの瞬間剛性
Dim ktan1 As Double '現ステップでの瞬間剛性
dt = dt0 / istep
End If
Dim omega1 As Double '初期固有円振動数
Dim c0 As Double '前ステップでの減衰係数
Dim c1 As Double '現ステップでの減衰係数
Dim
Dim
Dim
Dim
Dim
fd0
fd1
fr0
fr1
ffr
As
As
As
As
As
Double
Double
Double
Double
Double
'前ステップでの瞬間剛性に基づく減衰力
'現ステップでの減衰力
'前ステップでの復元力
'現ステップでの復元力
'現ステップでの増分復元力
Dim dfd As Double '減衰係数の変化による不釣り合い力
Dim dfr As Double '瞬間剛性の変化による不釣り合い力
Dim df0 As Double '現ステップでの不釣り合い力(=次ステップにて解除)
Dim dpn As Double '増分変位算定のための等価増分外力
Dim kn As Double '増分変位算定のための等価剛性
Dim
Dim
Dim
Dim
il
yr
yl
di
As
As
As
As
Integer '復元力ルールNo (1: 弾性，2:降伏，3:除荷)
Double '現ルールにおける変位の上限
Double '現ルールにおける変位の下限
Integer '現ルールにおける変位の向きを表す変数(1: 正方向，2:負方向)
Dim
Dim
Dim
Dim
dWk0
dWd0
dWs0
dEi0
As
As
As
As
Double
Double
Double
Double
'前ステップでの運動エネルギーの増分
'前ステップでの減衰吸収エネルギーの増分
'前ステップでのひずみエネルギーの増分
'前ステップでの入力エネルギーの増分
Dim
Dim
Dim
Dim
dWk1
dWd1
dWs1
dEi1
As
As
As
As
Double
Double
Double
Double
'現ステップでの運動エネルギーの増分
'現ステップでの減衰吸収エネルギーの増分
'現ステップでのひずみエネルギーの増分
'現ステップでの入力エネルギーの増分
Dim
Dim
Dim
Dim
Dim
Wk As Double '現ステップまでの運動エネルギー
Wd As Double '現ステップまでの減衰吸収エネルギー
Ws As Double '現ステップまでのひずみエネルギー
Ei As Double '現ステップまでの入力エネルギー
Err As Double '現ステップでのエネルギーバランスのエラー
Dim ymax As Double '最大相対変位
beta = ws1.Cells(10, 3)
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
1) = "Time"
2) = "地動加速度 G.Acc"
3) = "相対変位 y"
4) = "相対速度 v"
5) = "絶対加速度 a"
6) = "減衰力 Fd"
7) = "復元力 Fr"
8) = "il"
9) = "yr"
10) = "yl"
11) = "di"
12) = "ktan"
13) = "dFd"
14) = "dFr"
15) = "dF"
16) = "Wk"
17) = "Wd"
18) = "Ws"
19) = "Ei"
20) = "Err"
y0 = 0#
v0 = 0#
a0 = 0#
ag0 = 0#
df0 = 0#
omega1 = Sqr(stiff1 / mass)
c0 = 2 * (h1 / omega1) * stiff1
ktan0 = stiff1
il
yr
yl
di
=
=
=
=
1
qy / stiff1
-yr
1
Sheet3 - 3
Sheet3 - 4
dWk0
dWd0
dWs0
dEi0
=
=
=
=
0#
0#
0#
0#
dWk1
dWd1
dWs1
dEi1
=
=
=
=
0#
0#
0#
0#
Wk
Wd
Ws
Ei
0#
0#
0#
0#
GoTo 2
Else
ktan1 = stiff1
fr1 = stiff1 * y1
End If
ElseIf il = 2 Then
=
=
=
=
2
'ルール２：降伏の場合
ktan1 = alpha * stiff1
ymax = 0#
vmax = 0#
amax = 0#
frmax = 0#
If di = 1 Then
fr1 = qy + ktan1 * (y1 - qy / stiff1)
If y1 < yl Then
Wkmax
Wdmax
Wsmax
Eimax
=
=
=
=
0#
0#
0#
0#
il =
di =
yr =
yl =
GoTo
'****** 応答計算 ******
For i = 1 To nrecord
'
3
2
y0
y0 - 2 * (qy / stiff1)
3
Else
地動加速度の読み込み
yl = y1
ag0 = ag1
ag1 = amp * ws2.Cells(i, 1)
End If
Else
ws3.Cells(i + 1, 1) = i * dt0
ws3.Cells(i + 1, 2) = ag1
fr1 = -qy + ktan1 * (y1 + qy / stiff1)
dag = (ag1 - ag0) / istep
agt = ag0 + dag
If y1 > yr Then
il =
di =
yl =
yr =
GoTo
For ii = 1 To istep
'
増分変位の算定
dpn = -mass * dag + mass * (1# / (beta * dt) * v0 + (1# / (2 * beta)) * a0) + c0 * ((1# / (2 * beta)) * v0
+ (1# / (4 * beta) - 1) * a0) + df0
kn = ktan0 + c0 / (2 * beta * dt) + mass / (beta * dt ^ 2)
dy = dpn / kn
'
3
1
y0
y0 + 2 * (qy / stiff1)
3
Else
yr = y1
増分速度，加速度の算定
End If
dv = 1# / (2 * beta * dt) * dy - 1# / (2 * beta) * v0 - (1# / (4 * beta) - 1) * a0
da = 1# / (beta * dt ^ 2) * dy - 1# / (beta * dt) * v0 - 1# / (2 * beta) * a0
End If
Else
'
現ステップでの変位，速度，加速度の算定
3
y1 = y0 + dy
v1 = v0 + dv
a1 = a0 + da
'
'ルール３：除荷の場合
ktan1 = stiff1
fr1 = fr0 + ktan1 * dy
復元力特性による現ステップでの復元力の算定
If y1 > yr Then
ffr = ktan0 * dy
il
di
yr
yl
If il = 1 Then 'ルール１：弾性の場合
=
=
=
=
2
1
1000 * (qy / stiff1)
y1
If y1 > yr Then
GoTo 2
il
di
yr
yl
=
=
=
=
2
1
(qy / stiff1) * 1000
y1
ElseIf y1 < yl Then
il
di
yr
yl
GoTo 2
=
=
=
=
2
2
y1
-1000 * (qy / stiff1)
ElseIf y1 < yl Then
GoTo 2
il
di
yr
yl
=
=
=
=
2
2
y1
-(qy / stiff1) * 1000
End If
End If
Sheet3 - 5
'
Sheet3 - 6
現ステップでの減衰力の算定
Next i
c1 = (2 * h1 / omega1) * ktan1
fd1 = c1 * v1
'
'
duc = ymax / (qy / stiff1)
eta = Ws / (qy ^ 2 / stiff1)
現ステップでの不釣り合い力の算定
dfd = c0 * v1 - fd1
dfr = fr0 + ffr - fr1
ws1.Cells(14,
ws1.Cells(15,
ws1.Cells(16,
ws1.Cells(17,
df0 = dfd + dfr
'
現ステップでのエネルギーの算定
dWk1
dWd1
dWs1
dEi1
=
=
=
=
Wk
Wd
Ws
Ei
Wk
Wd
Ws
Ei
=
=
=
=
0.5
0.5
0.5
0.5
*
*
*
*
(dWk0
(dWd0
(dWs0
(dEi0
+
+
+
+
dWk1)
dWd1)
dWs1)
dEi1)
If Ei <> 0 Then
Err = Abs(Ei - Wk - Wd - Ws) / Ei
Else
Err = 0#
End If
dWk0
dWd0
dWs0
dEi0
'
=
=
=
=
dWk1
dWd1
dWs1
dEi1
データの更新
y0 = y1
v0 = v1
a0 = a1
agt = agt + dag
ktan0 = ktan1
c0 = c1
fr0 = fr1
If
If
If
If
Abs(y1) > ymax Then ymax = Abs(y1)
Abs(v1) > vmax Then vmax = Abs(v1)
Abs(a1 + agt) > amax Then amax = Abs(a1 + agt)
Abs(fr1) > frmax Then frmax = Abs(fr1)
If
If
If
If
Abs(Wk)
Abs(Wd)
Abs(Ws)
Abs(Ei)
>
>
>
>
Wkmax
Wdmax
Wsmax
Eimax
Then
Then
Then
Then
Wkmax
Wdmax
Wsmax
Eimax
Next ii
'
結果の出力
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
3)
3)
3)
3)
=
=
=
=
ymax
vmax
amax
frmax
ws1.Cells(14, 5) = duc
ws1.Cells(15, 5) = eta
mass * a1 * v1 * dt
fd1 * v1 * dt
fr1 * v1 * dt
-mass * agt * v1 * dt
+
+
+
+
解析結果の要約の出力
3) = y1
4) = v1
5) = a1 + ag1
6) = fd1
7) = fr1
8) = il
9) = yr
10) = yl
11) = di
12) = ktan1
13) = dfd
14) = dfr
15) = df0
16) = Wk
17) = Wd
18) = Ws
19) = Ei
20) = Err
=
=
=
=
Abs(Wk)
Abs(Wd)
Abs(Ws)
Abs(Ei)
ws1.Cells(20,
ws1.Cells(21,
ws1.Cells(22,
ws1.Cells(23,
3)
3)
3)
3)
=
=
=
=
Wkmax
Wdmax
Wsmax
Eimax
ws1.Cells(20,
ws1.Cells(21,
ws1.Cells(22,
ws1.Cells(23,
ws1.Cells(25,
5)
5)
5)
5)
5)
=
=
=
=
=
Wk
Wd
Ws
Ei
Err
End Sub
Sheet3 - 1
Sheet3 - 2
Sub SDF_Degrading_Bi_Linear()
'　１自由度系モデルの非線形地震応答解析（Degrading-Bi-Linear Model）（Excel用マクロ）Ver 1.1
'
'
Programmed by K. Fujii at 2014/05/29
'
Last Modified at 2014/06/27
'
'
本プログラムは，復元力特性がDegrading-Bi-Linearモデルである１自由度系モデルの非線形地震応答解析プログラムである．
'
'
'
'
Dim mass As Double '系の質量
Dim stiff1 As Double '系の初期剛性
Dim qy As Double '系の降伏耐力
Dim alpha As Double '系の降伏後剛性と初期剛性の比
Dim h1 As Double '系の減衰定数（瞬間剛性比例型）
Dim
Dim
Dim
Dim
Dim
dt0 As Double '地動加速度データの時間刻み
amp As Double '地動の入力倍率
nrecord As Integer '加速度の原データ数
istep As Integer '積分刻み
dt As Double '積分時間刻み（dt0＝dt/istep）
Dim
Dim
Dim
Dim
Dim
Wk As Double '現ステップまでの運動エネルギー
Wd As Double '現ステップまでの減衰吸収エネルギー
Ws As Double '現ステップまでのひずみエネルギー
Ei As Double '現ステップまでの入力エネルギー
Err As Double '現ステップでのエネルギーバランスのエラー
Dim
Dim
Dim
Dim
ymax As Double '最大相対変位
vmax As Double '最大相対速度
amax As Double '最大絶対加速度
frmas As Double '最大復元力
Dim duc As Double '最大塑性率
Dim eta As Double '累積塑性変形倍率
Dim
Dim
Dim
Dim
Wkmax
Wdmax
Wsmax
Eimax
As
As
As
As
Double
Double
Double
Double
'最大運動エネルギー
'最大減衰吸収エネルギー
'最大ひずみエネルギー
'最大入力エネルギー
Dim ypeak1 As Double '正側最大相対変位（ループ3-4開始点）
Dim ypeak2 As Double '負側最大相対変位（ループ3-4開始点）
Dim frpeak1 As Double '大ループの正側ピーク復元力（ループ3-4開始点）
Dim frpeak2 As Double '大ループの負側ピーク復元力（ループ3-4開始点）
Dim beta As Double 'Newmark-β法での定数
Dim
Dim
Dim
Dim
y0 As Double '前ステップでの相対変位
v0 As Double '前ステップでの相対速度
a0 As Double '前ステップでの相対加速度
ag0 As Double '前ステップでの地動加速度
Dim
Dim
Dim
Dim
y1 As Double '現ステップでの相対変位
v1 As Double '現ステップでの相対速度
a1 As Double '現ステップでの相対加速度
ag1 As Double '現ステップでの地動加速度
Dim
Dim
Dim
Dim
Dim
dy As Double '現ステップでの相対変位増分
dv As Double '現ステップでの相対速度増分
da As Double '現ステップでの相対加速度増分
dag As Double '現ステップでの地動加速度増分
agt As Double '分割ステップでの地動加速度
Dim ktan0 As Double '前ステップでの瞬間剛性
Dim ktan1 As Double '現ステップでの瞬間剛性
Dim
Dim
Dim
Dim
ktan04 As Double
ktan06 As Double
ktan08 As Double
ktan010 As Double
Dim omega1 As Double '初期固有円振動数
Dim c0 As Double '前ステップでの減衰係数
Dim c1 As Double '現ステップでの減衰係数
Dim
Dim
Dim
Dim
Dim
fd0
fd1
fr0
fr1
ffr
As
As
As
As
As
Double
Double
Double
Double
Double
'前ステップでの瞬間剛性に基づく減衰力
'現ステップでの減衰力
'前ステップでの復元力
'現ステップでの復元力
'現ステップでの増分復元力
Dim dfd As Double '減衰係数の変化による不釣り合い力
Dim dfr As Double '瞬間剛性の変化による不釣り合い力
Dim df0 As Double '現ステップでの不釣り合い力(=次ステップにて解除)
Dim dpn As Double '増分変位算定のための等価増分外力
Dim kn As Double '増分変位算定のための等価剛性
Dim
Dim
Dim
Dim
il
yr
yl
di
As
As
As
As
Integer '復元力ルールNo (1～11)
Double '現ルールにおける変位の上限
Double '現ルールにおける変位の下限
Integer '現ルールにおける変位の向きを表す変数(1: 正方向，2:負方向)
Dim
Dim
Dim
Dim
dWk0
dWd0
dWs0
dEi0
As
As
As
As
Double
Double
Double
Double
'前ステップでの運動エネルギーの増分
'前ステップでの減衰吸収エネルギーの増分
'前ステップでのひずみエネルギーの増分
'前ステップでの入力エネルギーの増分
Dim
Dim
Dim
Dim
dWk1
dWd1
dWs1
dEi1
As
As
As
As
Double
Double
Double
Double
'現ステップでの運動エネルギーの増分
'現ステップでの減衰吸収エネルギーの増分
'現ステップでのひずみエネルギーの増分
'現ステップでの入力エネルギーの増分
Dim ypeak01 As Double '小ループの正側ピーク変位（ループ5-6開始点）
Dim ypeak02 As Double '小ループの負側ピーク変位（ループ5-6開始点）
Dim frpeak01 As Double '小ループの正側ピーク復元力（ループ5-6開始点）
Dim frpeak02 As Double '小ループの負側ピーク復元力（ループ5-6開始点点）
Dim ypeak03 As Double '小ループの正側ピーク変位（ループ7-8開始点，ループ9-10開始点，ループ11-10開始点）
Dim ypeak04 As Double '小ループの負側ピーク変位（ループ7-8開始点，ループ9-10開始点，ループ11-10開始点）
Dim frpeak03 As Double '小ループの正側ピーク復元力（ループ7-8開始点，ループ9-10開始点，ループ11-10開始点）
Dim frpeak04 As Double '小ループの負側ピーク復元力（ループ7-8開始点，ループ9-10開始，ループ11-10開始点点）
Dim y01 As Double '大ループの正側復元力0点変位(ループ3-4中間点)
Dim y02 As Double '大ループの負側復元力0点変位(ループ3-4中間点)
Dim
Dim
Dim
Dim
y001
y002
y003
y004
As
As
As
As
Double
Double
Double
Double
'小ループの正側復元力0点変位(ループ5-6中間点)
'小ループの負側復元力0点変位(ループ5-6中間点)
'小ループの正側復元力0点変位（２周目，ループ7-8中間点，ループ9-10中間点，ループ11-10中間点)
'小ループの負側復元力0点変位（２周目，ループ7-8中間点，ループ9-10中間点，ループ11-10中間点)
'***** 初期設定 *****
Set ws1 = Worksheets("構造物データと解析結果の要約")
Set ws2 = Worksheets("加速度の原記録")
Set ws3 = Worksheets("運動方程式")
Worksheets("運動方程式").Cells.ClearContents
'***** データ読み込み *****
mass = ws1.Cells(2, 3)
stiff1 = ws1.Cells(4, 3)
qy = ws1.Cells(4, 5)
alpha = ws1.Cells(4, 7)
h1 = ws1.Cells(6, 3)
dt0 = ws1.Cells(8, 3)
amp = ws1.Cells(8, 5)
nrecord = ws1.Cells(8, 7)
istep = ws1.Cells(8, 9)
If istep <= 1 Then
istep = 1
dt = dt0
Else
dt = dt0 / istep
End If
beta = ws1.Cells(10, 3)
ws3.Cells(1, 1) = "Time"
ws3.Cells(1, 2) = "地動加速度 G.Acc"
Sheet3 - 3
Sheet3 - 4
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
3) = "相対変位 y"
4) = "相対速度 v"
5) = "絶対加速度 a"
6) = "減衰力 Fd"
7) = "復元力 Fr"
8) = "il"
9) = "yr"
10) = "yl"
11) = "di"
12) = "ktan"
13) = "dFd"
14) = "dFr"
15) = "dF"
16) = "Wk"
17) = "Wd"
18) = "Ws"
19) = "Ei"
20) = "Err"
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
21)
22)
23)
24)
=
=
=
=
"Ypeak1"
"Frpeak1"
"Ypeak2"
"Frpeak2"
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
25)
26)
27)
28)
=
=
=
=
"Ypeak01"
"Frpeak01"
"Ypeak02"
"Frpeak02"
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
29)
30)
31)
32)
=
=
=
=
"Ypeak03"
"Frpeak03"
"Ypeak04"
"Frpeak04"
Eimax = 0#
'****** 応答計算 ******
For i = 1 To nrecord
'
地動加速度の読み込み
ag0 = ag1
ag1 = amp * ws2.Cells(i, 1)
ws3.Cells(i + 1, 1) = i * dt0
ws3.Cells(i + 1, 2) = ag1
dag = (ag1 - ag0) / istep
agt = ag0 + dag
For ii = 1 To istep
'
増分変位の算定
dpn = -mass * dag + mass * (1# / (beta * dt) * v0 + (1# / (2 * beta)) * a0) + c0 * ((1# / (2 * beta)) * v0
+ (1# / (4 * beta) - 1) * a0) + df0
kn = ktan0 + c0 / (2 * beta * dt) + mass / (beta * dt ^ 2)
dy = dpn / kn
'
dv = 1# / (2 * beta * dt) * dy - 1# / (2 * beta) * v0 - (1# / (4 * beta) - 1) * a0
da = 1# / (beta * dt ^ 2) * dy - 1# / (beta * dt) * v0 - 1# / (2 * beta) * a0
'
ws3.Cells(1, 33) = "Y01"
ws3.Cells(1, 34) = "Y02"
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
ws3.Cells(1,
35)
36)
37)
38)
=
=
=
=
"Y001"
"Y002"
"Y003"
"Y004"
y0 = 0#
v0 = 0#
a0 = 0#
ag0 = 0#
ag1 = 0#
増分速度，加速度の算定
現ステップでの変位，速度，加速度の算定
y1 = y0 + dy
v1 = v0 + dv
a1 = a0 + da
'
復元力特性による現ステップでの復元力の算定
ffr = ktan0 * dy
If il = 1 Then 'ルール１：弾性の場合
If y1 > yr Then
il
di
yr
yl
df0 = 0#
omega1 = Sqr(stiff1 / mass)
c0 = 2 * (h1 / omega1) * stiff1
ktan0 = stiff1
il
yr
yl
di
=
=
=
=
=
=
=
=
0#
0#
0#
0#
dWk1
dWd1
dWs1
dEi1
=
=
=
=
0#
0#
0#
0#
Wk
Wd
Ws
Ei
0#
0#
0#
0#
=
=
=
=
2
1
(qy / stiff1) * 1000
y1
ktan1 = alpha * stiff1
ypeak2 = -qy / stiff1
frpeak2 = -qy
1
qy / stiff1
-yr
1
dWk0
dWd0
dWs0
dEi0
=
=
=
=
ypeak02 = -qy / stiff1
frpeak02 = -qy
GoTo 2
ElseIf y1 < yl Then
il
di
yr
yl
=
=
=
=
2
2
y1
-(qy / stiff1) * 1000
ypeak1 = qy / stiff1
frpeak1 = qy
ypeak01 = qy / stiff1
frpeak01 = qy
ymax = 0#
vmax = 0#
amax = 0#
frmax = 0#
ktan1 = alpha * stiff1
GoTo 2
Else
Wkmax = 0#
Wdmax = 0#
Wsmax = 0#
ktan1 = stiff1
fr1 = stiff1 * y1
Sheet3 - 5
Sheet3 - 6
ypeak04 = y1
frpeak04 = fr1
End If
ElseIf il = 2 Then
2
End If
'ルール２：降伏の場合
End If
If di = 1 Then
fr1 = qy + ktan1 * (y1 - qy / stiff1)
ElseIf il = 3 Then
3
If y1 < yl Then
'ルール３：除荷の場合（大ループ）
fr1 = fr0 + ktan1 * dy
il = 3
di = 2
If y1 > yr Then
If di = 2 Then
ypeak1 = y0
frpeak1 = fr0
il = 2
di = 1
ypeak01 = y0
frpeak01 = fr0
yr = 1000 * (qy / stiff1)
yl = y1
ktan1 = alpha * stiff1
ypeak03 = y0
frpeak03 = fr0
GoTo 2
ktan1 = stiff1 * Sqr(Abs((qy / stiff1) / ypeak1))
y01 = y0 - (fr0 / ktan1)
Else
yr = y0
yl = y01
il
di
yr
yl
GoTo 3
Else
=
=
=
=
4
1
ypeak1
y02
ktan1 = frpeak1 / (ypeak1 - y02)
fr1 = ktan1 * (y1 - y02)
yl = y1
ypeak1 = y1
frpeak1 = fr1
GoTo 4
End If
ypeak01 = y1
frpeak01 = fr1
ElseIf y1 < yl Then
ypeak03 = y1
frpeak03 = fr1
If di = 1 Then
il = 2
di = 2
yr = y1
yl = -1000 * (qy / stiff1)
ktan1 = alpha * stiff1
End If
Else
fr1 = -qy + ktan1 * (y1 + qy / stiff1)
GoTo 2
If y1 > yr Then
Else
il = 3
di = 1
ypeak2 = y0
frpeak2 = fr0
il
di
yr
yl
=
=
=
=
4
2
y01
ypeak2
ypeak02 = y0
frpeak02 = fr0
ktan1 = frpeak2 / (ypeak2 - y01)
fr1 = ktan1 * (y1 - y01)
ypeak04 = y0
frpeak04 = fr0
GoTo 4
End If
ktan1 = stiff1 * Sqr(Abs((qy / stiff1) / ypeak2))
y02 = y0 - (fr0 / ktan1)
End If
yl = y0
yr = y02
ElseIf il = 4 Then
4
'ルール４：逆方向載荷の場合（大ループ）
GoTo 3
If di = 1 Then
Else
fr1 = ktan1 * (y1 - y02)
yr = y1
ypeak2 = y1
frpeak2 = fr1
If y1 > yl Then
yl = y1
ypeak02 = y1
frpeak02 = fr1
If y1 > yr Then
Sheet3 - 7
Sheet3 - 8
il = 2
di = 1
yr = 1000 * (qy / stiff1)
yl = y1
ktan1 = alpha * stiff1
yl = y0
GoTo 5
End If
GoTo 2
End If
Else
ElseIf il = 5 Then
ypeak01 = y1
frpeak01 = fr1
5
'ルール５：除荷の場合（小ループ）
fr1 = fr0 + ktan1 * dy
ypeak03 = y1
frpeak03 = fr1
If y1 > yr Then
End If
If di = 2 Then
Else
il = 5
di = 2
il
di
yr
yl
=
=
=
=
4
1
ypeak1
y02
ypeak01 = y0
frpeak01 = fr0
ktan1 = frpeak1 / (ypeak1 - y02)
fr1 = ktan1 * (y1 - y02)
ypeak03 = y0
frpeak03 = fr0
GoTo 4
Else
ktan1 = stiff1 * Sqr(Abs((qy / stiff1) / ypeak1))
y001 = y0 - (fr0 / ktan1)
il
di
yr
yl
yr = y0
yl = y001
GoTo 5
=
=
=
=
6
1
ypeak1
y002
ktan1 = frpeak1 / (ypeak1 - y002)
fr1 = ktan1 * (y1 - y002)
End If
GoTo 6
Else
End If
fr1 = ktan1 * (y1 - y01)
ElseIf y1 < yl Then
If y1 < yr Then
If di = 1 Then
yr = y1
il
di
yr
yl
If y1 < yl Then
il = 2
di = 2
yl = -1000 * (qy / stiff1)
yr = y1
ktan1 = alpha * stiff1
=
=
=
=
4
2
y01
ypeak2
ktan1 = frpeak2 / (ypeak2 - y01)
fr1 = ktan1 * (y1 - y01)
GoTo 4
GoTo 2
Else
Else
il
di
yr
yl
ypeak02 = y1
frpeak02 = fr1
ypeak04 = y1
frpeak04 = fr1
=
=
=
=
6
2
y001
ypeak2
ktan1 = frpeak2 / (ypeak2 - y001)
fr1 = ktan1 * (y1 - y001)
End If
GoTo 6
Else
End If
il = 5
di = 1
End If
ypeak02 = y0
frpeak02 = fr0
ElseIf il = 6 Then
6
ypeak04 = y0
frpeak04 = fr0
ktan1 = stiff1 * Sqr(Abs((qy / stiff1) / ypeak2))
y002 = y0 - (fr0 / ktan1)
'ルール６：逆方向載荷の場合（小ループ）
If di = 1 Then
fr1 = ktan1 * (y1 - y002)
If y1 > yl Then
yr = y002
Sheet3 - 9
Sheet3 - 10
yl = y1
If y1 > yr Then
il = 2
di = 1
yr = 1000 * (qy / stiff1)
yl = y1
ktan1 = alpha * stiff1
ElseIf il = 7 Then
7
'ルール７：除荷の場合（小ループ2）
fr1 = fr0 + ktan1 * dy
If y1 > yr Then
If di = 2 Then
GoTo 2
il
di
yr
yl
Else
ypeak03 = y1
frpeak03 = fr1
=
=
=
=
6
1
ypeak1
y002
ktan1 = frpeak1 / (ypeak1 - y002)
fr1 = ktan1 * (y1 - y002)
End If
GoTo 6
Else
Else
il = 7
di = 2
ktan06 = frpeak1 / (ypeak1 - y004)
ktan08 = frpeak01 / (ypeak01 - y004)
ypeak03 = y0
frpeak03 = fr0
If ktan08 > ktan1 Or ktan08 < ktan06 Then
ktan1 = stiff1 * Sqr(Abs((qy / stiff1) / ypeak1))
y003 = y0 - (fr0 / ktan1)
y002 = y004
il
di
yr
yl
yr = y0
yl = y003
=
=
=
=
6
1
ypeak1
y002
GoTo 7
ktan1 = frpeak1 / (ypeak1 - y002)
fr1 = ktan1 * (y1 - y002)
End If
Else
GoTo 6
fr1 = ktan1 * (y1 - y001)
Else
If y1 < yr Then
il
di
yr
yl
yr = y1
=
=
=
=
8
1
ypeak01
y004
If y1 < yl Then
ktan1 = frpeak01 / (ypeak01 - y004)
fr1 = ktan1 * (y1 - y004)
il = 2
di = 2
yr = 1000 * (qy / stiff1)
yr = y1
ktan1 = alpha * stiff1
GoTo 2
Else
GoTo 8
End If
End If
ElseIf y1 < yl Then
ypeak04 = y1
frpeak04 = fr1
If di = 1 Then
il
di
yr
yl
End If
Else
il = 7
di = 1
=
=
=
=
6
2
y001
ypeak2
ktan1 = frpeak2 / (ypeak2 - y001)
fr1 = ktan1 * (y1 - y001)
ypeak04 = y0
frpeak04 = fr0
GoTo 6
Else
ktan1 = stiff1 * Sqr(Abs((qy / stiff1) / ypeak2))
y004 = y0 - (fr0 / ktan1)
yr = y004
yl = y0
GoTo 7
End If
End If
ktan06 = frpeak2 / (ypeak2 - y003)
ktan08 = frpeak02 / (ypeak02 - y003)
If ktan08 > ktan1 Or ktan08 < ktan06 Then
y001 = y003
il
di
yr
yl
=
=
=
=
6
2
y001
ypeak2
Sheet3 - 11
Sheet3 - 12
If y1 < yl Then
ktan1 = frpeak2 / (ypeak2 - y001)
fr1 = ktan1 * (y1 - y001)
il
di
yr
yl
GoTo 6
=
=
=
=
4
2
y01
ypeak2
Else
il
di
yr
yl
=
=
=
=
ktan1 = frpeak2 / (ypeak2 - y01)
fr1 = ktan1 * (y1 - y01)
8
2
y003
ypeak02
GoTo 4
Else
ktan1 = frpeak02 / (ypeak02 - y003)
fr1 = ktan1 * (y1 - y003)
ypeak04 = y1
frpeak04 = fr1
GoTo 8
End If
End If
Else
End If
il = 9
di = 1
End If
ElseIf il = 8 Then
8
ypeak04 = y0
frpeak04 = fr0
'ルール８：逆方向載荷の場合（小ループ2）
ktan1 = stiff1 * Sqr(Abs((qy / stiff1) / ypeak2))
If di = 1 Then
y004 = y0 - (fr0 / ktan1)
fr1 = ktan1 * (y1 - y004)
yr = y004
yl = y0
If y1 > yl Then
yl = y1
GoTo 9
If y1 > yr Then
il
di
yr
yl
=
=
=
=
End If
4
1
ypeak1
y02
End If
ElseIf il = 9 Then
9
'ルール９：ルール８からの除荷の場合（小ループ2）
ktan1 = frpeak1 / (ypeak1 - y02)
fr1 = ktan1 * (y1 - y02)
fr1 = fr0 + ktan1 * dy
GoTo 4
If y1 > yr Then
Else
If di = 2 Then
ypeak03 = y1
frpeak03 = fr1
il
di
yr
yl
End If
Else
=
=
=
=
8
1
ypeak01
y004
ktan1 = frpeak01 / (ypeak01 - y004)
fr1 = ktan1 * (y1 - y004)
il = 9
di = 2
ypeak03 = y0
frpeak03 = fr0
GoTo 8
Else
ktan1 = stiff1 * Sqr(Abs((qy / stiff1) / ypeak1))
ktan06 = frpeak1 / (ypeak1 - y004)
ktan010 = frpeak01 / (ypeak01 - y004)
y003 = y0 - (fr0 / ktan1)
If ktan010 > ktan1 Or ktan010 < ktan06 Then
yr = y0
yl = y003
GoTo 9
End If
Else
y002 = y004
il
di
yr
yl
=
=
=
=
6
1
ypeak1
y002
ktan1 = frpeak1 / (ypeak1 - y002)
fr1 = ktan1 * (y1 - y002)
fr1 = ktan1 * (y1 - y003)
GoTo 6
If y1 < yr Then
Else
yr = y1
il = 10
Sheet3 - 13
Sheet3 - 14
di = 1
yr = ypeak01
yl = y004
If ypeak1 > ypeak03 Then
il = 12
di = 1
ktan1 = frpeak01 / (ypeak01 - y004)
fr1 = ktan1 * (y1 - y004)
yr = ypeak1
yl = ypeak03
GoTo 10
ktan1 = (frpeak1 - frpeak03) / (ypeak1 - ypeak03)
fr1 = frpeak03 + ktan1 * (y1 - ypeak03)
End If
End If
GoTo 12
ElseIf y1 < yl Then
Else
If di = 1 Then
il
di
yr
yl
=
=
=
=
il = 2
di = 1
yr = 1000 * (qy / stiff1)
yl = y1
ktan1 = alpha * stiff1
8
2
y003
ypeak02
GoTo 2
ktan1 = frpeak2 / (ypeak02 - y003)
fr1 = ktan1 * (y1 - y001)
End If
GoTo 8
Else
Else
ypeak03 = y1
frpeak03 = fr1
ktan06 = frpeak2 / (ypeak2 - y003)
ktan010 = frpeak02 / (ypeak02 - y003)
End If
Else
If ktan010 > ktan1 Or ktan010 < ktan06 Then
il = 11
di = 2
y001 = y003
il
di
yr
yl
=
=
=
=
6
2
y001
ypeak2
ypeak03 = y0
frpeak03 = fr0
ktan1 = stiff1 * Sqr(Abs((qy / stiff1) / ypeak1))
ktan1 = frpeak2 / (ypeak2 - y001)
fr1 = ktan1 * (y1 - y001)
y003 = y0 - (fr0 / ktan1)
yr = y0
yl = y003
GoTo 6
Else
GoTo 11
il
di
yr
yl
=
=
=
=
10
2
y003
ypeak02
End If
Else
fr1 = ktan1 * (y1 - y003)
ktan1 = frpeak02 / (ypeak02 - y003)
fr1 = ktan1 * (y1 - y003)
GoTo 10
End If
End If
If y1 < yr Then
yr = y1
If y1 < yl Then
ypeak04 = yl
frpeak04 = ktan1 * (ypeak04 - y003)
End If
If ypeak2 < ypeak04 Then
ElseIf il = 10 Then
10
'ルール１０：逆方向載荷の場合（小ループ2）
If di = 1 Then
il = 12
di = 2
yl = ypeak2
yr = ypeak04
fr1 = ktan1 * (y1 - y004)
If y1 > yl Then
yl = y1
If y1 > yr Then
ypeak03 = yr
frpeak03 = ktan1 * (ypeak03 - y004)
ktan1 = (frpeak2 - frpeak04) / (ypeak2 - ypeak04)
fr1 = frpeak04 + ktan1 * (y1 - ypeak04)
GoTo 12
Else
il = 2
di = 2
yl = -1000 * (qy / stiff1)
Sheet3 - 15
Sheet3 - 16
yr = y1
ktan1 = alpha * stiff1
fr1 = ktan1 * (y1 - y004)
GoTo 10
GoTo 2
End If
End If
End If
Else
ElseIf y1 < yl Then
ypeak04 = y1
frpeak04 = fr1
If di = 1 Then
End If
il
di
yr
yl
Else
il = 11
di = 1
=
=
=
=
10
2
y003
ypeak02
ktan1 = frpeak2 / (ypeak02 - y003)
fr1 = ktan1 * (y1 - y003)
ypeak04 = y0
frpeak04 = fr0
GoTo 10
ktan1 = stiff1 * Sqr(Abs((qy / stiff1) / ypeak2))
Else
y004 = y0 - (fr0 / ktan1)
ktan06 = frpeak2 / (ypeak2 - y003)
ktan010 = frpeak02 / (ypeak02 - y003)
yr = y004
yl = y0
If ktan010 > ktan1 Or ktan010 < ktan06 Then
GoTo 11
y001 = y003
End If
il
di
yr
yl
End If
ElseIf il = 11 Then
11
'ルール１１：ルール１０からの除荷の場合（小ループ2）
=
=
=
=
6
2
y001
ypeak2
ktan1 = frpeak2 / (ypeak2 - y001)
fr1 = ktan1 * (y1 - y001)
fr1 = fr0 + ktan1 * dy
GoTo 6
If y1 > yr Then
Else
If di = 2 Then
il
di
yr
yl
=
=
=
=
il
di
yr
yl
10
1
ypeak01
y004
=
=
=
=
10
2
y003
ypeak02
ktan1 = frpeak02 / (ypeak02 - y003)
fr1 = ktan1 * (y1 - y003)
ktan1 = frpeak01 / (ypeak01 - y004)
fr1 = ktan1 * (y1 - y004)
GoTo 10
GoTo 10
End If
Else
End If
ktan06 = frpeak1 / (ypeak1 - y004)
ktan010 = frpeak01 / (ypeak01 - y004)
End If
If ktan010 > ktan1 Or ktan010 < ktan06 Then
y002 = y004
il
di
yr
yl
=
=
=
=
6
1
ypeak1
y002
ElseIf il = 12 Then
12
'ルール１２：ルール１０から骨格曲線への載荷（小ループ2）
If di = 1 Then
fr1 = frpeak03 + ktan1 * (y1 - ypeak03)
If y1 > yl Then
ktan1 = frpeak1 / (ypeak1 - y002)
fr1 = ktan1 * (y1 - y002)
yl = y1
GoTo 6
If y1 > yr Then
Else
il
di
yr
yl
=
=
=
=
10
1
ypeak01
y004
il = 2
di = 1
yr = 1000 * (qy / stiff1)
yl = y1
ktan1 = alpha * stiff1
GoTo 2
ktan1 = frpeak01 / (ypeak01 - y004)
Sheet3 - 17
Sheet3 - 18
Else
df0 = dfd + dfr
ypeak03 = y1
frpeak03 = fr1
'
End If
Else
il = 7
di = 2
ypeak03 = y0
frpeak03 = fr0
現ステップでのエネルギーの算定
dWk1
dWd1
dWs1
dEi1
=
=
=
=
Wk
Wd
Ws
Ei
Wk
Wd
Ws
Ei
=
=
=
=
mass * a1 * v1 * dt
fd1 * v1 * dt
fr1 * v1 * dt
-mass * agt * v1 * dt
+
+
+
+
0.5
0.5
0.5
0.5
*
*
*
*
(dWk0
(dWd0
(dWs0
(dEi0
+
+
+
+
dWk1)
dWd1)
dWs1)
dEi1)
If Ei <> 0 Then
ktan1 = stiff1 * Sqr(Abs((qy / stiff1) / ypeak1))
y003 = y0 - (fr0 / ktan1)
Err = Abs(Ei - Wk - Wd - Ws) / Ei
yr = y0
yl = y003
Else
Err = 0#
GoTo 7
End If
End If
dWk0
dWd0
dWs0
dEi0
Else
fr1 = frpeak04 + ktan1 * (y1 - ypeak04)
If y1 < yr Then
'
yr = y1
=
=
=
=
dWk1
dWd1
dWs1
dEi1
データの更新
y0 = y1
v0 = v1
a0 = a1
agt = agt + dag
ktan0 = ktan1
c0 = c1
fr0 = fr1
If y1 < yl Then
il = 2
di = 2
yr = 1000 * (qy / stiff1)
yr = y1
ktan1 = alpha * stiff1
GoTo 2
If
If
If
If
Abs(y1) > ymax Then ymax = Abs(y1)
Abs(v1) > vmax Then vmax = Abs(v1)
Abs(a1 + agt) > amax Then amax = Abs(a1 + agt)
Abs(fr1) > frmax Then frmax = Abs(fr1)
If
If
If
If
Abs(Wk)
Abs(Wd)
Abs(Ws)
Abs(Ei)
Else
ypeak04 = y1
frpeak04 = fr1
>
>
>
>
Wkmax
Wdmax
Wsmax
Eimax
Then
Then
Then
Then
Wkmax
Wdmax
Wsmax
Eimax
End If
Next ii
Else
'
il = 7
di = 1
ypeak04 = y0
frpeak04 = fr0
ktan1 = stiff1 * Sqr(Abs((qy / stiff1) / ypeak2))
y004 = y0 - (fr0 / ktan1)
yr = y004
yl = y0
GoTo 7
End If
End If
End If
'
現ステップでの減衰力の算定
c1 = (2 * h1 / omega1) * ktan1
fd1 = c1 * v1
'
現ステップでの不釣り合い力の算定
dfd = c0 * v1 - fd1
dfr = fr0 + ffr - fr1
結果の出力
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
3) = y1
4) = v1
5) = a1 + agt
6) = fd1
7) = fr1
8) = il
9) = yr
10) = yl
11) = di
12) = ktan1
13) = dfd
14) = dfr
15) = df0
16) = Wk
17) = Wd
18) = Ws
19) = Ei
20) = Err
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
+
+
+
+
1,
1,
1,
1,
21)
22)
23)
24)
=
=
=
=
ypeak1
frpeak1
ypeak2
frpeak2
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
+
+
+
+
1,
1,
1,
1,
25)
26)
27)
28)
=
=
=
=
ypeak01
frpeak01
ypeak02
frpeak02
=
=
=
=
Abs(Wk)
Abs(Wd)
Abs(Ws)
Abs(Ei)
Sheet3 - 19
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
+
+
+
+
1,
1,
1,
1,
29)
30)
31)
32)
=
=
=
=
ypeak03
frpeak03
ypeak04
frpeak04
ws3.Cells(i + 1, 33) = y01
ws3.Cells(i + 1, 34) = y02
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
ws3.Cells(i
+
+
+
+
1,
1,
1,
1,
35)
36)
37)
38)
=
=
=
=
y001
y002
y003
y004
Next i
'
解析結果の要約の出力
duc = ymax / (qy / stiff1)
eta = Ws / (qy ^ 2 / stiff1)
ws1.Cells(14,
ws1.Cells(15,
ws1.Cells(16,
ws1.Cells(17,
3)
3)
3)
3)
=
=
=
=
ymax
vmax
amax
frmax
ws1.Cells(14, 5) = duc
ws1.Cells(15, 5) = eta
ws1.Cells(20,
ws1.Cells(21,
ws1.Cells(22,
ws1.Cells(23,
3)
3)
3)
3)
=
=
=
=
Wkmax
Wdmax
Wsmax
Eimax
ws1.Cells(20,
ws1.Cells(21,
ws1.Cells(22,
ws1.Cells(23,
ws1.Cells(25,
5)
5)
5)
5)
5)
=
=
=
=
=
Wk
Wd
Ws
Ei
Err
End Sub