next up previous contents index
Next: 10.3 連続体の振動 曲げ剛性の無い構造要素の振動 Up: 10. 振動論の基礎 Previous: 10.1 1自由度系の振動

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


10.2 多自由度系の振動

10.2.1 2自由度系の振動

10.2.1.1 運動方程式

図 10.26: 2自由度系の振動

1自由度系の節では,振動問題の最も重要な固有振動数を誘導し, 設計に少しだけ関連付けながら,その特徴を示した。 もちろん,社会基盤構造を1質点系で近似するのがかなり乱暴であることは 容易に想像できる。読者の多くは多分,既に有限要素法を知っているから, 任意の初期値境界値問題が多自由度系の問題で近似できることも 知っているだろう。ここでは,その基礎として,複数の 質点から成る多自由度系を対象とし,もう一つの 重要な「振動モード」という概念について説明する。 そのためには,まず運動方程式を誘導する必要があるが, この方程式を求めるという訓練をする場は教育現場にはほとんど無いのが 現実ではないだろうか。文献[136]では,国語という科目が 義務教育以来あって多くの文章を書かされて採点されるにもかかわらず, その書き方を指導する場が一切存在しないことを問題視している。 運動方程式や支配方程式を誘導する技術についても,この国語における 文書作成技術の教育とほぼ同じ状況にある。 振動の運動方程式の誘導は,非常に簡単なNewtonの方程式の誘導に他ならないので, できるだけ答を見ないで各自努力してみて欲しい。 これは,東北大学の故北原道弘教授が振動の授業で学生に期待した事項でもある。 まずは2自由度系から始めよう。

図-10.26は,二つの質点がバネとダッシュポットを はさんで直列に並んでいる系で,二つの質点が連成(連立)しながら,しかし 自由に動くことができるので,2自由度系と 呼ばれる。同じ図の右に描いた力のダイアグラムから, 質点1に作用している右方向(変位$u_1$方向)の 力の総和は

\begin{displaymath}
f_1-k_1u_1-c_1\dot{u}_1+k_2(u_2-u_1)+c_2(\dot{u}_2-\dot{u}_1)
\end{displaymath}

である。一方,質点2の場合は

\begin{displaymath}
f_2-k_2(u_2-u_1)-c_2(\dot{u}_2-\dot{u}_1)
\end{displaymath}

である。したがって,2質点のNewtonの法則は

\begin{displaymath}
f_1-k_1u_1-c_1\dot{u}_1+k_2(u_2-u_1)+c_2(\dot{u}_2-\dot{u}_1...
... \quad
f_2-k_2(u_2-u_1)-c_2(\dot{u}_2-\dot{u}_1)=m_2\ddot{u}_2
\end{displaymath}

となる。1自由度系の運動方程式と同じ順番に整理すると

\begin{eqnarray*}
&& m_1 \ddot{u}_1(t)+\left(c_1+c_2\right) \dot{u}_1(t)-c_2 ...
...,\dot{u}_1(t)+c_2 \dot{u}_2(t)
-k_2 u_1(t)+k_2 u_2(t)=f_2(t)
\end{eqnarray*}

が2自由度系の運動方程式 である。これを1自由度系の運動方程式と同じ形式で行列表示すると

\begin{displaymath}
\fat{M} \ddot{\fat{u}}(t)+\fat{C} \dot{\fat{u}}(t)
+\fat{K} \fat{u}(t)=\fat{f}(t)
\end{displaymath} (10.73)

となる。ここに$\fat{u}$を変位ベクトル,$\fat{f}$を 外力ベクトル,$\fat{M}$質量行列$\fat{C}$減衰行列$\fat{K}$剛性行列 と呼び

    $\displaystyle \fat{u}(t)\equiv\lfloor u_1(t)    u_2(t) \rfloor\supersc{t}, \quad
\fat{f}(t)\equiv\lfloor f_1(t)    f_2(t) \rfloor\supersc{t},$  
    $\displaystyle \mbox{}\qquad \fat{M}\equiv\left(\begin{array}{cc}
m_1 & 0   0 ...
...}\equiv\left(\begin{array}{cc}
k_1+k_2 & -k_2   -k_2 & k_2
\end{array}\right)$ (10.74)

のように定義した。 なお簡単のためにこの章では,括弧無しの太字で行列を表している。 式(10.74)は,時間に関する二つの未知関数の2階の微分方程式なので, 初期条件が四つ必要になる。それは二つの質点の, それぞれの初期位置と初速の合計四つでいいから

\begin{displaymath}
\fat{u}(0)=\fat{u}_0, \quad \dot{\fat{u}}(0)=\fat{v}_0
\end{displaymath} (10.75)

と表すことができる。$\fat{u}_0$, $\fat{v}_0$は式(10.75)の 変位ベクトルに対応した初期位置ベクトルと初速ベクトルである。

10.2.1.2 非減衰自由振動

最も基本的な振動特性を調べるためには,最初は非減衰の自由振動を 考える必要があるので,$\fat{C}$$\fat{f}$を無視した運動方程式

\begin{displaymath}
\fat{M} \ddot{\fat{u}}(t)
+\fat{K} \fat{u}(t)=\fat{0}
\end{displaymath} (10.76)

を対象にしよう。1自由度系の振動解と同じように

\begin{displaymath}
\fat{u}(t)=\fat{A} \exp\left(\mbox{i}\omega t\right), \quad
\fat{A}\equiv\lfloor A_1    A_2 \rfloor\supersc{t}
\end{displaymath} (10.77)

と置こう。$A_1$, $A_2$がそれぞれの質点の 振動振幅で,$\fat{A}$を振幅ベクトル(この式では複素数と考えておく)と 呼ぶことにする。 式(10.78)を式(10.77)に代入すると

\begin{displaymath}
\left(\fat{K}-\omega^2 \fat{M}\right) 
\fat{A} \exp\left(\mbox{i}\omega t\right)=\fat{0}
\end{displaymath}

が任意の時刻$t$に成立するためには

\begin{displaymath}
\left(\fat{K}-\omega^2 \fat{M}\right) \fat{A}=\fat{0}
\end{displaymath} (10.78)

でなければならない。これは$A_1$, $A_2$に対する連立方程式であるが, 右辺が零であるから,もし係数の行列が正則であれば振幅ベクトルは零になる。 それは単なる静止状態を意味し,もちろんそれも答の一つだが, 自由振動解析としては欲しい答ではない。 したがって,この2元連立方程式の2式は1次従属の関係になければならないことから

\begin{displaymath}
\det\left\vert\fat{K}-\omega^2 \fat{M}\right\vert=0
\end{displaymath} (10.79)

が成立しなければならない。これが固有振動数$\omega$を決定する 振動数方程式である。つまり

\begin{displaymath}
\det\left\vert\begin{array}{cc}
\left(k_1+k_2\right)-m_1 \...
... & -k_2 \\
-k_2 & k_2-m_2 \omega^2
\end{array}\right\vert=0
\end{displaymath}

から

\begin{displaymath}
\left\{\left(k_1+k_2\right)-m_1 \omega^2\right\}
\left(k_2...
...eft\{m_1k_2+m_2\left(k_1+k_2\right)\right\} \omega^2+k_1k_2=0
\end{displaymath}

あるいは

\begin{displaymath}
\omega^4-\left\{\left(\dfrac{k_2}{m_2}\right)+\left(\dfrac{k...
...eft(\dfrac{k_1}{m_1}\right) \left(\dfrac{k_2}{m_2}\right)=0
\end{displaymath} (10.80)

が,固有振動数$\omega$を決定する。$\omega>0$の解しか意味がないから, この方程式からは二つの解が求められる。自由度が2になったため,2種類の 固有振動数が存在することを意味する。 そしてその解は

\begin{displaymath}
\omega^2=\dfrac12 \left[
\left\{\left(\dfrac{k_2}{m_2}\righ...
...rac{k_1}{m_1}\right) \left(\dfrac{k_2}{m_2}\right)}
\right]>0
\end{displaymath}

あるいは

\begin{displaymath}
\omega^2=\dfrac12 \left[
\left\{\omega_{02}^2+\omega_{01}^2...
...02}^2\right\}^2
-4 \omega_{01}^2 \omega_{02}^2}
\right]>0
\end{displaymath} (10.81)

となる。ここに

\begin{displaymath}
m_0\equiv\dfrac{m_2}{m_1}, \quad
\omega_{01}\equiv\sqrt{\dfrac{k_1}{m_1}}, \quad
\omega_{02}\equiv\sqrt{\dfrac{k_2}{m_2}}
\end{displaymath} (10.82)

と定義した。$\omega^2$の二つの解ともに正である。この二つの解を 小さい方から$\omega_1$, $\omega_2$ ( $0<\omega_1<\omega_2$)と記す。 よくある勘違いに, $\omega_1=\omega_{01}$, $\omega_2=\omega_{02}$と 思い込む人がいるが,そうではないので注意する。

それぞれの固有振動数に対して,式(10.79)を 満足する$A_1$$A_2$の比を求めることができる。 この方程式が1次従属なので,絶対値ではなく,比しか求めることはできない。 例えば2行目の式を用いて$A_1=0$ではないことを期待(いつも そうなるとは限らないので注意。特に定期試験では。)すると,その比が

\begin{displaymath}
d\equiv\dfrac{A_2}{A_1}=
\dfrac{k_2}{k_2-m_2 \omega^2}
\end{displaymath} (10.83)

となるので,それぞれの固有振動数に対して二つの比は

\begin{displaymath}
d_1=\dfrac{k_2}{k_2-m_2 \omega_1^2}, \quad
d_2=\dfrac{k_2}{k_2-m_2 \omega_2^2}
\end{displaymath}

と表すことができる。したがって,二つの質点の変位は

\begin{eqnarray*}
u_1(t)&=&D_1 \sin\omega_1 t+D_2 \cos\omega_1 t
+D_3 \sin\o...
...ght)
+d_2 \left(D_3 \sin\omega_2 t+D_4 \cos\omega_2 t\right)
\end{eqnarray*}

と求められる。積分定数$D_1$$D_4$は 初期条件式(10.76)で決める。つまり

$\displaystyle \fat{u}(t)$ $\textstyle =$ $\displaystyle \fat{\phi}_1 \left(D_1 \sin\omega_1 t+D_2 \cos\omega_1 t\right...
...left(D_3 \sin\omega_2 t+D_4 \cos\omega_2 t\right)
\qquad\qquad\mbox{あるいは}$  
$\displaystyle \fat{u}(t)$ $\textstyle =$ $\displaystyle a_1 \fat{\phi}_1 \sin\left(\omega_1 t-\alpha_1\right)
+a_2 \fa...
...ght\}, \quad
\fat{\phi}_2=\left\{\begin{array}{c} 1   d_2 \end{array}\right\}$ (10.84)

等と表すこともできる。$a_1$, $a_2$, $\alpha_1$, $\alpha_2$が積分定数で, 初期条件式(10.76)で決めればいい。この$\fat{\phi}_n$ ($n=1$, 2)は 二つの質点の各固有振動数に対応した振動成分の,任意時刻の位置関係を 表しているので,固有振動モード と呼ばれている。小さい方の振動数に対応したモードを1次モードと呼び, 大きい方が2次モードと呼ばれる。

図 10.27: 2自由度系の振動例
\begin{figure}\begin{center}
\unitlength=.25mm
\begin{picture}(156,91)(144,-5)...
...1,Legend(Title)
%,-1,Graphics End
%E,0,
%
\end{picture}\end{center}
\end{figure}

では一例として$k=k_1=k_2$, $m=m_1=m_2$の場合を解いてみよう。 ただし,初期条件式(10.76)は

\begin{displaymath}
\fat{u}(0)=\fat{u}_0=\lfloor u_0    0 \rfloor\supersc{t}, ...
...}}(0)=\fat{v}_0=\lfloor 0    0 \rfloor\supersc{t}
\eqno{(*)}
\end{displaymath}

としておこう。式(10.80)の特性方程式は

\begin{displaymath}
\omega^4-3 \omega_0^2 \omega^2+\omega_0^4=0, \quad
\omega_0\equiv \sqrt{\dfrac{k}{m}}
\end{displaymath}

となるから,固有振動数はそれぞれ

\begin{displaymath}
\omega_1=\sqrt{\dfrac{3-\sqrt{5}}{2}} \omega_0
=\dfrac{\sq...
...ega_0
=\dfrac{\sqrt{5}+1}{2} \omega_0 \simeq 1.618 \omega_0
\end{displaymath}

と求められる。それぞれのモード係数$d_n$は,式(10.84)より

\begin{displaymath}
d_1=\dfrac{2}{\sqrt{5}-1}\simeq 1.618, \quad
d_2=\dfrac{-2}{\sqrt{5}+1}\simeq -0.6180
\end{displaymath}

と求められる。$d_1$が正なのに対し,$d_2$は 負であることに注意しよう。図-10.27の左上に, その二つの質点の位置関係,つまり振動モードを簡単に表示した。$\omega_1$に 対応したモードは,二つの質点が同じ方向に常に移動していることから, 同位相モードと呼ばれる。一方$\omega_2$に対応したモードは逆位相モードと 呼ばれる。簡単のために$D_1$$D_4$を用いた解の表現を用い,それを 初期条件式($*$)に代入すると

\begin{displaymath}
D_2+D_4=u_0, \quad d_1 D_2+d_2 D_4=0, \quad
\omega_1 D_1+\omega_2 D_3=0, \quad
d_1 \omega_1 D_1+d_2 \omega_2 D_3=0
\end{displaymath}

となるから

\begin{displaymath}
D_1=0, \quad D_2=\dfrac{d_2 u_0}{d_2-d_1}, \quad
D_3=0, \quad D_4=\dfrac{-d_1 u_0}{d_2-d_1}
\end{displaymath}

と求められる。したがって,解が

\begin{displaymath}
\dfrac{u_1(t)}{u_0}=\dfrac{d_2}{d_2-d_1} \cos \omega_1 t
-...
...1 d_2}{d_2-d_1} 
\left(\cos\omega_1 t-\cos\omega_2 t\right)
\end{displaymath}

となる。図-10.27の左下の図が,この2質点の 変位の様子である。2種類の振動数で揺れているので,1自由度系より 複雑な動きに見えるかもしれない。そこで,それぞれの質点の変位を, 二つの振動数毎に分けて描いたのが右の二つの図である。 まず小さい振動数の$\omega_1$の図(右上)を見ると,二つの 質点が同じ符号で $\slfrac{1}{d_1}$の比で運動しているのが明らかである。 これが同位相モードである。一方,大きい振動数の$\omega_2$の図(右下)を 見ると,お互いがそっぽを向いて $\slfrac{1}{d_2}$の比で運動している。 これを逆位相モードと呼ぶのである。この例では,左の図から,あるいは 右2図の縦軸を比較すれば明らかなように, 初期条件の影響を強く受けたため逆位相の2次モードが卓越している。

10.2.1.3 減衰自由振動

では,減衰を含む自由振動の場合を考えてみよう。 式(10.74)で外力を無視すれば,運動方程式は

\begin{displaymath}
\fat{M} \ddot{\fat{u}}(t)+\fat{C} \dot{\fat{u}}(t)
+\fat{K} \fat{u}(t)=\fat{0}
\end{displaymath} (10.85)

になる。非減衰の場合と同じように

\begin{displaymath}
\fat{u}(t)=\fat{A} \exp\left(\mbox{i}\omega t\right), \quad
\fat{A}\equiv\lfloor A_1    A_2 \rfloor\supersc{t}
\end{displaymath} (10.86)

と置く。この場合は振幅ベクトル$\fat{A}$は複素数である。 というのも,式(10.87)を式(10.86)に代入すると

\begin{displaymath}
\left(\fat{K}-\omega^2 \fat{M}+\mbox{i} \omega \fat{C}\right) 
\fat{A} \exp\left(\mbox{i}\omega t\right)=\fat{0}
\end{displaymath}

となるからである。これが任意の時刻$t$に成立するためには

\begin{displaymath}
\left(\fat{K}-\omega^2 \fat{M}+\mbox{i} \omega \fat{C}\right) 
\fat{A}=\fat{0}
\end{displaymath} (10.87)

でなければならない。よって振動数方程式は

\begin{displaymath}
\det\left\vert\fat{K}-\omega^2 \fat{M}+\mbox{i} \omega \fat{C}\right\vert=0
\end{displaymath} (10.88)

となる。1自由度系の結果からも当然なように, この方程式の解は複素数である。 ただこの方程式は線形代数で習う標準的な固有値問題にはなっていないことから, 自由度が大きくなると固有振動数を求めることは困難になる。 つまり2自由度系の式(10.89)は$\omega$の複素係数4次方程式で あることから,電卓と手計算で解を求めることができる限界だろう。 したがって,このような減衰振動については別の解法を用いることにするが, それについては多自由度系の節-10.2.2で説明する。

10.2.1.4 非減衰強制振動--多自由度系の共振曲線

減衰系については,多自由度系の節-10.2.2で 取り扱うことにし,ここでは もう一つの基本的な状況として,正弦波状の外力が作用したときの強制振動を 解いて,2自由度系の共振曲線を 求めておく。図-10.26には,質点に 直接外力が作用しているものを描いたが,ここの例としては, 支点が正弦波で強制変位させられる図-10.28の 問題の方を解いておこう。図に示した力のダイアグラムから, 二つの質点に作用する右向きの力は,それぞれ

\begin{displaymath}
k_2 \left(u_2-u_1\right)-k_1 \left(u_1-w\right), \quad
-k_2 \left(u_2-u_1\right)
\end{displaymath}

なので,運動方程式は

\begin{displaymath}
m_1 \ddot{u}_1+k_1 \left(u_1-w\right)-k_2 \left(u_2-u_1\right)=0, \quad
m_2 \ddot{u}_2+k_2 \left(u_2-u_1\right)=0
\end{displaymath}

となる。1自由度系の場合と同様に,相対変位を

\begin{displaymath}
v_1(t)\equiv u_1(t)-w(t), \quad v_2(t)\equiv u_2(t)-w(t)
\end{displaymath} (10.89)

と定義して,上の運動方程式を書き直すと

\begin{displaymath}
m_1 \ddot{v}_1+\left(k_1+k_2\right) v_1-k_2 v_2=-m_1 \ddot{w}, \quad
m_2 \ddot{v}_2-k_2 v_1+k_2 v_2=-m_2 \ddot{w}
\end{displaymath}

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

\begin{displaymath}
\fat{M} \ddot{\fat{v}}(t)+\fat{K} \fat{v}(t)=\fat{f}_w(t),...
...\equiv -\ddot{w}(t) \lfloor m_1    m_2 \rfloor\supersc{t}
\end{displaymath} (10.90)

となる。$\fat{M}$$\fat{K}$は,前述の定義式(10.75)と 同じである。

図 10.28: 非減衰強制振動

ここでは支点の変位がsine関数で与えられ

\begin{displaymath}
w(t)=w_0 \sin pt
\end{displaymath} (10.91)

で動くものとすると,強制振動解も$\sin pt$で振動するだろうから

\begin{displaymath}
\fat{v}=\fat{A} \sin pt
\end{displaymath} (10.92)

と置く。式(10.92) (10.93)を式(10.91)に 代入すると

\begin{displaymath}
\left(\fat{K}-p^2 \fat{M}\right) \fat{A}=p^2 w_0 
\left\{\begin{array}{c}
m_1  m_2
\end{array}\right\}
\end{displaymath}

となる。あるいは

\begin{eqnarray*}
&& \left(\begin{array}{cc}
\dfrac{k_1+k_2}{m_1}-p^2 & -\dfrac...
...ht\}=p^2 w_0
\left\{\begin{array}{c}
1  1
\end{array}\right\}
\end{eqnarray*}

となる。ここでは式(10.83)の記号を用いた。 これを解くと

\begin{displaymath}
\left\{\begin{array}{c}
A_1  A_2
\end{array}\right\}=\dfr...
..._{01}^2+\omega_{02}^2+m_0\omega_{02}^2-p^2
\end{array}\right\}
\end{displaymath}

となる。ここに分母は

\begin{displaymath}
\Delta\left(p^2\right)\equiv
p^4-\left(\omega_{02}^2+\omega_{01}^2+m_0\omega_{02}^2\right) p^2
+\omega_{01}^2\omega_{02}^2
\end{displaymath}

と置いたが,この式は,式(10.81)の振動数方程式(特性方程式)の 左辺と全く同じなので,固有振動数を$\omega_1$, $\omega_2$とすると, これは

\begin{displaymath}
\Delta\left(p^2\right)=\left(p^2-\omega_1^2\right) \left(p^2-\omega_2^2\right)
\end{displaymath}

と置くことができる。したがって,振幅は

\begin{displaymath}
\left\{\begin{array}{c}
\dfrac{A_1}{w_0}  \noalign{\vskip...
...01}^2+\omega_{02}^2+m_0\omega_{02}^2-p^2
\end{array}\right\}
\end{displaymath} (10.93)

となる。分母から明らかなように,外力の振動数$p$が二つの固有振動数の いずれかに一致したときに振幅が無限大になる。つまり共振点が 二つ存在することになる。

図 10.29: 2自由度系の共振曲線
\begin{figure}\begin{center}
\unitlength=.01mm
\begin{picture}(6979,6500)(1250,-...
...1,Legend(Title)
%,-1,Graphics End
%E,0,
%
\end{picture}\end{center}
\end{figure}

前の例と同じように$m=m_1=m_2$, $k=k_1=k_2$の場合の 共振曲線を求めてみる。 固有振動数は, $\omega_0\equiv\sqrt{\slfrac{k}{m}}$と定義すると

\begin{displaymath}
\omega_1=\sqrt{\dfrac{3-\sqrt{5}}{2}} \omega_0\simeq 0.6180...
...2=\sqrt{\dfrac{3+\sqrt{5}}{2}} \omega_0\simeq 1.618 \omega_0
\end{displaymath}

であった。それを用いて,式(10.94)は

\begin{displaymath}
\left\{\begin{array}{c}
\dfrac{A_1}{w_0}  \noalign{\vskip...
...2ex}
3-\left(\dfrac{p}{\omega_0}\right)^2
\end{array}\right\}
\end{displaymath}

となる。これを図-10.29に示した。$p$が 二つの固有振動数に一致したときに,振幅は無限大になる。 また,上式右辺を見ればわかるように, 質点1は $\slfrac{p}{\omega_0}=\sqrt{2}\simeq 1.41$で, 質点2は $\slfrac{p}{\omega_0}=\sqrt{3}\simeq 1.73$で 動かなくなる(相対変位なので,支点と一緒に動く)ことも面白い。

  1. 図-10.30 (p.[*])に示したのは, 四角い剛体の重心Gの水平運動$u(t)$と回転運動$\theta(t)$の2自由度系と 考えていい。こういう運動をロッキング運動と呼んでいる。 剛体の質量を$M$とし,重心G回りの慣性モーメントを$J$とする。 剛体の下端中央には,回転角に比例してモーメント抵抗を発生させる 回転バネが付いており,その比例係数を$k\subsc{r}$とする。また, 剛体下端右端には,バネ定数$k\subsc{h}$の水平の線形バネが付いている。 これが自由振動をしているものとし, $u=U\exp(\mbox{i}\omega t)$, $\theta=\Theta\exp(\mbox{i}\omega t)$としたとき, 見かけ上の回転の中心Oまでの距離$x$は, $x=\slfrac{U}{\Theta}$と 考えていい。これがモードに相当する。重力$g$は無視し

    \begin{displaymath}
\omega\subsc{h}\equiv\sqrt{\dfrac{k\subsc{h}}{M}}, \quad
\om...
...mega\subsc{h}}\right)^2=4, \quad
\left(\dfrac{h}{r}\right)^2=3
\end{displaymath}

    という値をもつときの,固有振動数と固有モードを求めよ。

  2. 例えば図-10.26の質点2に正弦波外力 $f_2(t)=f_0\sin pt$が 作用したときの共振曲線を求めてみよ。
  3. 例えば$m_1$$k_1$の1自由度系を設計して,それができあがって しまったとする。しかし,使い始めてみたところ,その場所では 固有振動数 $\omega_0=\sqrt{\slfrac{k_1}{m_1}}$付近の共振に近い振動が 頻繁に発生することが判明した。こういう設計をしてはいけない。 そこで費用のことを考え,作り直す代わりに小さい質点と弱いバネをつけて 共振点を少しずらすことを考えた。例えば $\slfrac{m_2}{m_1}=0.1$, $\slfrac{k_2}{k_1}=0.1$の系を直列につないで2自由度系にした。 この2自由度系の固有振動数と,支点を式(10.92)で 揺らしたときの質点1の共振曲線を求め,その改善度合いを観察せよ。 これは,一種のパッシブダンパー による振動制御に相当する。

図 10.30: ロッキング運動
\begin{figure}\begin{center}
\unitlength=.25mm
\begin{picture}(160,220.276)(144,...
...48 (line)
\path (208,62.276)(208,2.276)
%
\end{picture}\end{center}
\end{figure}


10.2.2 多自由度系の振動

10.2.2.1 運動方程式

図 10.31: 多自由度系の振動

次に,図-10.31に示したような, もっと一般的に自由度が$N$の系を考えてみよう。 質点1と質点$N$の力のダイアグラムは2自由度系の2質点と 同じになるので,途中にある$n$番目の質点の力のダイアグラムだけを 同じ図の下に示した。これから$n$番目の質点の運動方程式

\begin{displaymath}
m_n \ddot{u}_n
-c_{n} \dot{u}_{n-1}+\left(c_n+c_{n+1}\righ...
... u_{n-1}+\left(k_n+k_{n+1}\right) u_n
-k_{n+1} u_{n+1}=f_n
\end{displaymath}

となることはわかると思う。したがって,行列表示すると

\begin{displaymath}
\fat{M} \ddot{\fat{u}}(t)+\fat{C} \dot{\fat{u}}(t)
+\fat{K} \fat{u}(t)=\fat{f}(t)
\end{displaymath} (10.94)

となる。ここに, 変位ベクトルや外力ベクトル・質量行列・減衰行列・剛性行列は

    $\displaystyle \fat{u}(t)\equiv\left\{\begin{array}{c}
u_1(t)   u_2(t)   \vd...
... m_{N-1} & 0 \\
\multicolumn{2}{l}{\mbox{Symm.}} & & & m_N
\end{array}\right),$ (10.95)
    $\displaystyle \fat{C}\equiv\left(\begin{array}{ccccc}
c_1+c_2 & -c_2 & 0 & \cdo...
...1}+k_N & -k_N \\
\multicolumn{2}{l}{\mbox{Symm.}} & & & k_N
\end{array}\right)$  

のように定義すればいい。質量行列 は対角行列になっていて,このような質量行列を集中質量行列 と呼ぶことがある。減衰行列と剛性行列は対称行列で, かつバンド幅3のバンド行列 (帯行列)でもあることは覚えておくこと。 初期条件については,これも2自由度系と同様

\begin{displaymath}
\fat{u}(0)=\fat{u}_0, \quad \dot{\fat{u}}(0)=\fat{v}_0
\end{displaymath} (10.96)

で与えればいい。もちろん$\fat{u}_0$$\fat{v}_0$$N$行のベクトルである。

10.2.2.2 非減衰自由振動

最初は,$\fat{f}$$\fat{C}$を無視した,非減衰の自由振動を考えておこう。 変位ベクトルを

\begin{displaymath}
\fat{u}=\fat{\phi} \exp\left(\mbox{i}\omega t\right)
\end{displaymath} (10.97)

と置いて,運動方程式に代入すると,任意時刻に

\begin{displaymath}
\left(\fat{K}-\omega^2 \fat{M}\right) \fat{\phi}=\fat{0}
\end{displaymath} (10.98)

でなければならない。これが非零の$\fat{\phi}$,つまり有意な解を持つための 条件が

\begin{displaymath}
\det\left\vert\fat{K}-\omega^2 \fat{M}\right\vert=0
\end{displaymath} (10.99)

であり,これが振動数方程式(特性方程式)になる。この式から 求められる固有振動数を, 小さい順に並べて$\omega_n$ ($n=1$, 2, $\cdots$ $N$)とする。 安定な系において固有振動数が$\omega_n^2>0$の実数である ことの証明は文献[103]にある。 このそれぞれの固有振動数に対応した固有振動モード$\fat{\phi}_n$とすると,それは

\begin{displaymath}
\left(\fat{K}-\omega_n^2 \fat{M}\right) \fat{\phi}_n=\fat{0}
\end{displaymath} (10.100)

を満足する。式(10.100)が成立することから, この式(10.101)の左辺の係数行列は特異行列である。したがって, この式は1次従属な関係式だから,$\fat{\phi}_n$は その成分を唯一には定めることはできず,各要素の比だけを求める ことができる。したがって,任意の初期値問題の自由振動解(一般解)は 式(10.98)の代わりに

\begin{displaymath}
\fat{u}(t)=
\sum_{n=1}^N \left\{\overline{A}_n \exp\left(\m...
...t)
+B_n  \cos\left(\omega_n t\right)\right\} \fat{\phi}_n
\end{displaymath} (10.101)

のような,2種類の未定定数を含む級数表示の形にならざるを得ない。 この複素数の $\overline{A}_n$ $\overline{B}_n$,あるいは実数の$A_n$$B_n$を 初期条件式(10.97)で決定すれば, 唯一の自由振動解が求められることになる。

ところで,固有振動数(固有周期)の大小関係に関しては

\begin{displaymath}
\omega_1<\omega_2<\cdots<\omega_N, \quad
T_1>T_2>\cdots>T_N
\end{displaymath} (10.102)

を仮定して定義したが,この不等号には等号が無いことを覚えておいて欲しい。 というのも,もしここに等号が成り立ち二つのモードの振動数が 同じになる場合には,それは系の「動的不安定 」を意味するからである。 振動数方程式は,微分方程式の特性方程式と同じものであり, その二つの特性根が一致する状況は特性根が重根であることを意味する。 その重根に対応した一般解は, $\exp\left(\cdot\right)$の 形ではなく $\left(c_1+c_2 t\right)\exp\left(\cdot\right)$の 形になることを意味し,それは時間と共に振幅が 無限大になる解,すなわち発散解になる。 動的に不安定な社会基盤構造があってはならないので, 正しく設計された構造系の上式(10.103)の 不等号には等号があってはならないのである。

10.2.2.3 モードの直交性

式(10.95)は連立常微分方程式で,外力が存在する場合に, そのまま解くのは困難であることはすぐに推測できると思う。 これを比較的容易にする手法に必要な,固有振動モードの重要な特性を ここで明らかにしておこう。例えば$i$次と$j$次の固有振動モード($i\neq j$)は それぞれ対応した式(10.101)を満足するから

\begin{displaymath}
\fat{K} \fat{\phi}_i=\omega_i^2 \fat{M} \fat{\phi}_i, \quad
\fat{K} \fat{\phi}_j=\omega_j^2 \fat{M} \fat{\phi}_j
\end{displaymath}

が成立する。この第1式に $\fat{\phi}_j\supersc{t}$を左から乗じ, 第2式に $\fat{\phi}_i\supersc{t}$を左から乗じると

\begin{displaymath}
\fat{\phi}_j\supersc{t} \fat{K} \fat{\phi}_i=
\omega_i^2\...
...^2 \fat{\phi}_i\supersc{t} \fat{M} \fat{\phi}_j
\eqno{(*)}
\end{displaymath}

となる。左の式は第$i$次の振動モードの運動方程式が第$j$次モードとする 仮想仕事と考えればいい。 質量行列は対角行列であり,剛性行列は対称行列であることから, 左の式の両辺の項は

\begin{displaymath}
\fat{\phi}_j\supersc{t} \fat{M} \fat{\phi}_i=
\fat{\phi}_...
...,\fat{\phi}_j
=\fat{\phi}_i\supersc{t} \fat{K} \fat{\phi}_j
\end{displaymath}

となる。これを上式($*$)に代入した上で, 二つの式を辺々引き算すると,左辺は同じであることから

\begin{displaymath}
0=\left(\omega_i^2-\omega_j^2\right) \fat{\phi}_i\supersc{t} 
\fat{M} \fat{\phi}_j
\end{displaymath}

となる。式(10.103)で解説したように,安定な構造系であれば 二つの異なる固有振動数は異なる値を持つことから,$i\neq j$の ときには $\omega_i\neq\omega_j$であり,上式は

\begin{displaymath}
\fat{\phi}_i\supersc{t} \fat{M} \fat{\phi}_j=0, \quad j\neq i
\end{displaymath} (10.103)

であることを示している。もし質量行列$\fat{M}$が 単位行列$\fat{I}$であれば,この式は, 二つのベクトル$\fat{\phi}_i$$\fat{\phi}_j$の内積であり, それが零になるということは,二つのベクトルが直交していることを 意味する。ここでは,内積を拡張して定義し,質量行列$\fat{M}$という重みを 含んだ内積 を用いて,二つの固有振動モードは「直交している 」と呼ぶことにしている。

ただし,$j=i$の時は零にはならないことから

\begin{displaymath}
\fat{\phi}_i\supersc{t} \fat{M} \fat{\phi}_j=
\left\{\begi...
...l}
0, & j\neq i \\
\overline{m}_i, & j=i
\end{array}\right.
\end{displaymath} (10.104)

のように,新しい量 $\overline{m}_i$を定義する。 この $\overline{m}_i$$i$次の固有質量のようなものであるが, 歴史的には一般化された質量 と呼ばれる。また式(10.101)を式(10.105)に 代入すると

\begin{displaymath}
\fat{\phi}_i\supersc{t} \fat{K} \fat{\phi}_j=
\omega_j^2 ...
... i \\
\overline{m}_i \omega_i^2, & j=i
\end{array} \right.
\end{displaymath} (10.105)

という関係も成立するため,剛性行列$\fat{K}$という重みを 含んだ内積に対しても直交性が成立する。$\fat{K}$は 対角行列ではないので,この内積はかなり一般化されたものであることもわかる。

10.2.2.4 自由振動の級数解--モード解析法

では,任意の初期条件式(10.97)に対して,自由振動解を 求めることを,再度考え直してみよう。どうやら解は,2自由度系の 振動の図-10.27からも予想されるように,$N$個の 振動モードで揺れて,それが重ね合わされた式(10.102)のような 級数解で表すことができそうだ。そこで解が,時間の未知関数$q_n(t)$を 係数とした

\begin{displaymath}
\fat{u}(t)=\sum_{n=1}^N q_n(t) \fat{\phi}_n
\end{displaymath} (10.106)

と表現できると仮定してみよう。未知関数$q_n(t)$が満足すべき 微分方程式は運動方程式から求められるだろうから, これを運動方程式(10.95)の 非減衰自由振動バージョンの $\fat{C}\equiv 0$, $\fat{f}(t)\equiv 0$の 式に代入すると

\begin{displaymath}
\sum_{n=1}^N \ddot{q}_n(t) \fat{M} \fat{\phi}_n
+\sum_{n=1}^N q_n(t) \fat{K} \fat{\phi}_n=\fat{0}
\end{displaymath}

となる。これに左から $\fat{\phi}_j\supersc{t}$を乗じ(つまり$j$次 モードとの仮想仕事を算定し), 式(10.105) (10.106)の直交条件を利用すると, 総和のうち$j\neq n$の項はすべて零になるので

\begin{displaymath}
\sum_{n=1}^N \ddot{q}_n(t) \fat{\phi}_j\supersc{t} \fat{M}...
...rline{m}_j \ddot{q}_j(t)+\overline{m}_j \omega_j^2 q_j(t)=0
\end{displaymath}

となり,結局,$q_j(t)$

\begin{displaymath}
\ddot{q}_j(t)+\omega_j^2 q_j(t)=0
\end{displaymath} (10.107)

を満足しないといけないことがわかる。 これは形式的に,1自由度系の非減衰自由振動の運動方程式そのものである。 すなわち$N$自由度系の問題は,固有振動数とモードさえ算定できていれば (しかも$N$個全部が求められている必要が無いのもミソ),$N$個の 独立した1自由度系の振動問題に置き換えられることを意味する。 したがって解は

\begin{displaymath}
q_j(t)=A_j \sin \omega_j t+B_j \cos\omega_j t
\end{displaymath} (10.108)

になる。この積分定数$A_j$, $B_j$を初期条件で決定すればいい。

式(10.97)の初期条件式は,級数解で表すと

\begin{displaymath}
\sum_{n=1}^N q_n(0) \fat{\phi}_n=\fat{u}_0, \quad
\sum_{n=1}^N \dot{q}_n(0) \fat{\phi}_n=\fat{v}_0
\end{displaymath}

であることから,例えば左から $\fat{\phi}_j\supersc{t} \fat{M}$を 乗じると

\begin{displaymath}
\sum_{n=1}^N q_n(0) \fat{\phi}_j\supersc{t} \fat{M} \fat{...
...M} \fat{\phi}_n=
\fat{\phi}_j\supersc{t} \fat{M} \fat{v}_0
\end{displaymath}

となり,再び直交性の式(10.105) (10.106)を 利用すると,結局

\begin{displaymath}
q_j(0)=\dfrac{\fat{\phi}_j\supersc{t} \fat{M} \fat{u}_0}{\...
...\fat{\phi}_j\supersc{t} \fat{M} \fat{v}_0}{\overline{m}_j}
\end{displaymath} (10.109)

が,それぞれの$q_j(t)$に対する初期条件になる。この式に 式(10.109)の一般解を代入して,積分定数の$A_j$, $B_j$を 決定すれば,唯一の解が求められる。再度,$j$次の$q_j(t)$に 対する初期値問題を整理して書いておくと

\begin{displaymath}
% latex2html id marker 9043\ddot{q}_j(t)+\omega_j^2 q_j(t)=0,
\eqno{(\ref{eq:fvib-105}) \mbox{再掲}}
\end{displaymath}


\begin{displaymath}
% latex2html id marker 9045\mbox{ただし,初期条件が}\quad
...
...v}_0}{\overline{m}_j}
\eqno{(\ref{eq:fvib-107}) \mbox{再掲}}
\end{displaymath}

を解けばいいことになる。これは1自由度系の振動問題である。 このように,非減衰自由振動モードを用いて式(10.107)のような 級数解を仮定して初期値問題を解く方法をモード解析法 と呼んでいる。このような問題の$q_j(t)$は 第$j$次モードの振幅であるが,歴史的には一般化された座標 (物理学での名称)と呼ばれている。 連続体の力学に限定すれば,著者としては一般化された変位 と呼びたいところではあるが。

ここでは一例として,3質点の自由振動を解いておこう。ただし, 簡単のために$m=m_1=m_2=m_3$, $k=k_1=k_2=k_3$とする。 また初期条件も簡単のために

\begin{displaymath}
\fat{u}_0=\lfloor u_0    0    0 \rfloor\supersc{t}, \quad
\fat{v}_0=\lfloor 0    0    0 \rfloor\supersc{t}
\end{displaymath}

としておこう。そのとき,式(10.99)は

\begin{displaymath}
\left(\begin{array}{ccc}
2k-m\omega^2 & -k & 0 \\
-k & 2k...
...i_1  \phi_2  \phi_3
\end{array}\right\}=\fat{0}
\eqno{(*)}
\end{displaymath}

となるから,振動数方程式はこの係数行列の行列式から

\begin{displaymath}
\omega^6-5\omega_0^2\omega^4+6\omega_0^4\omega^2-\omega_0^6=0, \quad
\omega_0\equiv\sqrt{\dfrac{k}{m}}
\end{displaymath}

となり,固有振動数が

\begin{displaymath}
\xi_1\equiv \dfrac{\omega_1}{\omega_0}\simeq 0.455, \quad
\x...
... 1.25, \quad
\xi_3\equiv \dfrac{\omega_3}{\omega_0}\simeq 1.80
\end{displaymath}

と求められる。3次方程式なので,公式を用いるか,例えば2分法 等で数値的に解ける。 再度上式($*$)にこれを代入すると,モードが例えば

\begin{displaymath}
\fat{\phi}_1=\left\{\begin{array}{c}
1  1.801  2.246
\e...
...ft\{\begin{array}{c}
1  -1.247  0.555
\end{array}\right\}
\end{displaymath}

と表現10.10できる。第1次モードはすべてが同じ向きの成分を持っているので 同位相のモードであることがわかる。これに対し,第2次と第3次の モードは,質点の一つが他の二つにそっぽを向いていることから 逆位相のモードになっている。式(10.105)に代入すれば,直交性も 確かめることができる。

図 10.32: 3自由度系の非減衰自由振動
\begin{figure}\begin{center}
\unitlength=.01mm
\begin{picture}(6845,4250)(1500,-...
...1,Legend(Title)
%,-1,Graphics End
%E,0,
%
\end{picture}\end{center}
\end{figure}

さて,この各モードを用い,それぞれの一般化された質量を 式(10.105)で求めると

\begin{eqnarray*}
&&\mu_1\equiv \dfrac{\overline{m}_1}{m}=1^2+1.801^2+2.246^2=9....
...}_2}{m}=1.841, \quad
\mu_3\equiv \dfrac{\overline{m}_3}{m}=2.863
\end{eqnarray*}

となる。 これを用いて式(10.110)に代入すると,たまたま$\fat{\phi}_i$の 第1成分をすべて1にしておいたため,各モードの 一般化された座標に対する初期条件が

\begin{eqnarray*}
q_1(0)&=&\dfrac{u_0}{\mu_1}, \quad \dot{q}_1(0)=0,  % \qquad...
...}_2(0)=0, \qquad
q_3(0)=\dfrac{u_0}{\mu_3}, \quad \dot{q}_3(0)=0
\end{eqnarray*}

と表される。各モードの一般解は式(10.109)であったから, この一般解を上の初期条件式に代入すると

\begin{displaymath}
A_1=0, \quad B_1=\dfrac{u_0}{\mu_1}, \qquad
A_2=0, \quad B_2=\dfrac{u_0}{\mu_2}, \qquad
A_3=0, \quad B_3=\dfrac{u_0}{\mu_3}
\end{displaymath}

と積分定数が求められる。したがって,解は

\begin{displaymath}
\fat{u}(t)=u_0\left\{
\dfrac{1}{\mu_1} \fat{\phi}_1 \cos\l...
...\mu_3} \fat{\phi}_3 \cos\left(\xi_3\omega_0 t\right)\right\}
\end{displaymath}

になる。図-10.32に3質点の変位履歴を示した。 三つの振動数の周期運動の重ね合わせなので,ちょっと見ただけでは 周期運動には見えないのも面白い。

10.2.2.5 非減衰強制振動

次に,外力が存在する強制振動を解くことを考えよう。 運動方程式は,式(10.95)の減衰項を無視したものになるため

\begin{displaymath}
\fat{M} \ddot{\fat{u}}(t)+\fat{K} \fat{u}(t)=\fat{f}(t)
\end{displaymath} (10.110)

である。ここでもモード解析法を用いてみよう。 つまり,解を自由振動の固有振動モードを用いた級数で

\begin{displaymath}
\fat{u}(t)=\sum_{n=1}^N q_n(t) \fat{\phi}_n
\end{displaymath} (10.111)

と表現できそうだと仮定して,これを上の運動方程式に代入すると

\begin{displaymath}
\sum_{n=1}^N \ddot{q}_n(t) \fat{M} \fat{\phi}_n
+\sum_{n=1}^N q_n(t) \fat{K} \fat{\phi}_n=\fat{f}(t)
\end{displaymath}

となる。これに左から $\fat{\phi}_j\supersc{t}$を乗じ(つまり$j$次モードとの 仮想仕事を算定し), 式(10.105) (10.106)の直交条件を利用すると, 総和のうち$j\neq n$の項はすべて零になるので

\begin{displaymath}
\sum_{n=1}^N \ddot{q}_n(t) \fat{\phi}_j\supersc{t} \fat{M}...
...ddot{q}_j(t)+\overline{m}_j \omega_j^2 q_j(t)
=\bar{f}_j(t)
\end{displaymath}

となる。ここに

\begin{displaymath}
\bar{f}_j(t)\equiv \fat{\phi}_j\supersc{t} \fat{f}(t)
\end{displaymath} (10.112)

と定義したものは,外力の第$j$次モード成分と考えられ,一般化された外力 と呼ばれている。結局

\begin{displaymath}
\ddot{q}_j(t)+\omega_j^2 q_j(t)=\dfrac{1}{\overline{m}_j} \bar{f}_j(t)
\end{displaymath} (10.113)

を解けば$q_j(t)$が求められる。これは1自由度系の非減衰強制振動の 運動方程式そのものである。この運動方程式の解は, 右辺の外力$\bar{f}_j(t)$に対する特解を$r_j(t)$と記すと

\begin{displaymath}
q_j(t)=A_j \sin\omega_j t+B_j \cos\omega_j t+r_j(t)
\end{displaymath} (10.114)

となるので,これを各モードの初期条件式(10.110)に代入して

\begin{eqnarray*}
q_j(0)&=&B_j+r_j(0)
=\dfrac{\fat{\phi}_j\supersc{t} \fat{M}\...
...t{\phi}_j\supersc{t} \fat{M} \fat{u}_0}{\overline{m}_j}-r_j(0)
\end{eqnarray*}

によって,積分定数$A_j$, $B_j$を決定すればいい。 あるいは1自由度系の式(10.48)の単位衝撃応答を利用した, 式(10.49)のDuhamel積分を用いて強制振動解を求めることにすると, 一般解は

\begin{displaymath}
q_j(t)=A_j \sin\omega_j t+B_j \cos\omega_j t
+\int_0^\inf...
...au)  \sin\left\{\omega_j\left(t-\tau\right)\right\} \dint\tau
\end{displaymath}

と表すこともできる。Duhamel積分は零初期条件の解なので, $A_j=
\dfrac{\fat{\phi}_j\supersc{t} \fat{M} \fat{v}_0}{\overline{m}_j\omega_j}$, $B_j=\dfrac{\fat{\phi}_j\supersc{t} \fat{M} \fat{u}_0}{\overline{m}_j}$と なる。

図 10.33: 2自由度系の非減衰強制振動
\begin{figure}\begin{center}
\unitlength=.01mm
\begin{picture}(6420,4500)(1800,-...
...1,Legend(Title)
%,-1,Graphics End
%E,0,
%
\end{picture}\end{center}
\end{figure}

3自由度系の計算は面倒なので演習問題に回すことにし, ここでは図-10.27で解いた$m=m_1=m_2$, $k=k_1=k_2$の2自由度系の 強制振動を,モード解析法で解いておこう。簡単のために 初期条件と外力条件は

\begin{displaymath}
\fat{u}_0=\left\{\begin{array}{c} u_0  0 \end{array}\right...
...=\left\{\begin{array}{c} 0  f_0 \sin pt \end{array}\right\}
\end{displaymath}

としておく。まず固有振動数は

\begin{displaymath}
\xi_1\equiv\dfrac{\omega_1}{\omega_0}\simeq 0.6180, \quad
\x...
...omega_0}\simeq 1.618, \quad
\omega_0\equiv\sqrt{\dfrac{k}{m}}
\end{displaymath}

であり,固有振動モードは

\begin{displaymath}
\fat{\phi}_1=\left\{\begin{array}{c} 1  d_1 \end{array}\ri...
...gin{array}{c} 1  d_2 \end{array}\right\},
\quad d_1=-0.6189
\end{displaymath}

であった。次に一般化された質量と$m$の比を求めると

\begin{displaymath}
\mu_1=\dfrac{\overline{m}_1}{m}=1^2+1.618^2=3.618, \quad
\mu_2=\dfrac{\overline{m}_2}{m}=1.382
\end{displaymath}

となる。初期条件式(10.110)に以上の情報を代入すると

\begin{displaymath}
q_1(0)=\dfrac{u_0}{\mu_1}, \quad
\dot{q}_1(0)=\dfrac{v_0 d_...
...\dfrac{u_0}{\mu_2}, \quad
\dot{q}_2(0)=\dfrac{v_0 d_2}{\mu_2}
\end{displaymath}

であり,一般化された外力は式(10.113)から

\begin{displaymath}
\bar{f}_1=d_1 f_0\sin pt, \quad
\bar{f}_2=d_2 f_0\sin pt
\end{displaymath}

となる。よって,Duhamel積分を用いない場合の一般解は

\begin{displaymath}
q_j(t)=A_j \sin \omega_j t + B_j \cos\omega_j t
+\dfrac{d_j f_0}{\overline{m}_j \left(\omega_j^2-p^2\right)} 
\sin pt
\end{displaymath}

と求められる。これを上の初期条件式に代入すると

\begin{displaymath}
A_j=\dfrac{d_j}{\mu_j \xi_j} \left(\dfrac{v_0}{\omega_0}\r...
...(\slfrac{p}{\omega_0}\right)^2},
\quad
B_j=\dfrac{u_0}{\mu_j}
\end{displaymath}

のように$A_j$, $B_j$が求められる。例として

\begin{displaymath}
\dfrac{v_0}{\omega_0}=0.5 u_0, \quad
\dfrac{f_0}{k}=0.8 u_0, \quad
\dfrac{p}{\omega_0}=1
\end{displaymath}

の場合の結果を図-10.33に示した。

  1. $m=m_1=m_2=m_3$, $k=k_1=k_2=k_3$の3自由度系で

    \begin{displaymath}
\fat{u_0}=\left\{\begin{array}{c} u_0  0  0 \end{array}\...
..., \quad
\dfrac{f_0}{k}=0.8 u_0, \quad
\dfrac{p}{\omega_0}=1
\end{displaymath}

    の場合の強制振動解を求めて図示せよ。


10.2.2.6 粘性減衰自由振動

では,粘性減衰も含めた式(10.95)を解くことを考えよう。 式(10.107)の形で解を仮定したモード解析法が, 減衰行列を含めた運動方程式にも使えるだろうか。 残念ながら答はもちろん否である。というのも, モード解析法に使うモードは非減衰自由振動モードであるからだ。 つまりそこには減衰の特性は含まれていないから, そんなものを使って正しい答が出るはずはないと考えるのが素直である。 それを知りながら,敢えて非減衰自由振動モードを用いようとしたとしても, 式(10.105) (10.106)のような直交性は, 減衰行列に対しては必ずしも成立しないはずである。 したがって,$N$個の連立微分方程式を,独立した$N$個の微分方程式に 分解することはできない。もしそのまま解きたい場合には, 次の節で説明する数値的な手法を用いざるを得ないだろう。 しかし機械部品とは異なり,特に社会基盤構造を対象とする場合には幸いなことに, 実は具体的に減衰係数$c_i$がわかっているわけではない。 したがって,ここではモード解析法が使えるのではないかと期待(近似)して, その解法をまず説明しよう。すなわち

\begin{displaymath}
\fat{\phi}_i\supersc{t} \fat{C} \fat{\phi}_j=
\left\{\begi...
...m}_i \overline{\beta}_i \omega_i, & j=i
\end{array}\right.
\end{displaymath} (10.115)

が成立することを前提にして定式化してみよう。 ここに $\overline{\beta}_i$はこの式から求められる定数で, 第$i$次の見かけ上の減衰定数,あるいは一般化された減衰定数 と呼ばれる。

最初は,外力$\fat{f}(t)$が無い自由振動を対象とする。 運動方程式(10.95)の減衰自由振動バージョンの $\fat{f}(t)\equiv 0$の 式に式(10.107)を代入すると

\begin{displaymath}
\sum_{n=1}^N \ddot{q}_n(t) \fat{M} \fat{\phi}_n
+\sum_{n=1...
...at{\phi}_n
+\sum_{n=1}^N q_n(t) \fat{K} \fat{\phi}_n=\fat{0}
\end{displaymath}

となる。これに左から $\fat{\phi}_j\supersc{t}$を乗じ(つまり$j$次モードとの 仮想仕事を算定し), 式(10.105) (10.106) (10.116)の 直交条件を利用すると,総和のうち$j\neq n$の項はすべて零になるので

\begin{displaymath}
\sum_{n=1}^N \ddot{q}_n(t) \fat{\phi}_j\supersc{t} \fat{M}...
... \omega_j \dot{q}_j(t)
+\overline{m}_j \omega_j^2 q_j(t)=0
\end{displaymath}

となり,結局

\begin{displaymath}
\ddot{q}_j(t)
+2 \overline{\beta}_j \omega_j \dot{q}_j(t)
+\omega_j^2 q_j(t)=0
\end{displaymath} (10.116)

を解けばいいことになる。これは1自由度系の減衰自由振動の運動方程式である。 したがって一般解は

\begin{displaymath}
q_j(t)=\exp\left(-\overline{\beta}_j \omega_j t\right)
\l...
...ega_d\right)_j\equiv \omega_j\sqrt{1-\overline{\beta}_j{}^2}
\end{displaymath} (10.117)

となる。初期条件式(10.110)を用いて積分定数$A_j$, $B_j$を 決定すれば解が決定される。

図-10.35 (p.[*])に 示した破線が,2自由度系で$m_1=m_2=m$, $c_1=c_2=c$, $k_1=k_2=k$とした上で

\begin{displaymath}
c=2 \beta_0 \omega_0, \quad \beta_0=0.02, \quad
\omega_0\equiv\sqrt{\dfrac{k}{m}}
\end{displaymath}

とし,初期条件を

\begin{displaymath}
\fat{u}_0=\left\{\begin{array}{c} u_0  0 \end{array}\right...
... v_0 \end{array}\right\}, \quad
\dfrac{v_0}{\omega_0}=0.5 u_0
\end{displaymath}

と与えたときの結果である。前節の例の記号を用いれば, 二つの一般化された減衰定数は

\begin{displaymath}
\overline{\beta}_i=\dfrac{\beta_0 \left(2-2d_i+d_i^2\right)}{\xi_i \mu_i}
\end{displaymath}

と表現できる。

10.2.2.7 Rayleigh減衰

しかし式(10.116)の「積極的に無視」というのは, いかに土木屋といえども, あまり好ましい状況とは思えない。また上述のように社会基盤構造の 場合には,具体的な減衰係数$c_i$については経験や実測から 求めなければならないのも事実である。 そこで,減衰係数が質量と剛性に比例したものとしてモデル化する 考え方[15]がある。これを比例減衰 あるいはRayleigh減衰 と呼ぶ。それは

\begin{displaymath}
\fat{C}=\zeta\supersc{m} \fat{M}+\zeta\supersc{k} \fat{K}
\end{displaymath} (10.118)

のようなモデル化である。このようにすれば, 式(10.105) (10.106)で表される 質量行列と剛性行列を通してのモードの直交性から, 減衰行列を通しても自動的にモードが直交することになる。 さて,粘性が弾性や質量に関係があるって本当だろうか。 後述のある種の梁モデルでは,剛性比例にはなることを示すことはできる。

図 10.34: 比例減衰
\begin{figure}\begin{center}
\unitlength=.01mm
\begin{picture}(6500,4675)(1500,-...
...1,Legend(Title)
%,-1,Graphics End
%E,0,
%
\end{picture}\end{center}
\end{figure}

式(10.119)に左右から$\fat{\phi}_i$を乗じると,

\begin{displaymath}
\fat{\phi}_i\supersc{t} \fat{C} \fat{\phi}_i
=\zeta\supers...
...\zeta\supersc{m}
+\overline{m}_i \omega_i^2 \zeta\supersc{k}
\end{displaymath}

という関係になる。ここでは 式(10.105) (10.106)の関係を用いた。これが, 形式的に式(10.18)の1自由度系の減衰項の 係数の質量倍の $2 \overline{m}_i \overline{\beta}_i \omega_i$という 表現になって欲しいことから,一般化された比例減衰定数を

\begin{displaymath}
\overline{\beta}_i=\overline{\beta}_i{}\supersc{m}
+\overli...
...{\beta}_i{}\supersc{k}=\dfrac{\zeta\supersc{k} \omega_i}{2}
\end{displaymath} (10.119)

のように定義すればいいことになる。したがって,減衰定数についても 直交性が成立し

\begin{displaymath}
\fat{\phi}_i\supersc{t} \fat{C} \fat{\phi}_j=
\left\{\begi...
...m}_i \overline{\beta}_i \omega_i, & j=i
\end{array}\right.
\end{displaymath} (10.120)

という関係が成り立つ。 この不思議な係数 $\zeta\supersc{m}$, $\zeta\supersc{k}$を例えば実測から 決定するには,例えば第$m$次モードで振動させたときの 減衰定数 $\overline{\beta}_m$と第$n$次モードの $\overline{\beta}_n$とを 測定しておく。これを式(10.120)に代入して連立させ

\begin{displaymath}
\left\{\begin{array}{c}
\zeta\supersc{m}  \zeta\supersc{k...
..._n}-\dfrac{\overline{\beta}_n}{\omega_m}
\end{array}\right\}
\end{displaymath} (10.121)

という関係から求めればいい。実際には,三つ以上の振動数に 対する減衰定数を用いて最小2乗法でも使うのではないだろうか。

10.2.2.8 粘性減衰強制振動

最後に,外力も含めた運動方程式を,モード解析法で解いておこう。 運動方程式(10.95)に式(10.107)を代入すると

\begin{displaymath}
\sum_{n=1}^N \ddot{q}_n(t) \fat{M} \fat{\phi}_n
+\sum_{n=1...
...\phi}_n
+\sum_{n=1}^N q_n(t) \fat{K} \fat{\phi}_n=\fat{f}(t)
\end{displaymath}

となる。これに左から $\fat{\phi}_j\supersc{t}$を乗じ(つまり$j$次モードとの 仮想仕事を算定し), 式(10.105) (10.106) (10.121)の 直交条件を利用すると,総和のうち$j\neq n$の項はすべて零になるので

\begin{eqnarray*}
&&\sum_{n=1}^N \ddot{q}_n(t) \fat{\phi}_j\supersc{t} \fat{M}...
...ot{q}_j(t)
+\overline{m}_j \omega_j^2 q_j(t)=\overline{f}_j(t)
\end{eqnarray*}

となる。ここに右辺の $\overline{f}_j(t)$は, 式(10.113)で定義した一般化された外力である。結局

\begin{displaymath}
\ddot{q}_j(t)
+2 \overline{\beta}_j \omega_j \dot{q}_j(t)...
...ega_j^2 q_j(t)=\dfrac{1}{\overline{m}_j} \overline{f}_j(t)
\end{displaymath} (10.122)

を解けばいいことになる。これは1自由度系の減衰強制振動の運動方程式である。

結局,モード解析法を用いさえすれば,本質的には1自由度系の運動方程式を 解くことができればいいことが明らかである。実は, 実際の振動過程を見ると,振動に含まれるすべての成分のうち, 低次の固有振動モードの成分の方が比較的大きくなる。 これは,高周波運動をさせるには大きなエネルギが必要となり, その逆の現象として,低周波での振動の方が発生し易いからである。 したがって,工学的な観点から社会基盤構造の設計を念頭に置くと, 例えば$N$自由度系の固有振動数とモードのすべてを知る必要があるわけでは ないことがわかる。工学的に重要と判断される低次の$M (\leq N)$個の 固有振動数と固有振動モードを求め,その一般化された 座標を求めてしまえば,必要な精度での近似解を得ることができる。 この$M$がいくつくらいなら適切かという判断は, 現場の技術者の技量によるのだろう。モードをじっくり考えながら, どのくらいまでのモードを対象にすればいいのか判断すべきである。 特に梁や柱・骨組のような構造振動の場合には, 振動の仕方についての力学的な予測・判断が非常に重要である。 そういう意味で力学的センスを身に着けることはとても大事である。 なお,減衰強制振動においても非減衰自由振動モードを用いて モード解析法ができることは, 初期値境界値問題がその固有関数によって解析できることに 通じている。


10.2.3 数値解析手法

図 10.35: 2自由度系の減衰自由振動
\begin{figure}\begin{center}
\unitlength=.01mm
\begin{picture}(6745,4250)(1500,-...
...1,Legend(Title)
%,-1,Graphics End
%E,0,
%
\end{picture}\end{center}
\end{figure}

ここでは一つの例として,1自由度系の 数値解析法の節でも紹介した中央差分法 を多自由度系に対して適応しておく。 表示を簡単にするために

\begin{displaymath}
\fat{u}^-\equiv \fat{u}(t-\Delta t), \quad
\fat{u}\equiv \fat{u}(t), \quad \fat{u}^+\equiv \fat{u}(t+\Delta t)
\end{displaymath}

と定義したとき,速度と加速度は

$\displaystyle \dot{\fat{u}}(t)$ $\textstyle =$ $\displaystyle \dfrac{1}{2\Delta t}\left(\fat{u}^+-\fat{u}^-\right)$ (10.123)
$\displaystyle \ddot{\fat{u}}(t)$ $\textstyle =$ $\displaystyle \dfrac{1}{\Delta t}\left\{
\dfrac{\fat{u}^+-\fat{u}}{\Delta t}-\d...
...\}=\dfrac{1}{\left(\Delta t\right)^2}
\left(\fat{u}^-2\fat{u}+\fat{u}^+\right)$  

と近似される。これを運動方程式(10.95)に代入すると

\begin{displaymath}
\dfrac{\fat{M}}{\left(\Delta t\right)^2}
\left(\fat{u}^-2\...
...left(-\fat{u}^-+\fat{u}^+\right)
+\fat{K} \fat{u}=\fat{f}(t)
\end{displaymath}

となるので,整理すると

\begin{displaymath}
\fat{u}^+=\fat{D} 
\left[\fat{f}(t)
-\left(\fat{K}-\dfrac...
...elta t^2} \fat{M}
+\dfrac{1}{2\Delta t} \fat{C}\right)^{-1}
\end{displaymath} (10.124)

という関係を得,時々刻々の解を求めることができる。 また,一番最初のステップの$\fat{u}^-$については, 式(10.124)の関係から$\fat{u}^+$を消去して得る

\begin{displaymath}
\fat{u}^-=\fat{u}-\Delta t \dot{\fat{u}}+
\dfrac{\Delta t^2}{2} \ddot{\fat{u}}
\end{displaymath}

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

\begin{displaymath}
\ddot{\fat{u}}(0)=\fat{M}^{-1} \left\{\fat{f}(0)
-\fat{C} \dot{\fat{u}}(0)-\fat{K} \fat{u}(0)\right\}
\end{displaymath}

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

\begin{displaymath}
\fat{u}(-\Delta t)=
\left(\fat{I}
-\dfrac{\Delta t^2}{2} ...
...{\fat{u}}(0)
+\dfrac{\Delta t^2}{2} \fat{M}^{-1} \fat{f}(0)
\end{displaymath} (10.125)

から,時刻$t=0$$\fat{u}^-$を 算定できる。$\fat{I}$$N\times N$の 単位行列である。図-10.53にフローチャートを 示した。図-10.35には,節-10.2.2 (6)で例題とした2自由度系を モード解析法で解いた厳密解(破線)との比較を示した。 この例は, $\omega_0\equiv\sqrt{\slfrac{k}{m}}$と したときの時間ステップを $\omega_0 \Delta t=0.4$とした場合の結果である。 これを $\omega_0 \Delta t=0.3$にすると,数値解は厳密解とほとんど 重なってしまう。数値解の安定性のためには,時間ステップは

\begin{displaymath}
\Delta t\leq \dfrac{T_N}{\pi}
\end{displaymath} (10.126)

である必要[8]がある。ここに$T_N$は最小固有周期である。 この例の場合の条件は $\omega_0\Delta t\leq \omega_0 T_2\simeq 1.24$程度で あった。つまり自由度$N$が大きくなると,時間ステップを非現実的なくらい小さく しなければならないことを示している。もっと良好な安定性や正確性を有する 他の手法には,例えばNewmark法やWilsonの$\theta$法等があるが, それについては文献[8]を参照のこと。


最新版を正確に読むためには pdf ファイル をどうぞ。これは web 検索のための簡易旧版です。
next up previous contents index
Next: 10.3 連続体の振動 曲げ剛性の無い構造要素の振動 Up: 10. 振動論の基礎 Previous: 10.1 1自由度系の振動
Iwakuma Tetsuo
Mon, 18 Feb 2013 12:50:17 +0900 : Stardate [-28]8120.79