振動問題をFreeMatで解いてみよう（その2）：無償ソフトで技術計算しよう【シミュレーション応用編】（3）（2/3 ページ）

[伊藤孝宏，MONOist]

減衰速度

　機械設計では、動いた後に速やかに停止してほしい場合があります。そこで、運動方程式のパラメータm、c、kを適当に変化させて、減衰の様子を見てみます。

　分かりやすく、外力を与えずに、初期値を変位1（m）とします。また、パラメータを任意に変えられるように関数mファイルでプログラムを作成します。コマンドウィンドウ上でex407（1,1,10）などと入力し、subplotコマンドで複数の結果をプロットすると、図4に示す減衰波形が得られます。

図4：パラメータの違いによる減衰の変化

　図を見ると、減衰係数cを大きくすると、速く減衰することが分かりますが、質量mやバネ定数kを小さくしても、減衰が速くなることが分かります。

　振動の減衰の速さは減衰比ζで表され、

　式から、減衰比は減衰係数に比例する他、質量mやバネ定数kのルートに反比例することが分かります。

　ただ、質量mは分かったとしても、バネ定数kや減衰係数cは分からないことが多くあります。そこで、停止時の減衰の様子が図5のようになったとして、これからバネ定数kと減衰係数cを計算してみます。

図5：減衰振動波形

　実際は測定結果のグラフから読み取ることになりますが、FreeMatで作成した図なので、ディスプレイ上にFigureウィンドウがでているうちに、コマンドウィンドウでt=pointと入力し、カーソルとFigureウィンドウに移動します。すると、カーソルが十字に変わるので、適当な山の位置でクリックすると、コマンドウィンドウ上に、時間t、変位yの順に数値が表示されます。同様に次の山でも繰り返し、2つの山の位置を測定します。

　2つの山の時間情報からバネ定数を求めます。

　質量mは測定済みで1（kg）として、2つの山の時間情報から、周期T＝4−2＝2（s）となり、バネ定数はk＝9.9（N/m）となります。これは、モデルで設定したk＝10（N/m）とほぼ一致します。

　次に、減衰係数を求めます。始めに、図5から対数減衰率を求めます。

　山と次の山あるいは、谷と次の谷の振幅の比の自然対数となります。

　次に、減衰比ζと対数減衰率δとは、

　対数減衰率δから、減衰係数cを求めることができます。では、減衰係数cを求めてみましょう。

　2つの山のy座標は、0.36と0.14であったので、

となります。一方、質量mは計測済みで1（kg）で、図5から計算したバネ定数kは9.9（N/m）でしたので、

　減衰係数cは、

となり、モデルで設定した減衰係数c＝1（Ns/m）とほぼ一致します。このように、実験結果から、バネ定数k、減衰係数cが推測できれば、減衰振動モデルから、各種の対策の効果を数値解析により求めることが可能になります。

　実は、説明の便宜上、区別をしていませんでしたが、上記で求めた振動の周期は、共振周波数であり、減衰がある場合、固有振動数とは異なります。従って、振動の周期から求めたバネ定数は、固有振動数から求めた値とは異なり、正確な値ではありません。では、どの程度異なるかを見積もってみましょう。共振周波数fsと固有振動数fとには、減衰比ζとの間に

　従って、バネ定数を求めるための振動数、すなわち固有振動数は、

となります。

　波形から求めた対数減衰率δ=0.94から減衰比ζを求めると、ζ＝δ/2π=0.15となり、固有振動数は、

　となります。

　すなわち、波形から求めた周波数はおおむね固有振動数と等しく、求められたバネ定数は、おおむね正しい値と考えられます。実際、求められたバネ定数は9.9（N/m）とモデルに設定された値10（N/m）とほぼ同じ結果となっています。

function ex407(m,c,k)
    dydt=@(t,y) [y(2);(-c*y(2)-k*y(1))/m];
    [t,y]=ode45(dydt,[0,10],[1;0]);
    plot(t,y(:,1));grid('on');
    title(['m=',num2str(m),' c=',num2str(c),' k=',num2str(k)]);
    ylabel('Displacement[m]');

ex407.m

＞＞「ex407.m」ダウンロード