MATLABで学ぶモード同定法最も簡便な方法 - MATLABパイセンが教える振動・騒音・音響・機械工学

はじめに
理論の説明
例題
MATLABプログラム（例題の解答）
1. 関連

はじめに

これまでに、モード円適合やプロニー法（ポリレファレンス法）、偏分反復法、直交多項式法（RFP法）を解説しました。これらの手法は専門的知識がないと自分でプログラムを組むことは難しいです。そのため、振動の専門家ではない設計者たちでは簡単に理解できる半値幅法がいまだに使われています。しかし、半値幅法は周波数分解能が大きい場合、減衰が小さい場合、低周波数に共振がある場合などで精度が低下するという課題があります。

つまり、私のようなモード同定マニア以外の方にはより簡便で精度の良いモード同定手法が求められていたことになります。2023年に中部大学から発表されたモード同定手法(リンクの109)は、半値幅法よりも簡単に理解でき、且つ半値幅法のデメリットが改善されたモード同定手法が提案されました。今回はそのモード同定手法について解説します。

理論の説明

本モード同定手法は1自由度法なので、1つの共振に対して、共振周波数とモード減衰を同定する手法です。隣り合う共振が近接しない場合（意味がわからない方はこちらを確認ください）、r次の共振周波数近傍のコンプライアンスG(変位／力の周波数応答関数)は式(1)で表せます。ここで、ωは角周波数、mrとcrおよびkrはr次のモード質量、モード減衰、モード剛性を意味し、φriとφriはi点およびj点のr次のモードベクトルの値を意味する。

$$
G(ω)= x_j / F_i e^{jωt} = ( \frac{ φ_{ri} φ_{rj} }{ -m_r ω^2 + j c_r ω + k_r })　 \tag{1}
$$

さらに式をまとめると式(2)となる。ここで、Ωはr次の角共振周波数を意味する

$$
G(ω) = \frac{ \frac{ φ_{ri} φ_{rj} }{k_r} }{ 1 – β^2 + 2jζβ } , β=\frac{ ω }{Ω} \tag{2}
$$

一般的に、周波数応答関数はコンプライアンス(変位／力の周波数応答関数)ではなく、アクセレランス(加速度／力の周波数応答関数)を用いることが多いので、式(2)をアクセラレンスLの式(3)に変更する。

$$
L(ω) = \frac{ -ω^2 \frac{ φ_{ri} φ_{rj} }{k_r} }{ 1 – β^2 + 2jζβ } , β=\frac{ ω }{Ω} \tag{3}
$$

式(3)を実部と虚部に分ける

式(3)を実部と虚部で表現すると式(4)および式(5)となる。

$$
Real[L(ω)] = \frac{ -ω^2 \frac{ φ_{ri} φ_{rj} }{k_r} (1 – β^2)^2 }{ (1 – β^2)^2 + (2ζβ)^2 } \tag{4}
$$

$$
Imag[L(ω)] = \frac{ ω^2 \frac{ φ_{ri} φ_{rj} }{k_r} (2ζβ) }{ (1 – β^2)^2 + (2ζβ)^2 } \tag{5}
$$

実部と虚部の比は式(6)となる。

$$
\frac{Imag[L(ω)]}{Real[L(ω)]} = \frac{-2ζβ }{1 – β^2} \tag{6}
$$

式(3)を振幅と位相で表現する

式(3)を振幅と位相で表現すると式(7)となる。

$$
L(ω) =|ω^2 G(ω)|e^{jθ} \tag{7}
$$

式(7)を実部と虚部に分けると式(8)および式(9)となる。

$$
Real[L(ω)] =|ω^2 G(ω)| cosθ \tag{8}
$$

$$
Imag[L(ω)] =|ω^2 G(ω)| sinθ \tag{9}
$$

実部と虚部の比は式(10)となる。

$$
\frac{Imag[L(ω)]}{Real[L(ω)]} = tanθ \tag{10}
$$

式(6)と式(10)から共振周波数とモード減衰比を求める

式(6)と式(10)は導出過程は異なりますが、アクセラレンスLの虚部と実部の比なので等価です。つまり式(6)＝式(10)となります。

実験モード解析では、アクセラレンスLの実測結果があります。r次の共振周波数近傍の2つの角周波数ω1とω2に対して、式(11)および式(12)が成立します。

$$
\frac{-2ζβ1 }{1 – β1^2}=\frac{-2ζ \frac{ ω1 }{Ω} }{1 – { \frac{ ω1 }{Ω} }^2} = tanθ1 \tag{11}
$$

$$
\frac{-2ζ \frac{ ω2 }{Ω} }{1 – { \frac{ ω2 }{Ω} }^2} = tanθ2 \tag{12}
$$

式(11)および式(12)をΩとζについて解くと式(13)と式(14)で表せます。

$$
Ω = \sqrt{ \frac{ ω1 ω2 (ω1 – ω2) \frac{tanθ2}{tanθ1} }{ω2 – ω1 \frac{tanθ2}{tanθ1}} } \tag{13}
$$

$$
ζ = \frac{1}{2} tanθ1 (\frac{Ω^2-ω1^2}{Ω ω1}) \tag{14}
$$

これまで数多くのモード同定手法を勉強してきていますが、高校生でもわかるような式展開で共振周波数（角共振周波数）とモード減衰比を同定できる手法はないです。

いやーすごいですね。

例題

では、図1を例題に解いてみましょう。図1に示すように1次共振周波数は8~10Hzの間に存在しますよね。そこで8Hzと10Hzのデータから1次の共振周波数とモード減衰比を求めてみましょう。

図1

式(13)と式(14)から求めた結果は表1です。高精度に同定できていることがわかりますね。

表1

	真値	同定結果	誤差 [%]
共振周波数[Hz]	9.4321	9.5061	-0.7851
モード減衰比	0.01	0.0093	7.0806

MATLABプログラム（例題の解答）

実行するとこんな図が出力されます。

実行コード

clear all;clc;close all

freq=1:1:100; %計算周波数の定義
w=2*pi*freq; %角振動数

m_vec=ones(1,4);
k_vec=(1:4)*2*10^4;
[M]=eval_Mmatrix(m_vec); %質量行列
[K]=eval_Kmatrix(k_vec); %剛性行列

[V,D]=eig(K,M); %固有値D、固有ベクトルV
wn=sqrt(diag(D)); %固有角振動数
fn=wn/(2*pi); %固有振動数

mr=diag(V.'*M*V); %モード質量行列
kr=diag(V.'*K*V); %モード剛性行列
c_cr=2*sqrt(mr.*kr);
modal_dampimg=0.01; %モード減衰比（とりあえず0.1%と設定）
cr=c_cr*modal_dampimg; %モード減衰係数

F=zeros(length(m_vec),1); F(1)=1; %外力ベクトル(全周波数帯を1で加振する)
Xj=zeros(length(m_vec),length(freq)); %加速度ベクトル（計算前にベクトルを定義してメモリを確保）

% % % % モード法
ii=1; %加振点 （加振点が複数ある場合はfor ii=1:kength(F)のループを追加すればよい ）
for ii1=1:length(m_vec) %応答点
    for ii2=1:length(wn)
        Xj(ii1,:)=Xj(ii1,:)+w.^2.*V(ii1,ii2)*V(ii,ii2)./( -mr(ii2).*w.^2 + 1i*cr(ii2)*w +kr(ii2) )*F(ii);
%       横ﾍﾞｸﾄﾙ　=横ﾍﾞｸﾄﾙ  +  ｽｶﾗｰ　　×ｽｶﾗｰ　　　./(　ｽｶﾗｰ×横ベクトル + ｽｶﾗｰ×横ベクトル +ｽｶﾗｰ）×ｽｶﾗｰ
    end
end


% % % % % ω1とω2の設定
id1=8;
id2=10;
w1=w(id1);
w2=w(id2);

% % % % % tanθ1とtanθ2を算出
tan1=tan(angle(Xj(1,id1)));
tan2=tan(angle(Xj(1,id2)));

wr=sqrt( (w1*w2*(w1-w2*tan2/tan1))/(w2-w1*tan2/tan1) ); %式13　角周波数
fn_re=wr/(2*pi);  %式13の周波数
zeta1=-1/2*tan1*((wr^2-w1^2)/(wr*w1));  %式14


figure(1)
semilogy(freq,abs(Xj(1,:)),'k')
hold on
semilogy([fn_re fn_re],[min(abs(Xj(1,:))) max(abs(Xj(1,:)))],'r--')
% legend('モード法')
xlabel('周波数 Hz')
ylabel('加速度/力 [(m/s^2)/N]')
message1=['共振周波数 ' num2str(fn_re) 'Hz   ,モード減衰比' num2str(zeta1)];
text(fn_re*0.5,max(abs(Xj(1,:)))*1.2,message1)
% % [fig1]=figure_EYEcatch_size(1);

disp('真値の共振周波数')
disp(fn(1))
disp('同定した共振周波数')
disp(fn_re)
disp('誤差 %')
disp((fn(1)-fn_re)/fn(1)*100)

disp('真値のモード減衰比')
disp(modal_dampimg)
disp('同定したモード減衰比')
disp(zeta1)
disp('誤差 %')
disp((modal_dampimg-zeta1)/modal_dampimg*100)

functionファイル

function [M]=eval_Mmatrix(m_vec)

M=diag(m_vec);

end

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