2013年8月11日
振動問題の基礎
目次
1. 1質点系での時刻歴計算法(2012.08.12) 1 1.1 運動方程式 . . . 1 1.2 線形加速度法. . . 1 1.3 平均加速度法. . . 2 1.4 Nermarkのβ法. . . 3 2. 応答スペクトル 5 2.1 加速度応答スペクトル(2012.08.12) . . . 5 2.2 3重応答スペクトル(2013.04.06) . . . 7 3. フーリエスペクトル (2012.08.13) 8 3.1 理論 . . . 8 3.2 簡単なフーリエ変換・逆変換の事例 . . . 9 3.3 フーリエスペクトル出力事例 . . . 13 4. 模擬地震波作成(2012.08.18) 15 4.1 模擬地震波作成処理の考え方 . . . 15 4.2 地震継続時間の影響 . . . 16 4.3 位相角について(Case 1の事例) (2013.04.07) . . . 20 5. 1質点系の周波数領域の計算による時刻歴計算法 (2012.10.05) 211.
1
質点系での時刻歴計算法
(2012.08.12)
運動方程式より変位・速度・加速度の時刻歴を計算する．1.1 運動方程式
時刻t + ∆tでの質点系の釣り合いを示す運動方程式は以下の通り．ここにドットは時間微分を示す．
m· ¨u(t + ∆t) + c · ˙u(t + ∆t) + k · u(t + ∆t) = f(t + ∆t) (1)
ここに， m ：質量 u¨ ：加速度 f ：外力 c ：減衰係数 u˙ ：速度 k ：バネ定数 u ：変位 1.2 線形加速度法 (1) 変位 時刻t + ∆tでの変位uをTaylor展開し，∆tの3次項まで考慮すると， u(t + ∆t) = u(t) +∆t 1! · ˙u(t) + (∆t)2 2! · ¨u(t) + (∆t)3 3! · ... u (t) +· · · (2) = u(t) + ∆t· ˙u(t) +(∆t) 2 3 · ¨u(t) + (∆t)2 6 · ¨u(t + ∆t) (3) 基本的な考え方として，全ての項を，時刻tと時刻t + ∆tでの加速度u，速度¨ u，変位˙ uで表現したい．この ため，線形加速度法では，∆tの3次項を以下のように線形化して上式を得ている． ... u (t) = u(t + ∆t)¨ − ¨u(t) ∆t (4) (2) 速度 加速度の積分を台形公式で近似することにより ˙ u(t + ∆t) = ˙u(t) + ∫ t+∆t t ¨ u(t)dt = ˙u(t) +∆t 2 · (¨u(t + ∆t) + ¨u(t)) (5) (3) 加速度 これまでにもとめた変位・速度を運動方程式に代入することにより加速度が以下の通り求まる． ¨ u(t + ∆t) = f (t + ∆t)− c · ( ˙ u(t) +∆t 2 · ¨u(t) ) − k · ( u(t) + ∆t· ˙u(t) +(∆t) 2 3 · ¨u(t) ) m + c·∆t 2 + k· (∆t)2 6 (6) (4) まとめ 線形加速度法 ¨ u(t + ∆t) = f (t + ∆t)− c · ( ˙ u(t) +∆t 2 · ¨u(t) ) − k · ( u(t) + ∆t· ˙u(t) +(∆t) 2 3 · ¨u(t) ) m + c·∆t 2 + k· (∆t)2 6 (7) ˙ u(t + ∆t) = ˙u(t) +∆t 2 · ¨u(t) + ∆t 2 · ¨u(t + ∆t) (8)
u(t + ∆t) =u(t) + ∆t· ˙u(t) +(∆t)
2
3 · ¨u(t) +
(∆t)2
1.3 平均加速度法 (1) 変位 時刻t + ∆tでの変位uをTaylor展開し，∆tの2次項まで考慮すると， u(t + ∆t) = u(t) + ∆t 1! · ˙u(t) + (∆t)2 2! · ¨u(t) + · · · (10) = u(t) + ∆t· ˙u(t) +(∆t) 2 4 · ¨u(t) + (∆t)2 4 · ¨u(t + ∆t) (11) 平均加速度法では，∆tの2次項を以下のように線形化して上式を得ている． ¨
u(t) = u(t + ∆t) + ¨¨ u(t)
2 (12) (2) 速度 線形加速度法と同様に ˙ u(t + ∆t) = ˙u(t) + ∫ t+∆t t ¨ u(t)dt = ˙u(t) +∆t 2 · (¨u(t + ∆t) + ¨u(t)) (13) (3) 加速度 これまでにもとめた変位・速度を運動方程式に代入することにより加速度が以下の通り求まる． ¨ u(t + ∆t) = f (t + ∆t)− c · ( ˙ u(t) +∆t 2 · ¨u(t) ) − k · ( u(t) + ∆t· ˙u(t) +(∆t) 2 4 · ¨u(t) ) m + c·∆t 2 + k· (∆t)2 4 (14) (4) まとめ 平均加速度法 ¨ u(t + ∆t) = f (t + ∆t)− c · ( ˙ u(t) +∆t 2 · ¨u(t) ) − k · ( u(t) + ∆t· ˙u(t) +(∆t) 2 4 · ¨u(t) ) m + c·∆t 2 + k· (∆t)2 4 (15) ˙ u(t + ∆t) = ˙u(t) +∆t 2 · ¨u(t) + ∆t 2 · ¨u(t + ∆t) (16)
u(t + ∆t) =u(t) + ∆t· ˙u(t) +(∆t)
2
4 · ¨u(t) +
(∆t)2
1.4 Nermarkのβ法 (1) 計算式 線形加速度法を一般化した方法である．時刻t + ∆tでの変位uおよび速度u˙ をTaylor展開したものを以 下に示す． u(t + ∆t) = u(t) +∆t 1! · ˙u(t) + (∆t)2 2! · ¨u(t) + (∆t)3 3! · ... u (t) +· · · (18) ˙ u(t + ∆t) = ˙u(t) +∆t 1! · ¨u(t) + (∆t)2 2! · ... u (t) +· · · (19) ここで，変位については右辺第4項において1/3! = βとおくことにより，速度については右辺第3項におい て1/2! = γとおき，双方で高次項を無視することにより以下の式を得る．
u(t + ∆t) = u(t) + ∆t· ˙u(t) +(∆t) 2
2 · ¨u(t) + β · (∆t) 3·...
u (t) (20)
˙
u(t + ∆t) = ˙u(t) + ∆t· ¨u(t) + γ · (∆t)2·...u (t) (21)
次に，線形加速度法と同様に，...uを以下のように線形化する．
...
u (t) = u(t + ∆t)¨ − ¨u(t)
∆t (22)
上記関係を用いて，運動方程式を整理することにより，以下の方程式が得られる．
u(t + ∆t) = u(t) + ∆t· ˙u(t) + ( 1 2− β ) (∆t)2· ¨u(t) + β(∆t)2· ¨u(t + ∆t) (23) ˙
u(t + ∆t) = ˙u(t) + (1− γ) · ∆t · ¨u(t) + γ · ∆t · ¨u(t + ∆t) (24)
ここで通常は，γ = 1/2が用いられることを考慮すると，Newmarkのβ法は，以下の式で表現できる． Newmarkのβ法 ¨ u(t + ∆t) = f (t + ∆t)− c · ( ˙ u(t) +∆t 2 · ¨u(t) ) − k · ( u(t) + ∆t· ˙u(t) + ( 1 2− β ) · (∆t)2· ¨u(t) ) m + c·∆t 2 + k· β · (∆t) 2 (25) ˙ u(t + ∆t) = ˙u(t) +∆t 2 · ¨u(t) + ∆t 2 · ¨u(t + ∆t) (26)
u(t + ∆t) =u(t) + ∆t· ˙u(t) +
( 1 2− β ) · (∆t)2· ¨u(t) + β · (∆t)2· ¨u(t + ∆t) (27) ここに，平均加速度法：β = 1 4 ，線形加速度法：β = 1 6 である．
(2) 解の安定性 非減衰振動に対し，Newmarkのβ法の安定条件は，以下のとおり表される． γ= 1 2 β5 1 2 ∆t Tmin 5 1 2π√γ/2− β (28) ここに，Tminは，構造物の最小固有周期である． また無条件安定条件は，以下のとおりである． 2β= γ = 1 2 (29) ここで，γが1/2を超えると，誤差が発生することが知られており，通常γ = 1/2が採用される． ■平均加速度法 平均加速度法では，γ = 1/2，β = 1/4なので，上式の条件を満たし，無条件安定となる． ■線形加速度法 線形加速度法においては，γ = 1/2，β = 1/6なので，条件付安定となる．その条件は， ∆t Tmin 5 1 2π√γ/2− β = 0.5513 (30)
よって，時間間隔を∆t = 0.01(sec)とすればTmin= 0.018(sec)が安定条件となる．通常，地震応答スペク
トルを求める場合は，地震動の観測時間間隔は0.01(sec)であり，固有周期の計算は0.02(sec)より開始する
ので，安定条件を満たす．
しかしFEMなど多自由度系の解析では，系の最小固有周期が0.018(sec)より小さくなることも考えられ
るため，線形加速度法では解が発散する可能性がある．このため，多自由度系の解析では，無条件安定である 平均加速度法を用いるほうが有利と思われる．
2.
応答スペクトル
2.1 加速度応答スペクトル (2012.08.12) φ u 地盤の上に構築された1質点系が，地盤から加速度を受ける場合の運動方程式は以下 のとおりとなる． m· ( ¨φ + ¨u) + c · ˙u + k · u = 0 (31) 上式において，質点系に作用する加速度は，地盤に作用する加速度φ¨と構造系内の基 準点に対する加速度u¨の合計となる．これに対し，質点の運動を減衰させようとする 力および変形したバネが元に戻ろうとする力は，構造物内の基準点に対する速度u˙ お よび変位uにより表され，地盤の変位とは無関係である． 応答スペクトルは，地震動の特性を把握するために用いられるものであり，ある地盤加速度を1質点系に与 えた時の質点系の絶対加速度φ + ¨¨ u，相対変位u˙あるいは相対変位uの最大値を，質点系の固有周期あるいは 固有振動数との関係で整理したグラフとして表現される．解くべき運動方程式は，基本式を質点の質量mで 除し，地盤加速度項を右辺に移項した以下の形のものが一般に用いられる． ¨ u + 2· h · ω0· ˙u + ω20· u = − ¨φ (32) ここに，h：系の減衰定数，ω0：系の固有円振動数であり，質点の質量mや減衰係数c，バネ定数kなどと以 下の関係を有する． c = 2· h√k· m ω0= √ k m T = 2π ω0 = 2π √ m k (33) 上記T は系の固有周期である． ここでは，地盤加速度を入力値として，加速度応答スペクトルを計算するプログラムを考える． 運動方程式としては一般的な構造物の動的応答解析に用いることを考慮して以下に示す一般形を用い，入力 値を地盤加速度φ，減衰定数¨ h，固有周期Tとするため，以下の関係を導入する． m· ¨u + c · ˙u + k · u = −m · ¨φ (34) m = 1 k = 4π 2· m T2 c = 2· h √ k· m (35) 地盤加速度を∆tピッチで与えるとすると，実際の適用を考え，固有周期はT = 2· ∆t ∼ 10(sec)の範囲と し，減衰定数はh = 0.05とする． なお，応答スペクトル図を作成する場合の注意点は，応答加速度としては「絶対加速度u + ¨¨ φ」の絶対値の 最大を検索・出力することである．すなわち，運動方程式を解いた場合，相対加速度としてu¨が求まるが，応 答スペクトルを作成する加速度としては，これに地盤加速度φ¨を加えた加速度を用いることに注意する． 用いたプログラムでは，運動方程式を線形加速度法で解いている．また，入力加速度の基線補正を組み込んでいる．加速度基線補正のための数値積分subroutine IACCおよび加速度基線補正subroutine CRAC は大
崎順彦先生の著書「新・地震動のスペクトル解析入門」に掲載のFORTRANソースをFortran90 に書き直
したものである．
計算に用いた地震加速度波形は，(独)防災科学技術研究所が運営しているKiK-netのデジタル波形デー
タを利用させていただきました．
入力データとした地盤の加速度時刻歴を図1に，線形加速度法により計算した加速度応答スペクトルを図2
に示す．赤線は，大崎先生の著書による地震応答スペクトル計算sabroutine ERESをFortran90に書き換え
−807.7 0.0 807.7 Acceleration (gal) 0 30 60 90 120 150 180 210 240 270 300 time (sec) TCGH16_UD2 −798.6 0.0 798.6 Acceleration (gal) 0 30 60 90 120 150 180 210 240 270 300 time (sec) TCGH16_NS2 −1196.7 0.0 1196.7 Acceleration (gal) 0 30 60 90 120 150 180 210 240 270 300 time (sec) TCGH16_EW2 図1 入力データとした地盤の加速度時刻歴 10 100 1000 10000
Response acceleration (gal)
0.01 0.1 1 10 Period (sec) 10 100 1000 10000
Response acceleration (gal)
0.01 0.1 1 10 Period (sec) 10 100 1000 10000
Response acceleration (gal)
0.01 0.1 1 10 Period (sec) 10 100 1000 10000
Response acceleration (gal)
0.01 0.1 1 10 Period (sec) 10 100 1000 10000
Response acceleration (gal)
0.01 0.1 1 10 Period (sec) 10 100 1000 10000
Response acceleration (gal)
0.01 0.1 1 10 Period (sec) 10 100 1000 10000
Response acceleration (gal)
0.01 0.1 1 10 Period (sec) Damping factor h=0.05 TCGH16_EW2 TCGH16_NS2 TCGH16_UD2 図2 線形加速度法により計算した加速度応答スペクトル
2.2 3重応答スペクトル (2013.04.06) 3重応答スペクトルは，速度応答スペクトル図を基本に，加速度軸および変位軸を加えたものである．近似 的にではあるが，1枚のグラフで加速度・速度・変位の応答地を示すことができる． 加速度軸および変位軸は，以下の関係を用いて定められている． Sa + 2π T · Sv Sd+ T 2π· Sv (36) ここに， Sa : 最大絶対加速度応答 Sv : 最大相対速度応答 Sd : 最大相対変位応答 1 2 5 10 20 50 100 200 500
Response velocity (kine)
0.02 0.05 0.1 0.2 0.5 1 2 5 10 Period (sec) 1 2 5 10 20 50 100 200 500
Response velocity (kine)
0.02 0.05 0.1 0.2 0.5 1 2 5 10 Period (sec) 1 2 5 10 20 50 100 200 500
Response velocity (kine)
0.02 0.05 0.1 0.2 0.5 1 2 5 10 Period (sec) 1 2 5 10 20 50 100 200 500
Response velocity (kine)
0.02 0.05 0.1 0.2 0.5 1 2 5 10 Period (sec) 1 2 5 10 20 50 100 200 500
Response velocity (kine)
0.02 0.05 0.1 0.2 0.5 1 2 5 10 Period (sec) 1 2 5 10 20 50 100 200 500
Response velocity (kine)
0.02 0.05 0.1 0.2 0.5 1 2 5 10 Period (sec) 1 2 5 10 20 50 100 200 500
Response velocity (kine)
0.02 0.05 0.1 0.2 0.5 1 2 5 10 Period (sec) 1 2 5 10 20 50 100 200 500
Response velocity (kine)
0.02 0.05 0.1 0.2 0.5 1 2 5 10 Period (sec) 1 2 5 10 20 50 100 200 500 1000 2000 5000 10000 20000 50000 100000 0.01 0.02 0.05 0.1 0.2 0.5 1 2 5 10 20 50 100 200 500 Damping factor h=0.05
Response acceleration (gal)
Response displacement (cm)
TCGH16_EW2 TCGH16_NS2 TCGH16_UD2
3.
フーリエスペクトル
(2012.08.13)
3.1 理論 (1) 有限フーリエ係数 時刻に関する連続関数x(t)について，以下の時刻歴データがあるものとする． 標本数 N 標本点間隔 ∆t 標本値 xm m = 0, 1, 2,· · · , N − 1 標本数N は高速フーリエ変換を適用することを前提にN = 2L（Lは正の整数）とし，標本値x m を k = N/2までの有限三角級数で表せるものとすれば xm= N/2 ∑ k=0 { Ak· cos ( 2πk N· ∆t · t ) + Bk· sin ( 2πk N· ∆t· t )} =A0 2 + N/2∑−1 k=1 ( Ak· cos 2πkm N + Bk· sin 2πkm N ) +An/2 2 · cos 2π(N/2)m N m = 0, 1, 2,· · · , N − 1 t = m· ∆t 上式を，関数x(t)の有限フーリエ近似，係数Ak, Bkを有限フーリエ係数という．また，有限フーリエ係数 を求める演算をフーリエ変換(Fourier transform)，有限フーリエ係数が既知で元の標本値を再現する演算をフーリエ逆変換(Fourier inverse transform)という．
Ak = 2 N N∑−1 m=0 xm· cos 2πkm N k = 0, 1, 2,· · · , N/2 − 1, N/2 (37) Bk = 2 N N∑−1 m=0 xm· sin 2πkm N k = 1, 2,· · · , N/2 − 1 (38) A0 2 = 1 N N∑−1 m=0 xm (全標本値の平均) (39) (2) 有限フーリエ変換の複素数表示
オイラーの公式（e±iθ = cos θ± i sin θ）により複素数表示を行う．
フーリエ変換 Ck= 1 N N∑−1 m=0 xm· e−i(2πkm/N) k =0, 1, 2,· · · , N − 1 (40) フーリエ逆変換 xm= N∑−1 k=0 Ck· e+i(2πkm/N ) m =0, 1, 2,· · · , N − 1 (41) (3) フーリエ係数と複素フーリエ係数の関係 Ck= Ak− iBk 2 Ak+ iBk= AN−k− iBN−k Ak= 2· Re(Ck) Bk=−2 · Im(Ck) k = 0, 1, 2,· · · , N/2
(4) フーリエ振幅スペクトル，フーリエ位相スペクトル 振動数 fk = k N· ∆t k = 0, 1, 2,· · · , N/2 フーリエ振幅 Fk = N· ∆t 2 √ Ak2+ Bk2= N· ∆t|Ck| = N · ∆t √ Re(Ck)2+ Im(Ck)2 位相角 φk= arctan ( −Bk Ak ) = arctan ( Im(Ck) Re(Ck) ) (−π < φk< π) 位相波 x(t) =ˆ N/2 ∑ k=1 cos(2π· fk· t + φk) フーリエ振幅を振動数に対してプロットしたものをフーリエ振幅スペクトル(単にフーリエスペクトルとも いう)，位相角を振動数に対してプロットしたものをフーリエ位相スペクトルという． 3.2 簡単なフーリエ変換・逆変換の事例 以下に簡単なFourier変換・逆変換の結果を示す．
データ（data）はFortran90の乱数発生サブルーチンrandom numberを用いて16個の一様乱数を発
生させたものである．データの並び順は発生させた乱数の順番である．
Re(Ck)およびIm(Ck)は，これをFFTにかけて複素フーリエ係数を求めたものである．
abs(Ck)は複素フーリエ係数の絶対値，fp(Ck)は位相角である．
Re(xm)およびIm(xm)は，上記で求めた複素フーリエ係数を逆変換して，入力データの再現を行った
ものである．
i data Re(Ck) Im(Ck) abs(Ck) fp(Ck) Re(xm) Im(xm)
0 0.998 0.478 0.000 0.478 0.000 0.998 0.000 1 0.567 0.154 -0.014 0.154 -5.171 0.567 0.000 2 0.966 -0.003 -0.053 0.053 -93.070 0.966 0.000 3 0.748 -0.018 -0.008 0.020 -155.386 0.748 0.000 4 0.367 0.057 -0.014 0.059 -14.125 0.367 0.000 5 0.481 0.000 0.092 0.092 89.861 0.481 0.000 6 0.074 0.012 0.030 0.033 67.645 0.074 0.000 7 0.005 0.027 -0.047 0.054 -60.520 0.005 0.000 8 0.347 0.062 0.000 0.062 0.000 0.347 0.000 9 0.342 0.027 0.047 0.054 60.520 0.342 0.000 10 0.218 0.012 -0.030 0.033 -67.645 0.218 0.000 11 0.133 0.000 -0.092 0.092 -89.861 0.133 0.000 12 0.901 0.057 0.014 0.059 14.125 0.901 0.000 13 0.387 -0.018 0.008 0.020 155.386 0.387 0.000 14 0.445 -0.003 0.053 0.053 93.070 0.445 0.000 15 0.662 0.154 0.014 0.154 5.171 0.662 0.000 Ave. 0.478 上表より以下のことがわかり，理屈どおりの結果となっている． フーリエ複素係数実数部は，最初の項(i=0)は入力全データの平均値となっている．また入力データ数 をn個とすればn/2+1(i=8)を中心に対称の分布をしている． フーリエ複素係数虚数部は，最初の項(i=0)は0となっている．また入力データ数をn個とすれば n/2+1(i=8)を中心に符号が逆転した対称の分布をしている． フーリエ逆変換の結果は，実数部は入力データを再現しており，虚数部は誤差の範囲で全て0となって
いる． ここでは南茂夫先生の「科学計測のための波形データ処理 計測システムにおけるマイコン／パソコン活用 技術」に掲載されている実数演算によるBASIC プログラムを，Fortran90に書き換えたプログラムを用いて いる． このプログラムの使用上の注意点は以下のとおり． 1 フーリエ変換時 式(40)の∑Nm=0−1xm· e−i(2πkm/N)の部分を出力するため，実際のフーリエ複素係数とす るためには出力値をデータ数Nで除す必要がある． 2 フーリエ逆変換時 与えられたフーリエ複素係数から，逆変換して加速度値を求めるためには，式(40)と 式(41)の指数部を比較して判るように，虚数部の符号を逆転させて入力する必要があることに注意 する． 本文にコードを示すのはなんとなく見づらくなる気がするが，非常に重要なプログラムなので，以下にコー ドを示しておく． module d e f p i i m p l i c i t none r e a l ( 8 ) , parameter : : p i = 3 . 1 4 1 5 9 2 6 5 3 5 8 9 7 9 3 2 3 8 4 6D0 end module d e f p i program f90 FFT !∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ ! FFT and I n v e r s e FFT !∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ use d e f p i i m p l i c i t none inte ger : : i !
inte ger : : ndata ! Number o f o r i g i n a l d a t a
inte ger : : nn ! Number o f c a l c u l a t e d d a t a ( p o w e r s o f 2 ) r e a l ( 8 ) : : dt ! Time i n c r e m e n t r e a l ( 8 ) , a l l o c a t a b l e : : a c c i n p ( : ) ! I n p u t t e d o r i g i n a l d a t a r e a l ( 8 ) , a l l o c a t a b l e : : x r ( : ) ! R e a l p a r t o f c a l c u l a t e d d a t a r e a l ( 8 ) , a l l o c a t a b l e : : x i ( : ) ! I m a g i n a r y p a r t o f c a l c u l a t e d d a t a r e a l ( 8 ) , a l l o c a t a b l e : : c r ( : ) ! Memory o f r e a l p a r t r e a l ( 8 ) , a l l o c a t a b l e : : c i ( : ) ! Remory o f i m a g i n a r y p a r t r e a l ( 8 ) , a l l o c a t a b l e : : ck ( : ) ! A b s o l u t e o f c o m p l e x F o u r i e r c o e f f i c i e n t r e a l ( 8 ) , a l l o c a t a b l e : : f p ( : ) ! F o u r i e r p h a s e s p e c t r u m
character : : fnameR∗50 ,fnameW∗50
character : : s t r c o m∗200 ,dummy∗50 , fmt1 ∗200 character : : l i n e b u f f∗1000
c a l l g e t a r g ( 1 , fnameR ) ! I n p u t d a t a f i l e name c a l l g e t a r g ( 2 , fnameW ) ! I n p u t o u t p u t f i l e name open ( 1 1 , f i l e =fnameR , status= ’ o l d ’ )
read ( 1 1 , ’ ( a ) ’ ) s t r c o m read ( 1 1 ,∗ ) dummy, dt read ( 1 1 ,∗ ) dummy, ndata ! To s e t number o f d a t a ! ( nn i s e q u a l t o 2 t o t h e power o f p o s i t i v e i n t e g e r number ! and nn i s g r e a t e r t h a n o r e q u a l t o number o f d a t a ) nn=2 do nn=nn∗2 i f ( ndata∗1<=nn ) exit end do a l l o c a t e ( a c c i n p ( 1 : nn ) ) a l l o c a t e ( x r ( 1 : nn ) ) a l l o c a t e ( x i ( 1 : nn ) ) a l l o c a t e ( c r ( 1 : nn ) ) a l l o c a t e ( c i ( 1 : nn ) ) a l l o c a t e ( ck ( 1 : nn ) ) a l l o c a t e ( f p ( 1 : nn ) )
do i =1 ,nn a c c i n p ( i ) = 0 . 0D0 x r ( i ) = 0 . 0D0 x i ( i ) = 0 . 0D0 end do do i =1 , ndata read ( 1 1 ,∗ ) accinp ( i ) x r ( i )= a c c i n p ( i ) end do c l o s e ( 1 1 ) c a l l FFT( nn , xr , x i ) ! F o u r i e r t r a n s f o r m do i =1 ,nn ! To d i v i d e r e t u r n e d v a l u e s by number o f d a t a x r ( i )= x r ( i ) / d b l e ( nn ) x i ( i )= x i ( i ) / d b l e ( nn ) c r ( i )= x r ( i ) ! To memory r e a l p a r t o f t r a n s f o r m e d v a l u e s c i ( i )= x i ( i ) ! To memory i m a g i n a l p a r t o f t r a n s f o r m e d v a l u e s ck ( i )= s q r t ( c r ( i )∗ ∗ 2 . 0D0+c i ( i ) ∗ ∗ 2 . 0D0) f p ( i ) = 0 . 0D0 i f ( 1 . 0 D−30<abs ( cr ( i ) ) ) then f p ( i )= a t a n ( c i ( i ) / c r ( i ) ) / p i∗180.0D0 i f ( c r ( i ) <0.0D0 . and . c i ( i ) <0.0D0) f p ( i )= f p ( i )−180.0D0 i f ( c r ( i ) <0.0D0 . and . 0 . 0 D0<c i ( i ) ) f p ( i ) = 1 8 0 . 0D0+f p ( i ) e l s e i f ( 0 . 0 D0<c i ( i ) . and . 0 . 0 D0<c r ( i ) ) f p ( i ) = 9 0 . 0D0 i f ( c i ( i ) <0.0D0 . and . c r ( i ) <0.0D0) f p ( i )=−90.0D0 i f ( 0 . 0 D0<c i ( i ) . and . c r ( i ) <0.0D0) f p ( i ) = 9 0 . 0D0 i f ( c i ( i ) <0.0D0 . and . 0 . 0 D0<c r ( i ) ) f p ( i )=−90.0D0 end i f end do do i =1 ,nn ! To c h a n g e t h e s i g n o f i m a g i n a r y p a r t s f o r i n v e r s e FFT x r ( i )= x r ( i ) x i ( i )=−x i ( i ) end do Call FFT( nn , xr , x i ) ! I n v e r s e F o u r i e r t r a n s f o r m fmt1=” ( i 5 , ’ , ’ , e 1 5 . 7 , 4 ( ’ , ’ , e 1 5 . 7 ) , 2 ( ’ , ’ , e 1 5 . 7 ) ) ” open ( 1 2 , f i l e =fnameW , status= ’ r e p l a c e ’ )
write ( 1 2 , ’ ( a ) ’ ) ’ i , data , Re ( Ck ) , Im ( Ck ) , a b s ( Ck ) , f p ( Ck ) , Re (xm) , Im (xm) ’ do i =1 ,nn write ( l i n e b u f f ,FMT=fmt1 ) i , a c c i n p ( i ) , c r ( i ) , c i ( i ) , ck ( i ) , f p ( i ) , x r ( i ) , x i ( i ) c a l l d e l s p a c e s ( l i n e b u f f ) write ( 1 2 , ’ ( a ) ’ ) t r i m ( l i n e b u f f ) end do c l o s e ( 1 2 ) stop contains subroutine FFT( nn , xr , x i ) !∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ ! F a s t F o u r i e r Transform , I n v e r s e t r a n s f o r m !∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ ! nn : Number o f d a t a ( p o w e r s o f 2 ) ! x r ( ) : R e a l p a r t o f i / o d a t a ! x i ( ) : I m a g i n a r y p a r t o f i / o d a t a integer , intent ( in ) : : nn r e a l ( 8 ) , intent ( inout ) : : x r ( 1 : nn ) r e a l ( 8 ) , intent ( inout ) : : x i ( 1 : nn ) integer : : g , h , i , j , k integer : : l ,m, n , p , q r e a l ( 8 ) : : a , b , xd r e a l ( 8 ) , a l l o c a t a b l e : : s ( : ) r e a l ( 8 ) , a l l o c a t a b l e : : c ( : ) n=nn a l l o c a t e ( s ( 1 : n /2+1)) a l l o c a t e ( c ( 1 : n /2+1)) i =0; j =0; k =0; l =0;p=0;h=0; g =0; q=0 m=i n t ( l o g ( d b l e ( n ) ) / l o g ( 2 . 0 D0) + 1 . 0D0) a =0.0D0
b=p i∗2.0D0/ dble ( n ) do i =0 ,n /2 s ( i +1)= s i n ( a ) c ( i +1)= c o s ( a ) a=a+b end do l=n h=1 do g =1 ,m l=l /2 k=0 do q=1 ,h p=0 do i=k , l+k−1 j=i+l a=x r ( i +1)−xr ( j +1) b=x i ( i +1)−x i ( j +1) x r ( i +1)=x r ( i +1)+x r ( j +1) x i ( i +1)= x i ( i +1)+ x i ( j +1) i f ( p==0)then x r ( j +1)=a x i ( j +1)=b e l s e x r ( j +1)=a∗ c ( p+1)+b∗ s ( p+1) x i ( j +1)=b∗ c ( p+1)−a∗ s ( p+1) end i f p=p+h end do k=k+l+l end do h=h+h end do j=n /2 do i =1 ,n−1 k=n i f ( j <i ) then xd=x r ( i +1) x r ( i +1)=x r ( j +1) x r ( j +1)=xd xd=x i ( i +1) x i ( i +1)= x i ( j +1) x i ( j +1)=xd end i f k=k /2 do while ( j>=k ) j=j−k k=k /2 i f ( k==0) e x i t end do j=j+k end do end subroutine FFT subroutine d e l s p a c e s ( s )
character (∗ ) , intent ( inout ) : : s character ( len=len ( s ) ) tmp integer i , j j = 1 do i = 1 , len ( s ) i f ( s ( i : i )== ’ ’ ) c y c l e tmp ( j : j ) = s ( i : i ) j = j + 1 end do s = tmp ( 1 : j−1) end subroutine d e l s p a c e s end program f90 FFT
3.3 フーリエスペクトル出力事例 フーリエスペクトルを描画する際は，以下の関係を描画する． 振動数 fk= k N· ∆t (k = 0∼ N/2) フーリエ振幅 Fk= N· ∆t √ Re(Ck)2+ Im(Ck)2 (k = 0∼ N/2) 図1に示した地盤の加速度時刻歴より計算した，フーリエスペクトルを図4に示す．図4では，スペクトル の計算値そのもの（細い赤線）では，変化が激しく，振動数のピークを読み取りにくいので，Parzen window による平滑化を行っている．
Parzen winndowによる平滑化subroutine SWINは大崎順彦先生の著書「新・地震動のスペクトル解析入
0.1 1 10 100 1000 10000
Fourier spectrum (gal*sec)
0 5 10 15 20 Frequency (Hz) 0.1 1 10 100 1000 10000
Fourier spectrum (gal*sec)
0 5 10 15 20 Frequency (Hz) 0.1 1 10 100 1000 10000
Fourier spectrum (gal*sec)
0 5 10 15 20 Frequency (Hz) TCGH16_EW2 Band=1.0Hz 0.1 1 10 100 1000 10000
Fourier spectrum (gal*sec)
0 5 10 15 20 Frequency (Hz) 0.1 1 10 100 1000 10000
Fourier spectrum (gal*sec)
0 5 10 15 20 Frequency (Hz) 0.1 1 10 100 1000 10000
Fourier spectrum (gal*sec)
0 5 10 15 20 Frequency (Hz) TCGH16_NS2 Band=1.0Hz 0.1 1 10 100 1000 10000
Fourier spectrum (gal*sec)
0 5 10 15 20 Frequency (Hz) 0.1 1 10 100 1000 10000
Fourier spectrum (gal*sec)
0 5 10 15 20 Frequency (Hz) 0.1 1 10 100 1000 10000
Fourier spectrum (gal*sec)
0 5 10 15 20 Frequency (Hz) TCGH16_UD2 Band=1.0Hz 図4 フーリエスペクトル出力事例 （黒線：平滑化後，赤線：平滑化前）
4.
模擬地震波作成
(2012.08.18)
4.1 模擬地震波作成処理の考え方 基本フローは，「国土技術政策総合研究所資料（国総研資料第244号：平成17年3月）：大規模地震に 対するダムの耐震性能照査に関する資料，資料1-6 加速度応答スペクトルに適合する時刻歴波形の作 成方法」を参考に作成しています． 機能は，目標とする加速度応答スペクトルと位相特性を参照する原種波形時刻歴を入力することによ り，目標スペクトルに適合させた模擬地震波を作成するものです． 適合の許容誤差は5%としていますが，収束しない場合もあるので，スペクトル比より計算した誤差が 減少傾向から増加に転じた時点で収束計算を打ち切っています． 模擬地震波の継続時間は，原種波形と同時刻に始まり同時刻に終了・打ち切っています． 目標スペクトルは，周期の関数として数字入力しますが，内部処理ではFourier変換を行う関係上，周 波数の関数として扱っています． 原種波形時刻歴入力 目標スペクトルSa入力 模擬地震波初期値セット (原種波形時刻歴) 模擬地震波の 応答スペクトルSa0 算定 応答スペクトル比算定 rk= Sa/S0a 適合性判定 er(rk)5 0 Fourier変換による Fourier係数Ck算出 Fourier係数補正 Ck= Ck× rk Fourier逆変換による 模擬地震波時刻歴作成 完了 ? ? ? ? ? ?No ? ? -Yes 適合性判定条件 er(rk) = v u u u u t n ∑ k=1 (1− rk)2 n 5 0 rk= Sa Sa0 0 ：許容値 error ：誤差 Sa ：目標スペクトル Sa0 ：模擬地震波スペクトル n ：スペクトル比較する周期の数 図5 模擬地震波作成手順4.2 地震継続時間の影響
前出の考え方で作成した模擬地震波作成プログラムを用いて，地震継続時間の影響をテストしてみる．
地震継続時間は，「設計用入力地震動作成手法技術指針(案)，建設省建築研究所・(財)日本建築センター，
平成4年3月」を参考に以下の考え方で設定した．
設定ケース 包絡関数
Case a(sec) b(sec) c(sec)
1 5 15 30 2 5 25 60 3 5 35 120 e(t) = (t/a)2 e(t) = 1.0 e(t) = exp{−α · (t − b)} α =−ln(0.1) c− b e(t) = 0.1 at t = c
0
1
e(t)
0
30
60
90
120
150
time t (sec)
0
1
e(t)
0
30
60
90
120
150
time t (sec)
0
1
e(t)
0
30
60
90
120
150
time t (sec)
0
1
e(t)
0
30
60
90
120
150
time t (sec)
0
1
e(t)
0
30
60
90
120
150
time t (sec)
e=0.1 a=5 b=15 b=25 b=35 c=30 c=60 c=120Envelop function of acceleration wave due to earthquake
図6 包絡関数
テスト用の波形は以下の要領で作成した．
Fortran90のrandom numberにより0∼ 1の範囲の一様乱数を必要個数発生させる．
上記で得た乱数列の範囲を−1 ∼ 1に拡張し，これに包絡関数と設定した最大加速度を乗じることによ りテスト波形を得る． サンプリングピッチは∆t = 0.01(sec)とする． 最大加速度は100galとする． 上記で得たテスト波形より，目標スペクトルに適合する模擬波形を作成する．目標スペクトルにおける最大 加速度は，300galである． 以降に各ケースにおける原種波形，模擬波形，加速度応答スペクトルを示す．なお，実際の耐震解析におい ては，乱数による位相の使用の是非の議論もあるようであるが，ここでは作成した模擬地震波の特性をつかむ 数値実験と割り切り，計算を行っている．
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
Original−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
Simulated0.01
0.1
1
10
100
1000
10000
Response acceleration (gal)
0.01
0.1
1
10
100
1000
Period (sec)
0.01
0.1
1
10
100
1000
10000
Response acceleration (gal)
0.01
0.1
1
10
100
1000
Period (sec)
0.01
0.1
1
10
100
1000
10000
Response acceleration (gal)
0.01
0.1
1
10
100
1000
Period (sec)
0.01
0.1
1
10
100
1000
10000
Response acceleration (gal)
0.01
0.1
1
10
100
1000
Period (sec)
Damping factor h=0.05 Target Simulated Original 図7 Case 1 (継続時間 C=30 sec)−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
Original−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
Simulated0.01
0.1
1
10
100
1000
10000
Response acceleration (gal)
0.01
0.1
1
10
100
1000
Period (sec)
0.01
0.1
1
10
100
1000
10000
Response acceleration (gal)
0.01
0.1
1
10
100
1000
Period (sec)
0.01
0.1
1
10
100
1000
10000
Response acceleration (gal)
0.01
0.1
1
10
100
1000
Period (sec)
0.01
0.1
1
10
100
1000
10000
Response acceleration (gal)
0.01
0.1
1
10
100
1000
Period (sec)
Damping factor h=0.05 Target Simulated Original 図8 Case 2 (継続時間 C=60 sec)−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
Original−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
−400
−300
−200
−100
0
100
200
300
400
Acceleration (gal)
0
30
60
90
120
time (sec)
Simulated0.01
0.1
1
10
100
1000
10000
Response acceleration (gal)
0.01
0.1
1
10
100
1000
Period (sec)
0.01
0.1
1
10
100
1000
10000
Response acceleration (gal)
0.01
0.1
1
10
100
1000
Period (sec)
0.01
0.1
1
10
100
1000
10000
Response acceleration (gal)
0.01
0.1
1
10
100
1000
Period (sec)
0.01
0.1
1
10
100
1000
10000
Response acceleration (gal)
0.01
0.1
1
10
100
1000
Period (sec)
Damping factor h=0.05 Target Simulated Original 図9 Case 3 (継続時間 C=120 sec)4.3 位相角について (Case 1 の事例) (2013.04.07) 位相角は，原種波形とこれより作成された模擬地震波で全く等しい． すなわち，この方法では原種波形の位相角の情報はそのまま保持され，フーリエ振幅のみを目標スペク トルに整合させるよう変化させていることがわかる． 模擬地震波の位相角より作成した位相波時刻歴は，原種波形の時刻歴と似通った形状となっている． 上記の理由は不明であるが，位相波には振幅情報は入っていないにもかかわらず，位相波の形状が，継 続時間や包絡曲線に大きく影響していると考えることができる．
−180
−90
0
90
180
Phase angle (degree)
0
10
20
30
40
50
Frequency (Hz)
−180
−90
0
90
180
Phase angle (degree)
0
10
20
30
40
50
Frequency (Hz)
Fourier Phase Spectrum (simulated wave)
−150
−100
−50
0
50
100
150
Acceleration (gal)
0
10
20
30
Time (sec)
−150
−100
−50
0
50
100
150
Acceleration (gal)
0
10
20
30
Time (sec)
Acceleration data (input wave)
−150
−100
−50
0
50
100
150
Acceleration (gal)
0
10
20
30
Time (sec)
−150
−100
−50
0
50
100
150
Acceleration (gal)
0
10
20
30
Time (sec)
Phase wave Acceleration data (simulated wave)
5.
1
質点系の周波数領域の計算による時刻歴計算法
(2012.10.05)
時刻歴計算法には，時間領域の数値積分によるものと，フーリエ変換・逆変換を活用した周波数領域の計算
によるものがある．ここでは，周波数領域の計算による1質点系の時刻歴計算法を述べる．
運動方程式は，以下の形式とする． ¨
u(t) + 2· h · ω0· ˙u(t) + ω02· u(t) = −a(t) (42)
ここに，u：質点の相対変位，a：地盤加速度，h：系の減衰定数，ω0：系の固有円振動数である． 地盤加速度および質点の変位・速度・加速度をフーリエ逆変換を用いて表すと， a(t) = 1 2π ∫ ∞ −∞ A(ω)eiωtdω (43) u(t) = 1 2π ∫ ∞ −∞ U (ω)eiωtdω (44) ˙ u(t) = 1 2π ∫ ∞ −∞
iω· U(ω)eiωtdω (45)
¨ u(t) = 1 2π ∫ ∞ −∞−ω 2· U(ω)eiωtdω (46) 式(43)(44)(45)(46)を式(42)に代入することにより， U (ω) = A(ω) (ω2− ω2 0)− 2hω0ω· i (47) u(t) = 1 2π ∫ ∞ −∞ A(ω) (ω2− ω2 0)− 2hω0ω· i eiωtdω (48) ˙ u(t) = 1 2π ∫ ∞ −∞ iωA(ω) (ω2− ω2 0)− 2hω0ω· i eiωtdω (49) ¨ u(t) = 1 2π ∫ ∞ −∞ −ω2A(ω) (ω2− ω2 0)− 2hω0ω· i eiωtdω (50) として，応答値である変位・速度・加速度が表現できる．ここで， H(ω) = 1 (ω2− ω2 0)− 2hω0ω· i (51) は伝達関数あるいは周波数応答関数と呼ばれている． 応答値の計算手順は以下のとおり． 1 入力加速度aの複素フーリエ係数A(ω)の算出 2 伝達関数H(ω)の算出 3 A(ω) · H(ω)のフーリエ逆変換による変位時刻歴u(t)の算出 4
iω · A(ω) · H(ω)のフーリエ逆変換による速度時刻歴u(t)˙ の算出 5 −ω2· A(ω) · H(ω)のフーリエ逆変換による加速度時刻歴u(t)¨ の算出 実際の計算では高速フーリエ変換を用いるため，以下に注意する． データ数は有限でN個（Nは2の累乗）とし，データの時間間隔を∆tとする． 計算される振動数fk および円振動数ωは以下のとおりとなる．ここでデータはN 個あるが実態とし ての振動数はN/2 + 1個しか定義できないため，N/2 + 1個目のデータを中心に対称形にN/2 + 2∼ N 個目の円振動数を定義しておく． fk = k N· ∆t ωk= 2πfk k = 0∼ N/2 (52) 複素フーリエ係数の演算も同様に，係数はN/2 + 1個目のデータを中心に対称形に格納するが，虚数 部は符号を逆転させる．
入力加速度時刻歴 a(t) −2000 −1500 −1000 −500 0 500 1000 1500 2000
Input acceleration (gal)
0 50 100 150 200 250 300 Time (sec) −2000 −1500 −1000 −500 0 500 1000 1500 2000
Input acceleration (gal)
0 50 100 150 200 250 300 Time (sec) 入力加速度フーリエ変換 |A(ω)| 0 1 2 3 4 5 |A| 0 10 20 30 40 50 Frequency (Hz) 0 1 2 3 4 5 |A| 0 10 20 30 40 50 Frequency (Hz) 伝達関数 |H(ω)| 0.000 0.001 0.002 0.003 0.004 0.005 |H| 0 10 20 30 40 50 Frequency (Hz) 0.000 0.001 0.002 0.003 0.004 0.005 |H| 0 10 20 30 40 50 Frequency (Hz) 応答加速度フーリエ変換 | − ω2· A(ω) · H(ω)| 0 1 2 3 4 5 |− ω 2*A*H| 0 10 20 30 40 50 Frequency (Hz) 0 1 2 3 4 5 |− ω 2*A*H| 0 10 20 30 40 50 Frequency (Hz) 応答絶対加速度時刻歴 a(t) + ¨u(t) −2000 −1500 −1000 −500 0 500 1000 1500 2000
Responce acceleration (gal)
0 50 100 150 200 250 300 Time (sec) −2000 −1500 −1000 −500 0 500 1000 1500 2000
Responce acceleration (gal)
0 50 100 150 200 250 300 Time (sec)
図11 1質点系の周波数領域の計算による応答加速度計算のステップ
10 100 1000 10000
Response acceleration (gal)
0.01 0.1 1 10 Period (sec) 10 100 1000 10000
Response acceleration (gal)
0.01 0.1 1 10 Period (sec) 10 100 1000 10000
Response acceleration (gal)
0.01 0.1 1 10 Period (sec) 10 100 1000 10000
Response acceleration (gal)
0.01 0.1 1 10 Period (sec) 10 100 1000 10000
Response acceleration (gal)
0.01 0.1 1 10 Period (sec) 10 100 1000 10000
Response acceleration (gal)
0.01 0.1 1 10 Period (sec) 10 100 1000 10000
Response acceleration (gal)
0.01 0.1 1 10 Period (sec) Damping factor h=0.05 TCGH16_EW2 TCGH16_NS2 TCGH16_UD2 図12 周波数領域の計算による加速度応答スペクトル 上図は，2章で用いた(独)防災科学技術研究所が運営しているKiK-netのデジタル波形データを用いて 描画した加速度応答スペクトルである．赤線は，大崎先生の著書による地震応答スペクトル計算sabroutine ERESをFortran90に書き換えたものを使用した場合のスペクトルであるが，両者はほとんど一致している．
