2012年10月6日
有限要素法による平面骨組振動解析
目次
1. マトリックス形式の運動方程式 (2012.08.14) 1 1.1 運動方程式 . . . 1 1.2 Newmarkのβ法により定式化した解くべき方程式 . . . 1 1.3 行列およびベクトルの要素 . . . 2 1.4 減衰マトリックス(2012.08.19) . . . 4 2. 有限要素法プログラムとしての入出力 5 2.1 固有値解析 . . . 5 2.2 時刻歴応答解析 . . . 6 3. 固有値解析(2012.08.20) 8 3.1 固有振動数と固有振動モードの算定方法の直感的解説 . . . 8 3.2 固有振動数の数値解析 . . . 8 4. 時刻歴応答解析(2012.08.23-2012.08.27) 11 4.1 非減衰振動 . . . 11 4.2 強制振動 . . . 23 4.3 減衰振動 . . . 261.
マトリックス形式の運動方程式
(2012.08.14)
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 Newmarkのβ法により定式化した解くべき方程式 (1) 加速度未知数 ( [m] +∆t 2 [c] + β(∆t) 2[k] ) {¨u(t + ∆t)} = {f(t + ∆t)} − [c]{va(t)} − [k]{vb(t)} (2) {va(t)} = { ˙u(t)} + ∆t 2 {¨u(t)} (3) {vb(t)} = {u(t)} + ∆t{ ˙u(t)} + ( 1 2− β ) (∆t)2{¨u(t)} (4) { ˙u(t + ∆t)} ={ ˙u(t)} +∆t 2 {¨u(t)} + ∆t 2 {¨u(t + ∆t)} (5)
{u(t + ∆t)} ={u(t)} + ∆t{ ˙u(t)} +
( 1 2− β ) (∆t)2{¨u(t)} + β(∆t)2{¨u(t + ∆t)} (6) (2) 変位未知数 ( 1 β(∆t)2[m] + 1 2β∆t[c] + [k] ) {u(t + ∆t)} = {f(t + ∆t)} + [m]{wa(t)} + [c]{wb(t)} (7) {wa(t)} = 1 β(∆t)2{u(t)} + 1 β∆t{ ˙u(t)} + ( 1 2β − 1 ) {¨u(t)} (8) {wb(t)} = 1 2β∆t{u(t)} + ( 1 2β − 1 ) { ˙u(t)} + ( 1 4β − 1 ) ∆t{¨u(t)} (9) { ˙u(t + ∆t)} = 1 2β∆t{u(t + ∆t)} − 1 2β∆t{u(t)} − ( 1 2β − 1 ) { ˙u(t)} − ( 1 4β − 1 ) ∆t{¨u(t)} (10) {¨u(t + ∆t)} = 1 β(∆t)2{u(t + ∆t)} − 1 2β(∆t)2{u(t)} − 1 β∆t{ ˙u(t)} − ( 1 2β − 1 ) {¨u(t)} (11) ここに，平均加速度法：β = 1 4 ，線形加速度法：β = 1 6 である． 加速度を未知数とした式は，1質点系での式をそのままマトリックス形式に書き換えたものであり，変位 未知数の式は，加速度未知数の式を，時刻t + ∆t時点での変位を未知数として解くように書き換えたもので ある． 1質点系や多質点系においては，加速度を未知数として解く式をよく見かけるが，有限要素法への適用にお いては，非線形を取り入れる場合に有利であるなどの理由から，変位を未知数とする場合も多いようである． また直感的には，m· a = f のイメージか，k· u = f のイメージかという気がするが，筆者の直感的感覚 は，有限要素法ではk· u = fであり，プログラムは変位を未知数として解くよう作成することにした．
1.3 行列およびベクトルの要素 平面骨組要素において，各種行列およびベクトルが以下のとおり定義されることは既知の事実とする． -6 x y -6 6 - Ni, ui Nj, uj Si, vi Sj, vj Mi, θi Mj, θj - 要素長さ ` 図1 要素概要図 -6 X Y U V * K x y u v φ-6 * K 図2 全体座標系(X,Y )と要素座標系(x,y)の関係図 {f} ：要素節点力ベクトル [k] ：要素剛性マトリックス [T ] ：座標変換マトリックス {u} ：要素節点変位ベクトル [m] ：要素質量マトリックス [k]T ：このT は転置を示す E ：要素弾性係数 γ ：要素単位体積重量 A ：要素断面積 g ：重力加速度 I ：要素断面２次モーメント φ ：要素座標系と全体座標系のなす角 ` ：要素長さ {f} ={Ni Si Mi Nj Sj Mj }T (12) {u} ={ui vi θi uj vj θj }T (13) [k] = EA/` 0 0 −EA/` 0 0
0 12EI/`3 6EI/`2 0 −12EI/`3 6EI/`2 0 6EI/`2 4EI/` 0 −6EI/`2 2EI/`
−EA/` 0 0 EA/` 0 0
0 −12EI/`3 −6EI/`2 0 12EI/`3 −6EI/`2
0 6EI/`2 2EI/` 0 −6EI/`2 4EI/` (14) [m] = γ· A · ` g 1/3 0 0 1/6 0 0 0 13/35 11`/210 0 9/70 −13`/420 0 11`/210 `2/105 0 13`/420 −`2/140 1/6 0 0 1/3 0 0 0 9/70 13`/420 0 13/35 −11`/210 0 −13`/420 −`2/140 0 −11`/210 `2/105 (15)
[T ] = cos φ sin φ 0 0 0 0 − sin φ cos φ 0 0 0 0 0 0 1 0 0 0 0 0 0 cos φ sin φ 0 0 0 0 − sin φ cos φ 0 0 0 0 0 0 1 (16) 座標変換マトリックス[T ]は，全体座標系変位を要素座標系変位に変換するものであり，全体座標系での剛 性マトリックス[K]と要素座標系での剛性マトリックス[k]には次の関係がある． [K] = [T ]T[k][T ] (17) 全体座標系での質量マトリックス[M ]と要素座標系での質量マトリックス[m]についても同様に次の関係 がある． [M ] = [T ]T[m][T ] (18) 減衰マトリックスについては，Rayleigh減衰を用いる場合，質量マトリックスと剛性マトリックスに，そ れぞれ定数を乗じたものの和となるため，以下のとおりとなる． [C] = ζm· [M] + ζk· [K] (19) ここに，[C]，[M ]，[K]はそれぞれ全体座標系での減衰マトリックス，質量マトリックス，剛性マトリック スであり，ζm，ζkは減衰を表す定数である．
1.4 減衰マトリックス(2012.08.19) 実際の計算で，構造物の減衰を考慮する場合，減衰マトリックス[c]を具体的に与える必要がある．ここで は，減衰マトリックスの与え方について考える． 地盤加速度を受ける1質点系の運動方程式について再掲する． m· ¨u + c · ˙u + k · u = −m · ¨φ ¨ u + c m· ˙u + k m· u = − ¨φ ¨ u + 2· h · ω · ˙u + ω2· u = − ¨φ h = c/m 2· ω ω = √ k m ここで，よく使われる減衰の形として，Rayleigh減衰を考える．これは，減衰係数は質量と剛性に比例す るという考え方であり，以下の式表現となる． c = ζm· m + ζk· k (20) 上式を考慮して減衰定数hを書き直してみると(ωは固有円振動数，f は振動数) h = c/m 2· ω= ζm· m + ζk· k 2· ω · m = ζm 2· ω+ ζk· ω 2 = 1 4πf· ζm+ πf · ζk よって，もし現地測定などで，i次モードとj次モードの減衰定数hi，hj および固有振動数fi，fjが判れ ば，上式を連立させることにより，比例定数ζmおよびζk を求めることができる． hi= 1 4πfi· ζm+ πfi· ζk hj = 1 4πfj· ζm+ πfj· ζk =⇒ ζm= 4π· fi· fj· (fj· hi− fi· hj) f2 j − fi2 ζk= fj· hj− fi· hi π(f2 j − f 2 i) 減衰係数をマトリックス表示すれば，上記で定めたスカラー量ζmおよびζkを用いて， [c] = ζm· [m] + ζk· [k] なお，実務上は以下のような手法が用いられる場合がある． 振動モードに応じた減衰定数がよくわからないので，減衰定数の値は1つとして，一般に使われている 値（例えばh1= h2= 0.05）などにセットする． 固有振動数については，固有値解析などで求めることができるため，1次および2次の固有振動数f1， f2を求める． 2つの固有振動数f1，f2と2つの減衰定数h1= h2= 0.05(実際には同一値)が定まれば，上式によ り，比例定数ζmおよびζk を求めることができる． babababababababababababababababababababab (注意)周期T (sec)，振動数f (Hz)，円振動数ω(rad/sec)の関係を混同しないこと！ T =2π ω = 1 f f = ω 2π
2.
有限要素法プログラムとしての入出力
独断と偏見が入る部分ではあるが，プログラム作成上は重要な部分なので，触れておこう． 2.1 固有値解析 固有値解析で得られる成果は，固有値と固有ベクトルであり，振動に関していえば，固有振動数と固有変形 モードである． 基本的にはこれらを出力すればいいだけであるが，現実問題としてはどこまで出力するか悩むところで ある． 例えば1節点3自由度の梁要素を用いて10要素11節点の片持梁の解析を行う場合，自由度などは以下の ようになる．（ここでは鉛直片持梁とする：y方向が長手方向） 1 全自由度11× 3 = 33 （1節点自由度と節点数の積） 2 固有方程式の次数11× 3 − 3 = 30 （マイナス3は基部で3自由度を拘束した分の自由度減分） 3 縦振動（軸方向振動）の固有値数10× 1 = 10 （基部を除いた10節点のy方向自由度の合計） 4 曲げ振動の固有値数10× 2 = 20 （基部を除いた10節点のx方向・回転方向自由度の合計) このため，片持梁の拘束条件である，基部のx・y・回転方向自由度の拘束のみを行えば，曲げ振動と縦振 動を含めた30個の固有値が得られる．したがって曲げ振動のみの固有値を求めたい場合は，全節点で軸方向 自由度を拘束する必要がある． なお全30個の固有値と30組の固有ベクトルを出力することは可能であるが，実務設計で全固有値および 固有ベクトルを必要とすることはまず無いし，ファイル出力するにしても全体の見通しが悪くなる． よってプログラムでは，小さいほうから並べていくつまで固有値を出力するかを指定できるように配慮 した． 出力する固有値の数および入出力ファイル名は，以下のように，組み込みサブルーチン：call getarg を用 い，コマンドラインより取得するようにしている． integer::negv !出力する固有値の数 character::fnameR*50,fnameW*50,fnameA*50 !入力ファイル・出力ファイル !コマンドラインからの入出力ファイル名読み込み call getarg(1,dummy) !出力する固有値の数 call getarg(2,fnameR) !構造入力ファイル call getarg(3,fnameW) !出力ファイルread(dummy,*) negv !文字列dummyを整数negvに変換
出力においては，出力する固有値の数がプログラム実行まで確定しない．このため，指定した固有値の数 (変数：negv)分，「コンマ＋数値」で構成される文字列を連結させて対応するようにしている．プログラミン
グ事例は以下のとおり．サブルーチンdel spaceは，コンマで連結された文字列から不要な空白を除去するサ
ブルーチンで，「Fortran入門」”http://www.nag-j.co.jp/fortran/index.html”に掲載されているものを使用 させていただいております．
!自由度番号 kk=mm if(negv<mm)kk=negv linebuf=’mode.’ do i=1,kk write(dummy,’(",",i0)’) i linebuf=trim(linebuf)//dummy end do call del_spaces(linebuf) write (12,’(a)’) trim(linebuf) !固有値 linebuf=’eigen value’ do i=1,kk write(dummy,’(",",e15.7)’) ev(i) linebuf=trim(linebuf)//dummy end do call del_spaces(linebuf) write (12,’(a)’) trim(linebuf)
図3 固有値出力時のプログラム事例 Impulse dt,2.00000000000000004E-003 ndata,10000 0.0000000000000000 0.0000000000000000 100.00000000000000 100.00000000000000 100.00000000000000 100.00000000000000 100.00000000000000 0.0000000000000000 0.0000000000000000 ・・・・・ 1行目 コメント 2行目 時間刻み 3行目 データ数 4行目以降 加速度値 図4 加速度入力ファイルの書式事例 2.2 時刻歴応答解析 (1) 出力 骨組構造の解析で得られるのは，基本的には各節点の変位・反力および要素の断面力である．時刻歴解析で は，各時刻で節点の加速度・速度変位および各要素の断面力が得られるため，これら全てを出力するとなると その量は膨大でる． そこで，何を出力するかを考えることになる．というか，本来は何が欲しいかでプログラムを書き換えれば いいのであるが，めんどうなのである程度必要と思われるものを出力できるようはじめから細工しておくほう が使い勝手が良い． ニーズから考えてみると，どの振動数で構造物が揺れやすいかを確認するためにはフーリエスペクトルを見 たい．そうなると，指定した節点あるいは指定した要素での断面力の時刻歴を出力できるようにしたい． また構造物の設計のことを考えると，変形図あるいは断面力図が欲しいため，指定した時刻での全節点変位 および全要素での断面力の出力が欲しい． そこで以下の考え方でプログラムを作成することにした． 指定した節点での加速度・速度・変位の全時刻での時刻歴が出力できる． 指定した要素の断面力の全時刻での時刻歴が出力できる． 指定した時刻での全節点変位および全要素断面力が出力できる． といってみたところで，指定した節点加速度等あるいは要素断面力の時刻歴についてはそのとおりプログラ ム化すればいいのだが，全節点変位および断面力を出力するタイミングが問題である．一番簡単なのは出力す る時刻を直接指定することである．規則的な波を入力する場合はそれでよいが，不規則波を入力する場合，指 定した時刻が必ずしも確認したい変位モードや応力状態になっているとは限らない．そこで，着目する節点と 自由度番号を入力し，その値が最大となった時の変位・断面力を記憶し，計算完了後に出力することにした． （実際にはプリント用の配列をつくり，最大値を検索しながら逐次書き換え，最終的に指定した自由度の変位 が最大となった時の値が残っているようプログラムする） (2) 入力 入力については問題を解くために必要な情報をいれるだけなのでそれほど問題はないが，時刻歴をどのよう に入力するかで少し考えた．
結局，熱伝導解析で，構造データと指定温度の時刻歴を別々に入力する方法をとっているので，同じ手法を 用いることにした． すなわち，入力データファイルを２つ準備する．１つは構造形状・材料特性などを入力するファイル，もう ひとつは加速度時刻歴を入力するファイルである． ２つの入力ファイルをもつことはある意味面倒であるが，同一構造で違う加速度を入力する場合，加速度時 刻歴ファイルを取り替えるだけなので，便利な場合も結構ある． 加速度入力書式は，フーリエスペクトルや応答スペクトル計算時の入力と同一とした．時刻歴の計算時間 は，加速度入力ファイルで指定するdt×ndataとなる．なお，時刻0における加速度・速度・変位・断面力は 0で固定している．
3.
固有値解析
(2012.08.20)
3.1 固有振動数と固有振動モードの算定方法の直感的解説 振動の分野では、固有値解析は，構造系の固有振動数と固有振動モードを求めるために利用される． 構造系の固有振動数と固有振動モードの算定方法を直感的に示そう． 1質点系の非減衰自由振動の運動方程式は以下のとおりとなる． m· ¨u + k · u = 0ここで，変位を正弦関数によりu(t) = A· sin(ω · t)とすれば，u =¨ −ω2· A · sin(ω · t)でるから，
(k− ω2· m) · u = 0 である．ここでu = 0，すなわちA = 0では振動現象は発生しないので，k− ω2· m = 0，すなわち ω = √ k m となり，1質点系の固有円振動数ωが求まった． これを梁のような連続体を解析するためのマトリックス表示にすれば， ([k]− ω2· [m]) · {u} = {0} ここで，左辺が行列であるため，振動が生じる条件で上式が成立するためには，{u} 6= 0であり，固有値問題 として，固有値ω2と固有ベクトル{u}を求めればよい． 上記を解くためには一般化ヤコビ法などが適用可能であり，固有値としての構造系の固有振動数(実際の固 有方程式の解は固有円振動数の2乗：ω2)，および固有ベクトルとして固有振動モードを求めることができる． なお，一般化ヤコビ法は，実対称行列[A]，[B]，固有値λとして，以下の固有方程式を解く数値解法である． ([A]− λ · [B]) · {x} = {0} 3.2 固有振動数の数値解析 片持梁および単純梁をモデルとして，一般化ヤコビ法を用いた有限要素法による固有値解析を行った． (1) モデル諸元 解析に用いるモデル諸元は以下のとおりとした． 弾性係数・材料の単位体積重量は，コンクリート相当とする． 片持梁は，鉛直方向を長手方向とし，長さは10mとする． 片持梁の断面は，幅0.5m，高さ0.5mの矩形断面とする． 単純梁は，鉛直方向を長手方向とし，長さは10mとする． 単純梁の断面は，幅0.3m，高さ0.3mの矩形断面とする． 片持梁モデルの諸元 弾性係数 E (t/m2) 2,000,000 断面積 A (m2) 0.25 断面二次モーメント I (m4) 0.00521 材料の単位体積重量 γ (tf/m3) 2.35 部材長さ ` (m) 10.0 解析モデル 梁要素FEM，10要素・11節点モデル 支持条件 下部固定端，上部無拘束(鉛直片持梁)
単純梁モデルの諸元 弾性係数 E (t/m2) 2,000,000 断面積 A (m2) 0.09 断面二次モーメント I (m4) 0.000675 材料の単位体積重量 γ (tf/m3) 2.35 部材長さ ` (m) 10.0 解析モデル 梁要素FEM，10要素・11節点モデル 支持条件 下部水平・鉛直変位拘束，上部水平変位拘束(鉛直単純梁) (2) 固有振動数の理論値 単純梁の曲げ振動の固有振動数は，次式で与えられる． fn = 1 2π ( nπ ` )2 · √ EI ρA 片持梁の曲げ振動の固有振動数は，次式で与えられる． fn= 1 2π ( α ` )2 · √ EI ρA αは右式を満たす値： cos(α)· cosh(α) + 1 = 0 また一端固定-一端自由の棒の縦振動の固有振動数は，次式で与えられる． fn = 2n− 1 4` √ E ρ 上式よりわかるように棒の長さ・弾性係数・材料密度が同じであれば同一値となる． (3) 固有値解析結果 固有値解析結果と理論固有値の比較結果は以下のとおりである．用いたモデルが10要素モデルと自由度が 小さいため，高次の振動数になるほど解析値は理論値から大きい側に乖離する．またこの解析結果によれば， 曲げ振動に比べ，縦振動の精度が悪いようである．これは縦振動の固有振動数が大きいことにも影響されてい る可能性がある． 表1 固有振動数の理論値と固有値解析結果の比較 片持梁の曲げ振動の固有振動数 モード 1 2 3 4 5 6 7 8 固有値解析 2.333 14.621 40.949 80.233 132.948 199.171 279.447 374.381 理論値 2.333 14.621 40.938 80.223 132.614 198.102 276.688 368.371 比 1.000 1.000 1.000 1.000 1.003 1.005 1.010 1.016 単純梁の曲げ振動の固有振動数 モード 1 2 3 4 5 6 7 8 固有値解析 3.929 15.716 35.377 62.963 98.604 142.555 195.234 257.227 理論値 3.929 15.715 35.358 62.859 98.217 141.432 192.594 251.434 比 1.000 1.000 1.001 1.002 1.004 1.008 1.014 1.023 棒の縦振動の固有振動数（片持梁・単純梁とも解析結果は同一値） モード 1 2 3 4 5 6 7 8 固有値解析 72.274 218.608 370.333 531.072 704.176 891.705 1091.854 1293.716 理論値 72.199 216.598 360.997 505.396 649.795 794.194 938.593 1082.992 比 1.001 1.009 1.026 1.051 1.084 1.123 1.163 1.195
babababababababababababababababababababab (参考事項) 固有値の理論値は以下のawkスクリプトで計算している．片持梁の固有振動数は，2分法で非線形 方程式を解いている．また解の範囲はcos関数の性質を利用して検索範囲を設定している． BEGIN{ pi=3.14159265358979323846 E=10000 I=1 m=1 al=10 for(iii=1;iii<=10;iii++){ x1=(iii-1)*pi x2=x1+pi do{ xi=0.5*(x1+x2) f1=cos(x1)*0.5*(exp(x1)+exp(-x1))+1 f2=cos(x2)*0.5*(exp(x2)+exp(-x2))+1 fi=cos(xi)*0.5*(exp(xi)+exp(-xi))+1 if(f1*fi<0)x2=xi if(f2*fi<0)x1=xi if(fi==0)exit if(f1*fi<0)eps=xi-x1 if(f2*fi<0)eps=x2-xi }while(1e-6<eps) f1=(iii*pi/al)^2*sqrt(E*I/m)/2/pi f2=(xi/al)^2*sqrt(E*I/m)/2/pi printf "%.3f,%.3f\n",f1,f2 } }
4.
時刻歴応答解析
(2012.08.23-2012.08.27)
4.1 非減衰振動 固有値解析に用いたものと同一のモデルを用い，片持梁と単純梁の固有振動数を，衝撃加速度に対する応答 から求めることを試みた． 振動方程式の数値解法としては，線形加速度法，平均加速度法が用いられることが多いと思われるが，それ ぞれ以下の特性がある． 線形加速度法は，解の精度は良いが，計算時間間隔が大きいと解が収束しない． 平均加速度法は，解の収束性は保証されるが，計算時間間隔が大きいと位相の遅れが顕著になる． そこで，ここでは収束性が保証される平均加速度法を用いることし，計算時間間隔を変化させ，その影響を 確認することにした． 具体的な固有振動数算定手順は以下のとおりである． a. 衝撃荷重を模擬して，短時間に一定加速度を構造物の支点に作用させる． b. 作用荷重に対する，予め選定した節点の変形あるいは断面力の時刻歴波形を追跡する． c. 得られた時刻歴波形のフーリエスペクトルから卓越振動数を把握し，これを系の固有値とみなす． (1) 片持梁の非減衰時刻歴応答解析時刻歴応答解析により鉛直片持梁の固有振動数を確認すべく，時刻0.01secに，0.01secの時間長さで100gal の水平加速度を下端部に与えた．鉛直片持梁の挙動特性は，天端自由端の水平変位と基部固定端モーメントに 特徴的に現れると考え，これらの数値の時刻歴を追跡した．構造物の減衰特性については考慮していない． 図6に時間刻み∆t = 0.01secの場合の支点の加速度・速度・変位および断面力のフーリエスペクトルを 示す． 図7・8に時間刻み∆tを変化させた場合の支点の加速度および断面力のフーリエスペクトルを示す． 図9に時間刻み∆tを変化させた場合の支点断面力の時刻歴波形を示す． これらより，時間刻み∆t = 0.01secでは，1次モードの周期以外は固有値解析で求めたモデルよりかなり 小さい卓越周期が計算されている．すなわち，平均加速度法の短所である位相遅れが顕著に現れているとみな すことができる．これに対し，時間刻み∆t = 0.01を小さくすることにより，卓越振動数は，固有値解析で得 られた系の固有値に近づくことが判る． (2) 単純梁の非減衰時刻歴応答解析
時刻歴応答解析により鉛直単純梁の固有振動数を確認すべく，時刻0.01secに，0.01secの時間長さで100gal の水平加速度を上端および下端支点に与えた．鉛直単純梁の挙動特性は，支点のたわみ角と支点のせん断力に 特徴的に現れると考え，これらの数値の時刻歴を追跡した．構造物の減衰特性については考慮していない． 図10に時間刻み∆t = 0.01secの場合の支点の加速度・速度・変位および断面力のフーリエスペクトルを 示す． 図11・12に時間刻み∆tを変化させた場合の支点の加速度および断面力のフーリエスペクトルを示す． 図13に時間刻み∆tを変化させた場合の支点断面力の時刻歴波形を示す． 単純梁においても，片持梁と同様，時間刻み∆t = 0.01secでは，系の卓越振動数と固有値解析で求めた固 有振動数の値は大きく異なるが，時間刻み時間刻み∆tを小さくすることにより，系の卓越振動数は，固有値 解析結果による固有振動数に近づくことが確認された． なお，片持梁では，1次モードから7次モード程度までのモードが比較的明確にフーリエスペクトルから読 み取れたが，単純梁では，対称変形モードとなる奇数モードしかフーリエスペクトルからは読み取れなかっ た．これは加速度の載荷方向に起因するものであり，上下の支点に逆方向の加速度を載荷すれば逆対称変形と
なる偶数モードの固有振動数が得られることは確認済みである． (3) 時間刻み変更時の加速度の入力方法 数値積分の時間刻みを変更した場合，荷重として入力する加速度をどのように入力すればいいのだろうか？ 入力は，定められた時間刻み毎に，その時々の時刻の加速度値を入力するようにしている． ここでは、時間刻みを小さくした場合の加速度値は以下の図のように入力している． これは，振動方程式の数値積分は，次の時刻の荷重が読み込まれたとき，その時刻の応答が決定されること を考えると，次の荷重が読み込まれるまでは前の値を保持しているとみなすのが妥当であるという考え方によ り設定したものである． この入力の結果として，図9・13に時間刻み∆tを変化させた場合の断面力が得られているが，期待通りの 出力となっている． （告白：下の図はgnuplotで作ってしまいました） 0 20 40 60 80 100 120 0.00 0.01 0.02 0.03 0.04 0.05 Acc. (gal) Time (sec) ∆t=0.01s 0 20 40 60 80 100 120 0.00 0.01 0.02 0.03 0.04 0.05 Acc. (gal) Time (sec) ∆t=0.005s 0 20 40 60 80 100 120 0.00 0.01 0.02 0.03 0.04 0.05 Acc. (gal) Time (sec) ∆t=0.002s 0 20 40 60 80 100 120 0.00 0.01 0.02 0.03 0.04 0.05 Acc. (gal) Time (sec) ∆t=0.001s 図5 加速度の入力方法
天端 水平加速度 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 13.70Hz ⇐ 28.96Hz ⇐ 37.99Hz ⇐ 2.32Hz ⇐ 42.53Hz ⇐ 44.95Hz ⇐ 46.39Hz Acceleration 天端 水平速度 0.0 0.1 0.2 0.3 0.4 0.5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0.0 0.1 0.2 0.3 0.4 0.5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 2.32Hz ⇐ 13.70Hz ⇐ 28.96Hz ⇐ 37.99Hz ⇐ 42.53Hz Velocity 天端 水平変位 0.000 0.005 0.010 0.015 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0.000 0.005 0.010 0.015 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 2.32Hz ⇐ 13.70Hz ⇐ 28.96Hz Displacement 基部 モーメント 0 1 2 3 4 5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0 1 2 3 4 5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 2.32Hz ⇐ 13.70Hz ⇐ 28.96Hz ⇐ 37.99Hz ⇐ 42.53Hz Moment 図6 片持梁の応答振動数特性（衝撃加速度100gal 入力，計算時間間隔 dt=0.01sec）
計算時間間隔 dt=0.01 sec 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 13.70Hz ⇐ 28.96Hz ⇐ 37.99Hz ⇐ 2.32Hz ⇐ 42.53Hz ⇐ 44.95Hz ⇐ 46.39Hz Acceleration 計算時間間隔 dt=0.005 sec 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 36.38Hz ⇐ 14.38Hz ⇐ 2.34Hz Acceleration 計算時間間隔 dt=0.002 sec 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 40.08Hz ⇐ 14.59Hz ⇐ 2.33Hz Acceleration 計算時間間隔 dt=0.001 sec 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 40.73Hz ⇐ 14.62Hz ⇐ 2.33Hz Acceleration 図7 片持梁応答振動数における時間刻みの影響（天端加速度，衝撃荷重100 gal， 平均加速度法）
計算時間間隔 dt=0.01 sec 0 1 2 3 4 5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0 1 2 3 4 5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 2.32Hz ⇐ 13.70Hz ⇐ 28.96Hz ⇐ 37.99Hz ⇐ 42.53Hz Moment 計算時間間隔 dt=0.005 sec 0 1 2 3 4 5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0 1 2 3 4 5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 2.34Hz ⇐ 14.38Hz ⇐ 36.38Hz Moment 計算時間間隔 dt=0.002 sec 0 1 2 3 4 5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0 1 2 3 4 5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 2.33Hz ⇐ 14.59Hz ⇐ 40.08Hz Moment 計算時間間隔 dt=0.001 sec 0 1 2 3 4 5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0 1 2 3 4 5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 2.33Hz ⇐ 14.62Hz ⇐ 40.73Hz Moment 図8 片持梁応答振動数における時間刻みの影響（基部モーメント，衝撃荷重 100 gal， 平均加速度法）
計算時間間隔 dt=0.01 sec −0.8 0.0 0.8 Moment (t−m) 0.0 0.5 1.0 1.5 2.0 time (sec) −0.8 0.0 0.8 Moment (t−m) 0.0 0.5 1.0 1.5 2.0 time (sec) 計算時間間隔 dt=0.005 sec −0.8 0.0 0.8 Moment (t−m) 0.0 0.5 1.0 1.5 2.0 time (sec) −0.8 0.0 0.8 Moment (t−m) 0.0 0.5 1.0 1.5 2.0 time (sec) 計算時間間隔 dt=0.002 sec −0.8 0.0 0.8 Moment (t−m) 0.0 0.5 1.0 1.5 2.0 time (sec) −0.8 0.0 0.8 Moment (t−m) 0.0 0.5 1.0 1.5 2.0 time (sec) 計算時間間隔 dt=0.001 sec −0.8 0.0 0.8 Moment (t−m) 0.0 0.5 1.0 1.5 2.0 time (sec) −0.8 0.0 0.8 Moment (t−m) 0.0 0.5 1.0 1.5 2.0 time (sec) 図9 片持梁基部モーメント時刻歴波形の比較（衝撃荷重 100 gal， 平均加速度法）
下端 回転加速度 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 26.68Hz ⇐ 44.85Hz ⇐ 48.05Hz ⇐ 3.91Hz ⇐ 40.06Hz ⇐ 49.07Hz Acceleration 天端 回転速度 0.00 0.02 0.04 0.06 0.08 0.10 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0.00 0.02 0.04 0.06 0.08 0.10 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 3.91Hz ⇐ 26.68Hz ⇐ 44.85Hz ⇐ 40.06Hz ⇐ 48.05Hz Velocity 下端 回転変位 0.000 0.001 0.002 0.003 0.004 0.005 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0.000 0.001 0.002 0.003 0.004 0.005 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 3.91Hz ⇐ 26.68Hz Displacement 下端 せん断力 0.0 0.1 0.2 0.3 0.4 0.5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0.0 0.1 0.2 0.3 0.4 0.5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 3.91Hz ⇐ 26.68Hz ⇐ 40.06Hz ⇐ 44.85Hz Shear force 図10 単純梁の応答振動数特性（衝撃加速度100gal 入力，dt=0.01sec）
dt=0.01 sec 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 26.68Hz ⇐ 44.85Hz ⇐ 48.05Hz ⇐ 3.91Hz ⇐ 40.06Hz ⇐ 49.07Hz Acceleration dt=0.005 sec 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 32.30Hz ⇐ 3.93Hz Acceleration dt=0.002 sec 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 34.81Hz ⇐ 3.92Hz Acceleration dt=0.001 sec 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0 5 10 15 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 35.23Hz ⇐ 3.92Hz Acceleration 図11 単純梁応答振動数における時間刻みの影響（下端回転加速度，衝撃荷重 100 gal， 平均加速度法）
dt=0.01 sec 0.0 0.1 0.2 0.3 0.4 0.5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0.0 0.1 0.2 0.3 0.4 0.5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 3.91Hz ⇐ 26.68Hz ⇐ 40.06Hz ⇐ 44.85Hz Shear force dt=0.005 sec 0.0 0.1 0.2 0.3 0.4 0.5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0.0 0.1 0.2 0.3 0.4 0.5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 3.93Hz ⇐ 32.30Hz Shear force dt=0.002 sec 0.0 0.1 0.2 0.3 0.4 0.5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0.0 0.1 0.2 0.3 0.4 0.5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 34.81Hz ⇐ 3.92Hz Shear force dt=0.001 sec 0.0 0.1 0.2 0.3 0.4 0.5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) 0.0 0.1 0.2 0.3 0.4 0.5 Fourier spectrum 0 10 20 30 40 50 Frequency (Hz) ⇐ 35.23Hz ⇐ 3.92Hz Shear force 図12 単純梁応答振動数における時間刻みの影響（下端せん断力，衝撃荷重 100 gal， 平均加速度法）
dt=0.01 sec −0.05 0.00 0.05 Shear force (t) 0.0 0.5 1.0 1.5 2.0 time (sec) −0.05 0.00 0.05 Shear force (t) 0.0 0.5 1.0 1.5 2.0 time (sec) dt=0.005 sec −0.05 0.00 0.05 Shear force (t) 0.0 0.5 1.0 1.5 2.0 time (sec) −0.05 0.00 0.05 Shear force (t) 0.0 0.5 1.0 1.5 2.0 time (sec) dt=0.002 sec −0.05 0.00 0.05 Shear force (t) 0.0 0.5 1.0 1.5 2.0 time (sec) −0.05 0.00 0.05 Shear force (t) 0.0 0.5 1.0 1.5 2.0 time (sec) dt=0.001 sec −0.05 0.00 0.05 Shear force (t) 0.0 0.5 1.0 1.5 2.0 time (sec) −0.05 0.00 0.05 Shear force (t) 0.0 0.5 1.0 1.5 2.0 time (sec) 図13 単純梁下端せん断力時刻歴波形の比較（衝撃荷重 100 gal， 平均加速度法）
(4) 時間刻みが固有振動数計算精度に及ぼす影響 これまでの解析結果に基づき，平均加速度法における時間刻みが固有振動数計算精度に及ぼす影響を整理し た．固有値解析結果による固有振動数は，連続体としての固有値とは異なるが，有限要素モデルがもつ固有の 値であるため，これを基準としている．時刻歴解析結果による固有振動数は，フーリエスペクトルより求めた 卓越振動数を固有値とみなしている． 表2 平均加速度法における時間刻みが固有振動数計算に及ぼす影響の整理 片持梁 1 固有値解析結果 2.33 14.62 40.95 80.23 132.95 199.17 279.45 2 fn(∆t=0.01s) 2.33 13.70 28.96 37.99 42.53 44.95 46.39 時刻歴解析による fn(∆t=0.005s) 2.34 14.38 36.38 57.32 71.56 80.30 85.74 固有振動数(Hz) fn(∆t=0.002s) 2.33 14.59 40.08 74.37 110.76 142.70 167.60 fn(∆t=0.001s) 2.33 14.62 40.73 76.66 125.93 177.98 229.34 ∆t=0.01s 0.023 0.137 0.290 0.380 0.425 0.450 0.464 無次元時間刻み ∆t=0.005s 0.012 0.072 0.182 0.287 0.358 0.402 0.429 fn∆t = ∆t/Tn ∆t=0.002s 0.005 0.029 0.080 0.149 0.222 0.285 0.335 ∆t=0.001s 0.002 0.015 0.041 0.077 0.126 0.178 0.229 ∆t=0.01s 1.000 0.937 0.707 0.474 0.320 0.226 0.166 固有振動数比 ∆t=0.005s 1.004 0.984 0.888 0.714 0.538 0.403 0.307 ( 2 / 1 ) ∆t=0.002s 1.000 0.998 0.979 0.927 0.833 0.716 0.600 ∆t=0.001s 1.000 1.000 0.995 0.956 0.947 0.894 0.821 単純梁 1 固有値解析結果 3.93 35.38 98.60 195.23 328.99 516.37 – 2 fn(∆t=0.01s) 3.91 26.68 40.06 44.85 – – – 時刻歴解析による fn(∆t=0.005s) 3.93 32.30 63.50 79.93 – – – 固有振動数(Hz) fn(∆t=0.002s) 3.92 34.81 88.27 141.14 – – – fn(∆t=0.001s) 3.92 35.23 95.63 175.13 255.25 324.16 – ∆t=0.01s 0.039 0.267 0.401 0.449 – – – 無次元時間刻み ∆t=0.005s 0.020 0.162 0.318 0.400 – – – fn∆t = ∆t/Tn ∆t=0.002s 0.008 0.070 0.177 0.282 – – – ∆t=0.001s 0.004 0.035 0.096 0.175 0.255 0.324 – ∆t=0.01s 0.995 0.754 0.406 0.230 – – – 固有振動数比 ∆t=0.005s 1.000 0.913 0.644 0.409 – – – ( 2 / 1 ) ∆t=0.002s 0.997 0.984 0.895 0.723 – – – ∆t=0.001s 0.997 0.996 0.970 0.897 0.776 0.628 – 図14に，平均加速度法における時間刻みと固有振動数の精度（時刻歴解析による固有振動数と固有値解析 による固有振動数の比）の関係を示す．図より，獲得したい固有周期Tnに対する時間刻み∆tの比が大きく なると固有振動数は小さく算定される（固有周期は大きく算定される）ことは明らかであり，固有振動数を基 準にして5%程度以内の精度の時刻歴を作成したい場合は，獲得したい最小の固有周期に対し1/10程度以下 の時間刻みを採用する必要があることがわかる． （告白：図14はgnuplotで作ってしまいました）
0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.0 0.1 0.2 0.3 0.4 0.5
Ratio of EV by THA to EV by EVA
Ratio of time step to natural period ∆t/Tn Cantilever Simple beam
図14 無次元時間刻み(時間刻み/固有周期)と固有振動数比(時刻歴解析/固有値解析)の関係 (EV：固有値，EVA：固有値解析，THA：時刻歴解析)
4.2 強制振動 強制振動による応答の計算事例として，時間刻みと加速度載荷方法に関するテストを行った． 構造物の持つ固有振動数と同じ振動数の外力を受ければ構造物は共振を起こすはずである． ここでは，片持梁を用いて，片持梁の3次固有振動数（40.73Hz：dt=0.001secの時刻歴解析結果による） と同じ振動数を持つ正弦波の水平加速度を支点に載荷することを考える．支点に入力する加速度は，振動数： 40.73Hz，加速度片振幅：10gal，波数：20波とし，時刻歴の計算方法は平均加速度法とする． 図15に加速度載荷方法の詳細を示す．本来，計算時間刻みdt=0.001sec単位で具体的な加速度値がわかっ ていれば，Originalとして示すように加速度を入力すればよい．しかし，計算時間刻みは計算精度を考え dt=0.001secとしたいが，何かしらの不都合で加速度値が0.01secピッチでしか知ることができなかったこと にしよう．この場合，以下のケースが考えられる． Case-1 ：計算時間刻みを細かくすることをあきらめ，既知加速度値と同じピッチである0.01secで加速度を 入力し計算も行う．（そもそも約40hzの波形をこの入力加速度で再現しようとすることに無理がある ことは直感的に納得できる）
Case-2 ：計算時間刻みはdt=0.001secとし，加速度は既知加速度値をdt=0.001secピッチでstep状にモデ ル化して入力する．
Case-3 ：計算時間刻みはdt=0.001secとし，加速度は既知加速度値をdt=0.001secピッチで線形内挿して 入力する． Original 計算時間間隔 dt=0.001 sec 0.001secピッチのsin波 −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 time (sec) −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 time (sec) −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 time (sec) Original (dt=0.001s) Case-1 計算時間間隔 dt=0.01 sec 0.01secピッチのsin波 −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 time (sec) −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 time (sec) −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 time (sec) Case−1 (dt=0.01s) Case-2 計算時間間隔 dt=0.001 sec Case-2のデータを step内挿 −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 time (sec) −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 time (sec) −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 time (sec) Case−2 (dt=0.001s) Case-3 計算時間間隔 dt=0.001 sec Case-2のデータを 線形内挿 −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 time (sec) −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 time (sec) −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 time (sec) Case−3 (dt=0.001s) 図15 加速度荷重詳細（初期 0.1sec，振動数 40.73Hz，加速度片振幅 10gal， 平均加速度法）
各入力ケースに対する片持梁の基部モーメントの応答波形を入力加速度波形とともに図16に示す． 図より以下のことがわかる． Case-1では時間刻みが粗いため挙動を追跡しきれず，共振現象は再現できない． Case-2では時間刻みを細かくしたため，加速度載荷中は基部モーメント値の振幅は増加していき，載 荷が完了した時点から概ね一定振幅の定常振動となる． Case-3においてもCase-2と同様，加速度載荷中は基部モーメント値の振幅は増加していき，載荷が完 了した時点から概ね一定振幅の定常振動となる．
Case-2とCase-3を比べると，応答値（モーメントの値）の振幅はCase-2のほうが大きく，Original
に近いことが判る．
以上より，このテストで得られた結論は以下のとおりとなる．
前述したとおり，計算時間刻みは挙動を再現したい固有周期の1/10程度以下にする必要がある．
もし荷重に設定した計算時間刻みの時間的精度がない場合は，へたに線形内挿するより，既知の値を
Original 計算時間間隔 dt=0.001 sec −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) Original (dt=0.001s) −0.4 −0.2 0.0 0.2 0.4 Moment (t−m) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) −0.4 −0.2 0.0 0.2 0.4 Moment (t−m) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) −0.4 −0.2 0.0 0.2 0.4 Moment (t−m) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) Original (dt=0.001s) Case-1 計算時間間隔 dt=0.01 sec −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) Case−1 (dt=0.01s) −0.4 −0.2 0.0 0.2 0.4 Moment (t−m) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) −0.4 −0.2 0.0 0.2 0.4 Moment (t−m) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) −0.4 −0.2 0.0 0.2 0.4 Moment (t−m) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) −0.4 −0.2 0.0 0.2 0.4 Moment (t−m) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) Case−1 (dt=0.01s) Case-2 計算時間間隔 dt=0.001 sec −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) Case−2 (dt=0.001s) −0.4 −0.2 0.0 0.2 0.4 Moment (t−m) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) −0.4 −0.2 0.0 0.2 0.4 Moment (t−m) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) −0.4 −0.2 0.0 0.2 0.4 Moment (t−m) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) Case−2 (dt=0.001s) Case-3 計算時間間隔 dt=0.001 sec −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 Acceleration (m/s) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) Case−3 (dt=0.001s) −0.4 −0.2 0.0 0.2 0.4 Moment (t−m) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) −0.4 −0.2 0.0 0.2 0.4 Moment (t−m) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) −0.4 −0.2 0.0 0.2 0.4 Moment (t−m) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time (sec) Case−3 (dt=0.001s) 図16 片持梁の共振数値実験（振動数 40.73Hz，加速度片振幅 10gal， 平均加速度法）
4.3 減衰振動 (1) モデル概要 片持梁モデルを用いて減衰振動における挙動を確認する．モデルはこれまで用いてきた鉛直片持梁モデル とし，荷重として，片振幅100galの水平加速度を下端固定部に0.01秒間作用させ，時刻歴解析を行った．計 算の時間刻みは，dt=0.001sec，計算時間は20秒間とした．着目点は，下端固定部モーメントとした． 減衰係数マトリックス[c]を支配するパラメータは以下の式に示すζmおよびζkである． [c] = ζm[m] + ζk[k] ζmおよびζkは，減衰定数h，振動数fを用いて以下の関係を有している． h = 1 4πf · ζm+ πf· ζk ここでは，固有値解析および減衰なしの場合の時刻歴解析で得られている片持梁の固有振動数に対し，以下 の係数を設定した． Case 固有振動数f 減衰定数h 係数ζm 係数ζm 1 (質量比例項のみ設定) 2.33Hz(1次モード) 0.05 1.4640 0 2 (剛性比例項のみ設定) 14.62Hz(2次モード) 0.02 0 0.0004354 (2) 解析結果 得られた下端固定部モ−メントの時刻歴を，減衰無しの場合を含めて，図17に示す．Case-1の減衰係数 は質量比例型であり，低次振動数の波の減衰が早く，高次振動数の波の減衰は遅いはずである．またCase-2 の減衰係数は剛性比例型であり，Case-1とは逆に，低次振動数の波の減衰は遅く，高次振動数の波の減衰は 早いはずである．Case-1・2とも上記特性を有していることは波形より認識することができる．しかし，定量 的に指定した減衰が放蕩に再現できているのであろうか？ 減衰無し −0.8 0.0 0.8 Moment (t−m) 0 1 2 3 4 time (sec) −0.8 0.0 0.8 Moment (t−m) 0 1 2 3 4 time (sec) No damping Case-1 ζm= 1.4649 ζk= 0 (f=2.33Hz, h=0.05) −0.8 0.0 0.8 Moment (t−m) 0 1 2 3 4 time (sec) −0.8 0.0 0.8 Moment (t−m) 0 1 2 3 4 time (sec) Zm=1.4640 Case-2 ζm= 0 ζk= 0.0004354 (f=14.62Hz, h=0.02) −0.8 0.0 0.8 Moment (t−m) 0 1 2 3 4 time (sec) −0.8 0.0 0.8 Moment (t−m) 0 1 2 3 4 time (sec) Zk=0.0004354
(3) 解析結果の確認
計算結果は指定した減衰を表現しているかの確認を行いたい．そこで，Simplex法を用いた回帰計算により
各モードの回帰式を求め，固有振動数および減衰定数の再現性の確認を行うことにした． 一般に減衰振動は次式で表される．
u(t) = C· exp(−h · 2πf · t) · sin(2πf · t + δ)
そこで，Simplex法により，以下の4つの未知パラメータp[i]を求めた．
u(t) = p[0]· exp(−p[1] · 2π · p[2] · t) · sin(2π · p[2] · t + p[3])
Simplex法による回帰結果を下表に，この結果に基づき，1・2次モードの振動波形を除去した時刻歴を，原
波形を含めて，図18に示す．
1次モード回帰結果
Case データ数 C h f δ
p[0] p[1] p[2] p[3]
Case-1 20000 -0.3942434E+00 0.4998969E-01 0.2323818E+01 -0.1581718E+00 Case-2 20000 -0.3912423E+00 0.3190515E-02 0.2332857E+01 -0.2065630E+00
2次モード回帰結果（回帰による1次モード除去後）
Case データ数 C h f δ
p[0] p[1] p[2] p[3]
Case-1 20000 0.2038415E+00 0.7693699E-02 0.1460409E+02 0.1865452E+01 Case-2 20000 0.2017669E+00 0.1894490E-01 0.1456833E+02 0.1951257E+01
3次モード回帰結果（回帰による2次モード除去後）
Case データ数 C h f δ
p[0] p[1] p[2] p[3]
Case-1 20000 0.9009480E-01 0.2681717E-02 0.4072046E+02 -0.5314064E+00 Case-2 20000 0.6462909E-01 0.3773392E-01 0.3938786E+02 0.1865192E+00
Case-1 ζm= 1.4649 ζk= 0 原波形 −0.8 0.0 0.8 Moment (t−m) 0 1 2 3 4 time (sec) −0.8 0.0 0.8 Moment (t−m) 0 1 2 3 4 time (sec) Zm=1.4640 Case-1 ζm= 1.4649 ζk= 0 1次モード除去波形 −0.8 0.0 0.8 Moment (t−m) 0 1 2 3 4 time (sec) −0.8 0.0 0.8 Moment (t−m) 0 1 2 3 4 time (sec) Zm=1.4640 Case-1 ζm= 1.4649 ζk= 0 1・2次モード除去波形 −0.8 0.0 0.8 Moment (t−m) 0 1 2 3 4 time (sec) −0.8 0.0 0.8 Moment (t−m) 0 1 2 3 4 time (sec) Zm=1.4640 Case-2 ζm= 0 ζk= 0.0004354 原波形 −0.8 0.0 0.8 Moment (t−m) 0 1 2 3 4 time (sec) −0.8 0.0 0.8 Moment (t−m) 0 1 2 3 4 time (sec) Zk=0.0004354 Case-2 ζm= 0 ζk= 0.0004354 1次モード除去波形 −0.8 0.0 0.8 Moment (t−m) 0 1 2 3 4 time (sec) −0.8 0.0 0.8 Moment (t−m) 0 1 2 3 4 time (sec) Zk=0.0004354 Case-2 ζm= 0 ζk= 0.0004354 1・2次モード除去波形 −0.8 0.0 0.8 Moment (t−m) 0 1 2 3 4 time (sec) −0.8 0.0 0.8 Moment (t−m) 0 1 2 3 4 time (sec) Zk=0.0004354
(4) 減衰振動のまとめ 減衰定数hは，比例定数ζm，ζkを用いて，振動数f の関数として以下のとおり現される． h = 1 4πf· ζm · · · for Case-1 h = πf· ζk · · · for Case-2 設定した比例定数に対し，時刻歴応答解析による挙動より回帰計算を用いて算定した固有振動数および 減衰定数は下表および図19のとおりである． 固有振動数が大きくなると回帰計算による算定値は理論値との乖離が大きくなる傾向にある． 特にCase-2のモード3については，図18からも判るように，1・2次モード除去後は高振動数波形の 減衰が大きく，きれいな減衰波形が残らず，減衰振動式への回帰の精度が著しく低くなったため，理論 値と計算値との乖離が大きくなったものと考えられる．
Case Case-1 Case-2
ζ ζm=1.4640 ζk=0.0004354 モード 1 2 3 1 2 3 回帰結果固有振動数(Hz) 2.324 14.604 40.720 2.333 14.568 39.388 回帰結果減衰定数 0.0500 0.0077 0.0027 0.0032 0.0189 0.0377 理論減衰定数 0.0501 0.0080 0.0029 0.0032 0.0199 0.0539 回帰結果と理論値の比 0.998 0.963 0.931 1.000 0.950 0.699 0.001 0.01 0.1 1 0 10 20 30 40 50 Damping factor h Frequency f (Hz) Case-1 Theory Case-1 Analysis Case-2 Theory Case-2 Analysis 図19 振動数と減衰定数の関係（理論値と計算値の比較） (告白：上図はgnuplotで作ってしまいました)
