MATLABで学ぶ振動工学 減衰系の強制振動

はじめに

振動工学やモード解析をご存じの方であれば、一度は読んだことがあると言われているバイブル的な書籍「モード解析入門」について再復習しております。今回はモード解析入門の2.5節の「減衰系の強制振動」のP.50から復習します。なお、振動工学の参考書はどの本も似たり寄ったりで内容が重複しています。例えば、「実用モード解析入門」は長松先生の息子さんが著者ということもあり、重複している部分が多いです。

各参考書で重複している部分は本記事で大まかに理解できると思いますので、他の参考書から本記事へアクセスした場合も最後まで読んでいただけると幸いです。

参考書解説一覧
MATLABとNVH(Noise, Vibration and Harshness)をこよなく愛する人生の先輩、それが私(MATLABパイセン)です。NVHを極めるにあたり、周辺の知識も網羅的に勉強することになり、その知識を共有すべく、本HPを運営しています。日本の製造業を応援すべく、機械エンジニアの「車輪の再開発」の防止と業務効率化の手助けをモットーに活動。専門である振動騒音工学や音響工学について、プログラムを示しながら解説。大学の授業よりもわかりやすく説明することを目指す。質問頂ければ、回答したいと思います。

応答

減衰系の強制振動の運動方程式は式(2.87)に示します。

$$
m\ddot{x(t)} + c\dot{x(t)} +kx(t) = F e^{jωt} \tag{2.87}
$$

式(2.87)を満たす解は式(2.76)です。

$$
x(t) =X e^{jωt} \tag{2.76}
$$

式(2.76)を式(2.87)に代入して、整理すると式(2.89)を得られる。

$$
(-ω^2 \frac{m}{k}+jω \frac{c}{k} +1) X =\frac{F}{k} \tag{2.89}
$$

振動工学ではおなじみですが、下記に置き換えます。Ω=√(k/m)は減衰のない角共振振動数ですよね。なので、Ωを不減衰固有角振動数と呼びます。

\[
\left\{
\begin{array}{l}
c_c=2 \sqrt{mk}=2m\sqrt{\frac{k}{m}}=2mΩ\\
ζ=c/c_c
\end{array}
\right.
\tag{2.47}
\]

式(2.47)を式(2.49)に代入します。

$$
(-(\frac{ω}{Ω})^2 +2jζ\frac{ω}{Ω} +1) X =\frac{F}{k} = X_{st} \tag{2.91}
$$

ここで、下記のβというものを導入する。

$$
β=\frac{ω}{Ω}
$$

ここからは、振動工学ではおなじみの式になります。

$$
\frac{X}{X_{st}} = \frac{1}{ 1 – β^2 + 2jζβ } \tag{2.92}
$$

Xstって何?

Xstは静的な変位です

静的な変位って………?

フックの法則はF=kXですよね?
簡単に説明するために、これを静的な変位と考えましょう
この時の変位をXstとすると、

Xst=F/kですよね

式(2.92)は下記を意味します。
(動的な力が作用したときの変位)/(静的な力が作用したときの変位)

つまり、静的な変位に対して、動的な変位がどの程度の倍率なのかを直感的に把握することができます。

さっそく、式(2.92)を計算してみましょう。
結果は下図です。

上記の結果からわかることは下記です。
・角周波数ωが小さいときは、静的な変位と動的な変位は等しい。
・共振周波数Ωよりも少し小さい角周波数ωのときに、応答が最大となる。
・共振周波数Ωと角周波数ωが一致するときが最大応答ではない。
・共振周波数Ωと角周波数ωが一致するときに位相角は-90°になる。
・共振周波数Ωよりも角周波数ωがはるかに大きくなると、応答は0に漸近する。
・共振周波数Ωよりも角周波数ωがはるかに大きくなると、応答が遅れ位相角が-180°になる。

 

簡単な式なので、MATLABプログラムは不要かもしれませんが一応載せておきます。

clear all;clc;close all;


k=1000; %ばね定数
m=1; %質量
% % c=1; %減衰c
% % cc=2*sqrt(m*k); %臨界減衰係数
Omega=sqrt(k/m); %固有値

freq=0.001:0.001:10; %周波数
w=2*pi*freq; %角周波数

beta=w/Omega;
% zeta=c/cc;
zeta=0:0.1:1;
j=sqrt(-1);

RGB=jet(length(zeta));

figure
for ii1=1:length(RGB)
X_Xst=1./(1-beta.^2+2*j*zeta(ii1)*beta);
subplot(121)
plot(beta,abs(X_Xst),'color',RGB(ii1,:))
hold on
grid on
subplot(122)
plot(beta,angle(X_Xst),'color',RGB(ii1,:))
hold on
grid on
end
subplot(121)
legend(split(num2str(zeta)))
ylim([0 4])
xlabel('β=ω/Ω')
ylabel('|X|/X_{st}')
subplot(122)
legend(split(num2str(zeta)))
xlabel('β=ω/Ω')
ylabel('angle( |X|/X_{st} )')

 

今日の解説はここまで。

Good luck

 

コメント