はじめに
共振周波数や振動モードを特定するときに、FRFのピークを見て、
「なんとなくここが共振周波数かな？」という感じで判断している方が多いのではないでしょうか？
ただ、必ずしも（FRFのピーク値の周波数）＝（共振周波数）が成立するわけではありません。
ハンマリング試験や加振試験を経験したことがある方はわかると思いますが、
低次のFRFピーク値は共振周波数と判断しやすいですが、
高次のFRFピーク値は狭帯域に密集するため、共振周波数なのか、ただのピーク値なのか、ノイズなのかの判断が難しいです。
本記事ではモード指示関数を用いた共振周波数の特定方法について説明します。
ちなみに、簡易的に共振周波数を求めたい場合は下記を参照ください。
前回紹介したCMIFも参照ください。
スポンサーリンク
モード指示関数
モード指示関数（Modal Indicator Functions, MIFs）についてみなさんどこまでご存じでしょうか？
振動の研究をしている方でも詳しく知らない方が多いと思います。
モード指示関数にはCMIF、MMIF、RMIFといった種類があります。
ここでは、MMIF(Multivariate Modal Indicator Functions)について紹介します。
ちなみに、Ewin著の”Modal Testing”ではモード指示関数について下記のような説明がありました。
複数入力－複数応答の伝達関数行列を特異値分解もしくは固有値分解することで、周波数依存する固有値を取り出すこと
伝達関数行列がわからない方は下記を参照ください。
MMIF(Multivariate Modal Indicator Functions)
MMIFの理論式を紹介します。理論式は下記文献を参考にしています。
MMIFは0～1の値をとり、共振周波数が存在する周波数で0に近くなるという特徴があります。なので、簡単に共振周波数を判別できます。
MMIFは下式の固有値βです。
ここで、H_Realは伝達関数行列H(ω)の実部、H_Imagは伝達関数行列H(ω)の虚部を意味します。式が長くなるのでωは省略しています。（Ewin著 “Modal Testing”）
$$
β([H_{Real}]^T[H_{Real}]+[H_{Imag}]^T[H_{Imag}]){F}=([H_{Imag}]^T[H_{Imag}]){F}
$$
上記のような式は最適化問題などで良く出てくる形です。
初見の方は どのように解くかわからない という方も多いと思います。
数式や導入式を理解するよりも、MATLABのコードを見た方が理解しやすいと思うので、理解できないという方はコードを見て頂けると幸いです。
ちなみに記憶が正しければ、たしか上式は間違ってたように思います。正しくは下式だったような気がします。（学生時代に検証したのですが、改めて再検証したいと思います。自信がないので文字を灰色にしました。）
$$
β([H_{Real}]^T[H_{Real}]+[H_{Imag}]^T[H_{Imag}]){F}=([H_{Real}]^T[H_{Real}]){F}
$$
スポンサーリンク
MATLABでの検証
実行ファイルでの検証
2点加振、4点応答の場合MMIFの結果を下図に示します。
図から、共振周波数とMMIFのピーク値が一致していることがわかります。
図1 実行ファイルの結果
ちなみに、市販の実験モード解析のソフトでは、共振周波数判定画面でMMIFを表示していることが多いです。それだけ、信用があり、直感的にわかりやすいのだと思います。ただ、ソフトの使用者はMMIFの意味を理解していないで使っていることが多いです。
私はMMIFとCMIFの両方を使用して、重根がないかを確認しながら共振周波数を判別・選定することを推奨します。理由は、MMIFとCMIFの両方を確認しておく癖があると、高周波数を分析する際に、共振周波数が密集していても共振周波数の判別が容易にできるようになるからです。
スポンサーリンク
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.2;
F1=zeros(length(m_vec),length(freq));
F1(1,:)=ones(1,length(freq));
F2=zeros(length(m_vec),length(freq));
F2(2,:)=ones(1,length(freq));
% % % % 変位X（ノイズなし）
X1=eval_direct_x_2ndedition(M,K,C,F1,freq);
X2=eval_direct_x_2ndedition(M,K,C,F2,freq);
% % % % 理想的なFRF
FRF_x4_f1=X1(4,:)./F1(1,:);
FRF_x3_f1=X1(3,:)./F1(1,:);
FRF_x2_f1=X1(2,:)./F1(1,:);
FRF_x1_f1=X1(1,:)./F1(1,:);
FRF_x4_f2=X2(4,:)./F2(2,:);
FRF_x3_f2=X2(3,:)./F2(2,:);
FRF_x2_f2=X2(2,:)./F2(2,:);
FRF_x1_f2=X2(1,:)./F2(2,:);
inputN=2;
MMIF=zeros(1,length(freq));
for ii1=1:length(freq)
H=[ FRF_x1_f1(ii1) FRF_x1_f2(ii1)
FRF_x2_f1(ii1) FRF_x2_f2(ii1)
FRF_x3_f1(ii1) FRF_x3_f2(ii1)
FRF_x4_f1(ii1) FRF_x4_f2(ii1)
];
HH=pinv(real(H).'*real(H)+imag(H).'*imag(H))*real(H).'*real(H);
[V,D]=eig(HH);
MMIF(1,ii1)=min(diag(D));
end
figure
subplot(211)
loglog(freq,abs([FRF_x1_f1;FRF_x2_f1;FRF_x3_f1;FRF_x4_f1]),'linewidth',3)
xlim([freq(1) freq(end)])
xlabel('周波数 Hz');ylabel('FRF');
legend('FRF_1_1','FRF_2_1','FRF_3_1','FRF_4_1')
subplot(212)
loglog(freq,MMIF,'--','linewidth',3)
xlim([freq(1) freq(end)])
xlabel('周波数 Hz');ylabel('MMIF');
legend('MMIF')
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
コメント