MATLABで学ぶ振動工学　固有値・共振周波数の特定方法1 (最大虚数部法)

2021.06.07

はじめに

共振周波数や振動モードを特定するときに、
FRFのピークを見て、
「なんとなくここが共振周波数かな？」という感じで判断している方が多いのではないでしょうか？

ただ、かならずしも（FRFのピーク値の周波数）＝（共振周波数）が成立するわけではありません。

ハンマリング試験や加振試験を経験したことがある方はわかると思いますが、
低次のFRFピーク値は共振周波数と判断しやすいですが、
高次のFRFピーク値は狭帯域に密集するため、共振周波数なのか、ただのピーク値なのか、ノイズなのかの判断が難しいです。

本記事では共振周波数の特定方法について説明します。
（もっと優れた方法があるのですが、今回は初級編ということでお願いいたします。）

スポンサーリンク

共振周波数でのFRFの特徴（ナイキスト線図）

モード同定（カーブフィット）という用語は聞いたことありますよね。モード同定はFRFの結果からモードパラメーターを抽出する方法を意味します。

モード同定の最も古い手法は、ナイキスト線図を図示してパラメーターを抽出する方法です。

これまでに教科書や参考書を読んで、下記のような印象を持っている方もいると思います。

ナイキスト線図自体はわかってるけど、実用上は使わないし、勉強する意味あるか？

わかりますよ。私も読み飛ばしていましたから。。。
ただ、何年勉強していると避けて通れないことに気が付くと思います。

 

下図のFRFからナイキスト線図を求めてみましょう。

図1

スポンサーリンク

ナイキスト線図は下のようになります。


図2

このままだと、少しわかりにくいので、ナイキスト線図上に周波数を表示させると以下になります。

図3

ナイキスト線図でFRFを図示すると以下の特徴があります。

・FRFの共振周波数近傍の帯域でナイキスト線図を図示すると円形になる
・原点からの距離が最大になるときに共振周波数になる

などなど（他にも特徴がありますが、本日内容には関係ないので割愛します）

つまり、この結果では1.75Hz近傍で共振周波数となることがわかります。

 

ただ、FRFがたくさんある場合や共振周波数がたくさんある場合は、1つのピーク値の周波数に対してナイキスト線図を図示して共振周波数か否かを判断するのは面倒ですよね。

ただ、共振周波数の特徴を利用して、簡便に共振周波数を特定する方法があるんです。

最大虚数部法

簡便に共振周波数を特定する方法として、最大虚数部法というものあります。参考文献は下記です。

最大虚数部法は、各FRFの虚部のパワー（絶対値）を足し合わせることで、
（足し合わせた値のピーク値の周波数）＝（共振周波数）
が成立する方法です。（正確には、成立するように頑張る方法です。）

 

図1のFRFで最大虚数部法を適用してみましょう。その結果が図4です。

図4

どうでしょうか？
共振周波数がわかりやすくなりましたよね。
特に、FRFの結果だけでは、5.03Hzは共振周波数なのかの判断が難しかったと思いますが、図4ではっきりと共振周波数であることがわかります。

スポンサーリンク

MATLABのコード

実行ファイル

clear all
close all
% % % % 初期設定
freq=0.5:0.01:20;
m_vec=[1 1 1 1];
[M]=eval_Mmatrix(m_vec);
k_vec=[1 1 1 1]*10^3;
[K]=eval_Kmatrix(k_vec);
C=ones(length(k_vec))*0.1;
F=zeros(length(m_vec),length(freq));
F(1,:)=ones(1,length(freq));

% % % % 変位X（ノイズなし）
X=eval_direct_x_2ndedition(M,K,C,F,freq);
% % % % 理想的なFRF
FRF_ideal_x3_f1=X(3,:)./F(1,:);

figure(1)
loglog(freq,abs(FRF_ideal_x3_f1),'k','linewidth',7)
hold on
grid on
xlim([freq(1) freq(end)])
xlabel('周波数 Hz');ylabel('FRF m/N');

figure(2)
plot(real(FRF_ideal_x3_f1(51:171)),imag(FRF_ideal_x3_f1(51:171)),'o-k')
for ii1=51:171
text(real(FRF_ideal_x3_f1(ii1))*1.2,imag(FRF_ideal_x3_f1(ii1)),[num2str(freq(ii1)) 'Hz'])
end
grid on
xlim([-0.025 0.025])

% % % 最大虚数部法
FRF_ideal_x1_f1=X(1,:)./F(1,:);
FRF_ideal_x2_f1=X(2,:)./F(1,:);
FRF_ideal_x4_f1=X(4,:)./F(1,:);

imag_power=..........
+abs(imag(FRF_ideal_x1_f1)).......
+abs(imag(FRF_ideal_x2_f1)).......
+abs(imag(FRF_ideal_x3_f1)).......
+abs(imag(FRF_ideal_x4_f1));

figure
semilogx(freq,imag_power,'k','linewidth',7)
hold on
grid on
xlim([freq(1) freq(end)])
xlabel('周波数 Hz');ylabel('最大虚数部法');

functionファイル

function [K]=eval_Kmatrix(k_vec)
K=zeros(length(k_vec)); 
for ii1=1:length(k_vec) 
    if ii1==1 K(ii1,ii1)=k_vec(ii1); 
    else 
        K(ii1-1:ii1,ii1-1:ii1)=K(ii1-1:ii1,ii1-1:ii1)+[1 -1;-1 1]*k_vec(ii1);
    end
end

end

function [M]=eval_Mmatrix(m_vec)
M=diag(m_vec);
end

function x=eval_direct_x_2ndedition(M,K,C,F,freq)
% % 2021.05.28 Fを修正 
x=zeros(size(F,1),length(freq));
j=sqrt(-1);
for ii1=1:1:length(freq)
    w=2*pi*freq(ii1);
    x(:,ii1)=inv(K+j*w*C-w^2*M)*F(:,ii1);
end

end

コメント