連載
» 2015年03月25日 09時00分 公開

振動問題を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」ダウンロード

Copyright © ITmedia, Inc. All Rights Reserved.