はじめに
共振周波数や振動モードを特定するときに、FRFのピークを見て、
「なんとなくここが共振周波数かな？」という感じで判断している方が多いのではないでしょうか？
ただ、必ずしも（FRFのピーク値の周波数）＝（共振周波数）が成立するわけではありません。
ハンマリング試験や加振試験を経験したことがある方はわかると思いますが、
低次のFRFピーク値は共振周波数と判断しやすいですが、
高次のFRFピーク値は狭帯域に密集するため、共振周波数なのか、ただのピーク値なのか、ノイズなのかの判断が難しいです。
また、重根（同じ周波数に共振周波数が2つ存在する場合）を判断することができません。
本記事ではモード指示関数を用いた共振周波数の特定方法について説明します。
ちなみに、簡易的に共振周波数を求めたい場合は下記を参照ください。
スポンサーリンク
モード指示関数
モード指示関数（Modal Indicator Functions, MIFs）についてみなさんどこまでご存じでしょうか？
振動の研究をしている方でも詳しく知らない方が多いと思います。
モード指示関数にはCMIF、MMIF、RMIFといった種類があります。
ここでは、CMIF(Complex Modal Indicator Functions)について紹介します。
ちなみに、Ewin著の”Modal Testing”ではモード指示関数について下記のような説明がありました。
複数入力－複数応答の伝達関数行列を特異値分解もしくは固有値分解することで、周波数依存する固有値を取り出すこと
伝達関数行列がわからない方は下記を参照ください。
私の認識では、モード（共振周波数）がどこに存在するのかを指し示す指標です。モード指示関数について詳しく説明がある日本語の参考書はないと思います。
おそらく、下記の早稲田大学の成田氏の論文が最も詳細な説明があると思います。
CMIF(Complex Modal Indicator Functions)
CMIFの理論式を紹介します。理論式は下記文献を参考にしています。
CMIFは伝達関数行列H(ω)を特異値分解して、下式のように表します。
$$
H(ω)=U∑V
$$
$$
CIMF(ω)=∑^{T}∑
$$
スポンサーリンク
MATLABでの検証
実行ファイル①での検証
1点加振、4点応答の場合CMIFの結果を下図に示します。ちなみに、1点加振の場合では1次のCMIFしか求められません。（最大特異値のみでの結果しか求められない）
図から、共振周波数とCMIFのピーク値が一致していることがわかります。
図1 実行ファイル①の結果
このままではCMIFの良さがあまりわかりませんね。CMIFの良さは重根の判別です。
実行ファイル②での検証
段階を踏むために、実行ファイル①を少しだけ変更して、2点加振、4点応答の場合CMIFを求めます。
図2 実行ファイル②の結果
2点加振の場合では1次と2次のCMIFが求まります。（1次は黒点線、2次は灰色点線）
実行ファイル③での検証
では、重根を持つ場合を検証してみましょう。5Hz近傍に共振周波数が2つ存在するように質量とバネ定数を調整しました。この条件で2点加振、4点応答の場合CMIFを求めます。
図3 実行ファイル③の結果
実行ファイル②と同様に、2点加振の場合では1次と2次のCMIFが求まります。（1次は黒点線、2次は灰色点線）
さて、FRFの結果だけをみると、共振周波数は3だけのように見えますね。
しかし、CMIFを確認してみるとどうでしょう？
1次のCMIFでは共振周波数が3つあることがわかります。
2次のCMIFでも5Hz近傍にピークがあり、5Hz近傍にもう一つの共振周波数があることがわかります。
CMIFの特徴・メリットのまとめ
つまり、CMIFの特徴は下記ということです。
- 共振周波数にCMIFのピークが表れる
- 狭帯域に共振周波数が複数あっても2次や3次（高次）の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.1;
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=1;
CMIF_vector=zeros(inputN,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)
];
[U,S,V] = svd(H);
CMIF_vector(:,ii1)=diag(S.'*S);
end
figure
loglog(freq,abs([FRF_x1_f1;FRF_x2_f1;FRF_x3_f1;FRF_x4_f1]),'linewidth',3)
hold on
loglog(freq,CMIF_vector,'k--','linewidth',3)
xlim([freq(1) freq(end)])
xlabel('周波数 Hz');ylabel('FRF');
legend('FRF_1_1','FRF_2_1','FRF_3_1','FRF_4_1','CMIF')
実行ファイル②
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;
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;
CMIF_vector=zeros(inputN,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)
];
[U,S,V] = svd(H);
CMIF_vector(:,ii1)=diag(S.'*S);
end
figure
loglog(freq,abs([FRF_x1_f1;FRF_x2_f1;FRF_x3_f1;FRF_x4_f1]),'linewidth',3)
hold on
loglog(freq,CMIF_vector,'k--','linewidth',3)
xlim([freq(1) freq(end)])
xlabel('周波数 Hz');ylabel('FRF');
legend('FRF_1_1','FRF_2_1','FRF_3_1','FRF_4_1','CMIF')
実行ファイル③
clear all
close all
% % % % 初期設定
freq=0.1:0.01:50;
m_vec=[0.80 0.98 30 0.95];
[M]=eval_Mmatrix(m_vec);
k_vec=[1 10 0.90 1]*10^3;
[K]=eval_Kmatrix(k_vec);
[V,D] = eig(K,M);
fn=sqrt(diag(D))/(2*pi);
C=ones(length(k_vec))*1;
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));
F3=zeros(length(m_vec),length(freq));
F3(3,:)=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);
X3=eval_direct_x_2ndedition(M,K,C,F3,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,:);
FRF_x4_f3=X3(4,:)./F3(3,:);
FRF_x3_f3=X3(3,:)./F3(3,:);
FRF_x2_f3=X3(2,:)./F3(3,:);
FRF_x1_f3=X3(1,:)./F3(3,:);
inputN=2;
CMIF_vector=zeros(inputN,length(freq));
for ii1=1:length(freq)
H=[ FRF_x1_f1(ii1) FRF_x1_f3(ii1)
FRF_x2_f1(ii1) FRF_x2_f3(ii1)
FRF_x3_f1(ii1) FRF_x3_f3(ii1)
FRF_x4_f1(ii1) FRF_x4_f3(ii1)
];
[U,S,V] = svd(H);
CMIF_vector(:,ii1)=diag(S.'*S);
end
figure
loglog(freq,abs([FRF_x1_f1;FRF_x2_f1;FRF_x3_f1;FRF_x4_f1]),'linewidth',3)
hold on
loglog(freq,CMIF_vector,'k--','linewidth',3)
xlim([freq(1) freq(end)])
xlabel('周波数 Hz');ylabel('FRF');
legend('FRF_1_1','FRF_2_1','FRF_3_1','FRF_4_1','CMIF')
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
コメント