ほぼ線形の動力学
Winnipeg にある歩道橋 Esplanade Riel 橋6
振動論の基礎
9.1
1 自由度系の振動
9.1.1 非減衰自由振動 (1) 運動方程式 この章は東京大学の伊藤学先生の 1970 年頃の講義ノートを元にした。当時の教科書は文献 [131] だったが， 多分元になったのは文献 [115] の方だと思われる。まず 1 自由度系の振動問題を用いて，振動論において最も 重要な固有振動数と共振という概念を説明する。多自由度系の場合にはベクトルとしての固有振動モードとい う概念を説明する。そして梁を例とした連続体の振動の場合には，その固有振動モードが関数になることを説 明する。最終的には有限要素法を用いることによって，連続体は見かけ上は多自由度系として近似的に取り扱 うことができることを説明する。しかし，共振や周波数応答関数の特性，さらに境界条件によって振動数やモー ドが異なること等が，数値解析をするまでもなく予想できるような技術者になって欲しい。数値解析手法は便 利なツールではあるが，技術者の力学的直感（と言い切ってはまずいか。正確な力学的知識と経験ということ になるかな）があってこそ有効に使えるものだ。なお，この文書では弾性範囲のみを対象とするので，強度設 計のためにはさらなる学習が必要である。またこの章でも簡単のために括弧無しの太字で行列を表している。 w(t) k m u(t) mg k f (t) mg kw(t) 図 9.1 1 自由度系の振動 さて，バネ定数 k の線形バネに取り付けた 1 個の質点の運動 を考えよう。バネの質量は無視できるものとする。図 9.1 のよ うに自然状態のバネの先端に質量 m の質点を付けたところ運 動し始めた。そのときの質点の下方への変位（バネの伸び）を w(t) とし，質点には外力 f (t) も作用することができるものとす る。このとき，質点の運動方向（図では下方）に作用している 力の総和は，図の右端にある力のダイアグラムから ∑ 【下向きの力】= f (t) + mg − k w(t) である。したがって Newton の運動方程式 (2.1a)1は ∑ 【下向きの力】= f (t) + mg − k w(t) = md 2w(t) dt2 → m ¨w(t) + k w(t) = mg + f (t) (9.1) 1米国のある雑誌（誌名は失念）の一コマ漫画。 Einstein のような人がチョークを持って黒板に向かって立っていて，黒板には 3 行の式 が上から順に ‘F= ma, E = mb×
2,E= mc2’ が書いてあった。式はうろおぼえだが，インターネット上にあるように 1 行目が E×
= ma2だ と次元が合わないからさすがに漫画家でも描かないだろう。第 1 著者にとっては運動方程式 ‘F= ma’ があったからこそ面白かったよう な記憶だが，でも確かに× 印が二つあったかもなぁ。 373となる。以下，上付きドットは時間による微分を表す。 もしw(t) の代わりに u(t)≡ w(t) − mg k (9.2) を用いると， ¨u= ¨w なので，運動方程式 (9.1) は m ¨u(t)+ k u(t) = f (t) (9.3) と書き直すことができる。この式 (9.2) の u(t) は，図 9.1 の右から二番目の図に示したように，質点をバネに そぉーっと取り付けて静的につり合った状態からの変位である。この形の方程式は，重力とは関係しない水平 方向にバネ質点系が摩擦の無い床の上で振動する場合の式に一致する。また社会基盤構造の振動問題では，そ の構造が供用され始めた時点の既に静的につり合っている状態からの動的な変動に興味があることの方が多い ので，以後，原則として自重（死荷重）の項は無視して，静的につり合っている状態からの運動方程式 (9.3) を 標準的な運動方程式とする。ここで，外力 f (t) が作用していない振動のことを自由振動と呼ぶ。節題目にある 「非減衰」については今はわからなくていい。 θ(t) ℓ m 図 9.2 振り子 上の例は質点の併進運動を対象としたが，社会基盤構造の例えばケーブル等は回転運 動もするだろう。そこで次に図 9.2 の振り子運動を考えてみよう。振り子は重力場でし か運動しないので重力加速度g を考慮しなければならない。そこで振り子の糸は質量が 無く伸び縮みもしないものとして，上の支点回りの反時計周りの回転をθ(t) とする。支 点回りにこの振り子に作用しているモーメントは重力の影響のみで，反時計回りを正と すると M(t)= −mg ℓ sin (θ (t)) ≃ −mg ℓ θ(t) だけだ。ここでは揺れは小さく|θ(t)| ≪ 1 が成立することを用いて近似した。質量は振り子の先端にしか無いの で，式 (2.4) で定義された慣性モーメントは J = m ℓ2となる。したがって回転についての Newton の運動方程 式 (2.3) は −mg ℓ θ(t) = m ℓ2θ(t) → ℓ ¨θ(t) + g θ(t) = 0¨ (9.4) と表される。これも形式的には式 (9.3) の標準的な運動方程式と同じ形 [(ℓ, g) ↔ (m, k)] になっていることがわ かる。したがってこの文書では，主に式 (9.3) を対象とした解法と振動特性を説明する。 (2) 非減衰自由振動 まず式 (9.3) で外力 f (t) が作用していない状態の自由振動解を求めてみよう。運動方程式は 2 階の常微分方 程式なので，唯一な解を得るためには時刻 t= 0 における状態，つまり初期条件を 2 個与える必要がある。その 最も単純な設定は，質点の初期の位置と初期の速度（以下，初速と呼ぶ）を与えることであり，例えば
u(0)= u0, ˙u(0) = v0 (9.5a, b)
のように与えられることにする。式 (9.3) を辺々 m で割り， m も k も正の定数なので ω ≡ √ k m (9.6) と定義すると，運動方程式は ¨u(t)+ ω2u(t)= 0 (9.7)
となる。この微分方程式の一般解を求めるために u(t)= exp (µ t) と置き，運動方程式 (9.7) に代入すると， µ に ついての特性方程式とその根は µ2+ ω2= 0 → µ = ±i ω と求められる。したがって解は u(t)= A eiω t+ B e−iω t と表現できるが A と B は複素数だ。そこで後述の Euler の公式 (9.62) を用いて実数解の形にすると
u(t)= A sin (ω t) + B cos (ω t) = C sin (ω t − ζ) , C ≡ √A2+ B2, ζ ≡ − tan−1
(B A ) (9.8a, b, c) と求められる。 C は振幅と呼ばれている。二つの積分定数 A と B あるいは C とζ を初期条件式 (9.5) で決定で きれば，それが解になる。ζ は， 0 から始まる sine 関数 sin(ω t) の初期 t = 0 の位相のずれである。 A と B を 使った方の表現を初期条件に代入すると u(0)= B = u0, ˙u(0) = ω A = v0 → A= v 0 ω, B = u0 と求められるので，結局，運動は u(t)= v0 ω sin (ω t) + u0 cos (ω t) (9.9) と表される。 O v0 u0 t T u(t) C 図 9.3 非減衰自由振動 図 9.3 にその運動を示した。図に示したように，元の位置 に同じ速度で戻るまでの時間 T を固有周期（単位は s）と 呼び，それは 1 s 間にω rad の速さの運動が 2π 回るのに必 要な時間なので T =2π ω (9.10) という関係になる。したがって式 (9.6) で定義されるω を固 有振動数2と呼ぶ。併進運動の場合の固有振動数が式 (9.6) なので，式 (9.3) と式 (9.4) の運動方程式同士を比較すれば 明らかなように，振り子運動の固有振動数は ω = √ g ℓ (9.11) で定義できる。つまり振り子の長さℓ を長くするとゆっくり動き，短くすると速く動くことを示している。時 計屋のいろいろな振り子時計（本当に振り子で動いているのかねぇ?）の動きを観察してみるといい。 (3) エネルギ さて，とりあえず理由は考えずに，式 (9.3) に ˙u(t) を乗じて時間積分してみよう。 ∫ m ¨u ˙u dt+ ∫ k u ˙u dt= ∫ f (t) ˙u dt+ const. ここで，第 1 項は K≡ ∫ m ¨u ˙u dt= ∫ 1 2m ˙ { ( ˙u)2}dt=1 2m ∫ d ( ˙u)2=1 2m( ˙u) 2 2正確には固有円振動数あるいは固有角振動数と呼ぶが，以下，省略して単に固有振動数と呼ぶことにする。固有振動数 f の正しい定義 は固有周期の逆数で f≡ 1 T = ω 2π（単位は 1/s あるいは Hz）である。ちなみに rad はもちろん無次元だ。
と変形でき，第 2 項も同様に U≡ ∫ k u ˙u dt= ∫ k u du= k ∫ u du= 1 2k u 2 となるので，上式は結局 K+ U = ∫ f (t) ˙u dt+ const., つまり 1 2m( ˙u) 2+1 2k u 2= ∫ f (t) ˙u dt+ const. (9.12) が成立する。これはエネルギ保存則であり， K は運動エネルギであり， U は位置エネルギである。右辺は外力 の持つエネルギである。どの状態を基準にするのか明記していないので，右辺には定数の項が残っている。自 由振動の場合の解の式 (9.9) で，例えば初速が零の場合の解を式 (9.12) の各項に代入すると K= 1 2mω 2u 0 sin2ωt = 1 2k u 2 0 sin 2ωt, U = 1 2k u 2 0 cos 2ωt となり， K+ U = const. であることがわかる。また K と U の位相が 90 度ずれているので
Kmax= Umax あるいは |K| = |U| (9.13)
であることもわかる。 ℓ x v(x, t) u(t) 図 9.4 質量のあるバネ この最後の結論を用いると，ちょっと面白いことを求めることができる。もしバ ネにも質量があり，単位長さ当たりの質量がµ0である場合の図 9.4 に示した質量 m の質点の振動を考えてみよう。バネの初期長さに比べると伸びは非常に小さいと考 えて，バネの長さはℓ のままだとする。簡単のために，質点の振動による変位を u(t)= A sin (ωt − ζ) とすると，図から明らかなように，バネの途中の任意点 x の変位v(x, t) は v(x, t) = xℓu(t)= Ax ℓ sin (ωt − ζ) になると考えていいだろう。これを用いて全体の運動エネルギを求めると K=1 2m( ˙u) 2+1 2 ∫ ℓ 0 µ0 (˙v)2 dx= 1 2ω 2A2 ( m+µ0ℓ 3 ) cos2(ωt − ζ) である。一方，位置エネルギ（バネの全変形が蓄えている変形のエネルギ）はバネの質量とは無関係で U=1 2k u 2=1 2k A 2 sin2(ωt − ζ) となる。この結果を式 (9.13) に代入すると 1 2ω 2A2 ( m+µ0ℓ 3 ) =1 2k A 2 → ω = v u u t k m+µ0ℓ 3 と固有振動数が求められる。なぜか，バネの全質量µ0ℓ の三分の一だけが振動数に貢献している。 (4) 自由落下と等速円運動 高校の物理学では関数と微分積分を使ってはいけないらしいが，その是非はともかく記憶科目として習って 覚えていた公式等の根拠を少し示しておこう。式 (9.1) でバネを取り去れば，それは自由落下の運動方程式にな るはずだ。簡単のために外力 f (t) も無いものとし，w(t) を g と逆向きに y(t) ≡ −w(t) と定義し直すと，それは m ¨y(t) = −mg → y(t) = −g¨ (∗)
である。つまり，加速度が一定g だから等加速度運動である。初期条件を質点の初期位置と初速で与えるため に，式 (9.5) に倣って y(0) = y0, ˙y(0) = v0 としておく。すると式 (∗) の一般解は w(t) = −1 2g t 2+ c 1t+ c2なので，上の初期条件に代入すると c1 = v0, c2= y0と求められる。したがって，自由落下の質点の位置は上向きを正とすると y(t) = −1 2g t 2+ v 0t+ y0 になる。もし重力場と直交する水平運動 u(t) の場合には式 (∗) に相当する運動方程式の右辺は零なので，上の y(t) の答の重力加速度 g の項を無視すればよく，それは等速運動 u(t)= u0t
である。ただし，初期条件は u(0)= 0, ˙u(0) = u0とした。この両式から t を消去して x(t) ≡ u(t) と置き換える
と，例えばボールを t = 0 に (x(0), y(0)) = (0, y0) の位置から，水平に速度 u0で，鉛直上向きに速度v0で投げ たときの，ある時刻 t のボールの位置 (x(t), y(t)) は y(t) = −1 2g ( x(t) u0 − v0 g )2 + y0+ v2 0 2g を満足する。これが t= v0 g に x= u0v0 g の位置で頂点y = y0+ v2 0 2g に達する放物運動であることは，すぐに理解 できると思う。 また式 (9.4) が水平面内の無重力場の回転運動の場合には， ¨θ = 0 → ˙θ = const. (θ = ωt + const.) の等速度運 動になり，それを等速円運動と呼ぶ。このときθ 方向の速さは vθ≡ r˙θ = rω なの3で，質点の速度ベクトルは
u = vθeθ, er = ex cosθ + ey sinθ, eθ= −exsinθ + ey cosθ
である。ここに四つの ed (d= x, y, r, θ) は各座標軸方向の単位基底ベクトルであるが， exと eyは空間に固定さ
れている。これから加速度ベクトルは
˙
u = vθ˙eθ= −vθθ˙(excosθ + eysinθ
) = −vθθ e˙ r と求められる。つまり加速度は erの負の方向（回転中心方向）に生じていて，その大きさ a は a= vθθ = v˙ θω (9.14) である。このことから，向心力が F= ma = m vθω = m r ω2 = mv 2 θ r となるのだ。高校の教科書に載っている誘 導過程はかなり難しいが，微分を使えば簡単だ。しかし，この向心力や慣性系・遠心力等の物理的な意味はよ くわからないし，流体力学の移流項から求められる式 (3.164a) の向心力もよくわからない。 演習問題 9-1 1. 有名な問題。木のてっぺんにいる猿が落ち始めた瞬間に，その木のてっぺんをねらって地上から投げた 石は，その落ちていく猿に当たることを証明せよ。もちろん，木に届くだけの十分な大きさの初速を与 えることを前提としている。方程式を使わずに，高校等ではどうやって教えてもらったんだろう。 3これも難しい。
9.1.2 減衰自由振動 (1) 減衰の原因 さて前節で求められた自由振動解は，一旦動かされると未来永劫4動き続けるという答になってしまってい る。しかし実際に振動する物体は，例えば空気抵抗等によって振動はそのうち止まってしまう。このような振 動を減衰振動と呼ぶ。前節では，その節題目が示すように減衰を無視していた。その減衰の原因としては，空 気抵抗（速度の 2 乗に比例）だけではなく，材料の中の何らかの非可逆変形でエネルギが消費されてしまった り，構造や材料のどこかに摩擦が発生していたりと，種々考えられる。ここではそれをモデル化したもののう ち，最も基本的な 2 種類，つまり速度に比例した粘性抵抗がある場合と，一定の力による固体摩擦抵抗がある 場合の振動を説明する。 (2) 粘性減衰 m c k f (t) f (t) k u c ˙u 図 9.5 粘性減衰 バネが変位に比例した抵抗をするのに対して，粘性抵抗というの は速度に比例した抵抗のことである。典型的なモデルはよく Kelvin-Voigt モデルと呼ばれる系で，図9.5 に示した。バネでない方の記号 はダッシュポットと呼ばれる装置を表し，これが粘性抵抗を受け持っ ている。図の中央に描いたように，ある種のオイルが密封されたシ リンダーに，穴の開いたピストンがついているものを思い描いて欲し い。ピストンをゆっくり動かすとオイルは穴をスムーズに通り抜け， ほとんど抵抗が発生しない。しかし速くピストンを動かそうとした場合には，オイルはそれほど速くは穴を通 り抜けられず，これが抵抗力5になる。このように，速度に比例した抵抗を発生しているのがダッシュポットで あり，その比例係数を c とする。図の右端に示した力のダイアグラムのように，質点にはダッシュポットから 上向きに c ˙u(t) の抵抗を受けるので，運動方程式は
f (t)− k u(t) − c ˙u(t) = m ¨u(t) → m ¨u(t) + c ˙u(t) + k u(t) = f (t) (9.15) となる。ここで ω ≡ √ k m, β ≡ c 2 mω (9.16a, b) と定義すると，外力 f (t) が作用していない自由振動の場合には，上の運動方程式 (9.15) は
¨u(t)+ 2 β ω ˙u(t) + ω2u(t)= 0 (9.17)
と書くことができる。ここに c は粘性減衰係数と呼ばれ，β は（粘性）減衰定数6と呼ばれる。 式 (9.17) の解を exp (µ t) と置き代入すると， µ についての特性方程式とその根が µ2+ 2 β ω µ + ω2= 0 → µ = ω ( −β ±√β2− 1 ) (9.18a, b) と求められる。この解はβ の大きさによって異なる振動挙動に対応したものになる。 4「一劫」という時間の長さについては諸説あるようだが，第 1 著者が（単位が取り易いという噂で）受けた笠原一男先生の「国史（中 身は浄土真宗の哲学）」の授業では「100 km 立方の岩を，年に一度だけ舞い降りてきた天女が，破れないようにそぉーっと羽衣で 1 回 だけ撫でてそれが磨り減るまで」と定義された。結構長いぞ，呵呵。 5ドアクローザ（ドアダンパー）がその具体例である。 6多分道路橋示方書 [172] で使っているからだと思うが，最近は h を使う書籍や論文が多いようである。しかしここでは，著者にとって 道路橋示方書が金科玉条でもないということはともかく，最も基本的で重要な物理パラメータの振動数ω と揃えるためだけの理由でギ リシャ文字β を用いている。
O t u(t) O t u(t) 図 9.6 過減衰と臨界減衰 O t u(t) un un+ j t 図 9.7 減衰振動 1. まずβ > 1 の場合の解は指数関数あるいは双曲線関数であり，振動応答はせず，図 9.6 の左図のような挙 動を示す。この現象を過減衰と呼ぶ。機械部品にとっては必要な特性かもしれないが，社会基盤構造の 振動ではこういう応答はほとんど期待できない。 2. 次にβ = 1 の場合には特性根が重根となるので，解は u = (A + Bt) exp (−ωt) となり，図 9.6 の右図のよ うな挙動を示す。指数関数の方が線形項より速く零に収束するため発散はしない。この現象を臨界減衰 と呼ぶ。式 (9.16b) の定義からβ = 1 のときの減衰係数 ccr≡ 2 m ω = 2 √ m k を臨界減衰係数と呼ぶ。 3. 最後に，一般的な社会基盤構造の場合にはβ < 1 と捉える7ことができる。実際にはさらにβ ≪ 1 と考え ていいくらい小さく，具体的にはせいぜいβ = 0.02 とか 0.05 程度である。したがって特性根は µ = −β ω ± i ωd, ωd≡ ω √ 1− β2 (9.19a, b) となるので，その実数解は
u(t)= exp (−β ω t) (A sin ωdt+ B cos ωdt)=
[ C e−β ω t]sin (ωdt− ζ) (9.20) と表される。つまり，振幅 C が指数関数的に小さくなりながら（上式の鉤括弧部分）振動することにな る。図 9.7 の左側の図にその例を示した。振幅のピークを結んだ一点鎖線が指数関数的に小さくなる様子 を示しており，振幅はその比率で小さくなっていくが，周期は一定のままの振動である。 以下，この最後のβ ≪ 1 の場合に限定して説明を続ける。このような減衰振動の場合の振動数は式 (9.19b) のωdだから，周期 T も T = 2π ωd = 2π ω√1− β2 (9.21) 7曖昧な表現をしたのは，実際の社会基盤構造の減衰が純粋に粘性減衰だけとは限らないからだ。現象の見かけ上の減衰がβ < 1 のよう に捉えることができるということだ。否むしろ粘性は実際の減衰の原因ではないだろうとも想像できるのだが。
となる。しかし一般には前述のように社会基盤構造の減衰定数β は非常に小さく，例えば鋼構造では 0.02 ∼ 0.03，コンクリート構造では 0.03 ∼ 0.05 と考えていいため，ωdとω の差は無視できるくらい小さく 社会基盤構造の場合は ω ≃ ωd, T ≃ 2π ω (9.22) と近似しても構わないのが普通である。なお実際の社会基盤構造の減衰が粘性減衰だとは限らないが，それで ほぼ近似でき，また近似できると便利なので，粘性減衰モデルを用いているのだろうと考えればいい。では， 実測した応答から粘性減衰定数を求める方法を考えてみよう。実測データが図 9.7 の右側の図のようだったと しよう。 n 番目の振幅と (n+ j) 番目の振幅（図では j = 3）をそれぞれ un, un+ jとすると，その振幅比（多くの 教科書では j= 1）を減衰比8と呼び，それは理論的には un un+ j = C exp (−βωt)
C exp{−βω (t + jT)} = exp ( jβωT) = exp √2π jβ 1− β2 (9.23) となるので，定数になる。したがって，この関係式を用いて構造物の振動試験の実測値を整理すれば，見かけ 上の減衰定数β を逆算することができる。また，この減衰比の対数を対数減衰率 と呼ぶ（これも j = 1 として 定義されるのが普通）が，それは δ ≡ ln un un+ j = 2π j β √ 1− β2 (9.24) となり，これも定数なので，この関係式から β = √ δ δ2+ 4π2j2 (9.25) と，減衰定数を求めることができる。あるいは，β ≪ 1 が成立するなら，この関係はもっと近似できて β ≃ δ 2π j (9.26) としてもいい。 2 4 6 −0.5 0 0.5 ¨u (m/s2) t (s) 図 9.8 実測データからのパラメータ同定 図 9.8 には，東北大学の学部の学生実験で得た 1 自由度 系の加速度記録を× 印で示した。実験では，鋼の薄片で板 バネを作りその先端に重りを取り付け，板バネを少し曲げ た状態で手を離し，そのあとの自由振動の重りの加速度を 記録した。手を放すときの乱れがノイズとして運動に影響 を及ぼすので，初期状態から数周期経たあとの記録を自由 振動していると判断して図示した。β は十分小さいと仮定 した上で，この記録の横軸をまたぐ 15 周期分の時間から周 期を T = 0.485 s と読み取り，近似式 (9.22) から固有振動数 をω = 12.95rad/sと同定した。さらに，この図の範囲の正方向の振幅データを用いて，いろいろな n や j で式 (9.24) と近似式 (9.26) を用いて算定したものを平均してβ = 0.00385 と同定できた。同定したパラメータを用 いて求めた応答が実線で，ほとんど重なっている。 E(t) L C R 図 9.9 単ループ回路 ところで唐突だが図 9.9 は電気回路の代表例である。抵抗器の電気抵抗を R と し，コンデンサーの電気容量を C，コイルの自己インダクタンスを L とする。起電 力が E(t) のとき，この回路内の電荷 Q(t) と電流 I(t) は L ¨Q(t)+ R ˙Q(t) + 1
CQ(t)= E(t), I(t) = ˙Q(t) (9.27a, b)
を満足するような変動をする。これはまさに粘性減衰振動の運動方程式と同じである。面白いですねぇ。
(3) 固体摩擦減衰 もう一つの減衰モデルの代表である固体摩擦を考える前に，式 (9.7) の運動方程式を変位と速度の平面上で表 すことを考えよう。このような平面を位相平面と呼び，第 9.5 節の安定問題等でもよく用いられる。まず運動 方程式の加速度項は ¨u=d ˙u dt = d ˙u du du dt = ˙u d ˙u du と書き直してもいいだろうから，運動方程式は ˙ud ˙u du+ ω 2u= 0 → ˙u d ˙u+ ω2u du= 0 となるので，これを項別に積分すれば ( ˙u)2+ ω2u2= const. あるいは (˙u ω )2 + u2= const. (9.28) と表すことができる。右辺の定数は初期条件で決まる値である。誘導過程と結果からも推測できるように，こ れは第 9.1.1 (3) 節で誘導したエネルギ保存則そのものだ。そしてこの式は， u と ˙u を 2 軸とする位相平面上で 楕円状の軌跡を描くことを意味している。あるいは，特にこの場合には u と ˙u ω を 2 軸にした円軌道になる。 k u(t) m F 図 9.10 固体摩擦減衰 さて図 9.10 に示したように，質点との間に一定の動摩擦力 F = µ′mg （µ′ は動摩擦係数）が発生する床の上の振動9を考えてみよう。速度と逆の向きに 摩擦が発生するので， ˙u> 0 のときの運動方程式は
−k u(t) − F = m ¨u(t) → m ¨u(t) + k u(t) + F = 0 となるのに対し， ˙u< 0 の場合には
−k u(t) + F = m ¨u(t) → m ¨u(t) + k u(t) − F = 0 となる。したがって，これをまとめると m ¨u(t)+ k u + Sgn (˙u) F = 0 (9.29) と書く10ことができる。ここに Sgn は引数の符号を表す記号で Sgn ( ˙u)= +1 もし ˙u ≥ 0 −1 もし ˙u < 0 (9.30) と定義した。初期条件は式 (9.5) で与えられる。この微分方程式を速度の符号で場合分けをして時々刻々解いて いけば解を求めることができる（図 9.11 の右側の図は実際そうやって求めた）が，ここでは上述の位相平面を 使ってみよう。なお，運動方程式 (9.29) の第 3 項には F があって斉次の微分方程式ではないので，次節で対象 とする強制振動として捉えないといけないが，次節以降では固体摩擦減衰を対象とはしないので，この節の最 後に扱ってしまうことにした。 この運動方程式 (9.29) に相当する位相平面上の軌道を求めようとすると，まず ˙u d ˙u+ ω2 { u+ Sgn (˙u)F k } du= 0 → 1 ω2˙u d ˙u+ { u+ Sgn (˙u)F k } d { u+ Sgn (˙u)F k } = 0 9実際には，バネが伸び切って振動が止まる瞬間の摩擦力は最大摩擦力程度以下の静止摩擦抵抗でいいわけだが，ここでは簡単のために 動摩擦力のみで近似する。 10あるいは，後述の式 (9.46) で定義する Heaviside 関数を用いて
m ¨u(t)+ k u + H(˙u) F − H(−˙u) F = 0 と書いてもいい。
−5 5 u(t) F/k ˙u(t) ωF/k 5 O −5 t= 0 2 4 −5 5 ωt π u(t) F/k O 図 9.11 固体摩擦減衰応答 としてもいいから，位相平面は (˙u ω )2 +{u+ Sgn (˙u)F k }2 = const. (9.31) と表すことができる。つまり図 9.11 の左側に示したように，速度の向きが変わる度に 2 箇所の円の中心位置 (±F/k, 0) が突然移動する。位相平面は速度軸を上向きに変位軸を右向きに描くのが普通だが，時刻歴応答を描き 易いように反時計回りに 90 度回転して描いた。初期条件の状態から時計回りに回転しながら，変位軸をまたぐ 半周毎に軌道の中心が突然移動し，それに伴い振幅も小さくなっていき，最後は k|u| < F = µ′mg の状態11に 入ったあとの速度が零になる時点（この位相平面の図の縦軸上）で止まることになる。 ちなみに式 (9.12) の誘導と同じようにして，固体摩擦がある場合のエネルギ保存則を誘導しておこう。式 (9.29) の運動方程式に ˙u(t) を乗じて積分すると，最初の 2 項は式 (9.12) と同じ結論になる。第 3 項は ∫ Sgn ( ˙u) F ˙u dt= F ∫ |˙u| dt = F ∫ |du| = F × (質点がたどった道のりの累積) になる。これも高校で習って覚えた公式であろう。したがって，エネルギ保存則は，外力項も含めると 1 2m( ˙u) 2+1 2k u 2+ µ′mg ∫ |˙u| dt = ∫ f (t) ˙u dt+ const. となる。 9.1.3 強制振動 (1) 正弦波外力による強制振動と共振曲線 前節までには振動の最も基本的なパラメータである固有振動数と減衰定数を導入したが，それがどうした? と思っている読者も少なくないだろう。実は我々にとって m と k と c は与えられるのではなく，何らかの設計 を通して設定する量なのだ。そのとき実際には，地震や風や活荷重等の外力が作用した状況における振動が問 題になり，その外力に対して安全な社会基盤構造になるように固有振動数（や減衰定数）を設計する必要があ るのだ。まず代表的かつ単純な外力が作用する場合の振動を対象として，この基本的な二つのパラメータと外 力特性の関係について考察する。外力 f (t) が作用した振動を強制振動と呼ぶ。一般的な問題として式 (9.15) を 解くが，それを m で除して式 (9.17) の形にすると
¨u(t)+ 2 β ω ˙u(t) + ω2u(t)= 1
m f (t) (9.32)
11正確には，静止摩擦係数をµ
s(> µ′) とした場合，速度が零になった瞬間には最大摩擦力程度以下の静止摩擦抵抗 k|u| ≃ µsmg のときに
止まるので，|u| ≤µsmg
と表現できる。
最も基本的な外力として正弦波状の外力を扱う。というのも，任意のほとんどの関数は Fourier 級数で表され ると考えていいので，その基本的な外力として f (t) = f0 sin pt について考えておけば，ほとんどすべての外力
に対する応答を求める基礎になると考えられるからだ。つまり，式 (9.32) は ¨u(t)+ 2 β ω ˙u(t) + ω2u(t)= g0 sin pt, g0≡
f0
m, f (t)= f0 sin pt (9.33a, b, c)
となる。一般解は，自由振動で求めた斉次解と外力に対する特解 up(t) の和で
u(t)= exp (−β ω t) (A sin ωdt+ B cos ωdt)+ up(t) (∗)
と表される。減衰が存在する場合には自由振動解は減衰解であり，時間と共に減少していくような過渡応答な ので，工学的に一番興味があるのは多分特解の方12だろう。そこでここでは特解 u p(t) だけを求めよう。 1 2 3 5 10 β = 0 0.05 0.08 0.2 β = √1 2 1.2 β = 0 0.05 0.08 0.2 0 1 2 3 π 2 π α Md p ω 0 O 図 9.12 動的増幅率と位相遅れ 数学で習った特解の求め方を用いるのが望ましいが， 多分 sin pt, cos pt が解の候補になると予想されるので up(t)= C sin pt + D cos pt
と置いて，これを式 (9.33a) に代入して sine, cosine 毎に 整理すると sin pt(−p2C− 2 β ω p D + ω2C− g0 ) + cos pt(−p2D+ 2 β ω p C + ω2D)= 0 となるので ω 2− p2 −2 β ω p 2β ω p ω2− p2 CD = g0 0 を解けばいい。つまり CD = ( 1 ω2− p2)2+ 4 β2ω2p2 ω 2− p2 2β ω p −2 β ω p ω2− p2 g00 = g0 ( ω2− p2)2+ 4 β2ω2p2 ω 2− p2 −2 β ω p = g0 ω2 { 1− (p ω )2}2 + 4 β2 (p ω )2 1− (p ω )2 −2 β p ω = f0 k { 1− (p ω )2}2 + 4 β2 (p ω )2 1− (p ω )2 −2 β p ω を得る。したがって特解は up(t)= f0 k { 1− (p ω )2}2 + 4 β2 (p ω )2 [{ 1− (p ω )2} sin pt− 2 β p ω cos pt ] (9.34) と求められる。あるいは ust≡ f0 k と定義すると up(t)= ustMd sin (p t− α) , Md≡ 1 √{ 1− (p ω )2}2 + 4 β2 (p ω )2 , tan α ≡ 2β p ω 1− (p ω )2 (9.35a, b, c) 12演習問題を解けばわかるはずだが，実際には，過渡応答を含めた振幅が動的増幅率を超えることがあることには注意するべきである。
と書くこともできる。 ustは外力の最大値 f0が静的に作用したときの静的な変位である。したがって Mdは， その静的な変位に対する動的な変位の割合を表しているので動的増幅率と呼ばれる。α は減衰による位相の遅 れを表している。図 9.13 の t = 0 付近の実線を拡大すれば，この位相の遅れを実感できる。 Mdとα を図 9.12 に示した。非減衰のβ = 0 の場合には p = ω で M d → ∞ となり発振してしまう。この現象を共振と呼び， Md のこの図の曲線を共振曲線と呼ぶことがある。そして p= ω の外力周波数の作用状態を共振点と呼ぶ。前述の ように社会基盤構造のβ は非常に小さいため，かなり大きい動的増幅率を持つ共振が生じる可能性がある。 20 40 −5 5 u(t) ust ωt Md p ω = 0.9 β = 0.05 O 図 9.13 減衰系の正弦波強制振動 では一例として，過渡応答も含めて，零初期条件つまり u0 = 0, v0 = 0 の場合の強制振動応答を求めておこう。式 (∗) に式 (9.34) を代入し，それを零初期条件に代入すると，二つ の積分定数を A ust = β √ 1− β2 B ust − p ωd { 1− (p ω )2}2 + 4 β2 (p ω )2 { 1− (p ω )2} , B ust = 2βp ω { 1− (p ω )2}2 + 4 β2 (p ω )2 と求めることができる。式が長くなるので，上式の A の右辺には B を残してあるので注意して欲しい。例とし て p ω = 0.9 で β = 0.05 のときの解を描いたのが図 9.13 である。実線がその応答で，零初期条件で始まってい ることがわかる。破線は応答の中のωdで振動する成分，つまり過渡応答（自由振動解）の成分を分離して取り 出したものである。減衰定数が大きめなので過渡応答の振幅はすぐに小さくなっていく。この成分が無くなっ てしまえば，実線の振幅は式 (9.35b) の Mdになり，周波数が p の振動成分のみ，つまり特解のみが残る。図 中の水平の 2 本の破線が Mdのレベルである。共振点の 90% 程度の外力周波数では静的変位の 5 倍もの大きな 動的な変位が生じていることが，この応答図からも明らかだ。また過渡応答も含めた場合の全振幅は，動的増 幅率よりときどき大きくなることにも十分注意しておく必要がある。 さて，社会基盤構造を設計したとしよう。それは構造の質量 m と剛性 k を設計した13ことに相当するから， 動的に重要なパラメータの一つの固有振動数ω を設計したことになる。もし，その構造に作用する外力の周波 数 p がその固有振動数に近い場合には，共振して応答の振幅は非常に大きくなり，構造は破壊するかもしれな い。動的な設計というのは，簡単に言うと，共振曲線や強制振動の時刻歴応答を確認しながら m と k を設計す るということである。 演習問題 9-2 2. 自分が使い慣れたソフトウェアを用いて減衰振動の応答曲線図 9.7 を描いてみよ。最近はソフトウェアが 非常に高性能になったため，スプライン関数等で滑らかな曲線を描くことが容易にできるが，実はそこ には大きな落とし穴がある。スプライン関数の性質を知らずに安易に使うと，とんでもない結果になっ てしまうので注意して欲しい。こういう周期的な曲線を描く場合，一番安全で確実なのは 1 周期当たり 20 個以上のデータ点を直線で結ぶことである。 3. 自分が使い慣れたソフトウェアを用いて共振曲線図 9.12 を描いてみよ。これも試しにスプライン関数で 描いてみよ。データ数が少ない場合には，とても変な曲線になることがある。 13動的制御の設計を減衰特性を用いて行う場合を除き，一般に社会基盤構造では減衰のパラメータ c を設計することはできない。
4. 零初期条件 u0 = 0, v0 = 0 の質点に正弦波外力 f (t) = f0 sin pt が作用したときの応答を求め，図 9.13 を 描け。その図では初速（線の傾き）が零になってますか? 傾きが零になるように工夫できますか? ソフト ウェアの結果を鵜呑みにするのは危険だ。例えばβ = 0.02 として， p ωが 0.2, 0.5, 0.97, 1.1, 1.5 等の場合 の応答を求め，±M dの線も描いた上で|u| > ustMdがどのような状況で発生するか比較してみよ。 5. 非減衰β ≡ 0 で外力が f (t) = f0 sinωt のとき，つまり p = ω の共振点では，式 (9.35b) の動的増幅率は Md→ ∞ になることから，特解が up(t)∼ sin pt = sin ωt ではないことが明らかである。では零初期条件 下の解を求めよ。数学が不得意でも，共振で発散することを念頭に置けば解は容易に推測できると思う が，例えば定数変化法を思い出そう。結果を図示してその発散の様子も確認せよ。 6. 固体摩擦による減衰の場合の運動方程式 (9.29) は速度が零になる毎に時々刻々特解の符号が変わるだけ なので，連続条件をその都度満足させれば容易に解くことができる。その解が以下のようになることを 確かめよ。例えば初期条件を u(0)= u0> 0, ˙u(0) = v0> 0 とし， τ ≡ ωt とした上で x≡ u0 F/k+ 1, y ≡ v0 ωF/k, τ0≡ 0, 0 < τ1≡ tan −1y x < π 2, τn = τn−1+ π (n ≥ 2) とすると，τn≤ τ < τn+1(n≥ 0) の解は次式になる。 u(τ) F/k = 1 − √ 2n x2+ y2 (x cos τ + y sin τ) + (−1)n+1 (2) 支点の正弦波状強制変位による応答 c k m u(t) w(t) 図 9.14 支点変位に対する 応答 次に地震に対する応答を思い描くと，構造自体に外力が作用している前節のよ うな状況ではなく，バネの支点が強制的に変位させられた状況における振動のよ うに見える。そこで図 9.14 に示したように，支点が水平方向に強制的にw(t) だ け変位させられた場合の質点の水平運動 u(t) を考察しよう。バネとダッシュポッ トは，支点の運動との相対変位と相対速度に比例した抵抗力しか生まないから， 運動方程式は
m ¨u(t)+ c {˙u(t) − ˙w(t)} + k {u(t) − w(t)} = 0 (9.36)
となる。 Newton の法則の慣性項は絶対加速度 ¨u(t) で表されているので，第 1 項はこれまでと同じである。さ てここで，支点の変位に対する相対的な質点の変位 v(t) ≡ u(t) − w(t) (9.37) を用いて，上の運動方程式 (9.36) を書き直すと m ¨v(t) + c ˙v(t) + k v(t) = −m ¨w(t) (9.38) と書くことができる。この式と式 (9.15) を比較すれば明らかなように，支点の強制変位に対する相対的な応答 v(t) は，見かけ上 f (t)= −m ¨w(t) = −(設計しようとしている構造の質量) × (地盤の地震入力加速度) (9.39) という外力が直接質点に作用したときの強制外力応答と同じであることがわかる。この式 (9.39) は非常に重要 な概念であり，耐震設計において用いられている震度法の根拠になっている。つまり，過去に観測された大地 震の地震加速度を参考にして設定された設計入力地震加速度に，設計されようとしている社会基盤構造の質量 を乗じることによってその構造に作用させるべき地震外力が算定でき，それを作用させたときの構造の抵抗を 照査することによって耐震設計ができることを示している。
1 2 3 5 10 β = 0 0.05 0.08 0.2 Mw p ω O 図 9.15 支点変位応答の動的増幅率 基本的な例として，支点変位が正弦波状14の w(t) = w0 sin pt (9.40) で与えられる場合を上式 (9.38) に代入すると，運動方程 式は m ¨v(t) + c ˙v(t) + k v(t) = m p2w0 sin pt と書ける。したがって前節の外力の振幅を f0 = m p2w0 と置き換えれば，そのまま前節の解を流用できるので， 特解は式 (9.35) とほぼ同様になり vp(t)= w0Mwsin (p t− α) , Mw= (p ω )2 Md (9.41a, b) となる。位相の遅れα は式 (9.35c) と同じである。式 (9.41b) の動的増幅率 Mwを図 9.15 に示した。つまり p ω ≫ 1 の場合には vp ≃ |w| となるが，位相差α はπ になる。これは長大橋のような柔構造の応 答に相当している。 p ω ≪ 1 の場合には相対変位はほぼ零になるが， Md ≃ 1, α ≃ 0 なので vp ≃ w0(p/ω)2 sin (pt) とな り vp ≃ 1/ω2| ¨w| である。したがってバネの抵抗は k vp = m| ¨w| となり，直接地震加速度に対 応する力で抵抗しなければならない。これが剛な構造の応答に相当している。 と15いうことがわかる。 ところで地震計というのは一体どうなっているんだろう。地盤と一緒に動いている地震計が，自分自身の運 動を記録できるのは不思議ではないだろうか。この疑問に答えてくれるのが，ここで求められた解である。つ まり地震計はここで定義した相対変位v(t) を記録紙に書き込んでいるのだから，相対変位 v(t) の挙動を考察す れば地震計の原理が理解できるはずだ。地震計の m と k を適切に設計して，観測したい地震動の振動数 p に対 して p ω ≫ 1 になるように製作すると上述のように 180 度の位相差で vp(t) ≃ |w(t)| となるので，地盤の変位応 答を記録できる。これと対照的に，地震計の m と k を調整して p ω ≪ 1 になるように製作すると，この場合は 上述のように vp(t) ≃1/ω2| ¨w(t)| となるので地盤の加速度応答が記録できるのである。 (3) 単位衝撃応答 ここまでは時間変化が滑らかな周期外力を考察してきたが，それとは全く異なる特性を持つと想像される衝 撃外力が作用したときの振動について考察する。まず，衝撃そのものをどのようにモデル化するかが問題にな るだろう。いま質点 m に，時刻 t = τ にある衝撃外力 f (t) を作用させたとする。衝撃なので，非常に短い時間 ϵ だけ作用していることにしよう。そのときの質点の運動方程式を t = τ −ϵ/2から t= τ +ϵ/2まで積分すると， まず慣性項は ∫ τ+ϵ/2 τ−ϵ/2 m ¨u dt= ∫ τ+ϵ/2 τ−ϵ/2 m d ˙u= m ˙u ( τ + ϵ 2 ) − m ˙u(τ − ϵ 2 ) 14残念なことに，地震波のような不規則波は Fourier 級数では表せないのだが，ま，共振曲線の特性を知るためという意味で。 15読者の多くは祭の屋台等で売っているヨーヨー風船を知っていると思う。これを極めてゆっくり (p≪ ω) 突いた場合の風船は手の平に は当たらず風船の変位と手の変位は同じ向きにほぼ同じ大きさ，つまり相対変位はほぼ零 (Mw≃ 0) で位相差 α は零である。逆に極め て速く (p≫ ω) 突くと風船はほとんど動かず手だけが上下運動する。つまり相対変位はほぼ手の変位に等しく (Mw≃ 1) 位相差が α ≃ π になる。そして共振点付近で手を動かすと風船の変位が増幅されてようやく手の平に到達して風船を突くことができるようになる。
となる。一方，バネの抵抗の項は位置エネルギの短時間の変化分であり，粘性抵抗も短時間の速度変化になる ので，次の式のようにϵ → 0 の極限では零になる。つまり ∫ τ+ϵ/2 τ−ϵ/2 k u dt→ 0, ∫ τ+ϵ/2 τ−ϵ/2 c ˙u dt= ∫ τ+ϵ/2 τ−ϵ/2 c du= c {u(τ + ϵ/2) − u(τ − ϵ/2)} → 0 である。したがって，運動方程式を積分したものは結局 m ˙u ( τ +ϵ 2 ) − m ˙u(τ −ϵ 2 ) = ∫ τ+ϵ/2 τ−ϵ/2 f (t) dt (9.42) となる。この式は実は運動量保存則であり，左辺は t= τ における運動量の変化量であり，右辺は衝撃外力 f (t) による力積である。つまり，ある瞬間 t= τ に衝撃力を与えるということは，その時点に短時間の力積を与える ことを意味する。このことは高校の物理学の教科書にバットでボールを打つ絵等で説明されていたと思う。 τ ϵ 1 ϵ 面積が 1 t f (t) O 図 9.16 単位衝撃外力 この力積が単位量，つまり 1 のときの外力 f (t) を単位衝撃外力 と呼ぶことにしよう。例として挙げたのが図 9.16 である。このよ うな関数 f (t) は，集中せん断外力をモデル化したときの式 (4.57) の Dirac のデルタ関数δ(x − a) と同じであり，それは ∫ c b ψ(x) δ(x − a) dx = ψ(a) もし b < a < c 0 もし a< b あるいは c < a (9.43) で定義される。図はちょっとだけ嘘だ。つまり f (t)= δ (t − τ) → 式 (9.42) 右辺の力積 = ∫ τ+ϵ/2 τ−ϵ/2 f (t) dt= ∫ τ+ϵ/2 τ−ϵ/2 δ (t − τ) dt = 1 となる。この質点が時刻 t = τ まで静止していたとすると， ϵ → 0 の極限において式 (9.42) の左辺第 2 項は零 になるので m ˙u(τ) − 0 = ∫ τ+ϵ/2 τ−ϵ/2 f (t) dt= ∫ τ+ϵ/2 τ−ϵ/2 δ (t − τ) dt = 1 → ˙u(τ) = 1 m を得る。つまり単位衝撃の作用は，その時点に強制的な速度1/mを与えた条件と等価なのである。したがって， 1 質点系に対して時刻 t= τ に単位衝撃外力を作用させたときの初期値問題は t ≥ τ において m ¨u(t)+ c ˙u(t) + k u(t) = 0, u(τ) = 0, ˙u(τ) = 1
m (9.44a, b, c)
で与えればいいことになる。もちろん t< τ における状態は u(t) ≡ 0 である。なお「単位」衝撃の「単位」は無 次元なので，振動問題のデルタ関数δ(t − τ) の次元は [時間]−1と考えればいい。
式 (9.44a) の一般解は
u(t)= exp (−β ω t) (A cos ωdt+ B sin ωdt)
なので，これを初期条件式 (9.44b) (9.44c) に代入すると
0= exp (−β ω τ) (A cos ωdτ + B sin ωdτ) ,
1
m= exp (−β ω τ) ωd (−A sin ωdτ + B cos ωdτ) − β ω exp (−β ω τ) (A cos ωdτ + B sin ωdτ) を満足しなければならない。これを整理すると
− sin ωcosωddττ cos ωsinωddττ AB = exp(βωτ) 01 mωd
となるので，これを解くと AB = exp(βωτ)
cossinωωddττ − sin ωcosωddττ 01 mωd = exp (mβ ω τ)ωd − sin ωcosωddττ のように積分定数が求められる。したがって解は u(t)= ui(t;τ) = ui(t− τ) ≡ 0 1 t< τ mωd exp{−β ω (t − τ)} sin {ωd(t− τ)} t> τ (9.45) と求められる。これを単位衝撃応答と呼んでいる。衝撃が与えられた時刻τ も関数の引数に明示するのが普通 であり，また式の形から明らかなように (t− τ) のみの関数なので，左から三つ目の表現がよく用いられる。な お ui(t− τ) の次元は [時間][質量]−1あるいは [長さ][力]−1[時間]−1である。 1 H(t− τ) O τ t 図 9.17 Heaviside 関数 さて，式 (9.45) に示したように，時刻 t をτ の前後に場合分けして二つの定義 をするのはとても面倒だし，使い難いので，これを一つの式にまとめてみたい。 そのために H(t− τ) = 0, t < τ 1, t > τ (9.46) 少し難しく，やや正しく書くと δ(t − τ)=d dH(t− τ) dt と定義される Heaviside 関数（階段関数）を導入する。なお上式のδ(t − τ) は式 (9.43) で定義した Dirac のデル タ関数で，いずれも正しくは関数ではなく超関数なので，最後の式の等号= は超関数のための特殊な等号であd る。超関数的な等号の意味については式 (10.10) 下の説明を参照のこと。さて，この Heaviside 関数の性質を利 用すれば，任意時刻 t> 0 の単位衝撃応答は ui(t− τ) = 1 mωd exp{−β ω (t − τ)} H(t − τ) sin {ωd(t− τ)} (9.47) と書くことができる。非減衰の場合はβ ≡ 0 なので ui(t− τ) = 1 mωH(t− τ) sin {ω (t − τ)} (9.48) と表される。この表現がとても便利だということは次節以降でわかる。なお初期値問題の初期条件式 (9.44b) (9.44c) から明らかなように，単位衝撃応答は時刻 t = τ における零初期条件の解ではないことには十分注意す る必要がある。 (4) 任意外力に対する応答 では大きさが単位量ではなく，図 9.18 の左上の図のように力積（面積）が F の衝撃外力に対する応答が u(t)= F ui(t− τ) になることは，誰も疑問には思わないだろう。それならば，その右側の図のように異なる時刻に三つの異なる 大きさ（面積）の衝撃外力（力積）が作用したときの応答が u(t)= 3 ∑ n=1 Fnui(t− τn)
になるのも納得してもらえると思う。否，任意時刻 t の応答は場合分けをしないといけないから，三つ全部を 足してしまってはいけないのではないかと思う人もいるかもしれないが，例えば上式の総和の第 2 項は F2 mωd exp{−β ω (t − τ2)} H(t − τ2) sin{ωd(t− τ2)} のように Heaviside 関数が含まれているから， t < τ2の間は無いのと同じになるのだ。したがって，任意の時 刻 t に対する解には三つを足し算しておいても大丈夫なのだ。 Heaviside 関数は便利でしょ。 f (t) f (t) f (t) O t O t O t F τ τ1 τ2 τ3 F1 F2 F3 τ ∆τ 図 9.18 任意外力 それでは図 9.18 の下の図のような任意の（図では連続関数だが不連 続でも構わない）関数 f (t) で表された外力に対する応答も，同様に， 単位衝撃応答を用いて表現できないだろうか。考え易くするために， 図のように十分小さい∆τ ずつの短冊に切って（そもそも Riemann 積 分をそのように定義したから）考えると，例えば網掛けした部分は， 面積∆τ f (τ) の衝撃として近似的に取り扱えばいいことがわかる。し たがって，この網掛け部分の外力によって生じる応答は ∆τ f (τ) ui(t− τ) でいいから，時刻 t= 0 以降にしか外力が作用しなかったとすれば u(t)= ∞ ∑ τ=0 ∆τ f (τ) ui(t− τ) がその応答になることは理解できると思う。 Heaviside 関数のおかげで未来は過去に影響しないから∞ まで加 算しても構わない。このあと∆τ → 0 の極限をとる16ことによって∑∆τ → ∫ dτ になるので，任意外力 f (t) に対する応答は u(t)→ ∫ ∞ 0 dτ f (τ) ui(t− τ) ⇒ u(t) = ∫ ∞ 0 ui(t− τ) f (τ) dτ (9.49) と表現できる。この表現を Duhamel 積分と呼んでいる。この解は t = 0 における零初期条件の解になってい る。このように単位衝撃応答 ui(t− τ) は， t = τ の瞬間の外力 f (τ) が全体応答に及ぼす影響を示していること から，静的な構造力学の影響線と数学的には同じものである。数学の，いわゆる Green 関数である。 例として，図 9.13 で示した f0 sin pt に対する零初期条件の応答を求めてみよう。減衰解は面倒なので非減衰 の場合を例にすると，同図の横にある解から，正解は u(t) ust = ω2 ω2− p2 sin pt− p ω ω2 ω2− p2 sinωt, ust≡ f0 k である。式 (9.49) に f (τ) = f0 sin pτ であることと式 (9.48) の単位衝撃応答とを代入すると u(t)= f0 mω ∫ ∞ 0 sin pτ H(t − τ) sin ω(t − τ) dτ = f0 mω ∫ t 0 sin pτ sin ω(t − τ) dτ となる。 Heaviside 関数から t> τ なので， τ の積分範囲は t までになる。さらに演算すると u(t)= f0 mω ( sinωt ∫ t 0
sin pτ cos ωτ dτ − cos ωt ∫ t 0 sin pτ sin ωτ dτ ) を得るので，この積分を実行すれば解が求められる。これも面倒なので割愛するが，最終的な結果は上の正解 と一致する。計算では t とτ の区別を決して間違わないように ! 16それが Riemann 積分つまり普通の積分の考え方に相当するだろうから。
(5) Duhamel 積分は正しいか
前節で外力の関数 f (t) を短冊に切ったことを怪しいと思った読者のために証明しておこう。なおこの節は North-western 大学の Olmstead 先生 (1980 年頃当時）の ‘Differential Equations of Mathematical Physics’ の講義ノー トを参考にした。解きたい問題は
m ¨u+ c ˙u + k u = f (t), u(0) = 0, ˙u(0) = 0 (a) である。一方， t= τ に単位衝撃外力が作用した場合の単位衝撃応答 ui(t− τ) は次式を満足する。
m ¨ui+ c ˙ui+ k ui= δ(t − τ), ui(0)= 0, ˙ui(0)= 0 (b) ここにδ(t − τ) は式 (9.43) で定義した Dirac のデルタ関数である。
さて目的は u(t) を ui(t−τ) で表すことだが，その補助のために，単位衝撃問題式 (b) に対する随伴問題と呼ば
れる問題を定義する。振動問題は実は自己随伴系ではなく，随伴問題は時刻 t = τ1に単位衝撃外力が作用した
時の応答 u∗i(t− τ1) に関する終局値問題になり，次式で与えられる。
m ¨u∗i − c ˙u∗i + k u∗i = δ(t − τ1), u∗i(∞) = 0, ˙u∗i(∞) = 0 (c)
ここで，運動方程式の減衰項の符号が負になっていることと， t = ∞ で終局条件が与えられていることに注意 して欲しい。この随判問題と元の問題の「関わりあい（仮想仕事）」を考えるために，式 (a)1の両辺に u∗i を乗 じて時間に関して 0 から∞ まで積分しよう。物理的には仮想仕事を求めていることになり ∫ ∞ 0 u∗ i(t− τ1) (m ¨u+ c ˙u + k u) dt = ∫ ∞ 0 u∗ i(t− τ1) f (t) dt (d) という関係が成立する。左辺を部分積分すると 左辺=(u∗ im ˙u− ˙u ∗ im u+ u ∗ ic u) ∞ 0 + ∫ ∞ 0 u(m ¨u∗ i − c ˙u ∗ i + k u ∗)dt となる。この第 1 項に初期条件式 (a)2(a)3と終局条件式 (c)2(c)3を代入し，第 2 項の被積分関数の中の括弧部 分に式 (c)1の運動方程式を代入すれば 左辺= ∫ ∞ 0 u(m ¨u∗ i − c ˙u∗i + k u∗ ) dt= ∫ ∞ 0 uδ(t − τ1) dt= u(τ1) となる。最後の等号では，式 (9.43) のデルタ関数の定義を用いた。結局上式 (d) からは 上式 (d) の左辺= u(τ1)= 上式 (d) の右辺 = ∫ ∞ 0 u∗i(t− τ1) f (t) dt という関係を得る。あるいは t とτ1とを入れ替えると u(t)= ∫ ∞ 0 u∗i(τ1− t) f (τ1) dτ1 (e)
と表される。一見目的に達したようだが，まだ “Close but no cigar!” である。
次に，同様の「関わりあい（仮想仕事）」を単位衝撃の問題と随判問題の間で考えると ∫ ∞ 0 u∗i(t− τ1) (m ¨ui+ c ˙ui+ k ui) dt= ∫ ∞ 0 u∗i(t− τ1)δ(t − τ) dt ( f ) であるから，上の演算と同様にして左辺を部分積分して，初期条件式 (a)2(a)3と終局条件式 (c)2(c)3と微分方 程式 (c)1を代入した上で，デルタ関数の定義を用いれば 左辺=(u∗ im ˙ui− ˙u∗im ui+ u∗i c ui) ∞ 0 + ∫ ∞ 0 ui (m ¨u∗ i − c ˙u ∗ i + k u ∗) dt=∫ ∞ 0 ui(t− τ) δ(t − τ1) dt= ui(τ1− τ)
となる。一方，上式 ( f ) の右辺もデルタ関数の定義を用いれば u∗i(τ − τ1) となるから，結局，式 ( f ) からは ui(τ1− τ) = u∗i(τ − τ1) (g) という関係を得る。つまり随判問題は，元の衝撃問題の観測時刻と載荷時刻を入れ替えた問題であることがわ かる。静力学の相反定理の式 (3.144) や式 (4.65) と同じである。 最終的にこの式 (g) を式 (e) に代入すれば u(t)= ∫ ∞ 0 ui(t− τ1) f (τ1) dτ1 を得るが，τ1はどんな記号を用いてもいい積分変数に過ぎないから添え字を省略して u(t)= ∫ ∞ 0 ui(t− τ) f (τ) dτ (9.50) となるので， Duhamel 積分が成立することを証明することができた。 (6) 数値的な解法 20 40 −5 5 u(t) ust ωt ω∆t = 1 ω∆t = 0.25 厳密解 O 図 9.19 差分法で強制振動を解く 社会基盤構造設計の実務上は，後述のように有限要素 法を用いてコンピュータによって動的解析をする場合が ある。そのとき微分方程式を時間軸方向に積分する際の 一つの数値的解析法に，微分を差分と呼ばれる方法で近 似して直接解く手法がある。ここでは中央差分 [8] の考 え方を簡単に述べておく。十分小さな時間ステップ∆t を設定して時刻 t 前後で
u−≡ u(t − ∆t), u ≡ u(t), u+≡ u(t + ∆t) と定義したとき，その時刻 t における速度と加速度を ˙u(t)= 1 2∆t ( u+− u−), ¨u(t) = 1 ∆t { u+− u ∆t − u− u− ∆t } = 1 (∆t)2 ( u−− 2u + u+) (9.51a, b) と近似するのだ。時刻 t の速度は，いわゆる中間値の定理を信じて t± ∆t 間の変位の線形勾配で定義してある。 また加速度は，時刻 t の直前直後の変位の線形勾配で定義した速度の変化率で定義してある。これを運動方程 式 (9.15) に代入すると m (∆t)2 ( u−− 2u + u+)+ c 2∆t (−u−+ u+)+ ku = f (t) となるので， u+を求める式として整理すると u+= ( 1 1 ∆τ2 + β ∆τ ) [ f (t) k − ( 1− 2 ∆τ2 ) u− ( 1 ∆τ2 − β ∆τ ) u− ] , τ ≡ ωt, ω ≡ √ k m (9.52a, b, c) という関係が求められる。つまり，時刻 t− ∆t と t の変位 u−と u がわかっていれば，この式から時刻 t の解 u+ を時々刻々と求めることができる。 ただし，一番最初のステップの u(−∆t) が定義されていない上に，初期条件のうちの初速が入力できていな い。そこで例えば時刻 t= 0 だけ後退差分で ˙u(0)=u(0)− u − ∆t = 初期位置− u− ∆t = 初速
start
?
ωとβを入力。∆tとNを設定。 f (t)/kを与える式を設定。
n= 0として，u := u(0), ˙u := ˙u(0)
? u(−∆t) = Eq.(9.53) ? t= n ∆t ? u+:= Eq.(9.52a) ?
˙u, ¨u := Eq.(9.51)
? tとu, ˙u, ¨uを出力 ? u−:= u, u := u+ n= n + 1 ? n= N no yes ? end -図 9.20 中央差分のフローチャート 10 20 −10 0 10 0 10 20 30 −4 −2 0 2 入力地震波 ¨ w(t) m/s2 t (s) u(t) w(t) ¨u(t) m/s2 実験結果 数値解析結果 t (s) 図 9.21 地震応答の例 から u−を概算することも可能だが，差分のとり方に整合性が無くなる。そこでこれについては，式 (9.51) の二 つの式から u+を消去して得られる関係 u−= u − ∆t ˙u +(∆t) 2 2 ¨u を時刻 t= 0 で表現した式の右辺第 3 項に，時刻 t = 0 の瞬間の運動方程式から得られる ¨u(0)= 1 m { f (0) − c ˙u(0) − k u(0)} という関係を代入すれば u(−∆t) = { 1−∆τ 2 2 } u(0)− ∆τ (1 + β ∆τ) ˙u(0) ω + 1 2∆τ 2 f (0) k (9.53) を得る。これを用いれば時刻 t= 0 の u−を算定でき，初速も入力できる。一例として，図 9.13 と同じ問題を解 いた結果を図 9.19 に示した。破線が厳密解なのでほとんど重なって区別できないが，∆τ = ω∆t = 0.25 にした 場合の太い実線は十分な精度を持っていることがわかる。これに対し，ω∆t = 1 とした結果の細い実線は厳密 解からのずれが時間と共に次第に大きくなる。 図 9.20 には計算のフローチャートを示した。また図 9.21 には，図 9.8 で同定した板バネの基部を地震加速度 で揺らしたときの応答に対する数値解析結果（実線）と実験結果（破線）のうち，強制外力が大きい t ∼ 20 s までの部分だけを示した。入力地震波は，兵庫県南部地震の神戸海洋気象台で観測された加速度の半分である が，この程度の精度で数値解析ができる。このような積分方法を陽解法と呼んでいる。一般に中央差分による 数値積分の「安定性」のためには，時間ステップは系の固有周期の1/ πより小さくないといけない [8] とされる。 図 9.19 の最初の例では無次元化した時間τ で固有周期は 7 程度なので， ω∆t = 1 は安定性のためには十分で解
析は継続できているが，「正確性」についてはもっと小さい時間ステップが必要だということを示している。 そもそも微分 dt は∆t → 0 の極限だから当たり前だ。 演習問題 9-3 7. 減衰系に正弦波外力 f0 sin pt が作用したときの応答を Duhamel 積分を用いて求め，図 9.13 の結果と一 致することを示せ。 8. 同じ応答を差分法でも求めてみよう。 9. f (t) O f0 t1 t2 t 図 9.22 一時的な一定外力 図 9.22 のような，ある時間だけ作用する一定外力が減衰系に作用した ときの応答を Duhamel 積分を用いて求め，縦軸に ust ≡ f0 k としたと きのu(t) ust をとり，横軸にωt をとったときの図を描け。ただし， β = 0.02, ωt1 = 10, ωt2 = 30 として， 0 ≤ ωt ≤ 100 程度までを図示せ よ。この外力は Heaviside 関数を用いれば f (t)= f0 {H(t − t1)− H(t − t2)} と表すことができることを用い，また Duhamel 積分も式 (9.47) のように Heaviside 関数を含んだままの 単位衝撃応答を用いて計算すると，場合分けをせずに楽に間違い無く答を得ることができる。 10. 図 9.22 の外力に対する応答を差分法で求め， Duhamel 積分の答と比較せよ。 11. 時刻 t= τ に単位衝撃加速度が支点に与えられたときの相対変位の条件が ˙v(τ) = −1 であることを示せ。 (7) 衝撃外力の近似 数値解析で落石のような衝撃外力を扱うことがあるだろう。また部材の破断による動的応答を対象とする場 合に，破断を衝撃外力でモデル化するかもしれない。そこで上の演習問題の 9 番で得た解析解（この節では数 値解析はしていない）を使って，図 9.16 のような有限な時間幅ϵ → δt を持つ一定外力で衝撃応答をどのくら い近似できるか検討しておこう。上の問題の図 9.22 で t1 = 0, t2 = δt として f0 = F δt とすれば，時刻 t = 0 における力積 F の衝撃力を近似したことになる。このδt → 0 の極限が真の衝撃外力であるが，種々の大きさ を持つ有限なδt に対する応答解を図 9.23 に示した。まず左側の図に一点鎖線で描いたのが真の衝撃応答なの で，それと比べることによって，もし時間方向の誤差（ずれ）に対しても十分な精度を確保するためには，こ の減衰定数β = 0.02 の場合には ωδt = 0.01 程度より短い時間の近似衝撃外力を用いる必要があることがわか る。あるいは少し精度を緩めて，応力等の応答の最大値（振幅）だけを精度よく求めたい場合には，右側の図 からわかるようにωδt = 0.2 程度以下であればよさそうだ。例えば m = 1, 000 kg の岩塊が h = 5 m の高さか ら落ちて突然止まるくらいの衝撃力だとすると，その力積は F = m√2hg = 9.9 kN·s 程度になる。したがっ て，もしωδt = 0.01 程度の力の作用時間で，例えば周期が 0.5 s の構造への衝撃力を近似するのであれば，結 局 f0 = 12 MN の力を δt = 0.00080 s 間与えなければならないことを意味する。なお多自由度系の衝撃応答を検 討する場合には高周波（大きなω）側の応答も重要になる場合もあり，それも精度よく求めるためにはさらに 大きな力のさらに短時間の作用にする必要があるので注意が必要になる。 (8) 複素応答 — 不規則外力応答のための準備 次節の不規則応答を取り扱うために必要なのだが， exp (ipt) という外力に対する応答を考える。これを複素 応答と呼ぶ。つまり
