MATLABで学ぶ振動工学 基礎加振による応答

はじめに

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

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

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

基礎加振とは

モード解析入門の解説をする前に「基礎加振」について説明します。

基礎加振とは下図の基礎(地面もしくは床)の部分が振動することで、質量mが加振されることを意味します。

図 1自由度系の振動モデル(斜線部分が基礎や床と呼ばれている部分)

基礎加振による応答

下図のように基礎が振動する場合を考えましょう。

ここで、絶対空間内の質量mの変位をx、振動する基礎から見た質量の変位をy、基礎の振動をxbとします。

$$
y = x-x_b \tag{2.110}
$$

質量の慣性力をfm、ばねの復元力をfk、粘性抵抗力をfcとすると式(2.111)の関係式で表せます。

$$
f_m = -m \ddot{x} , f_k = -ky , f_c = -c \dot{y}  \tag{2.111}
$$

この1自由度のバネマス系には、外力は作用していないので、力は釣り合っている状態になります。従って、式(2.112)が得られます。

$$
m \ddot{x} +c \dot{y} +ky =0 \tag{2.112}
$$

これまでの式で式(2.112)を整理すると式(2.113)が得られます。

$$
m \ddot{x} +c \dot{x} +kx =c \dot{x_b} +kx_b = B(k+jcω)e{jωt} \tag{2.113}
$$

以下省略

最終的に、xとyは下式で表せます。

$$
x = |X|e^{j(ωt+atan(cω/k)+atan(-2ζβ/(1-β^2)))} \tag{2.117}
$$

$$
y = \frac{Bβ^2}{(1-β^2+2jζβ)}e{jωt} \tag{2.120}
$$

参考文献のモード解析入門には、不減衰(c=0、ζ=0)のときなどの場合分けをして、考察をしています。

参考文献のような場合分けは、頭のいい人なら理解できるのかもしれません。
しかし、より簡単に理解するには、MATLABで式をコード化して、とりあえずパラメータスタディするのが良いと思っています。

なので、とりあえずMATLABでパラメータスタディしてみましょう。
まず、ζをいじってみましょう。

図 ζ=0~1 @31Hz

位相が少し変化したぐらいですかね。。。。

clear all;clc;close all;


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

% zeta=c/cc;
zeta=0:0.1:1;
j=sqrt(-1);
B=1;
RGB=jet(length(zeta));
t=0:0.0001:1;
freq=31; %周波数


figure
for ii1=1:length(RGB)
    w=2*pi*freq; %角周波数
    beta=w/Omega;
    y=B*beta.^2./(1-beta.^2 + 2*j.*zeta(ii1)*beta).*exp(j*w.*t);
    plot(t,real(y),'color',RGB(ii1,:))
hold on
grid on
end
legend(split(num2str(zeta)))
xlim([0 0.3])
xlabel('Time [s]')
ylabel('Ampli. [m]')

そもそもyは基礎から見た質量の変位なので、あまり変化が見られないのかもしれませんね。
振動特性を把握するには、周波数次元での応答倍率を確認するのが良いですよね。

ただ、この応答倍率を導出するのは難しいですよね。
導出過程は理解できていませんが、基礎の変位加振振幅に対する変位応答倍率λbは下式で表せるようです。

$$
λ_b = \frac{|X|}{B} = \frac{B\sqrt{1+(2ζβ)^2}}{B\sqrt{(1-β^2)^2+(2jζβ)^2}} \tag{2.123}
$$

この式は前回紹介した下記式と同じですよね。

$$
λ = sqrt{\frac{1+(2ζβ)^2}{ 1 – β^2 + 2jζβ }} \tag{2.107}
$$

従って、基礎の変位加振振幅に対する変位応答倍率λbは、質量mに作用する力が基礎に伝搬する応答倍率λと等しいことがわかります。

では、簡単にζを変化させた時のλbを図示してみましょう。結果は下図になります。

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:20; %周波数
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)
Lamda=sqrt( (1+(2*zeta(ii1).*beta).^2) ./ ((1-beta.^2).^2+(2*zeta(ii1).*beta).^2) );
plot(beta,Lamda,'color',RGB(ii1,:))
hold on
grid on
end
plot([0 sqrt(2)],[1 1],'k--','linewidth',3)
plot([sqrt(2) sqrt(2)],[0 1],'k--','linewidth',3)
legend(split(num2str(zeta)))
text(1.5,1,'(√2 , 1)')
ylim([0 3])
xlim([0 3])
xlabel('β=ω/Ω')
ylabel('λ')

 

今日の解説はここまで。

Good luck

 

コメント