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

• 10.4.1 梁の運動方程式
• 10.4.2 非減衰自由振動--正攻法
• 10.4.3 非減衰自由振動
  • 10.4.3.1 調和振動解
  • 10.4.3.2 両端単純支持梁
  • 10.4.3.3 片持ち梁
  
  
• 10.4.4 振動モードの直交性と自由振動解--モード解析法
• 10.4.5 粘性減衰自由振動
• 10.4.6 強制振動
  • 10.4.6.1 モード解析法
  • 10.4.6.2 単位衝撃応答
  • 10.4.6.3 減衰を含んだDuhamel積分
  • 10.4.6.4 例--走行荷重
  
  
• 10.4.7 有限要素と振動解析
  • 10.4.7.1 仮想仕事式と質量行列および減衰行列
  • 10.4.7.2 自由振動解析
  • 10.4.7.3 時刻歴応答解析
  • 10.4.7.4 平面骨組の振動解析
    • 10.4.7.4.1 軸方向変形についての運動:
    • 10.4.7.4.2 任意の平面骨組の運動:
    • 10.4.7.4.3 自由振動の数値解析例:
  
  
• 10.4.8 回転慣性とせん断変形
  • 10.4.8.1 回転慣性の影響
  • 10.4.8.2 せん断変形の影響

10.4 連続体の振動--梁の振動

10.4.1 梁の運動方程式

図 10.41: 梁の局所的な微分要素の運動

図-10.41には， 分布外力$p(x,t)$が作用し，合応力としての曲げモーメントと せん断力が発生して運動している状態の梁の微分要素$\dint x$を示した。 この微分要素の$z$方向の運動方程式は，Newtonの法則から

$$\D{V(x,t)}{x}+p(x,t)=m \D[2]{w(x,t)}{t}$$ (10.167)

となる。ここに$w(x,t)$は$z$方向のたわみであり，$m$は 梁の単位長さ当たりの質量で，密度を$\rho$として 断面積を$A$とすると，$m=\rho A$である。 次に点Aの反時計回りのモーメントのつり合いから

$$\D{M(x,t)}{x}-V(x,t)=0$$ (10.168)

が^10.11求められる。 式(10.169)を式(10.168)に代入することに よって，合応力で表した梁の運動方程式は

$$\D[2]{M(x,t)}{x}+p(x,t)=m \D[2]{w(x,t)}{t}$$ (10.169)

でいいことがわかる。さらに，$x$軸が断面の図心を通るように定義すれば， 式(4.13b)から

$$M(x,t)=-EI \D[2]{w(x,t)}{x}$$ (10.170)

なので，これを式(10.170)に代入すれば， 最終的な梁の運動方程式 が

$$-\D[2]{}{x}\left(EI \D[2]{w(x,t)}{x}\right)+p(x)=m \D[2]{w(x,t)}{t}$$

あるいは，以下，一様材料等断面の梁だけを対象とすれば

$$m \D[2]{w(x,t)}{t}+EI \D[4]{w(x,t)}{x}=p(x,t)$$ (10.171)

となる。場所$x$についての微係数が4階になっているが， そのことと符号に目をつぶれば， 形式的には式(10.130)の波動方程式によく似ている。

解が唯一であるためには，境界条件が式(4.28)のように， 両端$x_i, (i=1,2)$, $x_1=0$, $x_2=\ell$において適切に

(10.172)

と与えられている必要がある。 ここに$\theta$は反時計回りを正として定義した たわみ角であり，$n_i$は式(4.24)で 定義した断面の向きを表す符号で

$$% latex2html id marker 13120n_i\equiv \left\{ \begin{arra... ...array}\right. \eqno{(\ref{eq:boundarydirection}) \mbox{再掲}}$$

という値を持つ。 また図-10.42に示したように，$S_i(t)$と$C_i(t)$は 下向き正のせん断外力と反時計回りが正の外力モーメントである。 また，初期条件は，梁の初期位置と初速で与えればいいから

$$w(x,0)=w_0(x), \quad \D{w}{t}(x,0)=v_0(x)$$ (10.173)

とすればいい。

10.4.2 非減衰自由振動--正攻法

図 10.42: 両端の境界条件

最初は最も基本的な自由振動，つまり $p(x,t)\equiv 0$で端外力も 零の場合を解いておこう。ただしこの節では， 固有関数を用いたやや数学的なアプローチをしているので， 嫌いな人は読み飛ばしてもいい。 さて自由振動問題は， 振動解析で最も重要な固有振動数と固有振動モードを求める段階に当たる。 前の節-10.3.2を勉強した人は，この梁の 場合には，まず固有値問題

$$\D*[4]{\phi(x)}{x}=\Lambda \phi(x)$$

を解けばいいことがわかるかもしれないが，ここではゆっくり， まず変数分離で解いてみよう。 ただし，簡単のために，最初は長さ$\ell$の両端単純支持の問題を 対象とする。

まず，たわみを変数分離して

$$w(x,t)=\phi(x) q(t)$$

と置いて，運動方程式(10.172)に代入して，変数を見ながら整理すると

$$-\dfrac{\ddot{q}(t)}{q(t)}=\dfrac{EI}{m} \dfrac{\phi''''(x)}{\phi(x)}= \Lambda (\mbox{const.}) \eqno{[\star]}$$

のように，左辺右辺が定数でなければならないことがわかる。 ここにプライムは$x$による微分である。 境界条件は，両端でたわみ$w$とたわみの2階の微係数$\D[2]{w}{x}$が 零であればいいので，これも変数分離をすると，結局

$$\phi(0)=0, \quad \phi''(0)=0, \qquad \phi(\ell)=0, \quad \phi''(\ell)=0$$

となる。初期条件は分離できないので，まず解けるのは$\phi(x)$の方で， その微分方程式は上式から

$$\phi''''(x)=\Lambda \dfrac{m}{EI} \phi(x)$$

である。実は，上の変数分離した微分方程式のうち，$q(t)$の方を 見ると，$\Lambda$が正でない限り周期運動をしないことは明らかであるが， ここでは，固有関数を求める一般的なアプローチを続けておく。

まず $\Lambda \dfrac{m}{EI}=-4\mu^4<0$の場合は，一般解は

$$\phi(x)=\exp\left(\mu x\right) \left(a \sin\mu x+b \cos... ...p\left(-\mu x\right) \left(c \sin\mu x+d \cos\mu x\right)$$

となる。これを境界条件に代入すると，最終的に

$$\left(\begin{array}{cc} \sin\mu\ell\cosh\mu\ell & \cos\mu\e... ...y}\right\}=\left\{\begin{array}{c} 0 0 \end{array}\right\}$$

となるが，この係数行列は正則なので，$a$, $c$は零解しか存在しない。 もちろん$\Lambda=0$の場合も零解しか存在しないことはすぐにわかる。

したがって $\Lambda \dfrac{m}{EI}=\mu^4>0$とおくと，一般解は

$$\phi(x)=a \sin\mu x+b \cos\mu x+c \sinh\mu x+d \cosh\mu x$$ (10.174)

となる。これを境界条件に代入すると

$$\left(\begin{array}{cc} \sin\mu\ell & \sinh\mu\ell \\ -\s... ...y}\right\}=\left\{\begin{array}{c} 0 0 \end{array}\right\}$$

になる。この場合には，$a$または$c$が零にならない条件から， 係数行列は特異にならなければならないので， その行列式の零解を求めようとすると

$$\sin\mu\ell \sinh\mu\ell=0$$

から，任意の整数$n=1$, 2, $\cdots$に対して

$$\sin\mu\ell=0 \quad \to \quad \mu\ell=n\pi \quad \to \quad \mu_n=\dfrac{n\pi}{\ell}$$

で，$a$が不定になるので

$$\phi(x)= \phi_n(x)=\sin \mu_n x = \sin\left(\dfrac{n\pi x}{\ell}\right)$$

と求められる。この$\mu_n$が固有値で，$\phi_n(x)$が固有関数， つまり固有振動モードになる。

したがって

$$\Lambda_n=\dfrac{EI}{m} \mu_n^4 =\dfrac{EI}{m} \left(\dfrac{n\pi}{\ell}\right)^4$$

であるから，各$n$に対する$q(t)$についての微分方程式は

$$\ddot{q}_n(t)+\dfrac{EI}{m} \left(\dfrac{n\pi}{\ell}\right)^4 q_n(t)=0$$

となるので，ここで

$$\omega_n^2\equiv \dfrac{EI}{m} \left(\dfrac{n\pi}{\ell}\rig... ...t{\dfrac{EI}{m}} \left(\dfrac{n\pi}{\ell}\right)^2 \eqno{(*)}$$

と置けば，結局

$$q_n(t)=A_n \sin\omega_n t+B_n \cos\omega_n t$$

と求めることができ，自由振動解は

$$w(x,t)=\sum_{n=1}^\infty \phi_n(x) \left( A_n \sin\omega_n t+B_n \cos\omega_n t \right)$$ (10.175)

と表すことができる。 つまり，梁の振動は，式(10.175)のように三角関数と 双曲線関数の組み合わせで表される固有振動モードを持ち， 式($*$)のような固有振動数で単純に振動することになる。 固有振動数は，1自由度系のそれ $\sqrt{\slfrac{k}{m}}$と同様， 曲げ剛性を質量で除したものの平方根 $\sqrt{\slfrac{EI}{m}}$の 形を持っている。$A_n$, $B_n$を初期条件から 決める方法については，節-10.4.4で述べる。

10.4.3 非減衰自由振動

10.4.3.1 調和振動解

ここからは もう少し天下り的な工学的な方法を述べよう。非減衰の 振動問題なので，時間についてはある振動数$\omega$で振動する 解になるのは明らかかもしれない。したがって，解を調和振動解で

$$w(x,t)=\phi(x) \exp\left(\mbox{i}\omega t\right)$$ (10.176)

と仮定しよう。p.で変数分離した式[$\star$]から 振動解になる$q(t)$を先に求めたと考えればいい。$\phi(x)$が固有振動モード になる。この式(10.177)を， $p(x,t)\equiv 0$とした 運動方程式(10.172)に代入すると

$$\left\{EI \phi''''(x)-m \omega^2 \phi(x)\right\} \exp\left(\mbox{i}\omega t\right)=0$$

となる。これが任意の時刻$t$で成立しなければならないので，$\phi(x)$は

$$\phi''''(x)-\dfrac{m \omega^2}{EI} \phi(x)=0$$ (10.177)

を満足しなければならない。そこで

$$\mu^4\equiv \dfrac{m \omega^2}{EI}$$ (10.178)

と定義すると，上式(10.178)の一般解は

$$\phi(x)=c_1 \sin\mu x+c_2 \cos\mu x +c_3 \sinh\mu x+c_4 \cosh\mu x$$ (10.179)

となる。これは当然，前節で求めた式(10.175)と同じものである。

10.4.3.2 両端単純支持梁

まず両端単純支持梁を対象としよう。境界条件も式(10.177)を 代入すれば，任意の時刻に成立することから，$\phi(x)$のみで

$$\phi(0)=0, \quad \phi''(0)=0, \qquad \phi(\ell)=0, \quad \phi''(\ell)=0$$

と表される。この条件に式(10.180)の一般解を代入すると

$$\begin{eqnarray*} &&c_2+c_4=0, \quad \mu^2 \left(-c_2+c_4\right)=0, \quad c_1\... ...mu\ell-c_2\cos\mu\ell +c_3\sinh\mu\ell+c_4\cosh\mu\ell\right)=0 \end{eqnarray*}$$

となるので，行列表示すると

$$\left(\begin{array}{cccc} 1 & 0 & 1 & 0 \\ 1 & 0 & -1 & 0... ...gin{array}{c} 0 0 0 0 \end{array}\right\} \eqno{(a)}$$

となる。積分定数が有意な解を得るためには，この係数行列が特異であればいいので

$$\det\left\vert\mbox{係数行列}\right\vert=1\times \det\left\v... ...mu\ell \end{array}\right\vert =-4 \sin\mu\ell \sinh\mu\ell=0$$

が成立すればいいことになる。したがって

$$\sin\mu\ell=0 \quad \to \quad \mu\ell=n\pi \quad \to \quad \mu_n=\dfrac{n\pi}{\ell} \eqno{(b)}$$

となる。この結果を上式($a$)に代入すると，$c_1$のみが不定になるので

$$\phi_n(x)=\sin \mu_n x = \sin\left(\dfrac{n\pi x}{\ell}\right)$$ (10.180)

である。この$\mu_n$が固有値で，$\phi_n(x)$が 固有振動モードになる。式($b$)の 固有値を式(10.179)に代入すれば，固有振動数は

$$\mu^4=\dfrac{m\omega^2}{EI}=\left(\dfrac{n\pi}{\ell}\right)^... ...mega_n=\sqrt{\dfrac{EI}{m}} \left(\dfrac{n\pi}{\ell}\right)^2$$ (10.181)

になる。最終的に解は，すべてのモードの重ね合わせで

$$w(x,t)=\sum_{n=1}^\infty \sin\left(\dfrac{n\pi x}{\ell}\right) \left(A_n \sin\omega_n t+B_n \cos\omega_n t\right)$$ (10.182)

と表される。この未定係数$A_n$, $B_n$は初期条件から 求めなければならないが， それについては節-10.4.4で 説明する。図-10.43に第1〜3次モードと振動数を列挙した。

図 10.43: 単純支持梁の振動数と振動モード
図 10.44: 片持ち梁の振動モード

10.4.3.3 片持ち梁

左端が固定された片持ち梁の場合には，境界条件は

$$\phi(0)=0, \quad \phi'(0)=0, \quad \phi''(\ell)=0, \quad \phi'''(\ell)=0$$

になる。この条件に式(10.180)の一般解を代入すると

$$\begin{eqnarray*} &&c_2+c_4=0, \quad \mu \left(c_1+c_3\right)=0, \quad \mu^2 ... ...mu\ell+c_2\sin\mu\ell +c_3\cosh\mu\ell+c_4\sinh\mu\ell\right)=0 \end{eqnarray*}$$

となるので，行列表示すると

$$\left(\begin{array}{cccc} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 ... ...gin{array}{c} 0 0 0 0 \end{array}\right\} \eqno{(c)}$$

となる。積分定数が有意な解になるためには， この係数行列が特異であればいいので

$$\det\left\vert\mbox{係数行列}\right\vert=1\times \det\left\v... ...u\ell \end{array}\right\vert =2 \cos\mu\ell \cosh\mu\ell+2=0$$

が成立すればいいことになる。したがって

$$\cos\mu_n\ell \cosh\mu_n\ell=-1$$ (10.183)

を満足する$\mu_n$ ($n=1$, 2, $\cdots$)が固有値になる。 この結果を上式($c$)に代入すると

$$c_2=\Xi_n c_1, \quad c_3=-c_1, \quad c_4=-\Xi_n c_1, \quad... ...\cos\mu_n\ell+\cosh\mu_n\ell}% {\sin\mu_n\ell-\sinh\mu_n\ell}$$

と求められる。結局，振動モードは

$$\phi_n(x)=\left(\sin\mu_n x-\sinh\mu_n x\right) +\Xi_n \left(\cos\mu_n x-\cosh\mu_n x\right)$$ (10.184)

になる。式(10.184)は陽に解くことはできないので， 例えば2分法 等で数値的に解けばいい。 最初の三つを列挙しておくと

$$\mu_1\ell=1.875, \quad \omega_1=\dfrac{3.516}{\ell^2}\sqrt{\... ....855, \quad \omega_3=\dfrac{61.70}{\ell^2}\sqrt{\dfrac{EI}{m}}$$

くらいになる。図-10.44に対応する振動モードを描いた。

10.4.4 振動モードの直交性と自由振動解--モード解析法

前節までの範囲では，固有振動モードと固有振動数を求めることはできたが， 一般解を初期条件に代入して積分定数を求めたり， 強制振動解を得ることはまだできていなかった。 ここでは，多自由度系でも用いたモード解析法を梁に対して誘導しよう。 まず最も重要なのは振動モードの直交性である。 式(10.178)から第$n$次モードは

$$EI \phi_n''''(x)-m \omega_n^2 \phi_n(x)=0$$ (10.185)

を満足している。そこで，第$i$次モードを仮想変位とするような 仮想仕事を全スパンで計算すると

$$\int_0^\ell \phi_i(x) \left\{EI \phi_n''''(x)-m \omega_n^2 \phi_n(x)\right\}\dint x=0$$

が成立する。第1項を取り出し，それを部分積分して境界条件を考慮すると

$$\int_0^\ell EI \phi_i(x) \phi''''_n(x)\dint x =\left(EI\p... ..._n(x)\dint x =\int_0^\ell EI \phi''_i(x) \phi''_n(x)\dint x$$

となる。境界条件項が零になることについては，単純支持や片持ち梁では 明らかだろう。それ以外も確かめてみて欲しい。したがって，上の仮想仕事式は

$$\int_0^\ell EI \phi''_i(x) \phi''_n(x)\dint x - \omega_n^2 \int_0^\ell m \phi_i(x) \phi_n(x) \dint x=0$$

と表すことができる。これと同じプロセスを，第$i$次モードの運動方程式と 第$n$次モードとの仮想仕事の算定で行うと

$$\int_0^\ell EI \phi''_n(x) \phi''_i(x)\dint x - \omega_i^2 \int_0^\ell m \phi_n(x) \phi_i(x) \dint x=0$$

を得る。この二つの式を差し引きすると，結局，第1項は同じなので消えてしまい

$$\left(\omega_i^2-\omega_n^2\right) \int_0^\ell m \phi_n(x) \phi_i(x) \dint x=0$$

となる。多自由度系でも説明したように，安定した構造であれば， 異なる二つの固有振動数は異なる値を持っていることから，この式は

$$\left\langle \phi_n, \phi_i \right\rangle_m \equiv \int_0^\ell m \phi_n(x) \phi_i(x) \dint x=0$$ (10.186)

であることを示している。 この式は，二つの関数の（重み$m$を持つ）内積 を定義しており，その内積が零になることから，固有振動モードは お互いに直交している と呼んでいる。同じモード同士の内積は，線形代数からの類推で， 固有振動モードの「ノルム」と考えてもよく，それを

$$\left\vert\phi_n\right\vert^2_m\equiv \left\langle \phi_n,... ...\rangle_m = \int_0^\ell m \phi^2_n(x) \dint x=\overline{m}_n$$ (10.187)

と記すことにする。 $\overline{m}_n$は， 多自由度系のそれと同様一般化された質量 と呼ばれる。したがって，式(10.187)と一緒にすれば

$$\left\langle \phi_n, \phi_i \right\rangle_m = \int_0^\el... ...} 0, & i\neq n \\ \overline{m}_n, & i=n \end{array}\right.$$ (10.188)

と書くことができる。

この式(10.189)の$\phi_n$を，式(10.186)を 用いて$\phi''''_n$で置き換えると

$$\int_0^\ell m \phi_n(x) \phi_i(x) \dint x = \int_0^\ell \dfrac{EI}{\omega_n^2} \phi''''_n(x) \phi_i(x) \dint x$$

となることから，式(10.189)の直交性は別の内積の定義を用いて

$$\left\langle \phi_n, \phi_i \right\rangle_{EI} \equiv \i... ...q n \\ \overline{m}_n \omega_n^2, & i=n \end{array}\right.$$ (10.189)

と定義してもいいことになる。 あるいは，もう少しカッコ良くするために，2回部分積分をした上で境界条件を 考慮すると，結局，二つの関数の微係数の階数が同じになることから， もう一つ別の「対称な」（重み$EI$を持つ）内積が定義でき

$$\left\langle \phi_n, \phi_i \right\rangle_{EI}\super{(sym... ...q n \\ \overline{m}_n \omega_n^2, & i=n \end{array}\right.$$ (10.190)

と書いてもいい。仮想仕事的にはカッコいいが，以下では使わない。

では最終的に自由振動解を求めておこう。 前節までで，式(10.176)や式(10.183)のように， 無限個存在する固有振動モードを用いて，解は

$$w(x,t)=\sum_{n=1}^\infty \left( A_n \sin\omega_n t+B_n \cos\omega_n t \right) \phi_n(x)$$ (10.191)

と級数表現できることがわかった。 この未定係数を決定するために，初期条件式(10.174)に代入すると

$$w(x,0)=\sum_{n=1}^\infty B_n \phi_n(x)=w_0(x), \quad \D{w}{t}(x,0)=\sum_{n=1}^\infty A_n \omega_n \phi_n(x)=v_0(x)$$

となる。両式とある振動モード$\phi_i(x)$とを，式(10.189)で 定義した内積をとると

$$\sum_{n=1}^\infty B_n \left\langle \phi_n, \phi_i \righ... ...i \right\rangle_m= \left\langle v_0, \phi_i \right\rangle_m$$

となるが，モードの直交性から総和の項は第$i$項以外は零になり

$$B_i \left\langle \phi_i, \phi_i \right\rangle_m= \left\... ...i \right\rangle_m= \left\langle v_0, \phi_i \right\rangle_m$$

という関係を得る。したがって，未定係数が

$$A_i=\dfrac{1}{\overline{m}_i \omega_i} \int_0^\ell m v_0... ...c{1}{\overline{m}_i} \int_0^\ell m w_0(x) \phi_i(x) \dint x$$ (10.192)

と求められ，唯一の解を得ることができる。

1. 両端単純支持梁のモード(10.181)が 内積(10.187)で直交していることは容易にわかるだろう。 では，片持ち梁のモード(10.185)が直交することを確かめよ。
2. 両端固定梁と片端固定・片端単純支持梁の固有振動数，あるいは それを求める式を求めよ。
3. 左端固定・右端バネ支持梁の振動数方程式を求めよ。 バネは線形で抵抗係数は定数$k$である。

10.4.5 粘性減衰自由振動

もし，材料そのものに粘性抵抗がある場合には， 梁の応力ひずみ関係式(4.7)の曲げの部分だけ取り出して

$$\sigma_{xx}(x,z,t)=-E z \D[2]{w(x,t)}{x}-\eta z \D[3][2][x]{w(x,t)}{t}$$ (10.193)

といったモデルは物理的には容認できそうだ。 ここに$\eta$は材料の粘性係数である。 したがって，これを曲げモーメントの定義に代入すれば，式(10.171)の 代わりに

$$M(x,t)=-EI \D[2]{w(x,t)}{x} -\eta I \D[3][2][x]{w(x,t)}{t}$$ (10.194)

となる。したがって，運動方程式は，式(10.172)の代わりに

$$m \D[2]{w(x,t)}{t} +\eta I \D[5][4][x]{w(x,t)}{t} +EI \D[4]{w(x,t)}{x}=p(x,t)$$ (10.195)

とすればいいことになる。 また，境界条件も式(10.173)の代わりに，次のようになる。

(10.196)

では自由振動の場合( $p(x,t)\equiv 0$)に，解が

$$w(x,t)=\sum_{n=1}^\infty q_n(t) \phi_n(x)$$ (10.197)

の級数表示で表されるものとして，モード解析をしよう。 ここで$\phi_n(x)$は，非減衰の固有振動モードであることに 注意する。$q_n(t)$は梁の一般化された座標 である。この式を，運動方程式(10.196)に代入すると

$$\sum_{n=1}^\infty \left\{ m \ddot{q}_n \phi_n +\eta I ... ...rac{\eta}{E} \dot{q}_n+q_n\right) EI \phi''''_n \right\}=0$$

となる。この運動方程式と第$i$次モードとの内積をとるために

$$\int_0^\ell \phi_i \dint x \times$$

という作用を両辺にすると， 結局，式(10.189) (10.190)の直交条件を用いることによって， 総和の第$i$項のみが抽出でき，それが

$$\overline{m}_i \ddot{q}_i +\overline{m}_i \omega_i^2 \le... ...rac{\eta}{E} \omega_i^2 \dot{q}_i(t) + \omega_i^2 q_i(t)=0$$

という，1自由度系の粘性減衰自由振動の形の運動方程式になる。この 第2項の係数を形式的に，式(10.18)の1自由度系の減衰項の 係数 $2 \overline{\beta}\supersc{be}_i \omega_i$の表現に 置換したいので，Bernoulli-Euler梁の一般化された減衰定数 は

$$\overline{\beta}_i{}\supersc{be}=\dfrac{\eta \omega_i}{2E}$$ (10.198)

と定義できる。ちょうど，式(10.120)で定義したRayleigh減衰の うちの剛性比例の減衰定数 $\overline{\beta}_i{}\supersc{k}$と 同じ特性を持つことになる。 したがって，一般化された座標で表した運動方程式は

$$\ddot{q}_i(t)+2 \overline{\beta}_i{}\supersc{be} \omega_i \dot{q}_i(t) + \omega_i^2 q_i(t)=0$$ (10.199)

となる。形式的には1自由度系の減衰自由振動の運動方程式に一致する。 ただ実際には，有限要素法を用いた上で， 式(10.119)のRayleigh減衰を用いて解析するのが普通だと思われる。

10.4.6 強制振動

10.4.6.1 モード解析法

粘性減衰も考慮できることは前節で示したが，前節のような扱いは， 現場では一般的ではないかもしれないので，しばらくは無視し， 非減衰の強制振動を論じておく。 したがって，運動方程式は式(10.172)の

$$m \D[2]{w(x,t)}{t}+EI \D[4]{w(x,t)}{x}=p(x,t)$$

である。これに，自由振動モードで級数表示した式(10.198)の

$$w(x,t)=\sum_{n=1}^\infty q_n(t) \phi_n(x)$$

で解が表されることにして，強制振動問題を解こう。ここで$\phi_n(x)$は 自由振動の振動モードであることには注意する。上式をその上の運動方程式に 代入すると

$$\sum_{n=1}^\infty \left\{ m \ddot{q}_n \phi_n +EI q_n \phi''''_n \right\}=p(x,t)$$

となる。この運動方程式と第$i$次モードとの内積をとるために

$$\int_0^\ell \phi_i \dint x \times$$

という作用を両辺にすると

$$\sum_{n=1}^\infty \int_0^\ell \left\{ m \ddot{q}_n \phi_n... ...n \phi_i \right\}\dint x= \int_0^\ell p(x,t) \phi_i \dint x$$

となる。これに，式(10.189) (10.190)の直交条件を 用いることによって，総和の第$i$次項のみを抽出でき

$$\overline{m}_i \ddot{q}_i(t) +\overline{m}_i \omega_i^2 q_i(t) = \overline{p}_i(t)$$ (10.200)

となる。ここに $\overline{p}_i(t)$は

$$\overline{p}_i(t)\equiv \int_0^\ell p(x,t) \phi_i(x) \dint x$$ (10.201)

で定義した一般化された外力 である。

例えば，1自由度系の式(10.48)の単位衝撃応答を利用した， 式(10.49)のDuhamel積分を用いて強制振動解を 求めてみると，上の式(10.201)の解は

$$q_i(t)=A_i \sin\omega_i t+B_i \cos\omega_i t +\int_0^\inf... ...au) \sin\left\{\omega_i\left(t-\tau\right)\right\} \dint\tau$$

となるので，最終的にたわみは

$$w(x,t)=\sum_{n=1}^\infty \left(A_n \sin\omega_n t+B_n \co... ...au) \sin\left\{\omega_n\left(t-\tau\right)\right\} \dint\tau$$

となる。これに，一般化された外力の式(10.202)をさらに 代入すると

$$w(x,t) = \sum_{n=1}^\infty \left(A_n \sin\omega_n t+B_n \cos\omega_n t\right) \phi_n(x)$$ (10.202)

$$\mbox{}\quad +\sum_{n=1}^\infty \int_0^\ell \int_0^\infty \dfrac{... ... \sin\left\{\omega_n\left(t-\tau\right)\right\} p(\xi,\tau) \dint\tau \dint\xi$$

と表現できる。Duhamel積分解は零初期条件の解だから， 未定係数$A_n$, $B_n$は式(10.193)と同じである。

図 10.45: 偏った位置の集中外力に対する応答

図-10.45には，単純梁の左から1/4の位置に時刻$t=0$に一定の 集中外力$P(t)=P_0$を載せ，そのままにしたときのスパン中央の応答と， いくつか代表的な時刻のたわみ形状を示した。外力は

$$p(x,t)=P_0 \delta\left(x-\slfrac{\ell}{4}\right)$$

だから，式(10.203)のDuhamel積分の$\xi$についての積分部分については

$$\int_0^\ell \phi_n(\xi) p(\xi,\tau)\dint\xi=P_0 \phi_n(\slfrac{\ell}{4})$$

となる。あとの計算は簡単なので省略する。 なお図中の$\tau$と縦軸の$w_0$は

$$\tau\equiv \dfrac{1}{\ell^2}\sqrt{\dfrac{EI}{m}}t, \quad w_0\equiv \dfrac{11P_0\ell^3}{768EI}$$ (10.203)

と定義した。$w_0$は静的に$P_0$を載せたときの中央のたわみである。 この図の応答の1周期はほぼ0.65程度に見える。一方，無次元時刻$\tau$で 算定したときの1次($n=1$)の固有振動数は$\pi^2$なので， それに対応する1次の固有周期は $\slfrac{2}{\pi}\simeq 0.637$である。 このことから，静的な変位を中心にして，ほとんど1次振動数で 振動しているだけであるが，外力の位置が偏っているため， たわみ形状は1次モードだけとは限らないことを示している。

10.4.6.2 単位衝撃応答

さて，式(10.203)を

$$w(x,t)=\sum_{n=1}^\infty \left(A_n \sin\omega_n t+B_n \co... ...infty w\subsc{i}(x;\xi,t-\tau) p(\xi,\tau) \dint\tau \dint\xi$$ (10.204)

という形で表現してみる。ここに

$$w\subsc{i}(x;\xi,t-\tau)\equiv \sum_{n=1}^\infty \dfrac{1}{... ... H(t-\tau) \sin\left\{\omega_n\left(t-\tau\right)\right\}$$ (10.205)

と置いた。式(10.205)の右辺の最後の項は，「場所$\xi$に 時刻$\tau$に作用した外力$p(\xi,\tau)$の影響を$w\subsc{i}$倍して 足し算すれば，それが外力すべての貢献になる」というようには 読めないだろうか。もしそうなら，式(10.206)で 定義した$w\subsc{i}$は，その影響の度合いを表す影響因子であり， あるいは，場所$\xi$に時刻$\tau$に何か単位の衝撃が作用した 場合の応答であるとは読めないだろうか。 つまり， $w\subsc{i}(x;\xi,t-\tau)$は梁の「単位衝撃応答」では ないかと考えられる。そして，式(10.205)の外力に関する項は， 外力の貢献分を積分表現した，梁のDuhamel積分 だろう。 ここでは，式(10.206)が本当に単位衝撃応答かどうかを 確認しておこう。

図 10.46: 梁の単位衝撃応答

梁のある位置$x=\xi$に，時刻$t=\tau$に単位衝撃外力が 作用したとすると，それは図-10.46のような 状況になる。ここでは，1自由度系の単位衝撃応答を考えた際に 説明したDiracのデルタ関数で外力が表すことができることを利用している。 簡単のために両端単純支持梁を対象とする。 また，節-4.4.1で説明したように， 集中外力も同じくDiracのデルタ関数で表すことができることから， この図のような単位衝撃外力は

$$p(x,t)=\delta\left(x-\xi\right) \delta\left(t-\tau\right)$$

と表現できることがわかる。 このような外力に対する単位衝撃応答を $w\subsc{i}(x,t)$と 記すことにすると，その運動方程式は

$$m \D[2]{w\subsc{i}(x,t)}{t}+EI \D[4]{w\subsc{i}(x,t)}{x}= \delta\left(x-\xi\right) \delta\left(t-\tau\right)$$

と表される。また境界条件は

$$w\subsc{i}(0,t)=0, \quad \D[2]{w\subsc{i}}{x}(0,t)=0, \qquad... ...{i}(\ell,t)=0, \quad \D[2]{w\subsc{i}}{x}(\ell,t)=0 \eqno{(a)}$$

であり，初期条件は零初期条件とする。ただし，Fourier変換を使うので， 初期条件は拡張して，時刻$t=-\infty$に静止していたものとして 与えられるものとする。

まず，解の単位衝撃応答をFourier変換したものを

$$W\subsc{i}(x;\alpha) \equiv \mbox{\fontMathScript F}(w\subsc... ...}(x,t) \exp\left(-\mbox{i}\alpha t\right)\dint t \eqno{(b)}$$

と記すことにする。次に運動方程式をFourier変換しよう。 まず左辺第1項のFourier変換は，部分積分を2回すれば

$$\int_{-\infty}^\infty \D[2]{w\subsc{i}}{t} \exp\left(-\m... ...\infty w\subsc{i} \exp\left(-\mbox{i}\alpha t\right)\dint t$$

となるので，初期条件と式($b$)を代入すると

$$=\left(\D{w\subsc{i}}{t}+\mbox{i}\alpha w\subsc{i}\right)\... ...}\alpha t\right)\Bigr\vert _{t\to \infty}-\alpha^2 W\subsc{i}$$

となる。 最終的にFourier逆変換の積分をする際には，複素積分を使う必要が あるので，この$\alpha$も複素数に拡張して $\alpha=a+\mbox{i}b$と すると， $\exp\left(-\mbox{i}\alpha t\right)= \exp\left(bt\right)\exp\left(-\mbox{i}a t\right)$となるので， 上式第1項の$t\to\infty$の項を零にするためには

$$b<0 \quad \to \quad \Im(\alpha)<0 \eqno{(c)}$$

という条件，つまり，$\alpha$の虚数部$\Im(\alpha)$は 負である必要がある。この条件下で第1項のFourier変換は

$$\int_{-\infty}^\infty \D[2]{w\subsc{i}}{t} \exp\left(-\mbox{i}\alpha t\right)\dint t =-\alpha^2 W\subsc{i} \eqno{(d)}$$

になる。次に左辺第2項のFourier変換は，微分とFourier変換の 順番を入れ替えるだけなので，式($b$)を代入して

$$\int_{-\infty}^\infty \D[4]{w\subsc{i}}{x} \exp\left(-\mbox{i}\alpha t\right)\dint t =\D*[4]{W\subsc{i}}{x} \eqno{(e)}$$

となる。最後に右辺のFourier変換は，Diracのデルタ関数の定義から

$$\int_{-\infty}^\infty \delta(x-\xi) \delta(t-\tau) \exp... ...delta(x-\xi) \exp\left(-\mbox{i}\alpha \tau\right) \eqno{(f)}$$

となる。したがって，この式($d$) ($e$) ($f$)を運動方程式の 形に組みなおすと

$$EI \D*[4]{W\subsc{i}}{x}-m \alpha^2 W\subsc{i} =\delta(x-\xi) \exp\left(-\mbox{i}\alpha \tau\right)$$

となる。境界条件式($a$)も同様にFourier変換してしまうと， 結局$W\subsc{i}$は

$$\D*[4]{W\subsc{i}}{x}-\dfrac{m}{EI} \alpha^2 W\subsc{i} =\... ...delta(x-\xi) \exp\left(-\mbox{i}\alpha \tau\right) \eqno{(g)}$$

の微分方程式を満足し，境界条件が

$$W\subsc{i}(0;\alpha)=0, \quad \D*[2]{W\subsc{i}}{x}(0;\alpha... ...{i}(\ell;\alpha)=0, \quad \D*[2]{W\subsc{i}}{x}(\ell;\alpha)=0$$

で与えられる。

これをモード解析法で解くことにし，$W\subsc{i}$を式(10.181)の 固有振動モードの級数解として

$$W\subsc{i}(x;\alpha)=\sum_{n=1}^\infty c_n \phi_n(x) \eqno{(h)}$$

と仮定する。これを式($g$)に代入すると

$$\sum_{n=1}^\infty \left\{ c_n \phi_n'''' -c_n \dfrac{m}{E... ...c{1}{EI} \delta(x-\xi) \exp\left(-\mbox{i}\alpha \tau\right)$$

となる。プライムは$x$による微分である。 この左辺第1項に式(10.186)を代入すると

$$\sum_{n=1}^\infty c_n \left\{\dfrac{m}{EI} \omega_n^2 -\d... ...c{1}{EI} \delta(x-\xi) \exp\left(-\mbox{i}\alpha \tau\right)$$

あるいは

$$\sum_{n=1}^\infty c_n m \left(\omega_n^2-\alpha^2 \right) \phi_n(x) =\delta(x-\xi) \exp\left(-\mbox{i}\alpha \tau\right)$$

となる。これと$\phi_i(x)$との内積をとるために，両辺に

$$\int_0^\ell \phi_i \dint x \times$$

という作用をし，式(10.189)の直交条件を 用いれば，総和の第$i$次項のみを抽出でき

$$c_i \overline{m}_i \left(\omega_i^2-\alpha^2 \right) =\in... ...) \dint x =\phi_i(\xi) \exp\left(-\mbox{i}\alpha \tau\right)$$

となる。右辺の演算ではDiracのデルタ関数の定義を使った。したがって

$$c_i=\dfrac{1}{\overline{m}_i \left(\omega_i^2-\alpha^2 \right)} \phi_i(\xi) \exp\left(-\mbox{i}\alpha \tau\right)$$

と求められる。これを式($h$)に戻すことによって

$$W\subsc{i}(x;\alpha)=\sum_{n=1}^\infty \dfrac{\phi_n(x) ... ...n^2-\alpha^2 \right)} \exp\left(-\mbox{i}\alpha \tau\right)$$

を得る。1自由度系の式(10.55)の周波数応答関数の形や， 式(10.60)に示した単位衝撃応答と周波数応答関数の関係から， この$W\subsc{i}$は，梁の周波数応答関数に「相当」することがわかる。

これで，単位衝撃応答のFourier変換が求められたから， これを逆変換すれば解を得ることができる。つまり

$$w\subsc{i}(x,t)=\mbox{\fontMathScript F}^{-1}(W\subsc{i}) =\... ...c{i}(x;\alpha) \exp\left(\mbox{i}\alpha t\right)\dint\alpha$$

であるから，これに上式を代入すると

$$w\subsc{i}=\dfrac{1}{2\pi}\int_{-\infty}^\infty \exp\left(\... ...au\right)\right\}} {\alpha^2-\omega_n^2}\dint\alpha \eqno{(i)}$$

となる。この右辺の最後の積分

$$I\equiv \int_{-\infty}^\infty f(\alpha)\dint\alpha,\quad f(... ...ox{i}\alpha\left(t- \tau\right)\right\}} {\alpha^2-\omega_n^2}$$

が求められれば，単位衝撃応答を求めることができる。 この積分$I$は，$\alpha$を複素数と捉えた複素積分で求められる。 分母が $\alpha=\pm \omega_n$で零になり，複素平面上のそこは1位の極 である。また式($c$)の条件から， 積分経路は実軸の若干下（$\alpha$の虚数部が負でなければならないから）を 通ることに注意する必要がある。つまり図-10.47(a)に 示した$\Gamma$の経路を，$-\infty$から$\infty$まで 線積分した $\int_\Gamma f(\alpha)\dint\alpha$が$I$の 値になる。$\alpha$は複素数なので $\alpha=a+\mbox{i}b$と すると， $\exp\left\{\mbox{i}\alpha\left(t- \tau\right)\right\} =\exp\left\{\mbox{i}\left... ...ft(t- \tau\right)\right\} \uwave{\exp\left\{-b \left(t- \tau\right)\right\}}$となる。 したがって，この積分を複素平面で1周積分する場合， 無限遠点での$f(\alpha)$を考察すると

1. もし$t<\tau$のときは， $f(\alpha)\to 0$になるのは$b\to -\infty$の ときなので，図-10.47(b)のように 負の虚数領域で1周積分することを考える。すると
   
   $$\int_{\Gamma^-} f(\alpha)\dint\alpha \to 0 \eqno{(j)}$$
   
   

   を期待できる。

2. 逆に，もし$t>\tau$のときは， $f(\alpha)\to 0$になるのは$b\to +\infty$の ときなので，図-10.47(c)のように 正の虚数領域で1周積分することを考える。すると
   
   $$\int_{\Gamma^+} f(\alpha)\dint\alpha \to 0 \eqno{(k)}$$
   
   

   を期待できる。

これより，$t<\tau$のときは周積分( $-\Gamma\to\Gamma^-$)の 領域内に極は存在せず，$f(\alpha)$は正則関数であるから，Cauchyの積分定理 より

$$\ointctrclockwise_{-\Gamma+\Gamma^-} f(\alpha)\dint\alpha =... ...} f(\alpha)\dint\alpha +\int_{\Gamma^-} f(\alpha)\dint\alpha=0$$

となるので，式($j$)を考慮すると，結局

$$-I=\int_{-\Gamma} f(\alpha)\dint\alpha=0$$

と求められる。 つまり$w\subsc{i}=0$であり，式(10.206)と一致する。

図 10.47: 単位衝撃応答のFourier逆変換で用いる積分路

逆に$t>\tau$のときは周積分( $\Gamma\to\Gamma^+$)の 領域内に極が二箇所あることから，一周積分は留数 の定理で

$$\ointctrclockwise_{\Gamma+\Gamma^+} f(\alpha)\dint\alpha = \... ...a_n\right] + 2\pi\mbox{i} \mbox{Res}\left[f;-\omega_n\right]$$

のように求めることができる。 ここに $\mbox{Res}\left[f;z\right]$は，$\alpha=z$に おける関数$f$の留数を表している。したがって，上式($k$)を考慮すれば

$$I= \int_{\Gamma} f(\alpha)\dint\alpha = 2\pi\mbox{i} \mbox... ...a_n\right] + 2\pi\mbox{i} \mbox{Res}\left[f;-\omega_n\right]$$

となる。1位の極の留数は

$$\mbox{Res}\left[f;\pm\omega_n\right]=\lim_{\alpha\to\pm\omeg... ...\{\pm\mbox{i}\omega_n\left(t- \tau\right)\right\}} {2\omega_n}$$

と求められる。 したがって， $o_n\equiv\omega_n\left(t- \tau\right)$として，Eulerの 公式(10.62)を用いれば

$$I=\dfrac{2\pi\mbox{i}}{\omega_n} \dfrac{\exp\left\{\mbox{... ...i}{\omega_n} \sin\left\{\omega_n \left(t-\tau\right)\right\}$$

が求めたかった積分値になる。 これを式($i$)に代入して，$t>\tau$の単位衝撃応答は

$$w\subsc{i}(x;\xi,t-\tau) =\sum_{n=1}^\infty \dfrac{1}{\o... ...,\phi_n(\xi) \sin\left\{\omega_n \left(t-\tau\right)\right\}$$

と表すことができる。これは式(10.206)に一致する。

10.4.6.3 減衰を含んだDuhamel積分

ここでは粘性もある場合の定式化をしておく。 運動方程式は式(10.196)であり， 境界条件は式(10.197)で与えられる。 ここでは簡単のために両端単純支持梁を対象とする。 初期条件は零初期条件である。一方， 単位衝撃応答 $w\subsc{i}(x;x_0,t-t_0)$は 同じ境界条件・初期条件と

$$m \D[2]{w\subsc{i}(x;x_0,t-t_0)}{t} +\eta I \D[5][4][x]{... ...sc{i}(x;x_0,t-t_0)}{x}=\delta(t-t_0) \delta(x-x_0) \eqno{(l)}$$

で表される運動方程式を満足している。 これに対し，節-10.1.3 (5)で 用いたような随伴問題 も考える必要があり，たわみ $w^*\subsc{i}(x;x_1,t-t_1)$が

$$m \D[2]{w^*\subsc{i}(x;x_1,t-t_1)}{t} -\eta I \D[5][4][x... ...sc{i}(x;x_1,t-t_1)}{x}=\delta(t-t_1) \delta(x-x_1) \eqno{(m)}$$

で表される運動方程式を満足する。自己随伴系ではないことから， 随伴問題は上式の元の問題の減衰項の 符号が異なることに加えて，初期条件ではなく，$t=-\infty$の 終末条件が零条件として

$$w^*\subsc{i}(x,\infty)=0, \quad \D{w^*\subsc{i}}{t}(x,\infty)=0 \eqno{(n)}$$

と与えられるような終局値問題として定義される。境界条件も 減衰項の符号を変えなければならない( $\eta \to -\eta$)ことに注意する。

最初に元の問題と随伴問題の単位衝撃応答$w^*\subsc{i}$との内積をとると

$$\int_0^\infty\int_0^\ell w^*\subsc{i} \left\{m \D[2]{w}{... ... t=\int_0^\infty\int_0^\ell w^*\subsc{i} p(x,t)\dint x\dint t$$

であるが，まず左辺の慣性項は

$$\int_0^\infty w^*\subsc{i} m \D[2]{w}{t} \dint t =\left.\... ...} w\dint t =\int_0^\infty m \D[2]{w^*\subsc{i}}{t} w\dint t$$

となる。ここでは，$w$の零初期条件と式($n$)の零終末条件を用いた。 運動方程式の粘性項と剛性項はまず$x$で2回部分積分して， 両端の境界条件を考慮すると

$$\begin{eqnarray*} && \int_0^\ell w^*\subsc{i}\left( \eta I\D[5][4][x]{w}{t} +E... ...{x}\left( \eta I\D[3][2][x]{w}{t} +EI\D[2]{w}{x}\right)\dint x \end{eqnarray*}$$

と表すことができるので，一旦は

$$\int_0^\infty \int_0^\ell \D[2]{w^*\subsc{i}}{x}\left( \eta I\D[3][2][x]{w}{t} +EI\D[2]{w}{x}\right)\dint x\dint t$$

となる。この第1項を時間方向に1回部分積分すると

$$\mbox{(第1項)}= \left.\D[2]{w^*\subsc{i}}{x}\eta I \D[2]{w}{... ...0^\infty \eta I \D[3][2][x]{w^*\subsc{i}}{t}\D[2]{w}{x}\dint t$$

と書き直せる。 ここでは$w$の零初期条件から $\D[2]{w}{x}(x,0)=0$とし，$w^*\subsc{i}$の 零終末条件から $\D[2]{w^*\subsc{i}}{x}=0$も成立するとした。 したがって，粘性と剛性の項は

$$\int_0^\infty \int_0^\ell \D[2]{w}{x}\left( -\eta I\D[3][2... ...*\subsc{i}}{t} +EI\D[2]{w^*\subsc{i}}{x}\right)\dint x\dint t$$

となる。再度$x$について2回部分積分して境界条件を考慮すると

$$\begin{eqnarray*} && \int_0^\ell \D[2]{w}{x}\left( -\eta I\D[3][2][x]{w^*\subsc... ...[4][x]{w^*\subsc{i}}{t} +EI\D[4]{w^*\subsc{i}}{x}\right)\dint x \end{eqnarray*}$$

となる。この結果と慣性項を一緒にまとめると，最終的に

$$\int_0^\infty\int_0^\ell w(x,t) \left\{m \D[2]{w^*\subsc... ... t=\int_0^\infty\int_0^\ell w^*\subsc{i} p(x,t)\dint x\dint t$$

と表現できるので，式($n$)を左辺の中括弧に代入すると

$$\begin{array}{l} \displaystyle \int_0^\infty\int_0^\ell w(... ...sc{i}(x;x_1,t-t_1) p(x,t)\dint x\dint t \end{array}\eqno{(o)}$$

という関係が成り立つ。

全く同じプロセスを，単位衝撃問題の式($l$)と 随伴問題の単位衝撃応答$w^*\subsc{i}$との内積に対して実行すると

$$\int_0^\infty\int_0^\ell w\subsc{i} \left\{m \D[2]{w^*\s... ...ell w^*\subsc{i} \delta(x-x_0) \delta(t-t_0)\dint x\dint t$$

となるので，左辺の中括弧に式($m$)を代入すると

$$\begin{array}{l} \displaystyle \int_0^\infty\int_0^\ell w\... ...,t_1-t_0)= w^*\subsc{i}(x_0;x_1,t_0-t_1) \end{array}\eqno{(p)}$$

という関係を得る。 つまり相反定理 だ。 この式($p$)を式($o$)に 代入することによって

$$w(x_1,t_1)= \int_0^\infty\int_0^\ell w\subsc{i}(x_1;x,t_1-t... ..._0^\ell w\subsc{i}(x;\xi,t-\tau) p(\xi,\tau)\dint\xi\dint\tau$$ (10.206)

というDuhamel積分の表現を得る。 これは式(10.205)の右辺第2項に一致する。 式(10.199)の粘性減衰定数を用いて，1自由度系の 単位衝撃応答の解の式(10.47)を参考にすると， 減衰を含んだ梁の単位衝撃応答は

$$w\subsc{i}(x;\xi,t-\tau)\equiv \sum_{n=1}^\infty \dfrac{1}{... ...u) \sin\left\{\omega\super{d}_n \left(t-\tau\right)\right\}$$ (10.207)

になる。ここに

$$\omega\super{d}_n\equiv \omega_n \sqrt{ 1-\left(\overline{\beta}_n{}\supersc{be}\right)^2 }$$ (10.208)

と定義した。

10.4.6.4 例--走行荷重

図 10.48: 走行荷重

例として，一定の速度$V_0$で$x$の正方向に走る一輪車（集中外力）$P_0$による 単純支持梁の非減衰振動を解析してみよう。一輪車が梁の左端に載った瞬間を 時刻$t=0$とすると，図-10.48に示したように，$t$秒後には 一輪車は$x=V_0 t$の位置にあるので，この外力は

$$p(x,t)=P_0 \delta(x-V_0 t)$$

と表すことができる。$P_0$は単位長さ当たりの分布外力として表した 集中外力（力積のようなもの）なので，力の単位を持つことに注意する。

これをDuhamel積分の式(10.207)に 代入すると

$$w(x,t)= \int_0^\infty\int_0^\ell w\subsc{i}(x;\xi,t-\tau) P_0 \delta(\xi-V_0 \tau)\dint\xi\dint\tau$$

になる。外力の変数は$\xi$と$\tau$になっていることに注意する。 これに式(10.208)の単位衝撃応答の 減衰定数を $\overline{\beta}_n{}\supersc{be}=0$としたものを代入し， 単純支持梁の固有振動モードの式(10.181)も 代入すると， $\overline{m}_n=\dfrac{m\ell}{2}$なので

$$w(x,t)=\sum_{n=1}^\infty\dfrac{2P_0}{m\ell\omega_n} \sin... ...int_0^\infty H(t-\tau)\sin\omega_n\left(t-\tau\right)\dint\tau$$

となる。$\xi$についての積分項を取り出すと，Diracのデルタ関数の 定義によって，この項は積分範囲内に$V_0\tau$があるとき， つまり $0<V_0\tau<\ell$のときにしか値を持たず，それ以外では零に なることがわかる。つまり， $\tau<\slfrac{\ell}{V_0}$のときだけ値を 持つことができるので

$$\int_0^\ell \sin\left(\dfrac{n\pi \xi}{\ell}\right)\delta(\x... ...-\tau\right) \sin\left(\dfrac{n\pi V_0 \tau}{\ell}\right)$$

となる。これは当たり前だろう。というのも $\slfrac{\ell}{V_0}$は， 一輪車が梁の外に出る時刻だからである。ここで

$$k_n\equiv \dfrac{n\pi V_0}{\ell}, \quad t_0\equiv \dfrac{\ell}{V_0}$$ (10.209)

と置くと，上の解は

$$w(x,t)=\sum_{n=1}^\infty\dfrac{2P_0}{m\ell\omega_n} \sin... ...) \sin k_n \tau \sin\omega_n\left(t-\tau\right)\dint\tau$$

と書くことができる。二つのHeaviside関数からは

$$t>\tau, \quad t_0>\tau \quad \to \quad \min\left(t, t_0\right)>\tau$$

の範囲で積分すればいいことを示している。したがって上式は

$$w(x,t)=\sum_{n=1}^\infty\dfrac{2P_0}{m\ell\omega_n} \sin... ...)} \sin k_n \tau \sin\omega_n\left(t-\tau\right)\dint\tau$$

を計算すればいいことになる。$s$を$t$と$t_0$の 小さい方として

$$\int_0^s \sin k_n \tau \sin\omega_n\left(t-\tau\right)\dint\tau$$

の積分ができればいいだけであるが， 冗長になるので，その計算は省略して結果だけを示すと

(10.210)

となる。

図 10.49: 走行荷重による応答

図 10.50: 高速走行荷重に対する応答

つまり，一輪車が載っている間は固有振動数$\omega_n$と 外力の運動に起因する振動数$k_n$との 両方の特性で動いているが，通過したあとは固有振動数だけで 自由振動する。図-10.49には，梁の 代表的な3点の応答を示した。 この図は級数10項を用いた結果であるが，100項でも区別はつかない。 また図で用いた記号は，式(10.204)で定義した$\tau$の他には

$$\zeta\equiv \dfrac{1}{\ell V_0}\sqrt{\dfrac{EI}{m}}, \quad w_0\equiv \dfrac{P_0\ell^3}{48EI}$$

と定義した。$w_0$はスパン中央に集中外力$P_0$が作用したときの 中央のたわみである。また，$\tau=\zeta$に一輪車は梁から外に出る。 ここでは，かなり速く通過する場合を$\zeta=3$とし， ゆっくり通過する通常の橋梁等の場合を$\zeta=10$とした。 どちらも$\tau=\zeta$以降は自由振動をしている。 実はこの式(10.211)の解はとても面白い現象を示唆している。 それは，走行速度によっては分母を零にする可能性があるからである。 つまり$k_n=\omega_n$のとき，この解の振幅は無限大になる。 例えば，$n=1$のときには式(10.210)から

$$\dfrac{\pi V_0}{\ell}=\omega_1 \quad \to \quad \dfrac{\ell}... ...ga_1}=\dfrac12 T_1, \quad \left(\zeta=\dfrac{1}{\pi}\right)$$

となるから， つまり，一輪車が梁を通過するのに必要な時間が，梁の1次の 固有周期の半分のとき，解は発散解になることを示している。 これをZimmermann効果 と呼ぶ。ただし，現実的な問題の設定を考えると，これは 非現実的なくらい速い速度で一輪車が移動する場合に相当する。 例えば，狭い溝に鉄板を渡し，そこを新幹線くらいの高速車両が 通過する場合であろう。実際に $\zeta=\slfrac{1.001}{\pi}$の 解を図-10.50に示した。あっという間に 通過したあと，静的たわみ程度の振幅の自由振動になる。 文献[87]にも，振動が無限大になる前に一輪車が通過するので， 実際には振動が無限大になることはないと説明されている・・・が，面白い。

1. 図-10.45の問題を自ら解いてみて， 図の結果が求められることを確認せよ。
2. 図-10.45と同じ状況で， 外力を $P(t)=P_0 \sin pt$としたときの応答を求めよ。
3. 図-10.51に示したのは，速度$V_0$で 走行する2輪車による応答である。 ただし，スパン30mの橋梁に，いわゆるT荷重が載った情況を念頭に置いて
   
   $$\dfrac{d}{\ell}=0.13, \quad \dfrac{Q}{P}=4, \quad w_0\equiv \dfrac{\left(P+Q\right) \ell^3}{48EI}$$
   
   

   としてみた。Heaviside関数とDiracのデルタ関数の使い方の練習問題として， これを求めてみよ。 結果はあまり面白いものではなく，一輪車の結果を$d$だけ移動させたくらいの ものである。

図 10.51: 2輪の走行荷重による応答

10.4.7 有限要素と振動解析

10.4.7.1 仮想仕事式と質量行列および減衰行列

実際に複雑な社会基盤構造の設計時の振動解析で 偏微分方程式を直接用いることは， モデル化の点からもそれを解く観点からも無理がある。 したがって，現実的な方法は数値的に近似的に解く方法であろう。 そのためには，式(10.196)の運動方程式と 式(10.197)の境界条件を，何らかの方法で離散化しておく 必要がある。一つの代表的な手法は有限要素法の適用である。 ここでは梁を対象として，その要素の定式化と解析例を示す。

まず，長さ$\ell$の1要素を対象とし，式(10.196)に 対応する仮想仕事式を求めると

$$\int_0^\ell \delta w \left\{ m\D[2]{w}{t}+\D[2]{}{x}\left(\... ...ht) +\D[2]{}{x}\left(EI\D[2]{w}{x}\right)-p \right\}\dint x=0$$

となる。被積分関数の第2, 3項を2回部分積分すると，上式は

$$\begin{eqnarray*} &&\int_0^\ell \left\{ \delta w m\D[2]{w}{t}+\delta \D[2]{w}{... ...D[3][2][x]{w}{t} +EI \D[2]{w}{x} \right)\right\vert _0^\ell = 0 \end{eqnarray*}$$

となる。ここでは，たわみ角$\theta$が式(10.173)で， 反時計回りを正に定義されていたことを考慮してある。 この式に式(10.197)の境界条件を代入すると， 最終的に仮想仕事式は

$$\int_0^\ell \left\{ \delta w m\D[2]{w}{t}+\delta \D[2]{w}{... ...elta \theta_1 C_1 -\delta w_2 S_2 -\delta \theta_2 C_2 = 0$$ (10.211)

と表される。ここに，$w_1$, $\theta_1$は左端$x=0$のたわみと たわみ角であり，$w_2$, $\theta_2$は右端$x=\ell$の各量である。 また，$S_i$, $C_i$も両端の集中せん断外力と集中外力モーメントである。

非減衰振動の一般解は三角関数と双曲線関数であったが， 有限要素法の特徴は，変位を区分的多項式で近似することにあったから， ここでも静的な問題の式(5.21)と同様に，たわみの変位関数を

$$\left\{\begin{array}{ll} \psi_1(x)\equiv 1-\dfrac{3x^2}{\ell^2} +... ...3}, & \psi_4(x)\equiv \dfrac{x^2}{\ell} -\dfrac{x^3}{\ell^2} \end{array}\right.$$ (10.212)

と定義した上で

$$w(x,t)= w_1(t) \psi_1(x)+\theta_1(t) \psi_2(x) +w_2(t) \psi_3(x)+\theta_2(t) \psi_4(x)$$ (10.213)

と置くことにする。この式(10.214)を 上の仮想仕事式(10.212)に代入すると

$$\delta \fat{w}_b(t)\supersc{t} \left\{ \fat{m}_b \ddot{... ...{k}_b \fat{w}_b(t) - \fat{f}_b(t) - \fat{g}_b(t) \right\}=0$$

になるので，任意の仮想変位 $\delta \fat{w}_b(t)$に対してこれが 成立しなければならないことから，運動方程式は

$$\fat{m}_b \ddot{\fat{w}}_b(t) + \fat{c}_b \dot{\fat{w}}_b(t) + \fat{k}_b \fat{w}_b(t) = \fat{f}_b(t) + \fat{g}_b(t)$$ (10.214)

と表現できる。ここに，要素の変位ベクトルと二つの外力ベクトルを

$$\fat{w}_b(t) \equiv \lfloor w_1(t) \; \theta_1(t) \; w_2(t)... ... \lfloor q_1(t)\; q_2(t)\; q_3(t)\; q_4(t) \rfloor\supersc{t}$$ (10.215)

と定義した。また最後の分布外力の等価節点外力（成分の 下添え字$b$は略した）は，式(5.22b)と同様に

$$q_i(t)\equiv \int_0^\ell p(x,t) \psi_i(x)\dint x, \quad i=1,2,3,4$$ (10.216)

と定義した。この定式化で得ることができる質量行列$\fat{m}_b$は

$$\fat{m}_b\equiv \left(\int_0^\ell m \psi_i(x) \psi_j(x)\d... ...mn{3}{l}{\mbox{Symm.}}&\slfrac{\ell^3}{105} \end{array}\right)$$ (10.217)

と定義できる。多自由度系の集中質量行列（対角行列）とは 異なり，$4\times 4$のすべての成分が非零の行列になっている。 このような質量行列を整合質量行列 と呼ぶことがある。また剛性行列は，式(5.24)と同じで

$$\fat{k}_b\equiv\left(\int_0^\ell EI \psi_i''(x) \psi_j''(... ...column{3}{l}{\mbox{Symm.}}&\slfrac{4}{\ell} \end{array}\right)$$ (10.218)

と定義される。プライムは$x$による微分を表す。 また，粘性を材料そのものの粘性モデルで定式化した 有限要素の場合には，その減衰行列はRayleigh減衰の式(10.119)の うちの剛性比例減衰になり

$$\fat{c}_b\equiv \left(\int_0^\ell \eta I \psi_i''(x) \ps... ...column{3}{l}{\mbox{Symm.}}&\slfrac{4}{\ell} \end{array}\right)$$ (10.219)

と定義できる。すべて，1要素は一様断面で一様材料でできていることを 仮定して簡略化した。

式(10.215)は要素の運動方程式なので， あとは要素間の連続条件を用いて，節-5.3.2 (1)の直接剛性法を用いて， 全体構造の運動方程式を組み立てればいい。 最終的には多自由度系の運動方程式(10.95)と全く同じ

$$\fat{M} \ddot{\fat{u}}(t)+\fat{C} \dot{\fat{u}}(t) +\fat{K} \fat{u}(t)=\fat{f}(t)$$ (10.220)

の形の運動方程式になる。外力ベクトル$\fat{f}(t)$には， 集中外力$\fat{f}_b(t)$と分布外力による等価節点外力$\fat{g}_b(t)$の 両方をまとめた。 多自由度系と異なるのは，このようにして全体構造の運動方程式を 組み立てたあと，幾何学的な境界条件（自由度が固定されている 条件）をこの運動方程式に設定しないといけないことである。 例えば，上式(10.221)が見かけ上$N$自由度の 運動方程式になっていると考えて，その第$n$自由度が拘束されている場合には， それに対応する外力は反力になるので未知数である。 そこで，その行の方程式を，第$n$自由度の加速度を零にする条件式

$$\ddot{u}_n(t)=0$$

と置き換えることにする。つまり， 三つの行列と外力ベクトルの成分のうち， $i=1,2,\cdots N$に対して

$$M_{ni}=0, \quad M_{in}=0, \quad C_{ni}=0, \quad C_{in}=0, \q... ..., \quad K_{in}=0, \quad f_i=0 \quad\mbox{ただし}\quad M_{nn}=1$$ (10.221)

といった操作をすれば，対応する第$n$自由度を拘束したことになる。 ここに，例えば$M_{ni}$は行列$\fat{M}$の第$n$行第$i$列成分を表している。

10.4.7.2 自由振動解析

振動解析では，固有振動数と振動モードを知ることが最も重要であったから， まず自由振動解析をする必要がある。 外力を零にしておいて周期解を探せばいいだけなので

$$\fat{u}(t)=\fat{U} \exp(\mbox{i}\omega t)$$

と置けばいい。ここに$\fat{U}$は振幅列行列である。 これを非減衰で外力の無い式(10.221)に代入すると

$$\left( \fat{K} - \omega^2 \fat{M} \right) \fat{U} = \fat{0}$$

を満足していればいいことになる。したがって， 有意な$\fat{U}$が存在するための条件から

$$\det \left\vert\fat{K} - \omega^2 \fat{M} \right\vert = 0$$ (10.222)

でなければならず，この式が固有振動数$\omega$を決定する振動数方程式になる。 つまり，多自由度系の固有値問題と形式的には同じ形をしている。

例えば両端単純支持梁に1要素を用いた場合には，上式(10.223)は， 式(10.215)の第2, 4行だけの方程式

$$\left\{\begin{array}{c} 0 0 \end{array}\right\}= \left( E... ... \right)\left\{\begin{array}{c} U_2 U_4 \end{array}\right\}$$

となるので，振動数方程式(10.223)は

$$\det\left\vert \left(\begin{array}{cc} 4 & 2 \\ 2 & 4 \... ...{-1}{140} & \slfrac{1}{105} \end{array}\right) \right\vert=0$$

となる。$U_2$, $U_4$はそれぞれ$\theta_1(t)$, $\theta_2(t)$の振幅である。 これを解くと

$$\omega \ell^2\sqrt{\dfrac{m}{EI}}=2\sqrt{30}, \quad 6\sqrt{70}$$

と求めることができる。これが単純支持梁の第1, 2次の固有振動数に 相当するので，式(10.182)の厳密解と比較すると

$$\omega_1\ell^2\sqrt{\dfrac{m}{EI}}= 2\sqrt{30}\simeq 1.1099... ...\dfrac{m}{EI}}=6\sqrt{70} \simeq 1.27157 \left(2\pi\right)^2$$

程度になる。つまり1次の振動数が11%の誤差，2次の振動数が27%の 誤差で求められたことを意味する。

表 5.3: 固有振動数の近似度の改善

要素数	$\zeta_1$	$\zeta_2$	$\zeta_3$	$\zeta_4$	$\zeta_5$
1	1.10992	1.27157	--	--	--
2	1.00395	1.10992	1.23994	1.27157	--
4	1.00026	1.00395	1.01827	1.10992	1.12909
8	1.00002	1.00026	1.00129	1.00395	1.00927
16	1.00000	1.00002	1.00008	1.00026	1.00063

表-10.2には， 両端単純支持梁の固有振動数を，用いる要素数を増やしながら有限要素法で 求めたのものを示した。 ただし，表の$\zeta_k$は第$k$次振動数の厳密解と有限要素解の 比で， $\zeta_k \equiv \dfrac{\mbox{有限要素解}}{\mbox{厳密解}}= \dfrac{\omega_k}{(k\pi/\ell)^2 \sqrt{EI/m}}$と定義した。要素数が増えるにつれて 解の精度が上がっていくことが明らかである。 工学的に最も重要な最低次の固有振動数は，2要素（4自由度）で は1%未満になっていることは興味深い。 振動数の近似解が厳密解より大きくなるのは，近似することによって モデルの変形が真の変形より拘束されたものになって， その見かけ上の剛性が大きくなってしまうからである。

図 10.52: 両端単純支持梁の近似振動モード

もちろん，各固有振動数に対応する固有ベクトルは，固有振動モードである。 また当たり前であるが， ここの定式化では$\fat{u}(t)$あるいは$\fat{U}$には たわみとたわみ角が交互に並んでいるので， 実際に描画する場合には注意が必要である。 一例として，上述の両端単純支持梁の有限要素を用いた自由振動解析で 得られた固有ベクトルから，奇数行のたわみだけを直線で結んで 近似の固有振動モードを描いてみよう。図-10.52が 梁を40要素で分割したときの最初の三つの振動モードである。 固有ベクトルのうちのたわみだけ（図では$W_n$と記した）を用いた 近似であるが，図-10.43を見るまでもなく， 正しく振動モードが得られているのが明らかである。

図 10.53: 中央差分のフローチャート
図 10.54: 走行荷重による応答の有限要素解
図 10.55: 連続梁上の走行荷重による応答

10.4.7.3 時刻歴応答解析

次に，地震応答解析のように時刻歴応答解析をする場合には， 離散化された運動方程式(10.221)を多自由度系の場合と 同じ手法で解けばいい。 例えば，節-10.2.3で説明した中央差分を 用いる場合のフローチャートを図-10.53に示した。 いくつか例を解いておこう。まず非減衰の場合の例として， 単純支持梁上の走行荷重による図-10.49の応答のうち， スパン中央のたわみを2要素と4要素で 解析した結果を図-10.54に示した。 走行荷重はある時刻$t$に$V_0t$の位置にいるが， それが$n$番目の要素上にあると判定できた場合， その要素左端の$x$座標が$x_n$であるとすると， その要素内の一輪車の位置$a$は$a=V_0t-x_n$となる。したがって， 式(10.217)の等価節点外力の定義から， 走行荷重はDiracのデルタ関数を用いて

$$q_i=\int_0^\ell P_0 \delta(x-a) \psi_i(x)\dint x=P_0 \psi_i(a)$$

によって算定できる。ただし，この$\ell$は1要素の長さである。 モード解析法に よる厳密解の一点鎖線と4要素による実線の結果は ほとんど重なっていて，2要素による 破線の結果も，拡大図に示した程度の差（ほとんど重なっていると 思うが）しかない。工学的には2要素で 十分な結果が求められていると考えられる。図-10.55には， モード解析が面倒な不静定梁の例として，2径間連続梁を一輪車が 走行する場合の，左のスパン中央のたわみを示した。$\tau=\zeta$で 右のスパンに移動するので，たわみが負（上向き）に移行し，$\tau=1.5\zeta$で 走行荷重は梁の外に出て，そのあとは自由振動をしている。

図 10.56: 偏った位置の集中外力に対する応答

図-10.56(a)には，偏った固定集中外力に 対する応答を示した図-10.45の結果との 比較を，4要素で解析した結果を用いて示した。 モード解析結果が破線で描いてあるが，スパン中央の時刻歴応答も， 代表的な時刻のたわみ形状も，ほとんど重なっていて， 十分な精度で求められていることがわかる。 なお，たわみ形状は，変位関数の式(10.214)に 節点のたわみとたわみ角の両方を用いてより正確に描いたものである。 また図-10.56(b)には， 同じ構造と外力条件で，減衰を含めた場合の応答を示した。 減衰定数は

$$\beta\supersc{be}\equiv \dfrac12 \dfrac{\eta}{E} \dfrac{... ...frac{\eta}{E}=2 \beta\supersc{be} \ell^2\sqrt{\dfrac{m}{EI}}$$ (10.223)

と定義し，この例では $\beta\supersc{be}=0.02$とした。 この定義は，式(10.199)の $\overline{\beta}_i{}\supersc{be}$の 定義とは若干異なるので注意する。 破線がモード解析による非減衰解であるが，減衰を含めた数値 解は，静的なたわみ$w_0$（縦軸の1）に徐々に近づいていくことを示している。 時刻$t=0$に$P_0$の荷物を $x=\slfrac{\ell}{4}$に置いままにしただけの 問題なので，最終的にはその静的たわみに収束していく。

このように，数値解析によって，どんな外力に対する応答も計算できる わけだが，少し注意して欲しいのは，それが必ずしも直接設計に有用な 情報を与えてくれるとは限らないことである。 むしろ，必要に応じて卓越するモードを探し， その成分の応答を個々に比較する方が，設計には役立つのである。 したがって，前節までのモード解析についての理解も重要となる。 特に3次元の実構造物のように，曲げだけではなくねじれ等も含む 複雑な振動をする場合には，着目している設計断面の応力を 増大させる振動モードを探し出し，そのモードを制御するような 設計変更をする必要があるだろう。 このような場合，漫然とその断面の数値的な応答を眺めているより， モード解析法のような考察を試みることの方が有用であろう。 たいていの場合，数値解析というのは，あるケーススタディに なってしまっている可能性があるので注意して欲しい。

1. 図-10.45と同じ状況で， 外力を$P(t)=P_0$としたときの応答を，有限要素を用いて求めよ。
2. 図-10.45と同じ状況で， 外力を $P(t)=P_0 \sin pt$としたときの応答を，有限要素を用いて求めよ。

10.4.7.4 平面骨組の振動解析

10.4.7.4.1 軸方向変形についての運動:

前節までに示した梁と同様， トラス部材のような軸力で抵抗する棒部材も多自由度振動系にモデル化できる。 それを梁と一緒に用いることによって任意の平面骨組の振動解析も可能に なるだろう。 そこでまず，棒部材の有限要素運動方程式を求めておこう。 式(5.13)の仮想仕事式に慣性項を加えることによって， 長さ$\ell$の1要素の仮想仕事式は

$$\int_0^\ell EA \D{u(x,t)}{x} \delta \D{u}{x}\dint x +\int... ... \delta u \dint x -F_1(t) \delta u_1-F_2(t) \delta u_2 = 0$$ (10.224)

となる。ここに$A$は断面積であり，$p_a(x,t)$は軸方向の単位長さ当たりの 分布外力，$F_i(t)$, ($i=1,2$)は端部の軸外力である。 用いる変位関数は静的な問題の式(5.15)と同じ多項式でいいから

$$u(x,t)\sim u_1(t) \phi_1(x)+u_2(t) \phi_2(x), \qquad \phi... ...v \dfrac{\ell-x}{\ell}, \quad \phi_2(x)\equiv \dfrac{x}{\ell}$$ (10.225)

としよう。これを上の仮想仕事式(10.225)に代入して整理すると， 離散化された運動方程式は

$$\fat{m}_a \ddot{\fat{u}}_a(t)+\fat{k}_a \fat{u}_a(t)= \fat{f}_a(t)+\fat{g}_a(t)$$ (10.226)

と表すことができる。ここに，要素の変位ベクトルと二つの外力ベクトルを

$$\fat{u}_a(t) \equiv \lfloor u_1(t) \; u_2(t) \rfloor\supers... ...fat{g}_a(t) \equiv \lfloor p_1(t)\; p_2(t) \rfloor\supersc{t}$$ (10.227)

と定義した。また最後の分布外力の等価節点外力は，式(5.16)と同様に

$$p_i(t)\equiv \int_0^\ell p_a(x,t) \phi_i(x)\dint x, \quad i=1,2$$ (10.228)

と定義した。この定式化で得ることができる質量行列$\fat{m}_a$は

$$\fat{m}_a\equiv \left(\int_0^\ell m \phi_i(x) \phi_j(x)\d... ...ac{1}{6} \\ \slfrac{1}{6} & \slfrac{1}{3} \end{array}\right)$$ (10.229)

と定義でき，梁の場合と同様整合質量行列になっている。 また剛性行列は，式(5.17)と同じで

$$\fat{k}_a\equiv\left(\int_0^\ell EA \phi_i'(x) \phi_j'(x)... ...} \left(\begin{array}{rr} 1 & -1 -1 & 1 \end{array}\right)$$ (10.230)

と定義される。プライムは$x$による微分を表す。 すべて，1要素は一様断面で一様材料でできていることを仮定して簡略化した。 上式(10.227)は，節-11.1の1次元の 波動，つまり棒や筒の中の疎密波動の有限要素解析に使える。 また，$EA$を張力$T$で置き換えれば節-10.3.1の 弦の振動の有限要素運動方程式でもある。

粘性減衰については，例えばその材料モデルを 式(10.194)と同じように考えてもいいので

$$N(x,t)=EA \D{u(x,t)}{x} +\eta A \D[2][1][x]{u(x,t)}{t}$$ (10.231)

とすればいい。ここに$\eta$は材料の粘性係数である。 このモデルの場合には，その減衰行列は

$$\fat{c}_a\equiv \left(\int_0^\ell \eta A \phi_i'(x) \phi... ...} \left(\begin{array}{rr} 1 & -1 -1 & 1 \end{array}\right)$$ (10.232)

と定義できる。

図 5.25: 平面骨組有限要素の要素座標系と全体座標系

10.4.7.4.2 任意の平面骨組の運動:

平面骨組は，圧縮力も受ける曲げ部材が任意の方向に配置されたもので あるが，節-5.3.4にあるように考えれば， その有限要素運動方程式も定式化できる。 そこで，上の式(10.227)と 曲げの式(10.215)を合体させてみよう。 ただし簡単のために減衰と分布外力は無視する。 骨組の各部材はそれぞれ別々の方向を向いて 配置されているので，図-10.57のように 要素の$\xi$-$\zeta$座標系と全体座標系との向きの 違いを角度$\alpha$で定義しておく。 この要素座標系での質量行列と剛性行列はそれぞれ， 式(10.218)と式(10.230)を，そして 式(10.219)と式(10.231)を合体させて

$$\fat{m}_f \equiv \renewedcommand{arraystretch}{0.7} \left(\begin{array}{ccc\vert c... ...\\ (\fat{m}_f\supersc{b})\supersc{t} & \fat{m}_f\supersc{c} \end{array}\right)$$ (10.233)

$$\fat{k}_f \equiv \renewedcommand{arraystretch}{0.7} \left(\begin{array}{ccc\vert c... ...\\ (\fat{k}_f\supersc{b})\supersc{t} & \fat{k}_f\supersc{c} \end{array}\right)$$

と定義できる。

これに対して全体$x$-$z$座標系における骨組有限要素の運動方程式 は

$$\overline{\fat{m}}_f \ddot{\overline{\fat{u}}}_f(t)+ \overline{\fat{k}}_f \overline{\fat{u}}_f(t)=\overline{\fat{f}}_f(t)$$ (10.234)

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

$$\overline{\fat{u}}_f(t) \equiv \lfloor \overline{u}_1(t) \;... ... \; \overline{S}_2(t) \; \overline{C}_2(t) \rfloor\supersc{t}$$ (10.235)

と定義した。そして全体座標系の質量行列と剛性行列は， 式(10.234)の小行列を用いて

$$\overline{\fat{m}}_f= \left(\begin{array}{cc} \fat{T}\super... ...}\supersc{t} \fat{k}_f\supersc{c} \fat{T} \end{array}\right)$$ (10.236)

と定義できる。ここに，座標変換行列を

$$\fat{T}\equiv\left(\begin{array}{ccc} \cos\alpha & -\sin\al... ...icolumn{1}{r}{\cos\alpha} & 0\\ 0 & 0 & 1 \end{array}\right)$$ (10.237)

と定義した。

図 10.58: 門型骨組の固有振動モード

10.4.7.4.3 自由振動の数値解析例:

最終的な運動方程式は多自由度系のそれと同じなので， 前節でも示したような自由振動解析も時刻歴応答解析も同様に可能である。 ここでは自由振動だけに限定して，基本的な構造系の固有振動数と 振動モードを示すことによって，その特徴を説明しておこう。 最初に，門型骨組を解析したのが図-10.58である。 ここでも固有ベクトルのうちのたわみ角 $\overline{\theta}_i$を 無視し，変位成分だけを用いて直線で結んで描いた近似的な振動モードである。 また固有ベクトルの変位の値も増幅させて誇張した図になっている。 まず柱が同じ高さの左右対称な門型骨組を 解析したのが図-10.58(a)である。 対応する固有振動数を表-10.3の第2カラムに示しておいた。 柱と梁をそれぞれ40要素で近似しているので， 六つの境界条件のことを考慮すれば，117自由度の多自由度系で近似したことになる。 柱と梁は同じ長さ$L$であるとし， この長さで定義した細長さ$\slfrac{L}{r}$を100とした。 ここに$r$は断面の回転半径 $r\equiv\sqrt{\slfrac{I}{A}}$である。

両端単純支持梁の場合と異なり，第3次固有振動数はそれほど大きくなっていない。 これは，3次の振動モードが2次の振動モードの一部が逆向きになったくらいの 差しか無いためである。また第4次固有振動数と3次のそれも近い。 これも3次と4次の振動モードが類似しているからである。 このように，振動モードを境界条件等から予測して，どのような次数の 固有振動数が近付くのかを推測できるようになると， 共振し難い構造系を設計することができるようになるかもしれないよ。

次に，右の柱が左の半分の高さしかない非対称な門型骨組を 解析した結果が図-10.58(b)である。 固有振動数は表-10.3の第3カラムに示した。 今度は第3次と4次の振動数が離れている。これは4次の振動モードが かなり複雑なものになったことが原因である。 それによって，第5次の振動数も大きくなったと考えられる。

表: 固有振動数 $\omega L^2\sqrt{\slfrac{m}{EI}}$

	門型骨組		2ヒンジアーチ		固定アーチ
次数	対称	非対称	$\theta_0=30$度	$\theta_0=22.5$度	$\theta_0=45$度	$\theta_0=36$度
1	3.204	6.181	38.79	36.29	60.06	55.47
2	12.62	14.78	47.05	39.09	66.14	60.63
3	20.62	21.58	90.28	89.49	124.6	122.9
4	22.28	45.24	157.1	157.5	197.4	198.3
5	44.79	58.11	246.4	246.6	298.0	298.2

次に，非常に浅い2ヒンジアーチ状の構造を 対象としてみた。図-10.59の左上にあるのが， 内角$\theta_0$が30度の浅いアーチを示している。 この円弧の長さを$L$として，これを用いた細長比を $\slfrac{L}{r}=100$と設定し， この円弧を40個の弦で結んで近似してある。図-10.59(a)には この内角が30度の場合の振動モードを示した。興味深いことに，構造系が対称であるにもかかわらず，最低次の第1次の 振動モードは反対称である。これは，構造系そのものの形状が直線梁の 第1次対称振動モードと同じであることが原因である。表-10.3の 第4カラムにある固有振動数を見ると， 両端単純支持梁に比べて，それぞれがあまり離れていない。 つまり，共振振動数が近付いていることになり，対振動の 設計が少し難しくなることを示している。

これに対し，内角が22.5度になるような，もっと浅いアーチ状の構造の 場合には，図-10.59(b)にあるように， 最低次の第1次振動モードが対称になる。この二つの図を見てすぐにわかるように， この内角の違いによって， 第1次と2次のモードだけが逆転していることが明らかである。 実は， $\theta_0\simeq 24.385$度で小さい方の二つの固有振動数が 一致（ $\omega L^2\sqrt{\slfrac{m}{EI}}\simeq 39.024$）する。 そのため，その角度に近い22.5度の場合の固有振動数の 分布を表-10.3の第5カラムで見ると， 小さい方の二つの振動数が近付いているのも納得がいく。

図 10.59: 浅い2ヒンジ・アーチ状の棒の固有振動モード

図 10.60: 固定アーチの内角と固有振動数

そこで次に，固定アーチの場合の内角と小さい方の二つの固有振動数の 関係を求めてみたのが図-10.60である。 横軸がアーチの内角であり，約40.092度を境に， 最低次の振動モードの対称と反対称が入れ替わっている。 そこで$\theta_0$が45度の場合と36度の場合について， 固有振動モードを図-10.61に， 固有振動数を表-10.3の第6, 7カラムに示した。 明らかに，第1次と2次の振動モードが入れ替わっている。 内角が45度の場合には，第2, 3, 5次がほぼ対称で第1, 4次が反対称になる。 これに対し36度の場合には，第2, 4次が反対称である。 いずれにしても，表-10.3を横向きに比べれば 明らかなように，アーチ状の構造の小さい方の振動数は， 門型骨組や両端単純支持梁に比べてお互いに近付いていることに 注意する必要があるだろう。 さらに実際の社会基盤構造はすべて立体的なものであり， もっと複雑な変形の自由度を有していることから， 固有振動数同士が近付くことは普通に観察される状況である。 両端単純支持梁のようなかけ離れた固有振動数の分布をすることの 方が稀であることには十分注意する必要がある。

図 10.61: 固定アーチ状の棒の固有振動モード

10.4.8 回転慣性とせん断変形

10.4.8.1 回転慣性の影響

実は，式(10.169)のモーメントのつり合いを考えたときには， 動的な効果を無視していた。 つまり，回転運動の運動方程式は 式(10.4)のようなNewtonの法則が成立するはずなのに， 回転慣性を無視していた。回転慣性を導入すると， 断面の回転角がたわみ角と同じなので，式(10.169)の代わりに

$$\D{M(x,t)}{x}-V(x,t)=J \D[2]{\theta(x,t)}{t}, \quad \theta(x,t)=-\D*{w}{x}$$

が正しい運動方程式になる。ここの慣性モーメント$J$は 式(10.5)の定義を梁の断面に適用すれば，長さ$\dint x$当たりに

$$J=\int_A \rho z^2\dint A=\rho I$$

となる。ここに$\rho$は密度で，$I$は断面2次モーメントである。 したがって，モーメントの運動方程式は本当は

$$\D{M(x,t)}{x}-V(x,t)=\rho I \D[3][1][x]{\theta(x,t)}{t}$$ (10.238)

としなければならない。この右辺を梁の回転慣性 と呼んでいる。これを式(10.168)のせん断力に ついての運動方程式に代入すれば，回転慣性を考慮した 梁の運動方程式は，一様断面・一様材料の場合には

$$m \D[2]{w(x,t)}{t} +EI \D[4]{w(x,t)}{x} -\rho I \D[4][2][\xi]{w(x,t)}{t} =p(x,t)$$ (10.239)

となる。また境界条件式(10.197)も，せん断力の式を

$$n_i \left(-EI\D[3]{w}{x}(x_i,t) -\eta I \D[4][3][x]{w}{t}(x_i,t) +\rho I \D[3][1][x]{w}{t}(x_i,t)\right)=S_i(t)$$ (10.240)

で置き換えなければならない。

この回転慣性の定量的な影響は次の節で示す。 ここではその特性を明らかにするために，たわみと場所の独立変数を 梁の長さ$\ell$で無次元化しておく。つまり

$$\xi\equiv\dfrac{x}{\ell}, \quad v(\xi,t)\equiv\dfrac{w(x,t)}{\ell}$$ (10.241)

と定義して，上の運動方程式(10.240)を書き直すと

$$\left[ \D[2]{v(\xi,t)}{t} +\dfrac{EI}{m \ell^4} \D[4]{v(\x... ...\right] -\dfrac{1}{\lambda^2} \D[4][2][\xi]{v(\xi,t)}{t} =0$$ (10.242)

となる。ここに$\lambda$は式(6.37a)で定義した細長比

(10.243)

である。$r$は断面の回転半径 と呼ばれ，断面寸法を代表する長さである。 運動方程式の左辺の大括弧内にある最初の3項が 普通の梁理論の運動方程式の項であることを念頭に置くと， 第4項の回転慣性項は細長比$\lambda$が大きくなると 比較的小さくなって無視できることがわかる。 つまり，回転慣性の影響は，細長い梁では無視してもいいことになる。 あるいは，実質的に短い梁の振動現象に相当する高周波振動を 対象とする場合には，回転慣性の影響を無視できなくなる可能性が あることを示唆している。

10.4.8.2 せん断変形の影響

図 10.62: Timoshenko梁の運動場

もし細長比が小さい梁で回転慣性の影響が大きくなるのであれば， もう一つ，梁が短い場合に考慮しなければならなくなるせん断変形の 影響も無視できなくなることが予想される。 せん断変形を考慮した梁理論の代表はTimoshenko梁理論であろう。 ここでは，せん断変形の影響をTimoshenko梁のモデルで考慮し， 回転慣性も考慮した場合の運動方程式を導き， その影響を定量的に示しておく。Timoshenko梁の運動場を 記述したのが図-10.62である。 図から明らかなように，断面内任意点の$x$, $z$方向変位成分は 微小変位の範囲内で，Bernoulli-Euler梁の式(4.3)の代わりに

$$u_x=u+z \vartheta, \quad u_z=w$$ (10.244)

となる。ここに$\vartheta(x)$は，軸のたわみ角ではなく断面の 回転角を表していて，このたわみ角と断面の回転角の差

$$\gamma=\vartheta-\left(-w'\right)$$ (10.245)

がBernoulli-Euler梁では無視していたせん断変形である。 定式化の詳細は付録-Dに示した。図-10.62の うち，曲げについての成分のみをここでは考えることにすると， 回転慣性も含めた運動方程式は

$$\D{V}{x}+p=m \D[2]{w}{t}, \quad \D{M}{x}-V=\rho I\D[2]{\vartheta}{t}$$ (10.246)

となる。またせん断力と曲げモーメントは

$$V= Gk\subsc{t}A \gamma, \quad M= EI \left(\D{\gamma}{x}-\D[2]{w}{x}\right) = EI \D{\vartheta}{x}$$ (10.247)

という関係になっている。$G$はせん断弾性係数であり。$k\subsc{t}$は， 実際には図-4.68のように断面内で 非一様に分布するはずのせん断変形を，断面内一定の せん断変形$\gamma$で近似したことを補うための係数である。値そのものは 断面形状とPoisson比$\nu$に依存する定数[16]である。例えば

$$k\subsc{t}(\mbox{円形断面})=\dfrac{6(1+\nu)}{7+6\nu}, \qquad k\subsc{t}(\mbox{矩形断面})=\dfrac{10(1+\nu)}{12+11\nu}$$ (10.248)

といった値を持つ。境界条件は

$$w=w_i \quad\mbox{あるいは}\quad n_i V=S_i, \qquad \vartheta=\vartheta_i \quad\mbox{あるいは}\quad n_i M=C_i$$ (10.249)

となる。ここに$n_i$はp.に再掲した 式(4.24)で 定義した記号であり，$u_i$等は端部で与える変位量である。 回転角についての境界条件については，$\gamma$そのものや$-w'$を 与えるようなことはできないし， この2者が中間載荷点や中間支点で連続になる必要も無いことに 注意する。

式(10.247) (10.248)から$\vartheta$, $\gamma$を 消去すると，$w$を無次元化した式(10.242)の$v$で表した 運動方程式が

$$\left[ -\dfrac{EI}{m\ell^4}\D[4]{v}{\xi} - \D[2]{v}{t} +\df... ...m\ell^4}{EI}\dfrac{1}{\lambda^2}\D[2]{p}{t} \right\} \right]=0$$ (10.250)

となる。ここでは，$m=\rho A$と， 式(10.242)の無次元化された独立変数$\xi$を用いている。 また$\lambda$は式(10.244a)で定義した細長比である。 最初の大括弧がBernoulli-Euler梁の部分，最後の大括弧がTimoshenko梁の 部分である。また二つの大括弧の間の項は，前節のBernoulli-Euler梁に 対する回転慣性の項である。また $\alpha\subsc{t}$は， 式(4.85a)で定義された係数で

$$\alpha\subsc{t}\equiv\dfrac{E}{G \left(\lambda\subsc{t}\rig... ...a\subsc{t}\equiv\dfrac{\ell}{\sqrt{\slfrac{I}{(k\subsc{t}A)}}}$$ (10.251)

と定義されている。 $\lambda\subsc{t}$も細長比 に相当するから， 細長い梁では $\alpha\subsc{t}$も比較的小さくなるので， せん断変形の影響も回転慣性と同様に無視できることがわかる。

固有振動数を求めるために，分布外力が無いものとして，たわみを

$$v(\xi,t)=V(\xi) \exp(\mbox{i}\omega t)$$

と置いて，上の運動方程式(10.251)に代入して整理すると

$$V''''+\dfrac{\omega^2 \mu^4 \left(1+\Xi\right)}{\lambda^2}... ...i \mu^8 \omega^4-\mu^4 \lambda^4 \omega^2}{\lambda^4} V=0$$

となる。ここに，プライムは$\xi$に関する微分を表し

$$\Xi\equiv \dfrac{E}{k\subsc{t}G}, \quad \mu^2\equiv \ell^2... ...{m}{EI}} \quad\to\quad \alpha\subsc{t}=\dfrac{\Xi}{\lambda^2}$$ (10.252)

と定義した。さらに$V=\exp(p\xi)$と置くことによって，特性方程式が

$$p^4+\dfrac{\omega^2 \mu^4 \left(1+\Xi\right)}{\lambda^2} ... ...i \mu^8 \omega^4-\mu^4 \lambda^4 \omega^2}{\lambda^4} V=0$$

となるので，特性根が

$$p^2=\dfrac12 \left[ -\dfrac{\omega^2\mu^4\left(1+\Xi\right)... ...\omega^4\mu^4\left(1-\Xi\right)^2+4\omega^2\lambda^4 } \right]$$

と求められる。この根は $\omega^2<\dfrac{\lambda^4}{\Xi\mu^4}$のとき， 一つは実根$\pm p_2$でもう一つが虚根 $\pm\mbox{i}p_1$になる。 虚根は上式の平方根の前の符号が負の方

$$p_1^2=\dfrac12 \left[ \dfrac{\omega^2\mu^4\left(1+\Xi\right... ...^4\left(1-\Xi\right)^2+4\omega^2\lambda^4 } \right] \eqno{(a)}$$

である。したがって，一般解が

$$V(\xi)=c_1 \sin p_1\xi+c_2 \cos p_1\xi +c_3 \sinh p_2\xi+c_4 \cosh p_2\xi \eqno{(b)}$$

と求められる。

具体的な振動数を求めるために，簡単な単純支持梁を対象とする。 境界条件は式(10.250)から，両端で

$$w=0, \quad \D{\vartheta}{x}=0= -\D[2]{w}{x}+\dfrac{m}{Gk\subsc{t}A}\D[2]{w}{t}$$

となるので，式($b$)の一般解を代入して整理すると

$$\left(\begin{array}{cccc} 1 & 0 & 1 & 0 \\ p_2^2 & 0 & -p... ... \left\{\begin{array}{c} 0 0 0 0 \end{array}\right\}$$

となる。意味のある解になるために，係数行列の行列式が零になる条件は

$$-\left(p_1^2+p_2^2\right)^2 \sin p_1 \sinh p_2 =0$$

となる。したがって，固有振動数は

$$p_1=n\pi$$

を満足する解であり，$c_1$だけが不定になるので， 固有振動モードは$\sin p_1\xi$になる。 この式を式($a$)に代入して整理すると

$$\Xi \Omega_n^4 -\beta_n \left(\beta_n+\underline{1}+\Xi\right) \Omega_n^2+\beta_n^2=0$$

となる。ここに

$$\Omega_n\equiv \dfrac{\mu^2\omega}{\left(n\pi\right)^2}, \quad \beta_n\equiv \left(\dfrac{\lambda}{n\pi}\right)^2$$ (10.253)

と定義した。下線のある項が回転慣性の影響の項であり，$\Xi$の付いた項が せん断変形の影響の項である。つまり，$\Xi=0$として， 下線の項を無視することによって，回転慣性を無視したBernoulli-Euler梁の 固有振動数が$\Omega_n^2=1$になる。 具体例として矩形断面の梁を対象とし， 式(3.44)と式(10.249)から

$$\dfrac{E}{G}=2\left(1+\nu\right), \quad k\subsc{t}=\dfrac{10(1+\nu)}{12+11\nu}$$

と置いて，ポアソン比を$\nu=0.3$としたときの，1次から4次の 固有振動数を図-10.63に示した。 細長比が小さくなって短い梁になればなるほど，Bernoulli-Euler梁の 固有振動数より小さくなる。また高次の振動数ほど，小さくなる 影響は大きくなっている。 橋梁の場合は，細長比は小さくても30くらいではないかと思われるので，2次の 固有振動数までなら10%くらいの差しか無いが，太いコンクリート柱のような 場合には，こういった影響を考えないとまずいのではないかと思うが， どうだろう。

図 10.63: 固有振動数に対する回転慣性とせん断変形の影響

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