最新版を正確に読む場合には pdf ファイル をどうぞ。これは web 検索のための簡易旧版です。

• 10.1.1 非減衰自由振動
  • 10.1.1.1 運動方程式
  • 10.1.1.2 非減衰自由振動
  • 10.1.1.3 自由落下と等速円運動
  • 10.1.1.4 エネルギ
  
  
• 10.1.2 減衰自由振動
  • 10.1.2.1 減衰の原因
  • 10.1.2.2 粘性減衰
  • 10.1.2.3 固体摩擦減衰
  
  
• 10.1.3 強制振動
  • 10.1.3.1 正弦波外力による強制振動と共振曲線
  • 10.1.3.2 支点の正弦波状強制変位による応答
  • 10.1.3.3 単位衝撃応答
  • 10.1.3.4 任意外力に対する応答
  • 10.1.3.5 Duhamel積分は正しいか
  • 10.1.3.6 数値的な解法
  • 10.1.3.7 複素応答--不規則外力応答のための準備
  
  
• 10.1.4 不規則応答
  • 10.1.4.1 準備その1 -- Fourier変換
  • 10.1.4.2 準備その2 --波のスペクトル解析
  • 10.1.4.3 不規則外力に対する変位応答
  • 10.1.4.4 加速度応答と入力地震

10.1 1自由度系の振動

10.1.1 非減衰自由振動

10.1.1.1 運動方程式

この章は，東京大学の伊藤学名誉教授の1970年頃の 講義ノートを元にした。当時の教科書は文献[97]であったが， 多分元になったのは文献[87]の方だと思われる。まず1自由度系の 振動問題を用いて，振動論において最も重要な固有振動数と 共振という概念を説明する。多自由度系の場合には，さらに ベクトルとしての固有振動モードという概念を導入し，説明する。 そして梁を例とした連続体の振動の場合には，その固有振動モードが 関数になっていることを説明する。最終的には 有限要素法を用いることによって，連続体は見かけ上多自由度系として 近似的に取り扱うことができることを説明する。しかし，共振や 周波数応答関数の特性，さらに境界条件によってモードが変わること等が， 数値解析をするまでもなく予想できるような技術者になって欲しい。 数値解析は，単なるツールでしかなく，技術者の力学的直感の方が重要である。 なお，この文書では弾性範囲のみを対象とするので，強度設計のためには さらなる学習が必要である。

図 10.1: 1自由度系の振動

例えば，バネ定数$k$の線形バネに取り付けた1個の質点の運動を考えよう。 バネの質量は無視できるものする。図-10.1のように， 自然状態のバネの先端に質量$m$の質点を付けたところ運動し始めた。 そのときの質点の下方への変位（バネの伸び）を$w(t)$とし， 質点には外力$f(t)$も作用することができるものとする。 このとき，質点の運動方向（図では下方）に作用している力の 総和は，図の右端にある力のダイアグラムから

$$\sum\mbox{【力】} = f(t)+mg-k w(t)$$

である。したがってNewtonの法則（ $\mbox{力の総和}= \mbox{質量}\times\mbox{加速度}$）から，運動方程式は

$$\sum\mbox{【力】} = f(t)+mg-k w(t)=m \D*[2]{w(t)}{t} \quad\to\qquad\qquad m \ddot{w}(t)+k w(t)=mg+f(t)$$ (10.1)

となる。以下，上付きドットは時間による微分を表す。

もし$w(t)$の代わりに

$$u(t)\equiv w(t)-\dfrac{mg}{k}$$ (10.2)

を用いると， $\ddot{u}=\ddot{w}$なので，運動方程式(10.1)は

$$m \ddot{u}(t)+k u(t)=f(t)$$ (10.3)

となる。 この式(10.2)の$u(t)$は，図-10.1の 左から三番目の図に示したように，質点をバネにそぉーっと取り付けたあと 静的につり合った状態からの変位である。 この形の方程式は，このバネ質点系が摩擦の無い床の上で水平に 振動する場合の式に一致する。また，社会基盤構造の 振動問題では，既につり合っている状態からの動的な変動に興味が あることの方が多いので，以後，原則として自重（死荷重）の項は無視して， 式(10.3)を標準的な運動方程式とする。 ここで外力$f(t)$が作用していない振動のことを自由振動 と呼んでいる。 節題目にある「非減衰」については今はわからなくていい。

上の例は質点の併進運動を対象とした。次に振り子運動を考えてみよう。 振り子は重力場での運動なので重力加速度$g$を考慮することにする。 回転運動についてのNewtonの法則は，運動の回転角を$\theta(t)$としたとき， その回転角の向きに作用している外力モーメントの総和を$M(t)$とすれば

$$\sum\mbox{【モーメント】}= M(t)=J \ddot{\theta}(t)$$ (10.4)

である。ここに$J$は慣性モーメント と呼ばれ

$$J\equiv\int_V r^2 \rho \dint V$$ (10.5)

で定義される。ここに$\rho$は密度で，$r$は回転中心からの 距離，$V$は体積である。 これは併進運動の質量$m$の定義 $$\int_V\rho\dint V$$に対応した， 回転運動の物理パラメータである。

図 10.2: 振り子

ここで図-10.2の 振り子を考えよう。振り子の糸は質量が無く伸び縮みもしないものとする。 振り子の回転が反時計周りに$\theta(t)$だとすると， 支点回りにこの振り子に作用しているモーメントは重力の 影響のみで，それも反時計回りを正とすると

$$M(t)=-mg \ell \sin\left(\theta\left(t\right)\right)\simeq -mg \ell \theta(t)$$

だけである。ここでは揺れの大きさは小さく $\left\vert\theta(t)\right\vert\ll1$が 成立すると近似した。質量は振り子の先端にしか無いので， 慣性モーメントは$J=m \ell^2$となる。したがって，回転の 運動方程式(10.4)は

$$-mg \ell \theta(t)=m \ell^2 \ddot{\theta}(t) \quad\to\quad \ell \ddot{\theta}(t)+g \theta(t)=0$$ (10.6)

となる。これも形式的には，式(10.3)の標準的な運動方程式と 同じ形[($\ell$, $g$) $\leftrightarrow$ ($m$, $k$)]になっていることがわかる。

10.1.1.2 非減衰自由振動

では式(10.3)を用いて，外力$f(t)$が作用していない 状態の自由振動を求めてみよう。運動方程式は2階の常微分方程式なので， 唯一な解を得るためには，時刻$t=0$における状態， つまり初期条件を2式与える必要がある。それは， 最も単純な設定は，質点の初期の位置と 初期の速度（以下，初速と呼ぶ）を与えることであり，例えば

$$u(0)=u_0, \quad \dot{u}(0)=v_0$$ (10.7)

と与えられることにしよう。 式(10.3)を辺々$m$で割り，$m$も$k$も正の定数であることから

$$\omega\equiv\sqrt{\dfrac{k}{m}}$$ (10.8)

と定義すると，運動方程式は

$$\ddot{u}(t)+\omega^2 u(t)=0$$ (10.9)

となる。 この微分方程式の一般解を求めるために $u(t)=\exp\left(\mu t\right)$と置き， 運動方程式(10.9)に代入すると，$\mu$についての 特性方程式とその根は

$$\mu^2+\omega^2=0 \quad \to \quad \mu=\pm\mbox{i} \omega$$

と求められる。 したがって解は，Eulerの公式(10.62)（ $\exp \left(\pm\mbox{i}\omega t\right)\sim \sin\left(\omega t\right), \cos\left(\omega t\right)$）から

$$u(t)=A \sin\left(\omega t\right)+B \cos\left(\omega t\ri... ...^2+B^2}, \quad \zeta\equiv -\tan^{-1}\left(\dfrac{B}{A}\right)$$ (10.10)

と求められる。$C$は振幅 と呼ばれている。積分定数$A$と$B$，あるいは$C$と$\zeta$を 初期条件式(10.7)で 決定できれば，それが解になる。$A$と$B$を使った方の表現を初期条件に代入すると

$$u(0)=B=u_0, \quad \dot{u}(0)=\omega A=v_0 \qquad\to\qquad A=\dfrac{v_0}{\omega}, \quad B=u_0$$

と求められる。結局，運動は

$$u(t)=\dfrac{v_0}{\omega} \sin\left(\omega t\right)+ u_0 \cos\left(\omega t\right)$$ (10.11)

と表されることになる。$\zeta$は，$t=0$に0から 始まるsine関数 $\sin(\omega t)$の 初期位相のずれである。

図 10.3: 非減衰自由振動

図-10.3にその運動を示した。 図に示したように，元の位置に同じ速度で戻るまでの時間$T$を固有周期 （単位はs）と呼び，それは1秒間に$\omega$radの速さの運動が$2\pi$回るのに 必要な時間なので

$$T=\dfrac{2\pi}{\omega}$$ (10.12)

という関係になる。したがって式(10.8)で 定義される$\omega$を固有振動数 ^10.1と呼ぶ。 式(10.8)が併進運動の場合の固有振動数になるので， 式(10.3)と式(10.6)の運動方程式同士を 比較すれば明らかなように，振り子運動の固有振動数は

$$\omega=\sqrt{\dfrac{g}{\ell}}$$ (10.13)

となる。つまり振り子の長さ$\ell$を長くするとゆっくり動き， 短くすると速く動くことを示している。時計屋のいろいろな 振り子時計（本当に振り子で動いているのかねぇ?）の 動きを観察してみるといい。

10.1.1.3 自由落下と等速円運動

高校の物理で関数と微分積分を使ってはいけないことの是非はともかく， 記憶科目として習って覚えていたものの根拠を示しておこう。 式(10.1)でバネを取り去れば，それは自由落下 の運動方程式になるはずだ。 簡単のために外力$f(t)$も無いものとすると，それは

$$m \ddot{w}(t)=mg \quad\to\quad \ddot{w}(t)=g \eqno{(*)}$$

である。つまり，加速度が一定$g$であることから，等加速度運動である。 初期条件を質点の初期位置と初速で与えるために，式(10.7)に ならって

$$w(0)=w_0, \quad \dot{w}(0)=v_0$$

としておく。式($*$)を解くと

$$w(t)=\dfrac12 g t^2+c_1 t+c_2$$

が一般解なので，上の初期条件に代入すると$c_1=v_0$, $c_2=w_0$と 求められる。したがって，自由落下の質点の位置は下向きを正とすると

$$w(t)=\dfrac12 g t^2+v_0 t+w_0$$

になる。これが放物運動になっていることはすぐに理解できると思う。 もし重力場と直交する水平運動の場合には， この答の重力加速度$g$の項を無視すればいいので，それは等速運動になる。

また式(10.6)が無重力場（水平面内）の回転運動の 場合には， $\ddot\theta=0\to\dot\theta=\mbox{const.}$（ $\theta= \omega t+\mbox{const.}$）の等速度運動になり，それを等速円運動 と呼ぶ。 このとき質点の速度ベクトルは，$\theta$方向の 速さを $v_\theta\equiv r\omega$とすると

$$\fat{v}=v_\theta \fat{e}_\theta, \qquad \fat{e}_r=\fat{e}_... ...d \fat{e}_\theta=-\fat{e}_x \sin\theta+\fat{e}_y \cos\theta$$

である。ここに$\fat{e}_d$ ( $d=x,y,r,\theta$)は 各座標系の単位基底ベクトルである。 これから加速度ベクトルを求めると

$$\dot{\fat{v}}=v_\theta \dot{\fat{e}}_\theta =-v_\theta \d... ...+\fat{e}_y\sin\theta\right) =-v_\theta \dot\theta \fat{e}_r$$

となる。つまり加速度は$\fat{e}_r$の負の方向（回転中心方向）に生じていて， その大きさ$a$は

$$a=v_\theta \dot\theta=v_\theta \omega$$

である。このことから，向心力 が $F=ma=m v_\theta \omega=m r \omega^2 =m \dfrac{v_\theta^2}{r}$となるのである。 高校の教科書に載っている誘導過程は理解できないが，微分を使えば簡単だ。 しかし，第1著者は50歳になって初めてこれがわかったのだが， 物理的な意味は未だに説明できない。

10.1.1.4 エネルギ

さて，理由はとりあえず考えずに， 式(10.3)に$\dot{u}(t)$を乗じて時間積分してみよう。

$$\int m \ddot{u} \dot{u}\dint t+ \int k u \dot{u}\dint t=\int f(t) \dot{u}\dint t+\mbox{const.}$$

ここで，第1項は

$$K\equiv \int m \ddot{u} \dot{u}\dint t= \int\dfrac12 m \... ...t \left(\dot{u}\right)^2 =\dfrac12 m \left(\dot{u}\right)^2$$

と変形でき，第2項も同様に

$$U\equiv \int k u \dot{u}\dint t=\int k u\dint u=k \int u\dint u=\dfrac12 k u^2$$

となるので，上式は結局

$$K+U=\int f(t) \dot{u}\dint t+\mbox{const.}, \quad\mbox{つま... ...)^2+ \dfrac12 k u^2=\int f(t) \dot{u}\dint t+\mbox{const.}$$ (10.14)

が成立する。これはエネルギ保存則 であり，$K$は運動エネルギ であり，$U$は位置エネルギ である。右辺は外力の持つエネルギである。 どの状態を基準にするのか明記していないので，右辺には定数の項が残っている。 自由振動の場合の解の式(10.11)で，例えば初速が零の場合の 解を式(10.14)の各項に代入すると

$$K=\dfrac12 m \omega^2 u_0 \sin^2\omega t =\dfrac12 k u_0^2 \sin^2\omega t, \quad U=\dfrac12 k u_0^2 \cos^2\omega t$$

となり， $K+U=\mbox{const.}$であることがわかる。また

$$K\subsc{max}=U\subsc{max} \quad\mbox{あるいは}\quad \left\vert K\right\vert=\left\vert U\right\vert$$ (10.15)

であることもわかる。

図 10.4: 質量のあるバネ

この最後の結論を用いると，ちょっと 面白いことを求めることができる。 もしバネにも質量があり，単位長さ当たりの質量が$m_0$である 場合の，図-10.4に示した質量$m$の質点の振動を 考えてみよう。 伸びはバネの初期長さに比べて小さく無視できるとして， バネの長さを$\ell$とする。簡単のために，質点の振動による変位を

$$u(t)=A \sin\left(\omega t-\zeta\right)$$

とすると，図から明らかなように，バネの途中の任意点$x$の変位$v(x,t)$は

$$v(x,t)=\dfrac{x}{\ell} u(t)= A \dfrac{x}{\ell} \sin\left(\omega t-\zeta\right)$$

になると考えていいだろう。これを用いて全体の運動エネルギを求めると

$$K=\dfrac12 m\left(\dot{u}\right)^2+ \dfrac12\int_0^\ell m_... ...dfrac{m_0 \ell}{3}\right) \cos^2\left(\omega t-\zeta\right)$$

となる。一方，位置エネルギ（バネの全変形が蓄えている変形の エネルギ）はバネの質量とは無関係で

$$U=\dfrac12 k u^2=\dfrac12 k A^2 \sin^2\left(\omega t-\zeta\right)$$

となる。この結果を式(10.15)に代入すると

$$\dfrac12 \omega^2 A^2 \left(m+\dfrac{m_0 \ell}{3}\right)=... ...^2 \quad\to\quad \omega=\sqrt{\dfrac{k}{m+\dfrac{m_0\ell}{3}}}$$

と，固有振動数が求められる。 なぜか，バネの全質量$m_0\ell$の$\slfrac13$だけの貢献になる。

10.1.2 減衰自由振動

10.1.2.1 減衰の原因

さて，実際に振動する物体を観察すると，例えば空気抵抗等によって 振動はそのうち止まってしまうだろう。このような 振動を減衰振動と呼ぶ。減衰 の原因は，空気抵抗（速度の2乗に比例）だけではなく， 材料の中の何らかの非可逆変形でエネルギが 消費されてしまったり，構造と材料のどこかに摩擦が発生していたり， 種々考えられる。ここではそれをモデル化したもののうち， 最も基本的な2種類，つまり速度に比例した粘性抵抗がある場合と， 一定の力による固体摩擦抵抗がある場合の振動を説明する。

10.1.2.2 粘性減衰

図 10.5: 粘性減衰

バネが変位に比例した抵抗をするのに対して，粘性抵抗というのは速度に 比例した抵抗のことである。典型的なモデルは よくKelvin-Voigtモデル と呼ばれる系で，図-10.5に示した。 バネでない方の記号は，ダッシュポットと 呼ばれる装置で，これが粘性抵抗を受け持っている。 図の中央に描いたように，ある種のオイルが密封されたシリンダーに， 穴の開いたピストンがついているものを思い描いて欲しい。 ピストンをゆっくり動かすとオイルは穴をスムーズに 通り抜け，ほとんど抵抗が発生しない。しかし速くピストンを 動かそうとした場合には，オイルはそれほど速くは穴を通り抜けられず， これが抵抗力になる。このように，速度に比例した抵抗を 発生しているのがダッシュポットであり，その比例係数を$c$とする。 図の右端の力のダイアグラムのように，質点にはダッシュポットから 上向きに$c \dot{u}(t)$の抵抗を受けるので，運動方程式は

$$f(t)-k u(t)-c \dot{u}(t)=m \ddot{u}(t) \quad\to\quad m \ddot{u}(t)+c \dot{u}(t)+k u(t)=f(t)$$ (10.16)

となる。ここで

(10.17)

と定義すると，外力$f(t)$が作用していない自由振動の場合には， 上の運動方程式(10.16)は

$$\ddot{u}(t)+2 \beta \omega \dot{u}(t)+\omega^2 u(t)=0$$ (10.18)

と書くことができる。ここに$c$は粘性減衰係数 と呼ばれ，$\beta$は（粘性）減衰定数^10.2 と呼ばれる。

図 10.6: 過減衰と臨界減衰

式(10.18)の解を $\exp\left(\mu t\right)$と置き 代入すると，$\mu$についての特性方程式とその根が

$$\mu^2+2 \beta \omega \mu+\omega^2=0 \quad\to\quad \mu=\omega\left(-\beta\pm\sqrt{\beta^2-1}\right)$$ (10.19)

と求められる。 これに対応する解は$\beta$の大きさによって異なる挙動を 表現することになる。

1. まず$\beta>1$の場合には解は指数関数あるいは 双曲線関数であり，振動応答はせず，図-10.6の 左図のような挙動を示す。この現象を過減衰 と呼ぶ。機械部品にとっては必要な特性かもしれないが， 社会基盤構造ではこういう応答はほとんど期待できない。
2. 次に$\beta=1$の場合には，特性根が重根となるので， 解は $u=\left(A+Bt\right) \exp\left(-\omega t\right)$と なり，図-10.6の右図のような挙動を示す。 この現象を臨界減衰 と呼ぶ。式(10.17b)の定義から$\beta=1$のときの 減衰係数 $c\sub{cr}\equiv 2 m \omega=2\sqrt{m k}$を臨界減衰係数と 呼ぶ。
3. 最後に，一般的な社会基盤構造の場合には$\beta<1$である。 実際にはさらに$\beta\ll 1$と考えていいくらい小さく， 具体的には$\beta=0.02$とか$0.05$程度である。 そういう場合には，特性根は
   

   $$\mu=-\beta \omega\pm \mbox{i} \omega_d, \quad \omega_d\equiv \omega\sqrt{1-\beta^2}$$ (10.20)

   
   

   と置くことができる。この場合には，解は
   

   $$u(t)=\exp\left(-\beta \omega t\right) \left(A \sin\omega... ... \omega t\right)\right] \sin\left(\omega_d t-\zeta\right)$$ (10.21)

   
   

   となる。つまり，振幅$C$が指数関数的に小さくなりながら（上式の 鈎括弧部分）振動することに なる。図-10.7の左図にその例を示した。 振幅のピークを結んだ一点鎖線が指数関数的に小さくなる様子を 示しており，振幅はその比率で小さくなっていくが， 周期は一定のままの振動である。

図 10.7: 減衰振動

さて，この最後の$\beta\ll 1$の場合に限定して，以下の説明を 続けることにしよう。このような 減衰振動の場合の振動数は式(10.19)の$\omega_d$であるから， 周期$T$も

$$T=\dfrac{2\pi}{\omega_d}=\dfrac{2\pi}{\omega\sqrt{1-\beta^2}}$$ (10.22)

になる。 一般に前述のように，社会基盤構造の場合は減衰定数$\beta$は非常に小さく， 例えば鋼構造では0.02〜0.03，コンクリート構造では0.03〜0.05と 考えていいため，$\omega_d$と$\omega$の差は無視できるくらい小さく

$$\mbox{社会基盤構造物の場合は}\quad \omega\simeq\omega_d, \quad T\simeq\dfrac{2\pi}{\omega}$$ (10.23)

と近似してもいいのが普通である。 実際の社会基盤構造の減衰が粘性減衰であるとは限らないが， それでほぼ近似でき，また近似できると便利なので， 粘性減衰モデルを用いているのだろうと 思われる。では，実測した挙動から粘性減衰定数を求めることを 考えてみよう。実測データが図-10.7の右図のよう だったとしよう。$n$番目の振幅と$(n+j)$番目の振幅（図では$j=3$）を それぞれ$u_n$, $u_{n+j}$とすると， その振幅比（通常は$j=1$）を減衰比 ^10.3と呼び，それは理論的には

$$\dfrac{u_n}{u_{n+j}}= \dfrac{C \exp\left(-\beta\omega t\rig... ...right) =\exp\left(\dfrac{2\pi j\beta}{\sqrt{1-\beta^2}}\right)$$ (10.24)

となるので，定数になる。したがって，この関係式を用いて構造物の振動試験の 実測値を整理すれば，見かけ上の減衰定数$\beta$を逆算することができる。 また，この減衰比の対数を対数減衰率 と呼ぶ（$j=1$の場合で定義されているのが普通）が，それは

$$\delta\equiv \ln \dfrac{u_n}{u_{n+j}}=\dfrac{2\pi j \beta}{\sqrt{1-\beta^2}}$$ (10.25)

となり，これも定数なので，この関係式から

$$\beta=\dfrac{\delta}{\sqrt{\delta^2+4\pi^2 j^2}}$$ (10.26)

と，減衰定数を求めることができる。あるいは，$\beta\ll 1$が 成立するなら，この関係はもっと近似できて

$$\beta\simeq\dfrac{\delta}{2\pi j}$$ (10.27)

としてもいいことになる。

図 10.8: 実測データからのパラメータ同定

図-10.8には， 学生実験で得た1自由度系の加速度記録を×点で示した。 実験では，鋼の薄板で板バネを作り，その先端に重りを取り付け， 板バネを少し曲げた状態で手を離し，そのあとの自由振動の 重りの加速度を記録した。手を放すときの乱れがノイズとして 運動に影響を及ぼすので，初期状態から数周期経たあとの記録を， 自由振動していると判断して図示した。$\beta$は十分小さいと 仮定した上で，この記録の横軸をまたぐ15周期分の時間から 周期を$T=0.485$sと読み取り， 近似式(10.23)から固有振動数を$\omega= 12.95$ $\slfrac{\mbox{\scriptsize rad}}{\mbox{\scriptsize s}}$と同定した。 さらに，この図の範囲の正方向の振幅データを用いて， いろいろな$n$や$j$で式(10.25)と 近似式(10.27)を用いて算定したものを 平均して，$\beta=0.00385$と同定できた。 同定したパラメータを用いて求めた応答が実線で， ほとんど重なっていることがわかる。 実際の実験装置では，特に粘性減衰特性を有した装置を付けたわけでもなく， 減衰する理由は空気抵抗や板バネの取り付け部における何らかの エネルギ消失等だろうと予想されるが，巨視的な挙動を捉える限りは， この程度のシミュレーションが可能になっていることから，社会基盤構造の 振動モデルにおいても，粘性減衰モデルが用いられているのであろう。

10.1.2.3 固体摩擦減衰

固体摩擦のモデルを考える前に， 式(10.9)の運動方程式を変位と速度の平面上で 表すことを考えよう。このような平面を位相平面 と呼び，安定問題等でも用いられることがある。 まず運動方程式の加速度項は

$$\ddot{u}=\D*{\dot{u}}{t}=\D*{\dot{u}}{u}\D*{u}{t}=\dot{u}\D*{\dot{u}}{u}$$

と書き直してもいいから，運動方程式は

$$\dot{u}\D*{\dot{u}}{u}+\omega^2 u=0 \quad\to\quad \dot{u}\dint \dot{u}+\omega^2 u\dint u=0$$

となるので，これを項別に積分すれば

$$\left(\dot{u}\right)^2+\omega^2 u^2=\mbox{const.} \quad\mbo... ...quad \left(\dfrac{\dot{u}}{\omega}\right)^2+u^2=\mbox{const.}$$ (10.28)

と表すことができる。右辺の定数は初期条件で決まる値である。 誘導過程からも推測できるように，これはエネルギ保存則そのものでもある。 この式は$u$と$\dot{u}$を2軸とする位相平面上で楕円状の軌跡 を描くことを意味している。 あるいは，特にこの場合には$u$と $\dfrac{\dot{u}}{\omega}$を2軸にした 円軌道になる。

図 10.9: 固体摩擦減衰

さて，図-10.9に示したように，質点との 間に一定の動摩擦力$F=\mu' mg$（$\mu'$は動摩擦係数 ）が発生する床の上の振動^10.4を考えてみよう。 速度と逆の向きに摩擦が発生するので，$\dot{u}>0$のときの運動方程式は

$$-k u(t)-F=m \ddot{u}(t) \quad \to\quad m \ddot{u}(t)+k u(t)+F=0$$

となるのに対し，$\dot{u}<0$の場合には

$$-k u(t)+F=m \ddot{u}(t) \quad \to\quad m \ddot{u}(t)+k u(t)-F=0$$

となる。したがって，これをまとめると

$$m \ddot{u}(t)+k u + \mbox{Sgn}\left(\dot{u}\right) F=0$$ (10.29)

と書く^10.5ことができる。 ここにSgnは引数の符号を表す記号で

$$\mbox{Sgn}\left(\dot{u}\right)=\left\{\begin{array}{ll} +1 ... ...\ge 0 \\ -1 & \mbox{もし} \quad \dot{u}<0 \end{array}\right.$$ (10.30)

と定義した。初期条件は式(10.7)で与えられる。 この微分方程式を速度の符号で場合分けをして時々刻々解いていけば 解を求めることができるが，ここでは上述の位相平面を使ってみよう。

図 10.10: 固体摩擦減衰応答

この運動方程式(10.29)に相当する位相平面上の軌道を求めようとすると， まず

$$\dot{u}\dint \dot{u}+\omega^2 \left\{u +\mbox{Sgn}\left(\d... ...left\{u +\mbox{Sgn}\left(\dot{u}\right) \dfrac{F}{k}\right\}=0$$

としてもいいから，位相平面は

$$\left(\dfrac{\dot{u}}{\omega}\right)^2+ \left\{u +\mbox{Sgn}\left(\dot{u}\right) \dfrac{F}{k}\right\}^2= \mbox{const.}$$ (10.31)

と表すことができる。 つまり図-10.10に示したように， 速度の向きが変わる度に円の中心位置 $(\slfrac{\pm F}{k},0)$が突然移動する。 位相平面は，速度軸を上向きに変位軸を右向きに描くのが普通だが， 時刻歴応答を描き易いように，反時計回りに90度回転して描いた。 初期条件の状態から時計回りに回転しながら，変位軸をまたいだ半周毎に 軌道の中心が突然移動し，それに伴い振幅も小さくなっていき， 最後は $k \left\vert u\right\vert<F=\mu' mg$の状態^10.6で速度が零になる 時点（この位相平面の図の縦軸上）で止まることになる。

ちなみに，式(10.14)の誘導と同じアプローチで， 固体摩擦がある場合のエネルギ保存則を誘導しておこう。 式(10.29)の運動方程式に$\dot{u}(t)$を乗じて積分すると， 最初の2項は式(10.14)と同じ結論になる。第3項は

$$\int \mbox{Sgn}\left(\dot{u}\right) F \dot{u}\dint t =F \... ...t u\right\vert =F\times (\mbox{質点がたどった道のりの累積})$$

になる。これも高校で習って覚えた公式であろう。したがって， エネルギ保存則は，外力項も含めると

$$\dfrac12 m \left(\dot{u}\right)^2+ \dfrac12 k u^2+ \m... ...\right\vert\dint t= \int f(t) \dot{u}\dint t + \mbox{const.}$$

となる。

10.1.3 強制振動

10.1.3.1 正弦波外力による強制振動と共振曲線

ここまでは外力$f(t)$が作用していないものとして， 振動の最も基本的なパラメータである固有振動数と 減衰定数の意味について説明してきた。 実際には，地震や風や活荷重といった外力が作用した状況での振動が 問題になり，そのような外力に対して安全な社会基盤構造に なるように，固有振動数（や減衰定数）を設計する必要がある。 ここでは，いくつか代表的な外力が存在する場合の振動を対象として， この基本的な二つのパラメータと外力特性の関係について検討しよう。 外力$f(t)$が作用した振動を強制振動と呼ぶ。 一般的な問題として，式(10.16)を 解こう。あるいは$m$で除し，式(10.18)の形にすると

$$\ddot{u}(t)+2 \beta \omega \dot{u}(t)+\omega^2 u(t)=\dfrac{1}{m} f(t)$$ (10.32)

となる。

まず基本的な外力として正弦波状の外力を扱う。というのも，Fourier級数を 思い出せば，任意のほんんどの関数は三角関数で表示できると考えて いいだろうから， その基本的な外力として $f(t)=f_0 \sin pt$について考えておけば， ほとんどすべての外力に対する応答を求める基礎になると考えられるからだ。 つまり，式(10.32)は

$$\ddot{u}(t)+2 \beta \omega \dot{u}(t)+\omega^2 u(t)=g_0 \sin pt, \quad g_0\equiv \dfrac{f_0}{m}$$ (10.33)

となる。 一般解は，自由振動で求めた斉次解と外力に直接関係する特解との和になるから

$$u(t)=\exp\left(-\beta \omega t\right) \left(A \sin\omega_d t+B \cos\omega_d t\right)+u_p(t) \eqno{(*)}$$

と表される。減衰が存在する場合には，自由振動解は減衰解であり， 時間とともに減少していくような過渡応答 であるから，工学的に一番興味があるのは特解の方^10.7になる。 そこでここでは特解$u_p(t)$だけを求めてみよう。

数学で習った特解の求め方を用いるのが望ましいが，多分$\sin pt$, $\cos pt$が 解の候補になると予想されるので

$$u_p(t)=C \sin pt+D \cos pt$$

と置いて，式(10.33)を満足するような$C$, $D$を求めよう。 代入してsine, cosine毎に整理すると

$$\sin pt \left(-p^2 C-2 \beta \omega p D +\omega^2 C-... ...t \left(-p^2 D+2 \beta \omega p C +\omega^2 D\right)=0$$

となるので

$$\left(\begin{array}{cc} \omega^2-p^2 & -2 \beta \omega p... ...y}\right\}= \left\{\begin{array}{c} g_0 0\end{array}\right\}$$

を解けばいい。$2\times 2$の行列なので，簡単に逆行列が求められ

$$\begin{eqnarray*} \left\{\begin{array}{c} C D \end{array}\right\}&=& \dfrac{1}... ...{\vskip 1.2ex} -2 \beta \dfrac{p}{\omega} \end{array}\right\} \end{eqnarray*}$$

となる。したがって特解は

$$u_p(t)= \dfrac{\dfrac{f_0}{k}}{\left\{1-\left(\dfrac{p}{\ome... ...ight\} \sin pt -2 \beta \dfrac{p}{\omega} \cos pt \right]$$ (10.34)

と求められる。あるいは $u\sub{st}\equiv \dfrac{f_0}{k}$と定義すると

(10.35)

と書くこともできる。$u\sub{st}$は外力の最大値$f_0$が静的に作用したと 想定したときの静的な変位である。 したがって$M\subsc{d}$は，その静的な変位に対する動的な変位の割合を 表しているので，動的増幅率 と呼ばれる。$\alpha$は，減衰による位相の遅れを 表している。図-10.12の$t=0$付近を拡大すれば， この位相の遅れを実感できる。

図 10.11: 動的増幅率と位相遅れ

$M\subsc{d}$と$\alpha$を図-10.11に示した。 非減衰の$\beta=0$の場合には，$p=\omega$で $M\subsc{d}\to\infty$となり， つまり発振してしまう。この現象を共振 と呼び，$M\subsc{d}$のこの図の曲線を共振曲線と呼ぶことがある。 前述のように，社会基盤構造の$\beta$は非常に小さいため， このような共振に近い現象が生じる可能性がある。$\beta$が非常に小さい場合

$\dfrac{p}{\omega}\ll 1$の場合には，振幅は静的変位と同程度 であり，特解成分の位相のずれもほとんど無い。

$\dfrac{p}{\omega}\gg 1$の場合には，振幅は静的変位よりも 小さくなる一方，特解成分の位相のずれは$\pi$になる。

$\dfrac{p}{\omega}\simeq 1$の場合には，振幅は非常に大きくなり， 位相差はほぼ $\dfrac{\pi}{2}$程度になる。

となる。

図 10.12: 減衰系の正弦波強制振動

さて，社会基盤構造を設計したとしよう。 それは構造の質量$m$と剛性$k$を 設計したことに相当するから，動的に重要なパラメータの一つの 固有振動数$\omega$を設計したことになる。 もし，その構造に作用する外力の周波数$p$がその固有振動数に近い場合には， 共振して応答の振幅は非常に大きくなり，構造は破壊するかもしれない。 動的な設計というのは，簡単に言うと，そういう確認をすることである。

では一例として，過渡応答も含めて，零初期条件つまり$u_0=0$, $v_0=0$の場合の 強制振動応答を求めておこう。 式($*$)に式(10.34)を代入し，それを零初期条件に代入すると， 二つの積分定数を

$$\dfrac{A}{u\sub{st}}=\dfrac{\beta}{\sqrt{1-\beta^2}}\dfrac{B... ...t)^2\right\}^2 +4 \beta^2 \left(\dfrac{p}{\omega}\right)^2}$$

と求めることができる。上式の$A$の右辺に，$B$の結果を 代入すれば解が求められる。例として $\dfrac{p}{\omega}=0.9$で$\beta=0.05$の ときの解を描いたのが図-10.12である。 実線がその応答で，零初期条件で始まっていることがわかる。 破線は，その応答の中の$\omega_d$で振動する成分， つまり過渡応答の成分である。減衰定数が大きめなので過渡応答の 振幅は速く小さくなっていく。この成分が無くなってしまえば， 実線の振幅は，式(10.35b)の$M\subsc{d}$になり， 周波数が$p$の振動のみが残ることになる。 図中の水平の2本の破線が，この$M\subsc{d}$の値である。 共振点の90%程度の外力周波数でも， 静的変位の5倍もの大きな動的な変位が生じていることが， この応答図からもわかる。 また，過渡応答も含めた全振幅が動的増幅率より大きくなることも確認できる。

1. 自分が使い慣れたソフトウェアを用いて，減衰振動の応答 曲線図-10.7を描いてみよ。最近はソフトウェアが非常に 高性能になったため，スプライン関数等で滑らかな曲線を描くことが 容易にできるが，実はそこには大きな落とし穴がある。 スプライン関数の性質を知らずに安易に使うと，とんでもない結果に なってしまうので注意して欲しい。こういう周期的な曲線を描く場合， 一番安全で確実なのは，1周期当たり20個以上のデータ点を直線で結ぶことである。
2. 自分が使い慣れたソフトウェアを用いて，共振曲線図-10.11を 描いてみよ。これも試しにスプライン関数で描いてみよ。 データ数が少ない場合には，とても変な曲線になることがある。
3. 零初期条件$u_0=0$, $v_0=0$の質点に，正弦波外力 $f(t)=f_0 \sin pt$が 作用したときの応答を求め，図-10.12と 同じように図示せよ。$\beta=0.02$として， $\dfrac{p}{\omega}$が0.2, 0.5, 0.97, 1.1, 1.5等の場合の応答を求め， $\pm M\subsc{d}$の線も 描いた上で， $\left\vert u\right\vert>u\sub{st} M\subsc{d}$が どのような状況で発生するか比較してみよ。
4. 非減衰$\beta\equiv 0$で外力が $f(t)=f_0 \sin\omega t$のとき， つまり$p=\omega$の共振点では，式(10.35b)の動的 増幅率は $M\subsc{d}\to\infty$になることから， 特解が $u_p(t)\sim \sin pt=\sin\omega t$ではないことは明らかである。 零初期条件で，その解を求めよ。 数学が不得意でも，共振で発散することを念頭に置けば解は容易に求められるし， 図示してその発散の様子も確認せよ。

10.1.3.2 支点の正弦波状強制変位による応答

図 10.13: 支点変位に対する応答

さて，地震に対する応答のことを思い描くと， 構造物自体に外力が作用している前節のような状況ではなく， バネの支点が強制的に変位させられた状況での振動のように見える。 ここでは，そのような状況のモデルを考えてみよう。図-10.13に 示したように，支点が水平方向に強制的に$w(t)$だけ変位させられた 場合の質点の水平運動$u(t)$を考えよう。 バネとダッシュポットは，支点の運動との相対変位と相対速度に 比例した抵抗力しか生まないから，運動方程式は

$$m \ddot{u}(t)+c \left\{\dot{u}(t)-\dot{w}(t)\right\} +k \left\{u(t)-w(t)\right\}=0$$ (10.36)

となる。Newtonの法則の慣性項は絶対加速度$\ddot{u}(t)$で表されているので， 第1項はこれまでと同じである。さてここで，支点の変位に対する 相対的な質点の変位

$$v(t)\equiv u(t)-w(t)$$ (10.37)

を用いて，上の運動方程式(10.36)を書き直すと

$$m \ddot{v}(t)+c \dot{v}(t)+k v(t)=-m \ddot{w}(t)$$ (10.38)

と書くことができる。この式と式(10.16)を比較すれば 明らかなように，支点の強制変位に対する相対的な応答$v(t)$は， 見かけ上

$$f(t)=-m \ddot{w}(t) =-\mbox{(設計しようとしている構造の質量)}\times \mbox{(地盤の地震入力加速度)}$$ (10.39)

の外力が直接質点に作用したときの強制外力応答と同じであることがわかる。 この式(10.39)は非常に重要な概念であり， 耐震設計において用いられている震度法 の根拠になっている。つまり，過去に観測された大地震等 の地震加速度を参考にして設定された設計入力地震加速度に， 設計されようとしている社会基盤構造の質量を乗じることによって， その構造に作用させるべき地震外力を算定し，それを作用させたときの 構造の抵抗を照査することによって，耐震設計ができることを示している。

図 10.14: 支点変位応答の動的増幅率

基本的な例として，支点変位が正弦波状で

$$w(t)=w_0 \sin pt$$ (10.40)

となる場合を上式(10.38)に代入すると，運動方程式は

$$m \ddot{v}(t)+c \dot{v}(t)+k v(t)=m p^2 w_0 \sin pt$$

となる。したがって，前節の外力の振幅を $f_0=m p^2 w_0$と置き換えれば， そのまま前節の解を流用できるので， 特解は式(10.35)とほぼ同様になり

$$v_p(t)=w_0 M_w \sin\left(p t-\alpha\right), \quad M_w=\left(\dfrac{p}{\omega}\right)^2M\subsc{d}$$ (10.41)

となる。位相の遅れ$\alpha$は式(10.35c)と同じになる。 式(10.41)の動的増幅率$M_w$を図-10.14に 示した。つまり

$\dfrac{p}{\omega}\gg 1$の 場合には， $\left\vert v\right\vert\simeq\left\vert w\right\vert$となるが， 特解成分の位相差は$\pi$になる。これは柔構造の応答に相当している。

$\dfrac{p}{\omega}\ll 1$の場合には， 相対変位はほぼ零になっているが， $M\subsc{d}\simeq 1$なので $v_p\simeq w_0\left(\slfrac{p}{\omega}\right)^2 \sin\left(pt-\alpha\right)$と なり， $\left\vert v_p\right\vert\simeq \slfrac{1}{\omega^2}\left\vert\ddot{w}\right\vert$である。 したがって，バネの抵抗は $k\left\vert v_p\right\vert=m\left\vert\ddot{w}\right\vert$となり， 直接地震加速度に対応する力で抵抗しなければならない。これが剛な 構造の応答に相当している。

ということがわかる。

ところで，地震計 というのは一体どうなっているんだろう。 地面と一緒に動いている地震計が，自分自身の運動を記録できるのを 奇妙には感じないだろうか。 この疑問に答えてくれるのが，ここで求められた解である。 つまり，地震計はここで定義した相対変位$v(t)$を記録紙に書き込んでいる のだから， 相対変位$v(t)$の挙動を考えれば，地震計の原理が理解できるはずだ。 例えば，地震計の$m$と$k$を適切に設計して，観測したい地震動の 振動数$p$に対して $\dfrac{p}{\omega}\gg 1$になるように製作したとしよう。 すると，上で求められた解から

$$\left\vert v_p\right\vert\simeq w_0, \quad \alpha=\pi$$

になる。つまり，地面の変位応答を記録することができる。 これと対照的に地震計の$m$と$k$を調整して $\dfrac{p}{\omega}\ll 1$に なるように製作したとしよう。この場合には $M\subsc{d}\simeq 1$であることから， 式(10.41)から

$$\left\vert v_p\right\vert\simeq w_0 \left(\dfrac{p}{\omega}\right)^2, \quad \alpha=0$$

となっていることがわかる。つまり式(10.40)の微係数と比較して

$$\left\vert v_p\right\vert\simeq \left\vert\ddot{w}(t)\right\vert$$

と考えていいことになる。 したがって，地面の加速度応答が記録できるのである。ちょっと不思議。

10.1.3.3 単位衝撃応答

次に，周期外力とは正反対の特性を持つだろうと考えられる衝撃外力が 作用したときの振動について考えてみよう。 まず，衝撃そのものをどのようにモデル化するかが問題になるだろう。 いま質点$m$に，時刻$t=\tau$に，単位量の衝撃外力$f(t)$を作用させたとする。 衝撃なので，非常に短い時間$\epsilon$だけ作用しているとしよう。 そのときの質点の 運動方程式を $t=\tau-\slfrac{\epsilon}{2}$から $t=\tau+\slfrac{\epsilon}{2}$まで 積分すると，まず慣性項は

$$\int_{\tau-\epsilon/2}^{\tau+\epsilon/2} m \ddot{u} \dint t... ...on}{2}\right)- m \dot{u}\left(\tau-\dfrac{\epsilon}{2}\right)$$

となる。一方，バネの抵抗の項は位置エネルギの短時間の変化分であり， 粘性抵抗も短時間の速度変化になるので， 次の式のように$\epsilon\to 0$の極限では零，つまり

$$\int_{\tau-\epsilon/2}^{\tau+\epsilon/2} k u \dint t \to 0,... ...=c \left\{u(\tau+\epsilon/2)-u(\tau-\epsilon/2)\right\} \to 0$$

である。したがって，運動方程式を積分したものは結局

$$m \dot{u}\left(\tau+\dfrac{\epsilon}{2}\right)- m \dot{u}\... ...}\right)= \int_{\tau-\epsilon/2}^{\tau+\epsilon/2} f(t)\dint t$$ (10.42)

となる。この式は，実は運動量保存則 であり，右辺は衝撃外力$f(t)$による力積である。

図 10.15: 単位衝撃外力

この力積が単位量，つまり1のときの外力$f(t)$を単位衝撃外力 と呼ぶことにしよう。例として挙げたのが図-10.15である。 このような関数$f(t)$は，集中せん断外力をモデル化したときの 式(4.55)のDiracのデルタ関数 $\delta(x-a)$と同じであり， それは

$$\int_b^c \psi(x) \delta(x-a)\dint x = \left\{ \begin{array}... ...發��\quad a<b \quad\mbox{あるいは}\quad c<a \end{array}\right.$$ (10.43)

で定義される。つまり

$$f(t)=\delta\left(t-\tau\right) \quad \to \quad \mbox{右辺の... ...silon/2}^{\tau+\epsilon/2} \delta\left(t-\tau\right) \dint t=1$$

となる。この質点が時刻$t=\tau$まで静止していたとすると，$\epsilon\to 0$の 極限において式(10.42)の左辺第2項は零になるので

$$m \dot{u}(\tau) - 0 = \int_{\tau-\epsilon/2}^{\tau+\epsilon... ...tau\right) \dint t=1 \quad\to\quad \dot{u}(\tau)=\dfrac{1}{m}$$

となる。したがって，1質点系に時刻$t=\tau$に単位衝撃外力を作用させた ときの初期値問題は$t\ge \tau$において

$$m \ddot{u}(t)+c \dot{u}(t)+k u(t)=0, \quad u(\tau)=0, \quad \dot{u}(\tau)=\dfrac{1}{m}$$ (10.44)

で与えられることになる。もちろん$t<\tau$における状態は$u(t)\equiv 0$である。

式(10.44)の一般解は

$$u(t)=\exp\left(-\beta \omega t\right) \left( A \cos \omega_d t + B \sin \omega_d t\right)$$

なので，これを初期条件に代入すると

$$\begin{eqnarray*} 0&=&\exp\left(-\beta \omega \tau\right) \left( A \cos \o... ..., \left( A \cos \omega_d \tau + B \sin \omega_d \tau\right) \end{eqnarray*}$$

を満足しなければならない。これを整理すると

$$\left(\begin{array}{cc} \cos\omega_d \tau & \sin\omega_d ... ...begin{array}{c} 0 \dfrac{1}{m\omega_d} \end{array}\right\}$$

となるので，これを解くと

$$\left\{\begin{array}{c} A B \end{array}\right\}=\exp\lef... ... -\sin\omega_d \tau \cos\omega_d \tau \end{array}\right\}$$

のように積分定数が求められる。したがって解は

$$u(t)=u\subsc{i}(t;\tau)=u\subsc{i}(t-\tau)\equiv\left\{\begi... ...d\left(t-\tau\right)\right\} & \quad t>\tau \end{array}\right.$$ (10.45)

と求められる。これを単位衝撃応答 と呼んでいる。衝撃が与えられた時刻$\tau$も関数の引数に 明示するのが普通であり，また式の形から明らかなように， それは$(t-\tau)$のみの関数であることから， 左から三つ目の表現がよく用いられる。

図 10.16: Heaviside関数

さて，式(10.45)に示したように， 時刻$t$を$\tau$の前後に場合分けして二つの定義をするのはとても面倒だし， 使い難いので，これを一つの式にまとめてみたい。そのために

$$H(t-\tau) = \left\{\begin{array}{rl} 0, & t<\tau \\ 1, & t>\tau \end{array}\right.$$ (10.46)

$$\quad\mbox{少し難しく，やや正しく書くと}\quad \delta(t-\tau)\mathop{=}^{\mbox{\scriptsize d}}\D*{H(t-\tau)}{t}$$

と定義されるHeaviside関数（階段関数） を導入してみよう。 なお上式の $\delta(t-\tau)$は式(10.43)で定義したDiracの デルタ関数で，いずれも正しくは関数ではなく超関数であるから， 最後の式の等号 $$\mathop{=}^{\mbox{\scriptsize d}}$$は 超関数のための特殊な等号である。 さて，このHeaviside関数の性質を利用すれば，任意時刻$t>0$の 単位衝撃応答は

$$u\subsc{i}(t-\tau)= \dfrac{1}{m \omega_d} \exp\left\{-\b... ..., H(t-\tau) \sin\left\{\omega_d\left(t-\tau\right)\right\}$$ (10.47)

と書くことができる。非減衰の場合は$\beta\equiv 0$なので

$$u\subsc{i}(t-\tau)= \dfrac{1}{m \omega} H(t-\tau) \sin\left\{\omega\left(t-\tau\right)\right\}$$ (10.48)

となる。次の節以降で，この表現はとても便利だということがわかる。 なお，初期値問題の初期条件式(10.44)から 明らかなように，単位衝撃応答は時刻$t=\tau$における零初期条件の 解ではないことには十分注意する必要がある。

10.1.3.4 任意外力に対する応答

では，大きさが単位量ではなく，図-10.17の 左上の図のように，面積が$F$の衝撃外力に対する応答が

$$u(t)=F u\subsc{i}(t-\tau)$$

となることは，誰も疑問には思わないだろう。 それならば，その右の図のように，三つの異なる大きさの 衝撃外力に対する応答が

$$u(t)=\sum_{n=1}^3 F_n u\subsc{i}\left(t-\tau_n\right)$$

になるのも納得してもらえると思う。いや，任意時刻$t$の 応答は場合分けをしないといけないから，三つ共足してしまっては いけないのではないかと思う人も いるかもしれないが，例えば上式の総和の第2項は

$$\dfrac{F_2}{m \omega_d} \exp\left\{-\beta \omega \left... ...(t-\tau_2) \sin\left\{\omega_d\left(t-\tau_2\right)\right\}$$

であり，Heaviside関数が含まれているから$t<\tau_2$の間は無いのと同じである。 したがって，任意の時刻$t$に対する上の解の表現には三つ共 足しておいていいのである。Heaviside関数は便利でしょ。

図 10.17: 任意外力

それでは図-10.17の下の図のような 任意の関数$f(t)$で表された外力に対する応答も，同様に，単位衝撃応答を 用いて表現できないであろうか。 考え易くするために，図のように，$\Delta\tau$ずつの短冊に 切って考えると，例えば網掛けした部分は，面積 $\Delta\tau f(\tau)$の 衝撃として取り扱えばいいことがわかる。したがって， この網掛け部分による応答は

$$\Delta\tau f(\tau) u\subsc{i}(t-\tau)$$

でいいから，時刻$t=0$以降にしか外力が作用しなかったとすれば

$$u(t)=\sum_{\tau=0}^\infty \Delta\tau f(\tau) u\subsc{i}(t-\tau)$$

がその応答になることは理解できると思う。したがって $\Delta\tau\to 0$の 極限をとることによって， $$\sum \Delta\tau \to \int \dint\tau$$に なるので， 任意外力$f(t)$に対する応答は

$$u(t) \to \int_{0}^\infty \dint \tau f(\tau) u\subsc{i}(t-\... ...ad u(t)=\int_{0}^\infty u\subsc{i}(t-\tau) f(\tau) \dint \tau$$ (10.49)

と表現できる。この表現をDuhamel積分 と呼んでいる。この解は$t=0$における零初期条件の解になっている。 このように，単位衝撃応答 $u\subsc{i}(t-\tau)$は，$f(\tau)$が 全体応答に及ぼす影響を示していることから， 静的な構造力学の影響線と数学的には同じものである。 いわゆるGreen関数 である。

例として，図-10.12で示した$f_0 \sin pt$に 対する零初期条件の応答を求めてみよう。減衰解は面倒なので， 非減衰の場合を例にすると，同図の横にある解から，正解は

$$\dfrac{u(t)}{u\sub{st}}=\dfrac{\omega^2}{\omega^2-p^2} \sin... ...ga^2-p^2} \sin \omega t, \quad u\sub{st}\equiv\dfrac{f_0}{k}$$

である。式(10.49)に $f(\tau)=f_0 \sin p\tau$であることと 式(10.48)の単位衝撃応答とを代入すると

$$u(t)=\dfrac{f_0}{m \omega} \int_0^\infty \sin p\tau H(t... ... \omega} \int_0^t \sin p\tau \sin \omega(t-\tau) \dint \tau$$

となる。Heaviside関数から$t>\tau$なので，$\tau$の積分範囲は$t$までになる。 さらに演算すると

$$u(t)=\dfrac{f_0}{m \omega} \left( \sin\omega t \int_0^t \... ...mega t \int_0^t \sin p\tau \sin \omega\tau \dint \tau \right)$$

となるので，この積分を実行すれば解が求められる。 これも面倒なので割愛するが，最終的な結果は上の正解と一致する。 計算では$t$と$\tau$の区別を決して間違わないように!

10.1.3.5 Duhamel積分は正しいか

外力の連続関数$f(t)$を短冊に切ったことに不安を覚える人のために 証明しておこう。 なおこの節は，アメリカ合州国Illinois州Evanston, Northwestern大学Olmstead教授(1980年頃当時） の`Differential Equations of Mathematical Physics'の 講義ノートを参考にした。解きたい問題は

$$m \ddot{u}+c \dot{u}+k u=f(t), \quad u(0)=0, \quad \dot{u}(0)=0 \eqno{(a)}$$

である。一方，$t=\tau$に単位衝撃外力が作用した場合の 単位衝撃応答 $u\subsc{i}(t-\tau)$は 次式を満足する。

$$m \ddot{u}\subsc{i}+c \dot{u}\subsc{i}+k u\subsc{i}=\delt... ...quad u\subsc{i}(0)=0, \quad \dot{u}\subsc{i}(0)=0 \eqno{(b)}$$

ここに $\delta(t-\tau)$は式(10.43)で 定義したDiracのデルタ関数である。

さて目的は，$u(t)$を $u\subsc{i}(t-\tau)$で表すことだが， その補助のために，単位衝撃問題式($b$)に 対する随伴問題 と呼ばれる問題を定義する。 初期値問題は自己随伴系ではないことから， 随伴問題は時刻$t=\tau_1$に単位衝撃外力が 作用した時の応答 $u\subsc{i}^\ast(t-\tau_1)$に関する終局値問題になり， 次式で与えられる。

$$m \ddot{u}^\ast\subsc{i}-c \dot{u}^\ast\subsc{i}+ k u^\a... ...st(\infty)=0, \quad \dot{u}^\ast\subsc{i}(\infty)=0 \eqno{(c)}$$

ここで，運動方程式の減衰項の符号が負になっていることと，$t=\infty$での 終局条件が与えられていることに注意する。 この随判問題と元の問題の「関わりあい（仮想仕事）」を考えるために， 式($a$)の両辺に $u\subsc{i}^\ast$を乗じて時間に 関して$0$から$\infty$まで積分する。 物理的には仮想仕事を求めていることになる。

$$\int_0^\infty u\subsc{i}^\ast(t-\tau_1) \left( m \ddot{u... ...dint t = \int_0^\infty u\subsc{i}^\ast(t-\tau_1)\;f(t)\dint t$$

この左辺を部分積分すると

$$\mbox{左辺} = \left(u\subsc{i}^\ast m \dot{u}-\dot{u}\su... ...c{i}^\ast-c \dot{u}\subsc{i}^\ast+ k u^\ast \right) \dint t$$

となる。この第1項に初期条件式($a$)と 終局条件式($c$)を代入し，第2項の被積分関数の中の括弧部分に 式($c$)の運動方程式を代入すれば

$$= \int_0^\infty u \left( m \ddot{u}\subsc{i}^\ast-c \dot{... ...) \dint t = \int_0^\infty u \delta(t-\tau_1)\dint t=u(\tau_1)$$

となる。 最後の等号では，式(10.43)のデルタ関数の定義を用いた。 結局上式からは

$$u(\tau_1) = \int_0^\infty u\subsc{i}^\ast(t-\tau_1) f(t) \dint t$$

という関係を得る。あるいは$t$と$\tau_1$とを入れ替えると

$$u(t) = \int_0^\infty u\subsc{i}^\ast(\tau_1-t) f(\tau_1) \dint \tau_1 \eqno{(d)}$$

と表される。目的に一歩近付いているように見えるが， まだ``Close but no cigar!''である。

次に，同様の「関わりあい（仮想仕事）」を， 単位衝撃の問題と随判問題との間で考えると

$$\int_0^\infty u\subsc{i}^\ast(t-\tau_1) \left( m \ddot{... ...u\subsc{i}^\ast(t-\tau_1) \delta(t-\tau) \dint t \eqno{(e)}$$

となるから，上と同様に左辺を部分積分して，初期・終局条件と 微分方程式を代入した上で，デルタ関数の定義を用いれば

$$\mbox{左辺} = \left(u\subsc{i}^\ast m \dot{u}\subsc{i} -... ...t-\tau) \delta(t-\tau_1) \dint t = u\subsc{i}(\tau_1-\tau)$$

となる。一方，上式($e$)の右辺もデルタ関数の定義を 用いれば $u\subsc{i}^\ast(\tau-\tau_1)$となるから，結局，式($e$)からは

$$u\subsc{i}(\tau_1-\tau) = u\subsc{i}^\ast(\tau-\tau_1) \eqno{(f)}$$

という関係を得る。つまり随判問題は，元の衝撃問題の観測時刻と載荷時刻を 入れ替えただけの問題であることがわかる。 静力学の「相反定理 （式(3.73)や式(4.62)）」と同じである。

最終的にこの式($f$)を式($d$)に代入すれば

$$u(t) = \int_0^\infty u\subsc{i}(t-\tau_1) f(\tau_1) \dint \tau_1$$

となり，$\tau_1$はどんな記号を用いてもいい積分変数に過ぎないから

$$u(t) = \int_0^\infty u\subsc{i}(t-\tau) f(\tau) \dint\tau$$ (10.50)

と，Duhamel積分が成立することを示すことができた。

10.1.3.6 数値的な解法

社会基盤構造の設計に関する実務では，後述のように有限要素法を 用いて，コンピュータによって動的解析をする場合がある。 そのときの一つの解析法に，微分を差分 と呼ばれる方法で近似し，微分方程式を直接解く手法がある。 ここでは中央差分[8]の考え方を簡単に述べておく。 それは，十分小さな時間ステップ$\Delta t$毎に微分方程式を解く手法である。 表示を簡単にするために

$$u^-\equiv u(t-\Delta t), \quad u\equiv u(t), \quad u^+\equiv u(t+\Delta t)$$

と定義したとき，時刻$t$の速度と加速度は

$$\dot{u}(t)=\dfrac{1}{2\Delta t}\left(u^+-u^-\right), \quad \... ...}= \dfrac{1}{\left(\Delta t\right)^2} \left(u^-2u+u^+\right)$$ (10.51)

で近似できるとするのである。これを運動方程式(10.16)に代入すると

$$\dfrac{m}{\left(\Delta t\right)^2} \left(u^-2u+u^+\right) +\dfrac{c}{2\Delta t}\left(-u^-+u^+\right)+ku=f(t)$$

となるので，整理すると

$$u^+=\dfrac{1}{\left(\dfrac{1}{\Delta\tau^2}+\dfrac{\beta}{\D... ...ad \tau\equiv \omega t, \quad \omega\equiv\sqrt{\dfrac{k}{m}}$$ (10.52)

という関係が求められる。つまり，時刻$t-\Delta t$と$t$の変位$u^-$と$u$が わかっていれば，この式から時々刻々の解を順に求めることができる。

図 10.18: 差分法で強制振動を解く

ただし，一番最初のステップで$u^-$が何なのか定義できていない上に， 初期条件のうちの初速が入力できていない。そこで 例えば時刻$t=0$だけ後退差分で

$$\dot{u}(0)=\dfrac{u(0)-u^-}{\Delta t}= \dfrac{\mbox{初期位置}-u^-}{\Delta t}=\mbox{初速}$$

から$u^-$を概算することも可能だが，差分の取り方に整合性が無くなる。 そこでこれについては， 式(10.51)の二つの式から$u^+$を消去して得る

$$u^-=u-\Delta t \dot{u}+\dfrac{\left(\Delta t\right)^2}{2} \ddot{u}$$

を時刻$t=0$で表現し，さらに$t=0$の瞬間の運動方程式から得る

$$\ddot{u}(0)=\dfrac{1}{m} \left\{f(0)-c \dot{u}(0)-k u(0)\right\}$$

という関係を上の式の右辺第3項に代入すれば

$$u(-\Delta t)=\left\{1-\dfrac{\Delta\tau^2}{2}\right\} u(0) ... ...{\dot{u}(0)}{\omega} +\dfrac12 \Delta\tau^2 \dfrac{f(0)}{k}$$ (10.53)

となることから， 時刻$t=0$の$u^-$を算定でき，初速も入力できる。$f(0)=0$で零初期条件であれば， 結局 $u(-\Delta t)=0$でいい。 一例として，図-10.12と 同じ問題を解いた結果を図-10.18に示した。 ほとんど重なって見えないが， 図中の破線が厳密解なので， $\Delta\tau=\omega\Delta t=0.25$にした太い実線は 十分な精度を持っていることがわかる。これに対し， $\omega\Delta t=1$とした 結果の細い実線は厳密解からのずれが時間と共に次第に大きくなる。

図 10.19: 中央差分のフローチャート
図 10.20: 地震応答の例

図-10.19には計算のフローチャートを 示した。また図-10.20には，図-10.8で 同定した板バネの基部を地震加速度で揺らしたときの数値解析結果（実線）と 実験結果（破線）のうち，強制外力が大きい時間帯の部分だけを示した。 入力地震波は，兵庫県南部地震の神戸海洋気象台で観測された 加速度の半分であるが，この程度の精度で数値解析ができる。 このような積分方法を陽解法 と呼んでいる。一般に中央差分による数値積分の「安定性」のためには， 時間ステップは周期の $\slfrac{1}{\pi}$より小さくないと いけない[8]とされる。図-10.18の 最初の例では，無次元化した時間$\tau$で 周期は7程度なので， $\omega\Delta t=1$は安定性には十分で， 解析は継続できているが，「正確性」についてはもっと小さい 時間ステップが必要だということを示している。

1. 減衰系に正弦波外力$f_0 \sin pt$が作用したときの応答を，Duhamel積分を 用いて求め，図-10.12の結果と一致することを示せ。
2. 図-10.21のような， ある時間だけ作用する一定外力が減衰系に作用したときの応答を，Duhamel積分を 用いて求め，縦軸に $u\sub{st}\equiv \dfrac{f_0}{k}$と したときの $\dfrac{u(t)}{u\sub{st}}$をとり， 横軸に$\omega t$をとったときの図を描け。ただし，$\beta=0.02$, $\omega t_1=10$, $\omega t_2=30$として， $0\le \omega t \le 100$程度までを 図示せよ。この外力はHeaviside関数を用いれば
   
   $$f(t)=f_0 \left\{H(t-t_1)-H(t-t_2)\right\}$$
   
   

   と表すことができることを用い，またDuhamel積分も 式(10.47)のようにHeaviside関数を含んだままの 単位衝撃応答を用いて計算すると，場合分けをせずに，楽に 間違い無く答を得ることができる。

3. 図-10.21の外力に対する応答を差分法で 求め，Duhamel積分の答と比較せよ。

10.1.3.7 複素応答--不規則外力応答のための準備

図 10.21: 一時的な一定外力

次節の不規則応答を取り扱うために 必要なものだが， $\exp\left({\mbox{i}pt}\right)$と いう外力に対する応答（複素応答と呼ぶ）を考えてみる。あるいは

$$f(t)=f_0 \exp\left(\mbox{i}pt\right) =f_0 \left(\cos pt+\mbox{i} \sin pt\right)$$ (10.54)

であるから，正弦波応答と余弦波応答を同時に求めようとしていると考えてもいい。 これに対する応答も同様に $u(t)=\phi \exp\left(\mbox{i}pt\right)$と 置くことができそうだ。$\phi$は 複素数である。これを運動方程式(10.16)に代入すると

$$\left(-p^2+2\mbox{i} \beta \omega p+\omega^2 \right) \ph... ...\mbox{i}pt\right) =\dfrac{f_0}{m} \exp\left(\mbox{i}pt\right)$$

となるので

$$\phi= f_0 {\cal H}(\mbox{i}p), \quad {\cal H}(\mbox{i}p)\eq... ...\right\} +2 \mbox{i} \beta \left(\dfrac{p}{\omega}\right)}$$ (10.55)

と求められる。ここでは，$m \omega^2=k$を用いた。したがって複素応答は

$$u(t)=f_0 {\cal H}(\mbox{i}p) \exp\left(\mbox{i}pt\right) = {\cal H}(\mbox{i}p) \times f(t)$$ (10.56)

と表すことができる。したがって ${\cal H}(\mbox{i}p)$は，単位の周期外力に対する 応答の振幅であることから，周波数（振動数）応答関数 と呼ばれる。 この関数の虚数部が，節-10.1.3 (1)で 求めた正弦波応答になっている。 ${\cal H}(\mbox{i}p)$の負の偏角を$\alpha$とすると

$${\cal H}(\mbox{i}p)=\left\vert{\cal H}(\mbox{i}p)\right\vert... ...), \quad \alpha\equiv - \arg\left( {\cal H}(\mbox{i}p)\right)$$ (10.57)

と書くことができ

$$\left\vert{\cal H}(\mbox{i}p)\right\vert=\dfrac{1}{k} \dfra... ...,\beta \dfrac{p}{\omega}}{1-\left(\dfrac{p}{\omega}\right)^2}$$

となるから，式(10.35b)と比較すれば明らかなように， 動的増幅率$M\subsc{d}$との関係は

$$\left\vert{\cal H}(\mbox{i}p)\right\vert=\dfrac{M\subsc{d}}{... ...\subsc{d} \exp\left\{\mbox{i}\left(pt-\alpha\right)\right\}$$ (10.58)

であり，$\alpha$は式(10.35c)で定義された位相の遅れである。

節-10.1.3 (2)で示した支点を 強制的に変位させた場合も同様に

$$w(t)=w_0 \exp\left(\mbox{i}pt\right)$$

とすれば

$$\phi=w_0 \dfrac{p^2}{\left(\omega^2-p^2\right)+2\mbox{i} \beta \omega p}$$

と求められるので，支点との相対変位が

$$v(t)=w_0 {\cal H}_w(\mbox{i}p) \exp\left(\mbox{i}pt\right)... ...omega}\right)^2\right\} +2\mbox{i} \beta \dfrac{p}{\omega}}$$ (10.59)

となり， $\left\vert{\cal H}_w(\mbox{i}p)\right\vert=M_w$である。

さて，$f_0=1$の単位の周期外力に対する応答は，式(10.56)から

$$u(t)={\cal H}(\mbox{i}p) \exp\left(\mbox{i}pt\right) \eqno{(*)}$$

となる。同じ解は，Duhamel積分を用いても求めることができるはずだ。 ここでは，$t=-\infty$に何らかの初期条件が与えられて， それ以降継続して振動し続けているものと考えよう。 つまり $f(\tau)=\exp\left(\mbox{i}p\tau\right)$を代入したDuhamel積分は

$$u(t)=\int_{-\infty}^\infty u\subsc{i}(t-\tau) \exp\left(\mbox{i}p\tau\right)\dint \tau$$

と書くことができる。ここで$s=t-\tau$という変数変換を すると， $\dint s=-\dint \tau$, $s: \infty\to -\infty$なので上式は

$$=-\int_\infty^{-\infty} u\subsc{i}(s) \exp\left\{\mbox{i}... ...y}^\infty u\subsc{i}(s) \exp\left(-\mbox{i}ps\right)\dint s$$

と書き直すことができる。この最後の表現と上の式($*$)を比較すれば

$${\cal H}(\mbox{i}p)=\int_{-\infty}^\infty u\subsc{i}(t) \... ...mbox{i}p)=\mbox{\fontMathScript F}\left\{u\subsc{i}(t)\right\}$$ (10.60)

という関係にあることがわかる。すなわち，周波数応答関数 ${\cal H}(\mbox{i}p)$は 単位衝撃応答$u\subsc{i}$のFourier変換( F)なのだ。Fourier変換 については節-10.1.4 (1)で 説明する。したがって，特性が全く異なると考えられそうな二つの応答， つまり周期外力に対する応答と単位衝撃応答との間には 式(10.60)の関係，あるいは

$$M\subsc{d}=k \left\vert\mbox{\fontMathScript F}(u\subsc{i})\right\vert$$

という関係が存在することが明らかになった。 これはちょっと面白くありませんか。 というのも，衝撃応答には周期運動の要素は含まれていないように 感じるからだ。 しかし逆に，上の関係式およびその誘導からも明らかなように， 衝撃応答はすべての周波数の応答を均等に含んでいるのである。

10.1.4 不規則応答

10.1.4.1 準備その1 -- Fourier変換

前の節までは，種々の図の横軸には時刻$t$をとっていた。 しかし，例えば異なる日時に同じ活断層で発生した二つの地震波が， 同じ地盤を経て到達した場合，その時間軸上の加速度波形等の「時刻歴」が 同じにはならないことは容易に予想される。 それは，断層上の破壊の広さや向きが少々異なるからかもしれない。 しかし，同じ活断層から同じ地盤を経た二つの波の持つ特性の一つ， 「周波数成分の分布」には，あまり違いは無いのではないだろうかと， 期待している。 そのため，横軸に周波数をとった図で設計をすることを考えよう。 そのような，横軸に周波数をとったもので読者が既に知っていると 思われるものの代表はFourier係数だろう。 ただそのFourier級数は，最大周期$T$というものがある信号の処理を 前提とした離散的な周波数応答だった。 これに対し，例えば地震波の記録を見れば明らかなように， このような不規則信号には最大周期$T$というものが無い。 あるいは，それは$T=\infty$と考えなければならないということになる。 そこでFourier級数から議論を始め，その最大周期を無限大にとったときの 扱いを紹介し，それを用いた設計法について簡単に述べておこう。

さて，外力$f(t)$という関数のほとんどを近似的に， 最大周期$T$のFourier級数で

$$f(t)=\dfrac{a_0}{2}+\sum_{n=1}^\infty a_n \cos npt +\sum_{... ...\dfrac{2}{T}\int_{-T/2}^{T/2} f(t) \sin npt\dint t \nonumber$$

と表してもいいことは，多分知っていると思う。 ここに$p$は最小振動数 $p=\dfrac{2\pi}{T}$である。 これを式(10.56) (10.57)を用いた応答に 代入すると，$a_0$を省略して大雑把には

$$u(t)\sim \sum_{n=1}^\infty a_n \left\vert{\cal H}(\mbox{i... ...\cal H}(\mbox{i}np)\right\vert \sin \left(npt-\alpha\right)$$

となることから

$$\left\vert\mbox{応答のFourier係数}(np)\right\vert =\left\ve... ...t \times \left\vert\mbox{外力のFourier係数}(np)\right\vert$$

という関係が成り立つことになる。

図 10.22: 不規則応答の考え方

例えば外力のFourier係数の絶対値は， 外力の中の$np$ ( $n=1, 2,\cdots \infty$)の周波数成分の「強さ」を表している。 社会基盤構造を設計して${\cal H}$を設定し， 外力のFourier係数をそれに乗じて応答のFourier係数の絶対値が求められ， 応答の中の$np$の周波数成分の「強さ」が求められる。 もし，この応答の「強さ」が，材料や構造が抵抗できる「強さ」よりも 大きくなってしまうと，その構造は破壊してしまうことになる。 つまり，例えば地震の外力を念頭に置いたとき，図-10.22に 示したように，構造や地盤は外力信号を処理するフィルターの役目を していて，それを動的に設計する際には，この応答のFourier係数（一番 右の図）が所定の大きさを超えないように，材料や構造， つまり$\omega$と$\beta$あるいは${\cal H}$を適切に設計すればいいことになる。 以下，もう少し丁寧に，しかし天下り的に説明しておこう。

さて，式(10.61)にEulerの公式

$$\cos\theta=\dfrac{e^{\mbox{\scriptsize i}\theta} +e^{-\mbox... ...ptsize i}\theta} -e^{-\mbox{\scriptsize i}\theta}}{2\mbox{i}}$$ (10.61)

を用いて^10.8整理すると

$$f(t)=\dfrac{a_0}{2} +\sum_{n=1}^\infty \dfrac{a_n}{2} \left... ...}-\dfrac{b_n}{2\mbox{i}}\right) e^{-\mbox{\scriptsize i}np t}$$

と表すことができる。そこで式(10.61)を用いると，$n=0$も含めて

$$\begin{eqnarray*} \dfrac{a_n}{2}\pm\dfrac{b_n}{2\mbox{i}} &=&\dfrac{1}{T}\int_{-... ...{-T/2}^{T/2}f(\tau) e^{\mp\mbox{\scriptsize i}np \tau}\dint\tau \end{eqnarray*}$$

となるので，上式に代入し直すと

$$\begin{eqnarray*} f(t) &=&\sum_{n=0}^\infty e^{\mbox{\scriptsize i}np t} \dfrac... ...t_{-T/2}^{T/2}f(\tau) e^{-\mbox{\scriptsize i}np \tau}\dint\tau \end{eqnarray*}$$

という表現になる。これは複素Fourier級数 と呼ばれ，結局

(10.62)

と表すことができる。$c_n$は複素Fourier係数と呼ばれる。

続いて少し乱暴な演算をするが， ここで $\Delta p=\dfrac{2\pi}{T}$, $p_n=np$と置くと， 上式(10.63b)の$c_n$を式(10.63a)の 級数表現に代入したものは

$$f(t)=\dfrac{1}{2\pi} \sum_{n=-\infty}^{\infty} \exp\left(... .../2}^{T/2} f(\tau) \exp\left(-\mbox{i}p_n\tau\right)\dint\tau$$

と書くことができる。 ここで$T\to\infty$とすることによって，$\Delta p\to 0$, $$\sum_{n=-\infty}^{\infty} \Delta p \to \int_{-\infty}^\infty \dint p$$と 置き換えられるものとした上で，$p_n$を離散量から連続量に変更して$p_n\to p$と 置き直せば，上式は

$$f(t)=\dfrac{1}{2\pi} \int_{-\infty}^\infty \exp\left(\mbox... ...fty}^\infty f(\tau) \exp\left(-\mbox{i}p\tau\right)\dint\tau$$

となるので，結局

(10.63)

という関係を得る。式(10.64b)が 関数$f(t)$のFourier変換 であり，式(10.64a)がFourier逆変換である。 なお，指数関数の符号が逆に定義されている場合や， 逆変換の積分の前だけに $\slfrac{1}{2\pi}$とする代わりに， 両方の積分の前に $\slfrac{1}{\sqrt{2\pi}}$として定義する場合があるので 注意して欲しい。 いずれにしても物理的な意味は変わらない。

図 10.23: Fourier変換はちょうど光学プリズムのようなもの

複素数が出てきたりしてFourier変換と聞いただけで訳がわからなくなるかも しれないが，ちょうど図-10.22の左側の$f$のFourier係数の 分布のように，このFourier変換の「値」は， ある信号$f(t)$に含まれる成分の大きさを 連続的な周波数毎に表していると考えればいい。 そしてFourier変換という「操作」は，ちょうど光の分光のプリズムの役割と 同じだと考えればいい。 プリズムに白色光を当てると，波長が700nmくらいの赤（周波数の逆数の 周期に換算すると $2.5\times 10^{-15}$sくらい）から， 波長が300nmくらいの紫（周期が $1.3\times 10^{-15}$sくらい）の 可視光に分解して目で見ることができる。 地震加速度や構造の不規則応答も，ちょうど図-10.23の 右側に示したように，Fourier変換というプリズムを通して， 例えば社会基盤構造にとって重要な1sくらいまでの周期成分毎に分解して 取り扱おうとしていると考えれば，少しは理解が進むかもしれない。

10.1.4.2 準備その2 --波のスペクトル解析

不規則応答の代表例は地震応答である。地震波そのものも不規則信号であろう。 いずれも統計量としての平均値は零なので，あまり関心は無い。 興味があるのは力や応答の「強さ」であるから，振幅の2乗，あるいは 分散が重要な統計量であろう。そこで各種文献を見ると，例えば 外力$f(t)$の分散は

$$\sigma_f^2\equiv \left\{\mbox{$f(t)$の2乗平均}\right\} = \d... ...\dfrac{\left\vert F(\mbox{i}p)\right\vert^2}{T}\right\}\dint p$$ (10.64)

と定義してもいいらしい。 ここに$F(\mbox{i}p)$は$f(t)$のFourier変換である。 この被積分関数を

$$S_f(p)\equiv\lim_{T\to\infty}\dfrac{\left\vert F(\mbox{i}p)\right\vert^2}{T}$$ (10.65)

と記し，$f(t)$のパワースペクトル密度関数 と呼ぶ。例えば

1. もし，ある信号$f(t)$に，すべての周波数成分が均等に含まれている 場合には，$S_f(p)$は一定になる。このような波はwhite noise と呼ばれる。
2. また$f(t)=\sin qt$のように，1種類の周波数$q$しか 含まれていない場合には， $S_f(p)=\delta(p-q)$になる。$\delta$は 式(10.43)のDiracのデルタ関数である。

が一番わかり易い例である。 ちなみに $\left\vert F(\mbox{i}p)\right\vert$をFourierスペクトル と呼んでいる。

これに対し，ある外力$f(t)$の，時間$\tau$だけずれた同じ信号との 積の期待値

$$R_f(\tau)\equiv\lim_{T\to\infty}\dfrac{1}{T} \int_{-T/2}^{T/2} f(t) f(t+\tau)\dint t$$ (10.66)

を自己相関関数 と呼ぶ。$f(t)$が周期関数なら，$\tau$が周期の倍数のときに 最も相関が高くなることから，この自己相関関数も何らかの周期特性を 表現していることが予想される。ちなみに

$$R_f(0)=\sigma_f^2$$

である。この自己相関関数をFourier変換すると

$$\begin{eqnarray*} && \int_{-\infty}^\infty R_f(\tau) \exp\left(-\mbox{i}p\tau... ...o\infty} \dfrac{\left\vert F(\mbox{i}p)\right\vert^2}{T} =S_f(p) \end{eqnarray*}$$

のように，パワースペクトル密度関数になる。 つまり

$$S_f(p)=\mbox{\fontMathScript F}\left\{R_f(\tau)\right\}, \quad R_f(\tau)=\mbox{\fontMathScript F}^{-1}\left\{S_f(p)\right\}$$ (10.67)

という関係が成立し，これをWiener-Khintchineの関係^10.9 と呼んでいる。

10.1.4.3 不規則外力に対する変位応答

では，不規則な外力に対する変位応答を求めよう。 時間軸における時刻歴応答を求めたいならDuhamel積分や 数値手法を用いればいい。しかし前述のように，不規則外力に 対する応答も不規則であり，やはり周波数の関数として取り扱った方が 物理的には意味がありそうだ。 ここでは，その把握の方法を示しておく。1自由度系の 応答$u(t)$の自己相関関数を，外力の場合と同様

$$R_u(\tau)\equiv\lim_{T\to\infty}\dfrac{1}{T} \int_{-T/2}^{T/2} u(t) u(t+\tau)\dint t$$ (10.68)

と定義する。これにDuhamel積分の $$u(t)=\int_{-\infty}^\infty u\subsc{i}(t-\xi) f(\xi)\dint\xi$$を 代入すると

$$R_u(\tau) =\lim_{T\to\infty}\dfrac{1}{T} \int_{-T/2}^{T/2} ... ...fty u\subsc{i}(t+\tau-\eta) f(\eta)\dint\eta\right\} \dint t$$

と書くことができる。 ここで$s=t-\xi$, $\zeta=t+\tau-\eta$という変数変換をすると

$$\begin{eqnarray*} &=& \lim_{T\to\infty}\dfrac{1}{T} \int_{-T/2}^{T/2} \left\{... ...}{T} \int_{-T/2}^{T/2} f(t-s) f(t+\tau-\zeta)\dint t \right\} \end{eqnarray*}$$

さらに$\xi=t-s$の変数変換をすると，$f$の自己相関関数が現れ

$$= \int_{-\infty}^\infty u\subsc{i}(s)\dint s \int_{-\infty... ...{i}(s) u\subsc{i}(\zeta) R_f(\tau+s-\zeta)\dint s\dint\zeta$$

となる。

最終的に応答のパワースペクトル密度関数は， 上式の自己相関関数をFourier変換して

$$S_u(p)=\int_{-\infty}^\infty R_u(\tau) \exp\left(-\mbox{i... ...) \exp\left(-\mbox{i}p\tau\right)\dint s\dint\zeta\dint\tau$$

という表現になる。ここで変数変換を $\eta=\tau+s-\zeta$とすると

$$\begin{eqnarray*} &=&\int_{-\infty}^\infty\int_{-\infty}^\infty\int_{-\infty}^\i... ...nt\eta \\ &=&{\cal H}(\mbox{i}p) {\cal H}(-\mbox{i}p) S_f(p) \end{eqnarray*}$$

となるので，結局

$$S_u(p) = \left\vert{\cal H}(\mbox{i}p)\right\vert^2 S_f(p)$$ (10.69)

という関係が求められる。あるいは式(10.58)を用いると

$$S_u(p)=\left(\dfrac{M\subsc{d}}{k}\right)^2 S_f(p)$$

という関係にあることになる。したがって，$S_f(p)$が一定のwhite noiseを 構造に作用させて得ることができる応答のパワースペクトル密度関数$S_u(p)$は， そのまま共振曲線（の2乗）に相当することになる。 例えば社会基盤構造で観測される真夜中の微弱な振動（常時微動と呼ぶ）を 測定すれば，それはほとんどwhite noiseのような幅の広い周波数帯で 同じパワーを持った外力（周辺環境から伝わってくる種々雑多な動的作用） に対する応答だと近似していいので， そのデータが共振曲線つまり周波数応答関数の絶対値に相当する。

図 10.24: 連続スペクトルと設計

すなわち，外力のパワースペクトル密度関数が与えられたとき， この式(10.70)で求められる応答のパワースペクトル密度関数が 所定の条件を満足するように（構造が破壊しないように）， 構造や材料を選択して ${\cal H}(\mbox{i}p)$を 設計すればいいことになる。図-10.22ではFourier級数を 用いていたので，スペクトルが$p$刻みに表現された離散的なものであったのに対し， 外力そのものに最大周期というものが定義できない地震外力のような 不規則外力とその応答を検討する場合には，図-10.24の ように連続的な周波数で関数表示された$S_f(p)$, $S_u(p)$といった 連続スペクトルを用いるのである。 この図のような場合には，設計Bの方が望ましいと考えていいだろう。

10.1.4.4 加速度応答と入力地震

耐震設計のことを最後に考えておく。式(10.38)にあるように， 外力は地震加速度に質量を乗じたものである。 つまり，前節の$f(t)$には，標準的な地震加速度$\ddot{w}$を$-m \ddot{w}$と して入力すればいいように 思われるが，実際にはそうなっていない。文献[129]を 斜めに読むと，ある固有振動数（周期）を持つ1質点系に標準的な 地震加速度を入力して得ることができる「絶対加速度応答$\ddot{u}$」の 最大値を，基部に与える入力地震加速度にするとしている。 なんとなく逐次代入法の（否，屋上屋を架す， あるいは無闇に安全を確保しようとする）ようで第1著者には理解できないが， 構造の過渡応答等も含めて適切な入力加速度を求めるためとされている。

まずDuhamel積分表現に式(10.47)の単位衝撃応答を代入すると

$$u(t)=\int_0^\infty \dfrac{f(\tau)}{m \omega_d} \exp\left... ...au) \sin\left\{\omega_d\left(t-\tau\right)\right\}\dint\tau$$

と変位が表される。これを$t$で微分すると

$$\begin{eqnarray*} \dot{u}&=&\int_0^\infty \dfrac{f(\tau)}{m \omega_d}e^{-\beta... ..._d(t-\tau)\} +\omega_d \cos\{\omega_d(t-\tau)\} \right]\dint\tau \end{eqnarray*}$$

と表される。被積分関数の第2項は，デルタ関数の定義から零になる。 これを運動方程式に代入すれば加速度が求められるが， ここでは再度$t$で微分すると

$$\begin{eqnarray*} \ddot{u}&=&-\beta\omega\left(\int_0^\infty \dfrac{f(\tau)}{m\... ...} H(t-\tau) \cos\{\omega_d(t-\tau)\} \dint\tau +\dfrac{f(t)}{m} \end{eqnarray*}$$

となる。ここに式(10.37)から$u \to v$に変更し， 式(10.38)から $f(t)=-m\ddot{w}$を代入すると

$$\ddot{v}=\dfrac{\left(\omega_d^2-\beta^2\omega^2\right)}{\om... ...tau)} H(t-\tau) \cos\{\omega_d(t-\tau)\} \dint\tau -\ddot{w}$$

となる。再度式(10.37)を用いれば，絶対変位の加速度を

$$\ddot{u}(t)=\ddot{v}+\ddot{w}=\omega_d \int_0^\infty \widetilde{u}\subsc{i}\left(t-\tau\right) \ddot{w}(\tau) \dint\tau$$ (10.70)

と表すことができる。ここに

$$\widetilde{u}\subsc{i}(t-\tau)\equiv \exp\left\{-\beta \om... ...c{2\beta}{\sqrt{1-\beta^2}} \cos\{\omega_d(t-\tau)\} \right\}$$ (10.71)

と定義した。したがって，加速度応答のスペクトル密度関数は

$$S\subsc{a}(p) = \omega_d^2 \left\vert\widetilde{{\cal H}}... ...box{\fontMathScript F}\left\{\widetilde{u}\subsc{i}(t)\right\}$$ (10.72)

という関係になる。ここで$S\subsc{a}$は1質点系の加速度応答$\ddot{u}$の スペクトル密度関数で，$S\subsc{e}$はその1質点系に 実際に与えた入力地震加速度$\ddot{w}$のスペクトル密度関数である。

図 10.25: 加速度応答と入力地震

そこで現耐震設計法では，ある標準的な設計用の入力地震加速度$\ddot{w}$を， 種々の固有周期（振動数）$\omega_i$, $\omega_j$ etc.を 持つ1自由度系に入力して，式(10.73)で 加速度応答を求める。その最大値 $\left(S\subsc{a}\right)\sub{max}$を， その対応した固有周期（振動数）を 持つ構造に対する入力地震として，その構造の基部に与えることになっている。 考え方を簡単に図-10.25に示した。 文献[129]には，加速度の単位でこの 最大値 $\left(S\subsc{a}\right)\sub{max}$が与えられている。 例えば1自由度系を設計する場合には， 設計して得た$\omega$で固有周期$T$を求め， この右の図からその$T$の値に対応した 最大加速度応答 $\left(S\subsc{a}\right)\sub{max}$を求める。 そして， 設計して得た質量$m$をそれに乗じた $m \left(S\subsc{a}\right)\sub{max}$を その質点に外力として作用させたときのバネの抵抗力を算定し， それが破壊しないレベルにあるかどうかを確かめるのである。

最新版を正確に読むためには pdf ファイル をどうぞ。これは web 検索のための簡易旧版です。