19 drawnow; //画面更新
20 xxmax(i)=max(xx(1, 250:300)); //xx(250,1)〜xx(300,1)の最大値
21 end
22 subplot(2,1,2); plot(om1, xxmax, "o");
23 xlabel("om"); ylabel("max x(t)");
24 xtitle("Response Curve");
25 save("vib_res.dat","om1","xxmax"); // vib-res2.sceで使うデータの保存
Code 3を実行すると，次のような結果が得られる．なお，最後に表示されるグラ
フを，共振曲線(resonance curve)というが，これは14.3.4節p142で使う．
-4 -2 0 2 4 6 8
0 20 40 60 80 100
x(t)
t Waveform (om = 0.3)
f(t) x(t)
-4 -2 0 2 4 6 8
0 20 40 60 80 100
x(t)
t Waveform (om = 1.0)
f(t) x(t)
-4 -2 0 2 4 6 8
0 20 40 60 80 100
x(t)
t Waveform (om = 1.5)
f(t) x(t)
まず，(1)波形の細かさだが，入力f(t) =Pcosωtの振動数ωを増加させている のだから，当然，f(t)の山と谷の間隔は狭くなる．x(t)の間隔も似たようなもので，
見た感じ，x(t)とf(t)は同じ振動数に見える．次に，(2)応答x(t)の振幅Rの増 減だが，単調ではない．いったん5倍程度に逹してから，減少に転じている．最後 に，(3)入力f(t)に対するx(t)の遅れϕだが，ωが小さいとき(om=0.3)は，f(t) とx(t)は重なっている．ところが，ωが増えると，f(t)に対してx(t)が遅れ始める (x(t)が相対的に右にずれる)．さらにωを増やすと(om=1.5)，f(t)に対するx(t) の遅れは，半周期程度に逹する(f(t)の山とx(t)の谷が重なる)．
以上の観察をグラフにすると，だいたい次のようになる．この2枚のグラフを，周
図14.4 周波数応答のグラフ(概形)
波数応答(frequency response)またはボード線図(Bode diagram)という．共振現 象の特徴は，この2枚のグラフによってもれなく説明できる．
14.3. ハーモニックバランス 141
14.3.1 解析モデルの簡略化
調和入力を受ける強制振動系m¨x+cx˙+kx=Pcosωtは，章末A14節p147の 方法を使うと，振動波形を相似に保ったまま，次のように変形できる．
¨
x+ 2ζx˙+x=F(t) =AcosΩt, A= P mωn2 =P
k, Ω= ω ωn
(14.12) ζとωnは表13.2p129の減衰比と固有振動数，Ωは入力振動数ωの固有振動数ωn
の比である．以下，簡略化した(14.12)を使って，周波数応答を計算する．
14.3.2 定常応答の数式表現
実習14.2p139の(1)でも観察したが，一般に，次の物理法則が成立する15)．
力学法則14.1 (調和応答の性質) 入力が三角関数ならば，定常応答も三角関数16)．振
動数も同じ．
したがって，入力AcosΩtを受ける強制振動系(14.12)の定常応答β(t)は，振動 数Ωの三角関数となる．応答振幅Rと位相差ϕを含めて，式で書くと，
β(t) =Rcos(Ωt+ϕ) (14.13)
となる．ここで，ハーモニックバランスを楽にする工夫として，(14.13)を加法定理で，
β(t) = (Rcosϕ
| {z }
a1
) cosΩt+ (−Rsinϕ
| {z }
a2
) sinΩt
と展開した，
β(t) =a1cosΩt+a2sinΩt (14.14) を使うことにする．元に戻すには，
a21+a22=R2(cos2ϕ+ sin2ϕ) =R2 ∴R=√ a21+a22,
−a2
a1
= Rsinϕ
Rcosϕ = tanϕ ∴ϕ= tan−1−a2
a1
=−tan−1 a2
a1
(14.15) とすればよい．tan−1 はtanの逆関数である17)．
14.3.3 ハーモニックバランス法
(14.14)で仮定したβ(t) =a1cosΩt+a2sinΩtが，強制振動系(14.12)の解にな るように，a1, a2 を調整してやれば，定常応答β(t)の数式表現が求まる．
(14.14)のβ(t)とその時間微分，
15)この法則は数学の定理としても導ける．
17)正確には，(2.4)p13の atan を使ってϕ=−atan (a1, a2)である．
β(t) =˙ −a1ΩsinΩt+a2ΩcosΩt, β(t) =¨ −a1Ω2cosΩt−a2Ω2sinΩt を(14.12)のx,x,˙ ¨xに代入すると，
(−a1Ω2cosΩt−a2Ω2sinΩt) + 2ζ(−a1ΩsinΩt+a2ΩcosΩt)
+ (a1cosΩt+a2sinΩt) =AcosΩt となる．これに，ハーモニックバランス法(三角関数釣合い法)を適用する．全部左辺 へ移項して，cosΩtとsinΩtをくくり出すと
(−a1Ω2+ 2ζa2Ω+a1−A) cosΩt + (−a2Ω2−2ζa1Ω+a2) sinΩt= 0 となる．この式が任意の時刻tで成立する(釣合う)には，cosΩtとsinΩtの係数が 常に0でなければならない．ゆえに，
−a1Ω2+ 2ζa2Ω+a1−A= 0, −a2Ω2−2ζa1Ω+a2= 0
という2元2連立方程式が導かれる．以上の操作をハーモニックバランスという．
未知数a1, a2を求めるため，ベクトル表示すると，
1−Ω2 2ζΩ
−2ζΩ 1−Ω2
a1
a2
=
A 0
より
a1
a2
= 1
∆
1−Ω2 −2ζΩ 2ζΩ 1−Ω2
A 0
= A
∆
1−Ω2 2ζΩ
となる．ただし∆= (1−Ω2)2+ (2ζΩ)2は行列式である．
(14.15)p141で元に戻すと，
R=√
a21+a22=
√
A2(1−Ω2)2+A2(2ζΩ)2
∆2 =
√A2∆
∆2 = A
√∆,
ϕ=−tan−1a2
a1
=−tan−1 A(2ζΩ)
A(1−Ω2) =−tan−1 2ζΩ
1−Ω2 (14.16) となる．これらが，図14.4p140でスケッチした応答振幅Rと，位相差ϕのグラフの 数式表現である．(14.16)のRとϕを用いて，¨x+ 2ζx˙+x=AcosΩtの定常応答は，
β(t) =Rcos(Ωt+ϕ) = A
√(1−Ω2)2+ (2ζΩ)2cos (
Ωt−tan−1 2ζΩ 1−Ω2
)
(14.17) のように数式表現できる．目標達成！
14.3.4 共振曲線
Code 3で最後に出てくるグラフ，
14.3. ハーモニックバランス 143
1 1.5 2 2.5 3 3.5 4 4.5 5
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
max x(t)
om Response Curve
を，共振曲線と呼んだ．作り方を，Code 3の抜粋で説明しよう．
om1 = linspace (0.2,1.6,15); //外力周波数0.2〜1.6を15等分
for i = 1:15 //i=1,2,…,15を順に繰り返し
om = om1(i); //i番めの外力周波数を代入
xx = ode(x0, 0, tt , model); //振動波形を求める :
xxmax(i)=max(1,xx(250:300)); //xx(1,250)〜xx(1,300)の最大値 end //i<15なら9行目に戻る
plot(om1, xxmax); //横を周波数，縦を振幅としてプロット
抜粋の3行目で，ωに0.2,0.3,· · ·,1.6が順に代入される．4行目で，各ωでの振動 波形が計算される．6行目で，振動波形が定常応答になった後半部分を抜き出し，その 最大値を振幅xxmaxのi番目に代入．8行目で，(振動数，振幅)をプロットしてい る．これが共振曲線の実験的な求め方である18)．
次に，(14.17)p142からβ(t)の振幅を求めると，
max
t β(t) = max
t Rcos(Ωt+ϕ) =R= A
√(1−Ω2)2+ (2ζΩ)2 (14.18)
となる．これが共振曲線の数式表現である．
実習14.3 (共振曲線) Code 4を実行し，Code 3で実験的に求めた振幅maxx(t)
と，(14.18)で数式表現された振幅Rを，重ねて比較せよ．
Code 4 “vib-res2.sce” (Scilab)
1 clear; clf();
2 load("vib_res.dat","om1","xxmax"); // vib-res.sceで保存したデータ
3 m=1; c=0.2; k=1; P=1;
4 zeta = c/( 2*sqrt(m*k) ); // 表3.1
5 omn = sqrt(k/m); // 表3.1
6 function y = K(Om)
7 global zeta;
8 y=1/sqrt((1-Om^2)^2 +(2*zeta*Om)^2);
9 endfunction
10 om2 = linspace(0.2,1.6,100);
11 A = P/(m*omn^2);
12 for i=1:100
13 R(i) = A*K(om2(i)/omn);
14 end
15 plot(om1, xxmax,"o", om2, R,"-" );
16 xlabel("om"); ylabel("Amplitude");
17 xtitle("Response Curve (o max(x); - R)"); xgrid();
18 g=gca(); g.data_bounds=[0.2,0;1.6,6]; xgrid(); //座標軸の設定
Code 4を実行すると，図14.5p144の結果が得られる．
18)この例は，数値実験だが．
0 1 2 3 4 5 6 7
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
Amplitude
om Response Curve
max x(t) R
図14.5 共振曲線