振動工学 における線形多自由度系の振動 (せんけいたじゆうどけいのしんどう)は、線形 な特性を持ち、さらに2以上の自由度 を持つ系で起きる振動 である。運動方程式 は一般的に連立2階常微分方程式となり、行列 およびベクトル で表現される。
線形多自由度系の振動では、固有モードという多自由度系特有の概念が現れ、自由度の数だけ固有モードと固有振動数の組が存在する。固有モードの直交性によって、減衰 の無い系であれば固有モードごとの1自由度系の問題に帰着でき、振動解析を容易化できるのが特徴である。この手法を利用した振動解析手法はモード解析 と呼ばれる。減衰のある系でも、比例粘性減衰という仮定を導入することによって、同様なことが可能となる。
モード解析手法は、振動実験結果から振動特性を同定 するのにも使われる。有限要素法 による連続体 の振動計算においても、線形多自由度系の理論にもとづくモード解析手法が強力な効果を発揮し、振動解析を容易にする。
背景
振動 の問題を分類すると、あるいは振動現象を模した動力学モデルを分類すると、いくつかの視点が存在する。
質点 がばね とダッシュポット を介して基礎に固定されたモデル
物体の運動を表すのに必要な座標 あるいは変位 の数を、自由度 という。自由度が1の系を1自由度系 という。ある1つの質点 がばね とダッシュポット を介して地面に固定され、質点が上下方向のみ動く例を考える。三角関数 で表現される励振力がこの質点に加わるとき、この質点の運動方程式 は次のような形で与えられる。
m
x
¨
+
c
x
˙
+
k
x
=
f
0
cos
ω
t
{\displaystyle m{\ddot {x}}+c{\dot {x}}+kx=f_{0}\cos \omega t}
(1.1 )
ここで、m は質量 、c は粘性減衰係数 、k はばね定数 、f 0 は外力の振幅 、ω は外力の角振動数 、t は時間 である。x は静的釣り合い位置からの変位 で、上部の "˙" は時間微分 (d /dt ) を示す。式1.1 は、振動問題を考える際の最も基礎的な式となる。1自由度系の振動の問題は、振動現象を理解する上で不可欠な様々な概念を内包しており、1自由度系のモデルが振動問題を扱う基本モデルといえる。
一方、自由度が2以上の系は多自由度系 と呼ばれる。振動問題をモデル化する際には1自由度や2自由度にモデル化すれば解析が容易であり、モデル化時の自由度はできるだけ最小限とするのが大切である。しかし、現実の構造物は複雑で、すべてを1自由度系として扱うことができない。実際の機械や建物などでは、1自由度系のモデルではその特性を説明しきれず、多自由度系としての取り扱いが必要になることも多い。
また、系の構成要素が全て線形 であるとき、その系を線形系 という。系が線形であれば、出力と入力の間には重ね合わせの原理 が成り立ち、出力は入力に単純に比例 する。振動系が線形であれば、式1.1 のように慣性力、減衰力、復元力が、それぞれ加速度 、速度 、変位 に比例する。比例係数が時間変化する場合を除き、線形系の振動問題については理論的手法がほぼ確立できている。
一方で、線形ではない系を非線形系 という。摩擦、がたつき、大変位時の材料特性などが非線形の要因となる。厳密に考えると、実際の機械や構造物はなんらかの非線形特性を持ち、実際の振動現象のほとんどは非線形系といえる。非線形系の厳密解を得ることはほとんどできず、一般的には非線形系の問題に対しては近似解法や数値計算に頼らざるを得ない。しかし、線形系の振動として取り扱うことによって、十分に多くの問題を解決することもできる。非線形系であっても、安定な平衡状態周りの微小変位運動であれば線形の理論が当てはまる。特に線形系であれば、多自由度系であっても後述のように1自由度系の問題に帰着できるという顕著な特性がある。
運動方程式
一般の多自由度系
多自由度系では多数の変数や係数を扱うため、行列 とベクトル を使って運動方程式 を記す。自由度が n の線形多自由度系の一般的な運動方程式は、次のように書き表される。
M
x
¨
+
C
x
˙
+
K
x
=
f
{\displaystyle {\boldsymbol {M{\ddot {x}}}}+{\boldsymbol {C{\dot {x}}}}+{\boldsymbol {Kx}}={\boldsymbol {f}}}
(2.1 )
ここで、ẍ , ẋ , x , f は n 次元縦ベクトル、M , C , K は n 次正方行列 で、以下のように表される。
x
¨
=
(
x
¨
1
x
¨
2
⋮
x
¨
n
)
{\displaystyle {\boldsymbol {\ddot {x}}}={\begin{pmatrix}{\ddot {x}}_{1}\\{\ddot {x}}_{2}\\\vdots \\{\ddot {x}}_{n}\end{pmatrix}}}
,
x
˙
=
(
x
˙
1
x
˙
2
⋮
x
˙
n
)
{\displaystyle {\boldsymbol {\dot {x}}}={\begin{pmatrix}{\dot {x}}_{1}\\{\dot {x}}_{2}\\\vdots \\{\dot {x}}_{n}\end{pmatrix}}}
,
x
=
(
x
1
x
2
⋮
x
n
)
{\displaystyle {\boldsymbol {x}}={\begin{pmatrix}x_{1}\\x_{2}\\\vdots \\x_{n}\end{pmatrix}}}
,
f
=
(
f
1
f
2
⋮
f
n
)
{\displaystyle {\boldsymbol {f}}={\begin{pmatrix}f_{1}\\f_{2}\\\vdots \\f_{n}\end{pmatrix}}}
M
=
(
m
11
m
12
⋯
m
1
n
m
21
m
22
⋯
m
2
n
⋮
⋮
⋱
⋮
m
n
1
m
n
2
⋯
m
n
n
)
{\displaystyle {\boldsymbol {M}}={\begin{pmatrix}m_{11}&m_{12}&\cdots &m_{1n}\\m_{21}&m_{22}&\cdots &m_{2n}\\\vdots &\vdots &\ddots &\vdots \\m_{n1}&m_{n2}&\cdots &m_{nn}\end{pmatrix}}}
(2.2 )
C
=
(
c
11
c
12
⋯
c
1
n
c
21
c
22
⋯
c
2
n
⋮
⋮
⋱
⋮
c
n
1
c
n
2
⋯
c
n
n
)
{\displaystyle {\boldsymbol {C}}={\begin{pmatrix}c_{11}&c_{12}&\cdots &c_{1n}\\c_{21}&c_{22}&\cdots &c_{2n}\\\vdots &\vdots &\ddots &\vdots \\c_{n1}&c_{n2}&\cdots &c_{nn}\end{pmatrix}}}
(2.3 )
K
=
(
k
11
k
12
⋯
k
1
n
k
21
k
22
⋯
k
2
n
⋮
⋮
⋱
⋮
k
n
1
k
n
2
⋯
k
n
n
)
{\displaystyle {\boldsymbol {K}}={\begin{pmatrix}k_{11}&k_{12}&\cdots &k_{1n}\\k_{21}&k_{22}&\cdots &k_{2n}\\\vdots &\vdots &\ddots &\vdots \\k_{n1}&k_{n2}&\cdots &k_{nn}\end{pmatrix}}}
(2.4 )
x は変位ベクトル と呼ばれる。その成分の x は、静的な釣り合いの位置を原点とした各自由度の変位 あるいは一般化座標 を表している。ẋ は速度ベクトル 、ẍ は加速度ベクトル と呼ばれる。成分の ẋ は x の時間1階微分すなわち各自由度の速度 を意味し、ẍ は x の時間2階微分すなわち各自由度 の加速度を意味する。
f は外力ベクトル と呼ばれる。成分の f は各一般化座標に対応して作用する外力で、一般的には時間の関数である。
M は質量行列 や慣性行列 、質量マトリックス や慣性マトリックス と呼ばれる。K は剛性行列 や剛性マトリックス と呼ばれる。C は減衰行列 や減衰マトリックス と呼ばれる。各成分の m , c , k は質量、粘性減衰係数、剛性 (ばね定数)を表す。ここでの m は慣性係数 ともいい、いわゆる質量 だけでなく慣性モーメント なども含む。また、k は復元係数 ともいい、通常のばね定数の他に回転ばね定数なども含む。
線形多自由度系の一般基礎式2.1 は、線形1自由度系の基礎式1.1 と形式は同じで、変位と外力がベクトルに、質量とばね定数と粘性減衰係数が行列に置き換わった式となる。ただし、対象の系の規模が大きくなり、複雑なものとなると、ニュートンの運動法則 から運動方程式を導出するのは容易ではない。実際に多自由度系の運動方程式を立てる際は、エネルギーのスカラー量から形式的に運動方程式を導けるラグランジュの運動方程式 が便利で、間違いを犯しにくく、多用されている。
2自由度系の例
2層構造物の揺れのモデリング
建築物の例では、複数の階を持つ多層構造物が多自由度系の問題となる。柱と床から構成されるラーメン構造 の建物が振動する場合を考える。柱は床に比べて柔らかいので、構造物が揺れるとき各階の床は水平方向に揺れ、柱は水平方向のばねとして働くと見なせる。これは、水平方向のみに動く2つの質点をばねで連結した2自由度系モデルと等価となる。2層構造物における、1層目の床質量を m 1 、2層目の床質量を m 2 、1層目の水平変位を x 1 、2層目の水平変位を x 2 、基礎と1層目の間の等価ばね定数を k 1 、1層目と2層目の間の等価ばね定数を k 2 とすると、このモデルの運動方程式は
{
m
1
x
¨
1
+
(
k
1
+
k
2
)
x
1
−
k
2
x
2
=
0
m
2
x
¨
2
−
k
2
x
1
+
k
2
x
2
=
0
{\displaystyle {\begin{cases}m_{1}{\ddot {x}}_{1}+(k_{1}+k_{2})x_{1}-k_{2}x_{2}=0\\m_{2}{\ddot {x}}_{2}-k_{2}x_{1}+k_{2}x_{2}=0\end{cases}}}
(2.5 )
となる。行列表示すれば、
(
m
1
0
0
m
2
)
(
x
¨
1
x
¨
2
)
+
(
k
1
+
k
2
−
k
2
−
k
2
k
2
)
(
x
1
x
2
)
=
(
0
0
)
{\displaystyle {\begin{pmatrix}m_{1}&0\\0&m_{2}\end{pmatrix}}{\begin{pmatrix}{\ddot {x}}_{1}\\{\ddot {x}}_{2}\end{pmatrix}}+{\begin{pmatrix}k_{1}+k_{2}&-k_{2}\\-k_{2}&k_{2}\end{pmatrix}}{\begin{pmatrix}x_{1}\\x_{2}\end{pmatrix}}={\begin{pmatrix}0\\0\end{pmatrix}}}
(2.6 )
となり、この式では、
(
m
1
0
0
m
2
)
{\displaystyle {\begin{pmatrix}m_{1}&0\\0&m_{2}\end{pmatrix}}}
が質量行列、
(
k
1
+
k
2
−
k
2
−
k
2
k
2
)
{\displaystyle {\begin{pmatrix}k_{1}+k_{2}&-k_{2}\\-k_{2}&k_{2}\end{pmatrix}}}
が剛性行列である。
自動車の上下振動とピッチング振動の簡易モデル
並進運動と回転運動が組み合わさった2自由度系の例として、自動車の上下振動とピッチング振動の簡易モデルがある。自動車を側面から見て、上下方向変位とピッチング回転だけの動きを考える。前位側のタイヤとサスペンションを1つのばねと見なして、後位側のタイヤとサスペンションも同様に1つのばねと見なし、剛体 を前位と後位を2つのばねが支えているモデルを考える。前位側と後位側のばねのばね定数をそれぞれ k 1 、k 2 とする。剛体の重心位置から前位側のばねまでの距離を l 1 、前位側のばねまでの距離を l 2 とする。剛体の質量を m 、重心位置周りのピッチング方向の慣性モーメント を IG とする。この剛体の重心の上下運動 x とその周りのピッチング運動 θ の運動方程式は次のようになる。
{
m
x
¨
+
(
k
1
+
k
2
)
x
+
(
k
2
l
2
−
k
1
l
1
)
θ
=
0
I
G
θ
¨
+
(
k
2
l
2
−
k
1
l
1
)
x
+
(
k
2
l
2
2
+
k
1
l
1
2
)
θ
=
0
{\displaystyle {\begin{cases}m{\ddot {x}}+(k_{1}+k_{2})x+(k_{2}l_{2}-k_{1}l_{1})\theta =0\\I_{G}{\ddot {\theta }}+(k_{2}l_{2}-k_{1}l_{1})x+(k_{2}l_{2}^{2}+k_{1}l_{1}^{2})\theta =0\end{cases}}}
(2.7 )
連成項
自動車の上下振動とピッチング振動の簡易モデルである式2.7 において、もし k 1 l 1 = k 2 l 2 であれば、x と θ は互いに独立した運動方程式となる。この条件では運動方程式は
{
m
x
¨
+
(
k
1
+
k
2
)
x
=
0
I
G
θ
¨
+
(
k
2
l
2
2
+
k
1
l
1
2
)
θ
=
0
{\displaystyle {\begin{cases}m{\ddot {x}}+(k_{1}+k_{2})x=0\\I_{G}{\ddot {\theta }}+(k_{2}l_{2}^{2}+k_{1}l_{1}^{2})\theta =0\end{cases}}}
(2.8 )
となり、2つの式はそれぞれ単独で解くことができる。すなわち、x の振動と θ の振動が、互いに干渉すること無しに独立して起こる状態になっている。このような状態を非連成 という。運動方程式が単独では解けない連立方程式になっているとき、方程式は連成 しているといい、式2.7 の (k 1 l 1 − k 2 l 2 )x と (k 1 l 1 − k 2 l 2 )θ は連成項 と呼ばれる[ 47] 。k 1 l 1 = k 2 l 2 のときに、これら連成項が 0 になり、x と θ が独立した運動になる。連成状態にある系を連成系 、逆に各自由度の運動が完全に独立している系を非連成系 と呼ぶ。
一般に、質量行列 M (式2.2 )、減衰行列 C (式2.3 )、剛性行列 K (式2.4 )における非対角成分の存在は、振動が連成していることを示している。もし系が非連成系であれば、これら全ての行列は、非対角成分が全て 0 の対角行列 となる。特に、質量行列が非対角成分によって振動が連成していることを動連成 や動的連成 といい、質量行列の非対角成分を動連成項 という[ 50] 。また、剛性行列が非対角行列であれば静連成 や静的連成 といい、剛性行列の非対角成分を静連成項 という[ 50] 。多自由度系の問題は概して複雑で、個々のパラメータが結果に与える影響の見通しを立てることが難しい。多自由度系では一般的に連成が存在し、これが多自由度系の解析を難しいものにしている。動力学的な設計を行う上では、系を非連成化して影響をわかりやすくすることが有効となる。
不減衰自由振動
振動数方程式
継続的な外力が作用せず、外乱だけが与えられて起こる振動を自由振動 という。自由振動の振動は、系自体が持つ特性によって決まり、系の動特性を知る上で重要な振動の形態である。自由振動は減衰要素が存在しない場合と存在する場合に分かれる。減衰要素と励振力が存在しない場合の式2.1 は
M
x
¨
+
K
x
=
0
{\displaystyle {\boldsymbol {M{\ddot {x}}}}+{\boldsymbol {Kx}}=\mathbf {0} }
(3.1 )
となる。ここで、0 は下記のようなゼロベクトル である。
0
=
(
0
0
⋮
0
)
{\displaystyle \mathbf {0} ={\begin{pmatrix}0\\0\\\vdots \\0\end{pmatrix}}}
さらに、M と K は正定値 であると仮定する。すなわち、M と K は成分が全て実数 の対称行列 で、それらの二次形式は常に正となる。実際に、質量行列と剛性行列が正定値行列であることは線形多自由度系の一般的な特徴であり、多くの振動系でこの仮定は満たされる。この条件を満たす M と K を前提にすれば、振動系の具体的な構成に依存しない一般性の高い議論を展開できる。
式3.1 は定数係数の線形常微分方程式 であるため、その解法に従って解を
x
=
u
e
λ
t
{\displaystyle {\boldsymbol {x}}={\boldsymbol {u}}e^{\lambda t}}
(3.2 )
という形式で仮定できる。ここで e はネイピア数 、λ は未知定数、 u は n 個の未知定数 u から成る下記のような縦ベクトルである。
u
=
(
u
1
u
2
⋮
u
n
)
{\displaystyle {\boldsymbol {u}}={\begin{pmatrix}u_{1}\\u_{2}\\\vdots \\u_{n}\end{pmatrix}}}
式3.1 に対して仮定として与えられる解には、他に、最初から単振動を仮定して複素指数関数 や三角関数 の形式もある。
式3.2 を式3.1 に代入して整理すると下記のような式になる。
(
λ
2
M
+
K
)
u
e
λ
t
=
0
{\displaystyle (\lambda ^{2}{\boldsymbol {M}}+{\boldsymbol {K}}){\boldsymbol {u}}e^{\lambda t}=\mathbf {0} }
(3.3 )
この式が恒等的に成り立つには下記の条件が満たされる必要がある。
(
λ
2
M
+
K
)
u
=
0
{\displaystyle (\lambda ^{2}{\boldsymbol {M}}+{\boldsymbol {K}}){\boldsymbol {u}}=\mathbf {0} }
(3.4 )
数学的には、この形式の方程式は一般化固有値問題 として知られ、λ 2 は固有値 、u は固有ベクトルと呼ばれる。u = 0 であれば式3.4 の条件は満たされるが、これは最初の釣り合いの位置からそのまま静止しているだけ状態を意味する解である。よって、ここでは u = 0 は興味の対象外で、それ以外の式3.4 を満たす λ と u について知りたい。u ≠ 0 で、なおかつ式3.4 が満たされる条件は、U の係数行列 λ 2 M + K の逆行列 が存在しないことである。したがって、係数行列の行列式 が 0 であればこの条件が満たされる。したがって、λ が、
det
(
λ
2
M
+
K
)
=
0
{\displaystyle {\mbox{det}}(\lambda ^{2}{\boldsymbol {M}}+{\boldsymbol {K}})=0}
(3.5 )
を満たすとき、式3.4 が u = 0 以外の解を持つ。ここで、det( ) は行列式を表す。式3.5 の行列式を展開すると、 λ 2 についての n 次多項式になる。したがって、原理的には λ 2 の値を求めることができる。ただし、この多項式を解析的に解くことができるのはせいぜい2自由度あるいは3自由度までで、それ以上の自由度の系になると数値解析 で固有値を計算する。ヤコビ法 などを使って数値計算するときは、コレスキー分解 を使い、一般化固有値問題形式の3.4 を標準的な固有値問題 の形へ変換する。
上記のとおり、M と K は正定値であると仮定した。このとき、式3.5 の解は全て負の実数となる(解を複素指数関数と仮定した場合は全て正の実数)。したがって、n 個の λ 2 の解を
λ
2
=
−
ω
1
2
,
−
ω
2
2
,
⋯
,
−
ω
n
2
{\displaystyle \lambda ^{2}=-\omega _{1}^{2},\ -\omega _{2}^{2},\ \cdots ,\ -\omega _{n}^{2}}
(3.6 )
とおくことができ、平方根を取って 2n 個のλ の値
λ
=
±
j
ω
1
,
±
j
ω
2
,
⋯
,
±
j
ω
n
{\displaystyle \lambda =\pm j\omega _{1},\ \pm j\omega _{2},\ \cdots ,\ \pm j\omega _{n}}
(3.7 )
が得られる。ここで j は虚数単位 である。n 個の ω は固有角振動数 または固有円振動数 と呼ばれ、値が小さいものから順に1次 、2次 、…、n 次 の固有角振動数と呼ぶ。特に、最も値が小さい1次の固有角振動数は基本振動数 と呼ばれる。式3.5 は角振動数を求める式であるため振動数方程式 と呼ばれる。あるいは、式3.5 は M や K といった系の特性によって構成される式であることから特性方程式 とも呼ぶ。
固有モード
得られた固有値あるいは固有角振動数の値を式3.4 に代入すれば、0 以外の u の解が得られる。ここで得られる u の各成分は一意な値を持たず、定まるのは各成分の互いの比 u 1 : u 2 : … : un だけである。一つの固有値あるいは固有角振動数に対応して一つの u が定まり、u は n 個存在する。このような固有値と固有ベクトルの組は固有ペアと呼ばれる。r 次の固有角振動数 ωr に対応する u を u r を表現すれば、 u r は
u
r
=
(
u
1
r
u
2
r
⋮
u
n
r
)
{\displaystyle {\boldsymbol {u}}_{r}={\begin{pmatrix}u_{1r}\\u_{2r}\\\vdots \\u_{nr}\end{pmatrix}}}
(3.8 )
というベクトルである。ただし、上記のとおり、各成分 u 1r , u 2r , …, unr は互いの比の値を表している。
式2.5 で表現される2層構造物振動モデルの固有モードの例。k 1 = k 2 , m 1 = m 2 の場合を示しており、各固有モードの成分の内の一番大きな絶対値が1となるようにして値を定めている。
式3.8 で表されるベクトルが式3.4 における固有ベクトルであり、振動工学では固有モード 、振動モード 、固有振動モード 、基準振動モード 、モードベクトル などと呼ぶ。n 自由度系には n 個の固有角振動数があり、固有角振動数それぞれに対応する形で n 個の固有モードが存在している。各自由度の振幅比を決める固有モードは、固有角振動数が「振動の速さ」を表しているのに対して、「振動の形」を表していると言える。自由度の数の分だけ固有角振動数が存在し、それら固有角振動数に対応して固有モードが存在していることが、線形多自由度系の特有な性質といえる。
各々の固有モードを、対応する固有角振動数が小さい順に1次固有モード、2次固有モード、…、r 次固有モード、…、n 次固有モードと呼ぶ。特に、基本振動数(1次の固有角振動数)に対応する固有モードは基本モード と呼ばれる。固有モード u r を次数が低い順に並べて作る、下記のような行列をモード行列 やモードマトリックス という。
U
=
(
u
1
,
u
2
,
⋯
,
u
n
)
=
(
u
11
u
21
⋯
u
1
n
u
21
u
22
⋯
u
2
n
⋮
⋮
⋱
⋮
u
n
1
u
n
2
⋯
u
n
n
)
{\displaystyle {\boldsymbol {U}}=({\boldsymbol {u}}_{1},\ {\boldsymbol {u}}_{2},\ \cdots ,\ {\boldsymbol {u}}_{n})={\begin{pmatrix}u_{11}&u_{21}&\cdots &u_{1n}\\u_{21}&u_{22}&\cdots &u_{2n}\\\vdots &\vdots &\ddots &\vdots \\u_{n1}&u_{n2}&\cdots &u_{nn}\end{pmatrix}}}
(3.9 )
固有モードの直交性と正規化
多自由度系の固有モードには直交性 という重要な性質がある。r 次の固有モード u r と s 次の固有モード u s について考える。ここで、r と s は任意だが、r ≠ s である。M と K は上記のとおり対称行列とする。このとき、 u r , u s , M , K には次のような関係がある。
u
r
⊤
M
u
s
=
0
{\displaystyle {\boldsymbol {u}}_{r}^{\top }{\boldsymbol {M}}{\boldsymbol {u}}_{s}=0}
(3.10 )
u
r
⊤
K
u
s
=
0
{\displaystyle {\boldsymbol {u}}_{r}^{\top }{\boldsymbol {K}}{\boldsymbol {u}}_{s}=0}
(3.11 )
ここで、u r ⊤ は u r の転置行列 を表している。これらの式は u r と Mu s の内積 および u r と Mu s の内積が零であることを示しており、M と K に関して u r と u s が直交 していることを意味している。これが質量行列および剛性行列を介した固有モードの直交性である。別の言い回しでは u r と u s が一般直交性を有しているともいう。
固有モードの直交性の物理的な意味は、次数の異なる固有モードの振動の間で力学的エネルギー の移動が全く起きないことを表している。ある固有モードの振動において、運動エネルギー と復元力によるポテンシャルエネルギー は時間的に変化しているが、それらの和の力学的エネルギーは時間に依らず一定に保たれている。そのため、自由振動中にある固有モードの振動が他の固有モードに移り変わったり、他の固有モードが新たに誘起されたりすることはない。
一方で、r = s 、すなわち同じ次数同士の固有モードの場合は、上記の式でも左辺は 0 とはならず、ある定数となる。これらの定数を Mr と Kr と表せば、
u
r
⊤
M
u
r
=
M
r
{\displaystyle {\boldsymbol {u}}_{r}^{\top }{\boldsymbol {M}}{\boldsymbol {u}}_{r}=M_{r}}
(3.12 )
u
r
⊤
K
u
r
=
K
r
{\displaystyle {\boldsymbol {u}}_{r}^{\top }{\boldsymbol {K}}{\boldsymbol {u}}_{r}=K_{r}}
(3.13 )
と表される。固有モード u r は成分間の比だけを持ち、定まった値を持たないので、Mr と Kr の値もこの段階では定まらない。固有モードの絶対値が決まった後に、Mr と Kr の値も定まる。
定数 Mr と Kr は正の値であり、それぞれを(r 次の)モード質量 およびモード剛性 という。モード行列を使って式3.12 と式3.13 を表せば、
U
⊤
M
U
=
(
M
1
0
⋯
0
0
M
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
M
n
)
{\displaystyle {\boldsymbol {U}}^{\top }{\boldsymbol {M}}{\boldsymbol {U}}={\begin{pmatrix}M_{1}&0&\cdots &0\\0&M_{2}&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &M_{n}\end{pmatrix}}}
(3.14 )
U
⊤
K
U
=
(
K
1
0
⋯
0
0
K
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
K
n
)
{\displaystyle {\boldsymbol {U}}^{\top }{\boldsymbol {K}}{\boldsymbol {U}}={\begin{pmatrix}K_{1}&0&\cdots &0\\0&K_{2}&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &K_{n}\end{pmatrix}}}
(3.15 )
となる。モード質量とモード剛性は、同じ次数の固有角振動数 ωr と下記のような関係がある。
ω
r
=
K
r
M
r
{\displaystyle \omega _{r}={\sqrt {\frac {K_{r}}{M_{r}}}}}
(3.16 )
固有モードは比の関係を定めるベクトルであった。絶対的な大きさを持たない固有モードに大きさを定める方法として、次のような方法がある。
下記のように、式3.12 の右辺 Mr の値が 1 となるように定める。ここで、ū r は、等式を満たすに大きさが決定された固有モードを意味している。
u
¯
r
⊤
M
u
¯
r
=
1
{\displaystyle {\boldsymbol {\overline {u}}}_{r}^{\top }{\boldsymbol {M}}{\boldsymbol {\overline {u}}}_{r}=1}
(3.17 )
下記のように、同じ次数の固有モード同士の内積 が 1 となるように定める。
u
¯
r
⊤
u
¯
r
=
1
{\displaystyle {\boldsymbol {\overline {u}}}_{r}^{\top }{\boldsymbol {\overline {u}}}_{r}=1}
(3.18 )
固有モードの成分の内、絶対値 が最大のものを 1 とおいて定める。
1番目の自由度に対応する固有モード成分 (u 1r ) を 1 とおいて定める。
以上のように固有モードの大きさを一意に定めることを正規化 と呼び、正規化された固有モード ū r を正規固有モード や正規化モード と呼ぶ。1番目の手法による正規固有モードは、特に質量正規固有モード やM -正規固有モード という名で呼ばれる。式3.17 を満たすようにしておくと利点が多く、質量正規固有モードはよく使われる。正規固有モード ū r を並べて作る下記のような行列を、正規モード行列 という。
U
¯
=
(
u
¯
1
,
u
¯
2
,
⋯
,
u
¯
n
)
=
(
u
¯
1
,
1
u
¯
1
,
2
⋯
u
¯
1
,
n
u
¯
1
,
2
u
¯
2
,
2
⋯
u
¯
2
,
n
⋮
⋮
⋱
⋮
u
¯
1
,
n
u
¯
n
,
2
⋯
u
¯
n
,
n
)
{\displaystyle {\boldsymbol {\overline {U}}}=({\boldsymbol {\overline {u}}}_{1},\ {\boldsymbol {\overline {u}}}_{2},\ \cdots ,\ {\boldsymbol {\overline {u}}}_{n})={\begin{pmatrix}{\overline {u}}_{1,1}&{\overline {u}}_{1,2}&\cdots &{\overline {u}}_{1,n}\\{\overline {u}}_{1,2}&{\overline {u}}_{2,2}&\cdots &{\overline {u}}_{2,n}\\\vdots &\vdots &\ddots &\vdots \\{\overline {u}}_{1,n}&{\overline {u}}_{n,2}&\cdots &{\overline {u}}_{n,n}\end{pmatrix}}}
(3.19 )
質量正規固有モードを採用すると、式3.17 に対して剛性行列は固有角振動数と
u
¯
r
⊤
K
u
¯
r
=
ω
r
2
{\displaystyle {\boldsymbol {\overline {u}}}_{r}^{\top }{\boldsymbol {K}}{\boldsymbol {\overline {u}}}_{r}=\omega _{r}^{2}}
(3.20 )
という関係を持つ。
モード座標と一般解
固有モード(固有ベクトル)は互いに直交なので、これらを使った線形結合 で変位ベクトル x を表すことができる。つまり、固有モード U を用いて、x を
x
=
q
1
u
¯
1
+
q
2
u
¯
2
+
⋯
+
q
n
u
¯
n
=
U
¯
q
{\displaystyle {\boldsymbol {x}}=q_{1}{\boldsymbol {\overline {u}}}_{1}+q_{2}{\boldsymbol {\overline {u}}}_{2}+\dots +q_{n}{\boldsymbol {\overline {u}}}_{n}={\boldsymbol {{\overline {U}}q}}}
(3.21 )
と表すことができる。ここで、q は時間の未知関数で、モード座標 、基準座標 、正規座標 、規準座標 、主座標 などと呼ばれる。特に正規化した固有モードを使っているモード座標を指して、正規座標と呼ぶこともある。q は qr を 1 次から n 次まで並べた次のような縦ベクトルである。
q
=
(
q
1
q
2
⋮
q
n
)
{\displaystyle {\boldsymbol {q}}={\begin{pmatrix}q_{1}\\q_{2}\\\vdots \\q_{n}\end{pmatrix}}}
モード座標 q に対して、元の座標 x を物理座標 と呼ぶ。3.21 は、正規固有モードを介して物理座標をモード座標に変換していることを意味する。モード座標 q r は、x に対して r 次固有モードū r が寄与する度合いを表しているとも言える。式3.21 は展開定理とも呼ばれる。
式3.21 を運動方程式3.1 へ代入して、左から転置したモード行列 U ⊤ を掛けると次のようになる。
U
⊤
M
U
q
¨
+
U
⊤
K
U
q
=
0
{\displaystyle {\boldsymbol {U^{\top }MU{\ddot {q}}}}+{\boldsymbol {U^{\top }KUq}}=\mathbf {0} }
(3.22 )
固有モードの直交性(式3.14 と式3.15 )を上式に当てはめると、
{
M
1
q
¨
1
+
K
1
q
1
=
0
M
2
q
¨
2
+
K
2
q
2
=
0
⋮
M
n
q
¨
n
+
K
n
q
n
=
0
{\displaystyle {\begin{cases}M_{1}{\ddot {q}}_{1}+K_{1}q_{1}=0\\M_{2}{\ddot {q}}_{2}+K_{2}q_{2}=0\\\vdots \\M_{n}{\ddot {q}}_{n}+K_{n}q_{n}=0\\\end{cases}}}
(3.23 )
というような q に関する n 個の運動方程式が得られる。それぞれの式の両辺を Mr で割ると、下記のような形になる。
{
q
¨
1
+
ω
1
2
q
1
=
0
q
¨
2
+
ω
2
2
q
2
=
0
⋮
q
¨
n
+
ω
n
2
q
n
=
0
{\displaystyle {\begin{cases}{\ddot {q}}_{1}+\omega _{1}^{2}q_{1}=0\\{\ddot {q}}_{2}+\omega _{2}^{2}q_{2}=0\\\vdots \\{\ddot {q}}_{n}+\omega _{n}^{2}q_{n}=0\\\end{cases}}}
(3.24 )
式2.7 で表現される自動車の上下・ピッチング振動モデルの振動の例。図はピッチング θ の振動を示している。パラメータは m = 1600 kg, IG = 2500 kg-m2 , k 1 = 3500 N/m, k 2 = 4100 N/m, l 1 = 1.4 m, l 2 = 1.6 m で、初期条件は上下変位速度が ẋ = 1 m/s で他は全て 0 という条件の例。
式3.24 は非連成化されており、各式はそれぞれ独立している。そのため、一つ一つの式は1自由度系の不減衰自由振動と同じであるから、q の各解は以下のようになる。
{
q
1
=
c
1
sin
(
ω
1
t
+
θ
1
)
q
2
=
c
n
sin
(
ω
2
t
+
θ
2
)
⋮
q
n
=
c
n
sin
(
ω
n
t
+
θ
n
)
{\displaystyle {\begin{cases}q_{1}=c_{1}\sin(\omega _{1}t+\theta _{1})\\q_{2}=c_{n}\sin(\omega _{2}t+\theta _{2})\\\vdots \\q_{n}=c_{n}\sin(\omega _{n}t+\theta _{n})\\\end{cases}}}
(3.25 )
ここで、cr と θr (r = 1, 2, …, n ) は初期条件で決まる定数で、cr が振幅 、θr が初期位相 を意味する。x の一般解は、式3.25 を式3.21 に代入して
x
=
∑
r
=
1
n
u
¯
r
C
r
sin
(
ω
r
t
+
θ
r
)
{\displaystyle {\boldsymbol {x}}=\sum _{r=1}^{n}{\boldsymbol {\overline {u}}}_{r}C_{r}\sin(\omega _{r}t+\theta _{r})}
(3.26 )
と得られる。すなわち、n 自由度系における各自由度の自由振動は、n 個の調和振動 の重ね合わせ(和)となっており、それら調和振動のそれぞれの振動数は 1 次から n 次までの固有角振動数となっている。固有モードは、対応する固有角振動数を持つ調和振動成分の各自由度間の振幅比を定めている。位相角を陽に表さずに、余弦関数と正弦関数の和
x
=
∑
r
=
1
n
u
¯
r
(
a
r
cos
ω
r
t
+
b
r
sin
ω
r
t
)
{\displaystyle {\boldsymbol {x}}=\sum _{r=1}^{n}{\boldsymbol {\overline {u}}}_{r}(a_{r}\cos \omega _{r}t+b_{r}\sin \omega _{r}t)}
(3.27 )
や、指数の正負が異なる複素指数関数の和
x
=
∑
r
=
1
n
u
¯
r
(
d
r
e
i
ω
r
t
+
g
r
e
−
i
ω
r
t
)
{\displaystyle {\boldsymbol {x}}=\sum _{r=1}^{n}{\boldsymbol {\overline {u}}}_{r}(d_{r}e^{i\omega _{r}t}+g_{r}e^{-i\omega _{r}t})}
(3.28 )
といった形で一般解を表すこともできる。ar , br , dr , gr も初期条件で決まる定数で、式3.27 と式3.27 は式3.26 へ式変形可能な同値な式である。
以上のような、物理座標をモード座標へ変換し、固有モードの直交性を利用して多自由度系の問題を1自由度系の重ね合わせの問題に帰着させ、解析を行う手法をモード解析 (特にモード重畳法)という[ 115] 。後述のように、モード解析は特に有限要素法 へ適用することで有効性を発揮する。実物の振動特性を求める上でもモード解析の理論を適用することで各特性を同定 でき、この手法は実験モード解析 と呼ばれる。
減衰自由振動
一般の粘性減衰
基礎の上でばねと減衰器が一組となって質点と連結し、直列に連なった3自由度減衰系の例
減衰行列 C が存在する線形多自由度系の振動について考える。外力が無い場合の n 自由度系の運動方程式は以下のようになる。
M
x
¨
+
C
x
˙
+
K
x
=
0
{\displaystyle {\boldsymbol {M{\ddot {x}}}}+{\boldsymbol {C{\dot {x}}}}+{\boldsymbol {Kx}}={\boldsymbol {0}}}
(4.1 )
基礎の上でばねと減衰器が一組となって質点と連結し、直列に連なった3自由度減衰系の典型的な例では、M , C , K は次のような行列となる。
M
=
(
m
1
0
0
0
m
2
0
0
0
m
3
)
{\displaystyle {\boldsymbol {M}}={\begin{pmatrix}m_{1}&0&0\\0&m_{2}&0\\0&0&m_{3}\end{pmatrix}}}
(4.2 )
C
=
(
c
1
−
c
1
0
−
c
1
c
1
+
c
2
−
c
2
0
−
c
2
c
2
+
c
3
)
{\displaystyle {\boldsymbol {C}}={\begin{pmatrix}c_{1}&-c_{1}&0\\-c_{1}&c_{1}+c_{2}&-c_{2}\\0&-c_{2}&c_{2}+c_{3}\end{pmatrix}}}
(4.3 )
K
=
(
k
1
−
k
1
0
−
k
1
k
1
+
k
2
−
k
2
0
−
k
2
k
2
+
k
3
)
{\displaystyle {\boldsymbol {K}}={\begin{pmatrix}k_{1}&-k_{1}&0\\-k_{1}&k_{1}+k_{2}&-k_{2}\\0&-k_{2}&k_{2}+k_{3}\end{pmatrix}}}
(4.4 )
式4.1 の解を式3.2 と同じように仮定し、式3.2 を式4.1 に代入して整理すると、
(
λ
2
M
+
λ
C
+
K
)
u
=
0
{\displaystyle (\lambda ^{2}{\boldsymbol {M}}+\lambda {\boldsymbol {C}}+{\boldsymbol {K}}){\boldsymbol {u}}={\boldsymbol {0}}}
(4.5 )
となり、u = 0 以外の解を持つという条件から
det
(
λ
2
M
+
λ
C
+
K
)
{\displaystyle {\mbox{det}}(\lambda ^{2}{\boldsymbol {M}}+\lambda {\boldsymbol {C}}+{\boldsymbol {K}})}
(4.6 )
という特性方程式が得られる[ 120] 。式4.6 を展開すると、λ に関する 2n 次多項式となる[ 120] 。この多項式を解いて得た解 λ 1 , λ 2 , …, λ 2n を式4.5 に代入し、u を定め、u に適当な正規化を行う[ 120] 。正規化された各ベクトルを ū 1 , ū 2 , …, ū 2n と表すと、式4.1 の一般解は
x
=
u
¯
1
a
1
e
λ
1
t
+
u
¯
2
a
2
e
λ
2
t
+
⋯
+
u
¯
2
n
a
2
n
e
λ
2
n
t
{\displaystyle {\boldsymbol {x}}={\boldsymbol {\overline {u}}}_{1}a_{1}e^{\lambda _{1}t}+{\boldsymbol {\overline {u}}}_{2}a_{2}e^{\lambda _{2}t}+\cdots +{\boldsymbol {\overline {u}}}_{2n}a_{2n}e^{\lambda _{2n}t}}
(4.7 )
となる[ 120] 。ここで、a 1 , a 2 , …, a 2n は初期条件で決まる定数である[ 120] 。2n 個の λ は、解が減衰振動 であれば n 組の互いに共役な複素数 になる[ 120] 。このとき、固有モードも複素数となり、複素固有モード と呼ばれる。
比例粘性減衰
式4.1 のように一般的な減衰行列 C が運動方程式に存在する場合、正規モード行列によって M と K は対角化できるが、C も同時に対角化することはできない。そのため、不減衰自由振動で可能だったモード座標に変換しての非連成化が不可能となり、モード解析の利点を活かすことができなくなる。そこで、C が下記のような比例粘性減衰 として与えられると仮定し、実際の振動解析が行われることも多い。
C
=
α
M
+
β
K
{\displaystyle {\boldsymbol {C}}=\alpha {\boldsymbol {M}}+\beta {\boldsymbol {K}}}
(4.8 )
ここで、α と β は定数で、比例減衰定数 と呼ばれる[ 123] 。比例粘性減衰を有する系を比例粘性減衰系 などと呼ぶ。比例粘性減衰が成り立つと仮定すれば減衰行列を対角化できる。
式4.8 の形で与えられる比例粘性減衰は、特にレイリー減衰 やレイリー型減衰 と呼ばれる。αM のみで仮定されるものは質量比例型減衰 、βK のみで仮定されるものは剛性比例型減衰 などと呼ぶ。
実際の減衰が比例粘性減衰になっていることはまれであり、比例粘性減衰はあくまでも近似的なものである。しかし、比例粘性減衰の仮定を導入することで、不減衰系と同じ取り扱いが可能となり、モード解析の手法が適用可能になる。もし減衰が全体に分布しているような構造であれば、適当な比例減衰定数を設定すれば、実際の現象を実用問題ないレベルで再現できるという一定の妥当性もある。減衰は摩擦・材料減衰・流体粘性など様々な要因で起こるため、そもそも減衰の適切な定式化自体が難しいといった事情もある。式2.3 のように速度の比例定数として与えられる一般の粘性減衰も、多種多様な発生機構によって減衰が起きるという実情に起因して厳密には成立しない。
式4.8 を式4.1 に代入した場合、
M
x
¨
+
(
α
M
+
β
K
)
x
˙
+
K
x
=
0
{\displaystyle {\boldsymbol {M{\ddot {x}}}}+(\alpha {\boldsymbol {M}}+\beta {\boldsymbol {K}}){\boldsymbol {\dot {x}}}+{\boldsymbol {Kx}}={\boldsymbol {0}}}
(4.9 )
という運動方程式になる。解を式3.2 のように仮定して上式に代入し、整理すると、
{
(
λ
2
+
α
λ
)
M
+
(
β
λ
+
1
)
K
}
u
=
0
{\displaystyle \{(\lambda ^{2}+\alpha \lambda ){\boldsymbol {M}}+(\beta \lambda +1){\boldsymbol {K}}\}{\boldsymbol {u}}={\boldsymbol {0}}}
(4.10 )
となる。さらに、
γ
2
=
λ
2
+
α
λ
β
λ
+
1
{\displaystyle \gamma ^{2}={\frac {\lambda ^{2}+\alpha \lambda }{\beta \lambda +1}}}
(4.11 )
とおけば、式4.10 は
(
γ
2
M
+
K
)
u
=
0
{\displaystyle (\gamma ^{2}{\boldsymbol {M}}+{\boldsymbol {K}}){\boldsymbol {u}}={\boldsymbol {0}}}
(4.12 )
となる。式4.12 は、(複素指数関数や三角関数で解を仮定したときの)不減衰振動の固有値問題における λ を γ に置き換えただけの式になる。したがって、比例粘性減衰を仮定した減衰振動の固有モードは、同一の質量行列と剛性行列を有する不減衰振動の固有モードと同じである。
一方、比例粘性減衰を仮定した減衰振動の固有角振動数は、同一の質量行列と剛性行列の不減衰振動の固有角振動数よりも小さくなる。モード行列 U で減衰行列を対角化すると
U
⊤
C
U
=
(
α
M
1
+
β
K
1
0
⋯
0
0
α
M
2
+
β
K
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
α
M
n
+
β
K
n
)
{\displaystyle {\boldsymbol {U}}^{\top }{\boldsymbol {C}}{\boldsymbol {U}}={\begin{pmatrix}\alpha M_{1}+\beta K_{1}&0&\cdots &0\\0&\alpha M_{2}+\beta K_{2}&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &\alpha M_{n}+\beta K_{n}\end{pmatrix}}}
(4.13 )
となる。この行列の成分 αMr + βKr (r = 1, 2, …, n ) をモード減衰 やモード減衰係数 と呼び、 Cr などで表す。ここで、Mr はモード質量、Kr はモード剛性である。さらに、
C
c
,
r
=
2
M
r
K
r
{\displaystyle C_{c,r}=2{\sqrt {M_{r}K_{r}}}}
(4.14 )
という量を導入して、これでモード減衰を割った量
ζ
r
=
C
r
C
c
,
r
=
α
M
r
+
β
K
r
2
M
r
K
r
{\displaystyle \zeta _{r}={\frac {C_{r}}{C_{c,r}}}={\frac {\alpha M_{r}+\beta K_{r}}{2{\sqrt {M_{r}K_{r}}}}}}
(4.15 )
をモード減衰比 と呼ぶ。多自由度系の減衰系では固有モードごとに減衰の効果が異なっており、モード減衰比が固有モードごとの減衰の効果の程度を表している[ 137] 。ζr ≥ 1 ならば過減衰の状態であり、その固有モードの振動は起こらない[ 137] 。ζr < 1 ならば減衰振動となり、その固有角振動数は
ω
d
,
r
=
ω
r
1
−
ζ
r
2
{\displaystyle \omega _{d,r}=\omega _{r}{\sqrt {1-\zeta _{r}^{2}}}}
(4.16 )
で与えられる。ωd,r を減衰固有角振動数と呼ぶ。以上のように、線形1自由度系の減衰振動の考え方が固有モードごとの振動にも当てはまる。正規座標へ変換を行い、モード減衰比と固有角振動数を用いて運動方程式を表すと、下記のように表現できる。
{
q
¨
1
+
2
ζ
1
ω
1
q
˙
1
+
ω
1
2
q
1
=
0
q
¨
2
+
2
ζ
2
ω
2
q
˙
2
+
ω
2
2
q
2
=
0
⋮
q
¨
n
+
2
ζ
r
ω
r
q
˙
r
+
ω
n
2
q
n
=
0
{\displaystyle {\begin{cases}{\ddot {q}}_{1}+2\zeta _{1}\omega _{1}{\dot {q}}_{1}+\omega _{1}^{2}q_{1}=0\\{\ddot {q}}_{2}+2\zeta _{2}\omega _{2}{\dot {q}}_{2}+\omega _{2}^{2}q_{2}=0\\\vdots \\{\ddot {q}}_{n}+2\zeta _{r}\omega _{r}{\dot {q}}_{r}+\omega _{n}^{2}q_{n}=0\\\end{cases}}}
(4.17 )
ζr < 1 であれば、qr の一般解は積分定数を ar と br として、
q
r
=
e
−
ζ
r
ω
r
t
(
a
r
cos
ω
d
,
r
t
+
b
r
sin
ω
d
,
r
t
)
{\displaystyle q_{r}=e^{-\zeta _{r}\omega _{r}t}(a_{r}\cos \omega _{d,r}t+b_{r}\sin \omega _{d,r}t)}
(4.18 )
となる。物理座標 x へ逆変換すれば、比例粘性減衰系の自由振動の一般解は下記のようになる。
x
=
∑
r
=
1
n
u
¯
r
e
−
ζ
r
ω
r
t
(
a
r
cos
ω
d
,
r
t
+
b
r
sin
ω
d
,
r
t
)
{\displaystyle {\boldsymbol {x}}=\sum _{r=1}^{n}{\boldsymbol {\overline {u}}}_{r}e^{-\zeta _{r}\omega _{r}t}(a_{r}\cos \omega _{d,r}t+b_{r}\sin \omega _{d,r}t)}
(4.19 )
強制振動
不減衰系
任意の励振力を受ける減衰の無い系の運動方程式は以下のように表される。
M
x
¨
+
K
x
=
f
{\displaystyle {\boldsymbol {M{\ddot {x}}}}+{\boldsymbol {Kx}}={\boldsymbol {f}}}
(5.1 )
例えば、励振力が角振動数 Ω の余弦関数で与えられるとすれば、
M
x
¨
+
K
x
=
f
a
cos
Ω
t
{\displaystyle {\boldsymbol {M{\ddot {x}}}}+{\boldsymbol {Kx}}={\boldsymbol {f_{a}}}\cos {\mathit {\Omega }}t}
(5.2 )
となる。ここで、fa は、以下のような各自由度に対する励振力の振幅値の縦ベクトルである。
f
a
=
(
f
a
,
1
f
a
,
2
⋮
f
a
,
n
)
{\displaystyle {\boldsymbol {f_{a}}}={\begin{pmatrix}f_{a,1}\\f_{a,2}\\\vdots \\f_{a,n}\end{pmatrix}}}
式5.2 に対して特解を x = u cosΩt と仮定し、式5.2 に代入すれば
(
M
−
Ω
2
K
)
u
=
f
a
{\displaystyle ({\boldsymbol {M}}-{\mathit {\Omega }}^{2}{\boldsymbol {K}}){\boldsymbol {u}}={\boldsymbol {f_{a}}}}
(5.3 )
となるので、u および x は、
u
=
(
M
−
Ω
2
K
)
−
1
f
a
{\displaystyle {\boldsymbol {u}}=({\boldsymbol {M}}-{\mathit {\Omega }}^{2}{\boldsymbol {K}})^{-1}{\boldsymbol {f_{a}}}}
(5.4 )
x
=
(
M
−
Ω
2
K
)
−
1
f
a
cos
Ω
t
{\displaystyle {\boldsymbol {x}}=({\boldsymbol {M}}-{\mathit {\Omega }}^{2}{\boldsymbol {K}})^{-1}{\boldsymbol {f_{a}}}\cos {\mathit {\Omega }}t}
(5.5 )
となる。ここで "−1 " は逆行列 を意味する。したがって、逆行列 (M − Ω 2 K )−1 を計算すれば式5.5 から x の値が分かる。しかし、この逆行列を解析的に解くことは困難で、数値計算を行うにしても自由度の数が増えると膨大な計算量になる。もし励振力の角振動数 Ω を変えると、そのたびに逆行列を計算する必要がある。そのため、実際に x の解を求めるために行われるのは、下記のようなモード解析による手法である。
解を得るために、x を複素数に拡張し、励振力 f を複素指数関数 fa e jΩt の形で与えるとする。この場合、計算して解が得られた後に実部あるいは虚部を取ることで、実際の解が得られる。運動方程式5.1 の右辺を 0 としたときのモード行列 U が、事前に求められているとする。モード座標への変換式3.21 を運動方程式5.1 へ適用して、左から U ⊤ を掛け、式3.14 と式3.15 の対角化を適用する。すると、励振力を受ける不減衰系の運動方程式は
{
M
1
q
¨
1
+
K
1
q
1
=
F
1
e
j
ω
t
M
2
q
¨
2
+
K
2
q
2
=
F
2
e
j
ω
t
⋮
M
n
q
¨
n
+
K
n
q
n
=
F
n
e
j
ω
t
{\displaystyle {\begin{cases}M_{1}{\ddot {q}}_{1}+K_{1}q_{1}=F_{1}e^{j\omega t}\\M_{2}{\ddot {q}}_{2}+K_{2}q_{2}=F_{2}e^{j\omega t}\\\vdots \\M_{n}{\ddot {q}}_{n}+K_{n}q_{n}=F_{n}e^{j\omega t}\\\end{cases}}}
(5.6 )
という独立・非連成の n 個の運動方程式に帰着する。ただし、右辺の Fr は次のような値である。
F
r
=
u
r
⊤
f
a
{\displaystyle F_{r}={\boldsymbol {u}}_{r}^{\top }{\boldsymbol {f_{a}}}}
(5.7 )
式5.6 は線形1自由度系と同じなので、強制振動 を表す特解は、
q
r
=
F
r
K
r
−
M
r
Ω
2
e
j
Ω
t
{\displaystyle q_{r}={\frac {F_{r}}{K_{r}-M_{r}{\mathit {\Omega }}^{2}}}e^{j{\mathit {\Omega }}t}}
(5.8 )
となる。式3.21 を使ってモード座標を物理座標へ逆変換し、固有角振動数を使って整理すれば強制振動の解は次のようになる。
x
=
∑
r
=
1
n
F
r
K
r
−
M
r
Ω
2
u
r
e
j
Ω
t
=
∑
r
=
1
n
F
r
M
r
(
ω
r
2
−
Ω
2
)
u
r
e
j
Ω
t
{\displaystyle {\boldsymbol {x}}=\sum _{r=1}^{n}{\frac {F_{r}}{K_{r}-M_{r}{\mathit {\Omega }}^{2}}}{\boldsymbol {u}}_{r}e^{j{\mathit {\Omega }}t}=\sum _{r=1}^{n}{\frac {F_{r}}{M_{r}(\omega _{r}^{2}-{\mathit {\Omega }}^{2})}}{\boldsymbol {u}}_{r}e^{j{\mathit {\Omega }}t}}
(5.9 )
式5.9 から分かることは、励振力の角振動数 Ω が系の固有角振動数 ω 1 , ω 2 , …, ωn のどれかに近いとき、係数の極限 が 0 となってその振動成分が極めて大きくなり、共振 が起こるという点である。Ω が固有角振動数のいずれかに一致するとき、x の振幅は無限大に発散する。
比例粘性減衰系
系が比例粘性減衰系であれば、励振力を受ける場合でもモード座標変換によって独立した1自由度系に分解でき、大自由度の問題も扱いが容易になる。一般の運動方程式2.1 に、レイリー減衰の式4.8 およびモード座標への変換式3.21 を代入する。左から転置した質量正規固有モードの正規モード行列 Ū ⊤ を掛けて対角化する。さらに、2ζr ωr = α + βωr 2 という関係を用いれば、運動方程式は下記のような独立・非連成の n 個の方程式に帰着する。
{
q
¨
1
+
2
ζ
1
ω
1
q
˙
1
+
ω
1
2
q
1
=
u
¯
1
⊤
f
q
¨
2
+
2
ζ
2
ω
2
q
˙
2
+
ω
2
2
q
2
=
u
¯
2
⊤
f
⋮
q
¨
n
+
2
ζ
n
ω
n
q
˙
n
+
ω
n
2
q
n
=
u
¯
n
⊤
f
{\displaystyle {\begin{cases}{\ddot {q}}_{1}+2\zeta _{1}\omega _{1}{\dot {q}}_{1}+\omega _{1}^{2}q_{1}={\boldsymbol {\overline {u}}}_{1}^{\top }{\boldsymbol {f}}\\{\ddot {q}}_{2}+2\zeta _{2}\omega _{2}{\dot {q}}_{2}+\omega _{2}^{2}q_{2}={\boldsymbol {\overline {u}}}_{2}^{\top }{\boldsymbol {f}}\\\vdots \\{\ddot {q}}_{n}+2\zeta _{n}\omega _{n}{\dot {q}}_{n}+\omega _{n}^{2}q_{n}={\boldsymbol {\overline {u}}}_{n}^{\top }{\boldsymbol {f}}\\\end{cases}}}
(5.10 )
励振力 f が fa cosΩt で与えられるとする。線形1自由度系と同様の解法によって、正規座標 qr の強制振動の解が下記のように求まる。
q
r
=
u
¯
r
⊤
f
a
ω
r
2
{
1
−
(
Ω
ω
r
)
2
}
2
+
{
2
ζ
r
Ω
ω
r
}
2
cos
(
Ω
t
−
ϕ
r
)
{\displaystyle q_{r}={\dfrac {{\boldsymbol {\overline {u}}}_{r}^{\top }{\boldsymbol {f_{a}}}}{\omega _{r}^{2}{\sqrt {\left\{1-\left({\dfrac {\mathit {\Omega }}{\omega _{r}}}\right)^{2}\right\}^{2}+\left\{2\zeta _{r}{\dfrac {\mathit {\Omega }}{\omega _{r}}}\right\}^{2}}}}}\cos({\mathit {\Omega }}t-\phi _{r})}
(5.11 )
ϕ
r
=
tan
−
1
(
2
ζ
r
ω
r
Ω
ω
r
2
−
Ω
2
)
{\displaystyle \phi _{r}=\tan ^{-1}\left({\frac {2\zeta _{r}\omega _{r}{\mathit {\Omega }}}{\omega _{r}^{2}-{\mathit {\Omega }}^{2}}}\right)}
(5.12 )
解を複素数で表現した場合は、
q
r
=
u
¯
r
⊤
f
a
ω
r
2
{
1
−
(
ω
r
Ω
)
2
+
j
2
ζ
r
ω
r
Ω
}
e
j
Ω
t
{\displaystyle q_{r}={\dfrac {{\boldsymbol {\overline {u}}}_{r}^{\top }{\boldsymbol {f_{a}}}}{\omega _{r}^{2}\left\{1-\left({\dfrac {\omega _{r}}{\mathit {\Omega }}}\right)^{2}+j2\zeta _{r}{\dfrac {\omega _{r}}{\mathit {\Omega }}}\right\}}}e^{j{\mathit {\Omega }}t}}
(5.13 )
となり、物理座標の解は
x
=
∑
r
=
1
n
u
¯
r
⊤
f
a
u
¯
r
ω
r
2
{
1
−
(
Ω
ω
r
)
2
+
j
2
ζ
r
Ω
ω
r
}
e
j
Ω
t
{\displaystyle {\boldsymbol {x}}=\sum _{r=1}^{n}{\frac {\boldsymbol {{\overline {u}}_{r}^{\top }f_{a}{\overline {u}}_{r}}}{\omega _{r}^{2}\left\{1-\left({\dfrac {\mathit {\Omega }}{\omega _{r}}}\right)^{2}+j2\zeta _{r}{\dfrac {\mathit {\Omega }}{\omega _{r}}}\right\}}}e^{j{\mathit {\Omega }}t}}
(5.14 )
となる。Ω /ωr は強制振動比などと呼ばれる。通常の固有モードの場合は次式となる。
x
=
∑
r
=
1
n
u
r
⊤
f
a
u
r
K
r
{
1
−
(
Ω
ω
r
)
2
+
j
2
ζ
r
Ω
ω
r
}
e
j
Ω
t
{\displaystyle {\boldsymbol {x}}=\sum _{r=1}^{n}{\frac {\boldsymbol {u_{r}^{\top }f_{a}u_{r}}}{K_{r}\left\{1-\left({\dfrac {\mathit {\Omega }}{\omega _{r}}}\right)^{2}+j2\zeta _{r}{\dfrac {\mathit {\Omega }}{\omega _{r}}}\right\}}}e^{j{\mathit {\Omega }}t}}
(5.15 )
周波数応答関数
強制振動の特解から、周波数応答関数 あるいは伝達関数 が求められる。周波数応答関数は、振動系への入力が調和振動の形で与えられるときの、出力と入力の振幅比および位相差を周波数の関数として表す。比例減衰系において、角振動数 Ω の調和外力が q 番目の自由度のみに加わるとする。このときの p 番目の自由度の振動応答は
G
p
q
(
Ω
)
=
∑
r
=
1
n
u
p
r
u
q
r
K
r
{
1
−
(
Ω
ω
r
)
2
+
j
2
ζ
r
Ω
ω
r
}
{\displaystyle G_{pq}(\Omega )=\sum _{r=1}^{n}{\frac {u_{pr}u_{qr}}{K_{r}\left\{1-\left({\dfrac {\mathit {\Omega }}{\omega _{r}}}\right)^{2}+j2\zeta _{r}{\dfrac {\mathit {\Omega }}{\omega _{r}}}\right\}}}}
(5.16 )
と表される。この場合の周波数応答関数 G (Ω ) は[変位]/[力]の単位を持ち、コンプライアンス とも呼ばれる。q 番目加振・p 番目応答の周波数応答関数と p 番目加振・q 番目応答の周波数応答関数は、相反定理 によって互いに等しく、Gpq (Ω ) = Gqp (Ω ) の関係である。不減衰系の場合は、周波数応答関数は次式となる。
G
p
q
(
ω
)
=
∑
r
=
1
n
u
p
r
u
q
r
K
r
−
M
r
Ω
2
{\displaystyle G_{pq}(\omega )=\sum _{r=1}^{n}{\frac {u_{pr}u_{qr}}{K_{r}-M_{r}{\mathit {\Omega }}^{2}}}}
(5.17 )
式4.2 ,4.3 ,4.4 によって与えられた3自由度減衰系の周波数応答関数の応答曲線の例。図は質点 m 3 が加振されたときの質量 m 2 の応答の例。山になっている点が共振点で、谷になっている点が反共振点。
これらの周波数応答関数によって、系の周波数特性が把握できる。線形多自由度系の場合は、固有角振動数と固有モードを求めれば、調和外力による強制振動に対する周波数特性も同時に理解することができる。横軸に振動数や強制振動比を取り、縦軸に調和外力に対する振幅比や位相差を図示したものを応答曲線や共振曲線と呼ぶ。縦軸をデシベル にして図示したものは特にボード線図 と呼ばれる。系の振動特性を振動実験から同定する実験モード解析では、実験値に上記の周波数応答関数でカーブフィッティング して系の各種パラメータを推定する。
上述のように、励振力の振動数が系の固有角振動数に一致すると、振動は極大化する。一方、多自由度系では共振して振幅が発散する現象だけでなく、励振力が特定の角振動数のときに周波数応答関数が極小になる現象もある。このような現象を反共振 と呼ぶ。ただし、共振が系全体で起こるのに対して、反共振は一組の加振点・応答点ごとにしか起こらない。周波数応答関数の絶対値を縦軸にした応答曲線上では、共振点は曲線の鋭い山のように現れ、反共振は曲線の鋭い谷のように現れる[ 166] 。反共振は、代表的な制振器 である動吸振器 の原理として活用される。振動を抑えたい対象にばね・ダンパを介して付加質量を取り付けることによって振動を抑制でき、地震や強風に対する建築構造物の防振や回転機械の防振などに使われる。
刺激係数
刺激係数 は、地震 のように構造物の基礎が揺れている場合に、その励振が各モードに対してどのぐらい強く寄与するのかを表す。構造物を多質点系でモデル化して、基礎の変位を y として、各質点の基礎からの相対変位を x とする。各質点は ẍi + ÿ の加速度を受けるため、運動方程式は
M
x
¨
+
C
x
˙
+
K
x
=
−
M
1
y
¨
{\displaystyle {\boldsymbol {M{\ddot {x}}}}+{\boldsymbol {C{\dot {x}}}}+{\boldsymbol {Kx}}=-{\boldsymbol {M1}}{\ddot {y}}}
(5.18 )
となる。ここで 1 は成分が全て 1 の縦ベクトルである。式5.18 において、C をレイリー減衰で表し、さらに x をモード座標に置き換え、左から転置したモード行列を掛けて整理すれば、下記のような各モード座標についての式になる。
{
q
¨
1
+
2
ζ
1
ω
1
q
˙
1
+
ω
1
2
q
1
=
−
u
1
⊤
M
1
M
1
y
¨
q
¨
2
+
2
ζ
2
ω
2
q
˙
2
+
ω
2
2
q
2
=
−
u
2
⊤
M
1
M
2
y
¨
⋮
q
¨
n
+
2
ζ
n
ω
n
q
˙
n
+
ω
n
2
q
n
=
−
u
n
⊤
M
1
M
n
y
¨
{\displaystyle {\begin{cases}{\ddot {q}}_{1}+2\zeta _{1}\omega _{1}{\dot {q}}_{1}+\omega _{1}^{2}q_{1}=-{\dfrac {\boldsymbol {u_{1}^{\top }M1}}{M_{1}}}{\ddot {y}}\\{\ddot {q}}_{2}+2\zeta _{2}\omega _{2}{\dot {q}}_{2}+\omega _{2}^{2}q_{2}=-{\dfrac {\boldsymbol {u_{2}^{\top }M1}}{M_{2}}}{\ddot {y}}\\\vdots \\{\ddot {q}}_{n}+2\zeta _{n}\omega _{n}{\dot {q}}_{n}+\omega _{n}^{2}q_{n}=-{\dfrac {\boldsymbol {u_{n}^{\top }M1}}{M_{n}}}{\ddot {y}}\\\end{cases}}}
(5.19 )
さらに、上式の ÿ の係数を次のように変形する。
u
r
⊤
M
1
M
r
=
u
r
⊤
M
1
u
r
⊤
M
u
r
=
β
r
{\displaystyle {\dfrac {\boldsymbol {u_{r}^{\top }M1}}{M_{r}}}={\dfrac {\boldsymbol {u_{r}^{\top }M1}}{\boldsymbol {u_{r}^{\top }Mu_{r}}}}=\beta _{r}}
(5.20 )
この βr を r 次の刺激係数と呼び、r 次モードの運動方程式は次のようになる。
q
¨
r
+
2
ζ
r
ω
r
q
˙
r
+
ω
r
2
q
n
=
−
β
r
y
¨
{\displaystyle {\ddot {q}}_{r}+2\zeta _{r}\omega _{r}{\dot {q}}_{r}+\omega _{r}^{2}q_{n}=-\beta _{r}{\ddot {y}}}
(5.21 )
つまり、βr は励振加速度が各モードに対してどのぐらい寄与しているか、あるいは各モードは基礎励振に対してどのぐらい影響を受けやすいかを表している。さらに、r 次の刺激係数と固有モードの積を r 次の刺激関数 と呼び、刺激関数の総和は次のように 1 に等しい。
β
1
u
1
+
β
2
u
2
+
⋯
+
β
n
u
n
=
1
{\displaystyle \beta _{1}{\boldsymbol {u}}_{1}+\beta _{2}{\boldsymbol {u}}_{2}+\ \cdots \ +\beta _{n}{\boldsymbol {u}}_{n}={\boldsymbol {1}}}
(5.22 )
すなわち、対称の振動系に 1 という外力ベクトルが加わったときに各モードが受ける度合いを、刺激関数は表している。
有限要素法による連続体の振動への応用
多自由度系では、質点または剛体から成る系を想定し、剛性、減衰、慣性などの振動特性も局所的に集中して系内で点在しているというモデルを考えていた。一方で、実際の物体は、物体自体が変形する。実際の物体は連続体 としての性質を有しており、質量、剛性などは連続的に系内に分布しているモデルとなる。実際の問題では、多自由度系に近似して取り扱っても十分な場合も多いが、振動時の機械・構造物の各部の変形や応力といったものを知るには連続体として取り扱う必要ある。
連続体の振動の運動方程式は時間と空間に関する偏微分方程式 で記述され、厳密解が得られることは限られる。実際の複雑な形状の構造物で連続体の振動を扱うには、実用的には有限要素法 という手法が用いられる。有限要素法では、対象の連続体を小さな有限要素に分割し、連続体を多自由度系に置き換えて解を計算する。有限要素法においても、線形多自由度系の理論にもとづくモード解析手法が強力な効果を発揮し、振動解析が精度良くかつ容易にできるようになる。
有限要素法による振動解析では、立てられた振動数方程式の数値計算を行い、まず固有振動数と固有モードを得る[ 115] 。次いで、固有モードの直交性を利用して周波数応答関数を得て、応答解析を行う[ 115] 。もし、定式化された運動方程式が1万自由度だとしたら、解くべき方程式は1万元の2次連立微分方程式となり、コンピュータを用いても計算に長時間を要する。しかし、モード解析手法を用いれば、1自由度系の解を1万回解いて、重ね合わせるだけで解が得られる。さらに、対象が大規模自由度になったとしても、自由度の分だけ現れるモードを全て計算する必要性もない。実用的に興味のある外力振動数を含む次数までモードの重ね合わせでも、十分な精度の振動応答解析が可能となる。上記の1万自由度の例えで言えば、1自由度系の解を1万回解く必要もなく、もっと少ない回数の計算で事足りるようになる。これらの長所によって、モード解析手法は有限要素法による振動解析で絶大な威力を発揮し、数十万規模の自由度を扱うような有限要素法計算であっても、モード解析手法の適用によって特段の支障なく計算が可能となる。
係数行列が非対称行列の場合
質量行列 M (式2.2 )、減衰行列 C (式2.3 )、あるいは剛性行列 K (式2.4 )が正定値 の条件を満たさない場合、すなわち実対称行列ではなく、非対称行列であるとき、その系では不安定振動が起こることがある。このような条件では、式4.6 で表される特性方程式の固有値 λ に、実部が正の固有値が含まれることがありえる。固有値に実部を正とする複素数が含まれるとき、時間とともに振幅が大きくなっていく振動が起こる[ 120] 。このようなメカニズムは、自励振動 が起こりえる系で平衡点から振動が成長するか否かを考察する上で基礎となる。自励振動は1自由度系でも起きる現象だが、係数行列が非対称であることによって引きこされる種類の自励振動は、多自由度系特有のものである。
例として、次のような2自由度不減衰系を考える[ 188] 。
(
m
1
0
0
m
2
)
(
x
¨
1
x
¨
2
)
+
(
k
11
k
12
k
21
k
22
)
(
x
1
x
2
)
=
(
0
0
)
{\displaystyle {\begin{pmatrix}m_{1}&0\\0&m_{2}\end{pmatrix}}{\begin{pmatrix}{\ddot {x}}_{1}\\{\ddot {x}}_{2}\end{pmatrix}}+{\begin{pmatrix}k_{11}&k_{12}\\k_{21}&k_{22}\end{pmatrix}}{\begin{pmatrix}x_{1}\\x_{2}\end{pmatrix}}={\begin{pmatrix}0\\0\end{pmatrix}}}
(7.1 )
ただし、k 12 ≠ k 21 で、剛性行列は非対称行列である[ 188] 。さらに、k 12 と k 21 のどちらかが正でどちらかが負であるような異符号の関係にあるとき、固有値は
λ
=
±
λ
r
±
i
λ
i
{\displaystyle \lambda =\pm \lambda _{r}\pm i\lambda _{i}}
(7.2 )
という形の複素数となる[ 188] 。λr とλi は
λ
r
=
1
2
−
ω
t
r
+
ω
t
r
2
−
ω
d
i
f
f
2
−
4
ω
s
k
{\displaystyle \lambda _{r}={\frac {1}{2}}{\sqrt {-\omega _{tr}+{\sqrt {\omega _{tr}^{2}-\omega _{diff}^{2}-4\omega _{sk}}}}}}
(7.3 )
λ
i
=
1
2
ω
t
r
+
ω
t
r
2
−
ω
d
i
f
f
2
−
4
ω
s
k
{\displaystyle \lambda _{i}={\frac {1}{2}}{\sqrt {\omega _{tr}+{\sqrt {\omega _{tr}^{2}-\omega _{diff}^{2}-4\omega _{sk}}}}}}
(7.4 )
で与えられ、ここで、ωtr = k 11 /m 1 + k 22 /m 2 , ωdiff = k 11 /m 1 − k 22 /m 2 , ωsk = (k 12 /m 1 )(k 21 /m 2 ) である[ 188] 。λr は発散率と呼ばれ、自励振動の強さを表す[ 188] 。
このような係数行列の非対称性によって起きる自励振動の事例は機械振動の中で多く見られ、クーロン摩擦 による摩擦振動や滑り軸受 で起こるオイルホイップなどがある。
出典
参照文献