はじめに
振動工学やモード解析をご存じの方であれば、一度は読んだことがあると言われているバイブル的な書籍「モード解析入門」について再復習しております。今回はモード解析入門の6.2.1項の「周波数応答関数の大きさを用いる方法」P.344から復習します。なお、振動工学の参考書はどの本も似たり寄ったりで内容が重複しています。例えば、「実用モード解析入門」は長松先生の息子さんが著者ということもあり、重複している部分が多いです。
各参考書で重複している部分は本記事で大まかに理解できると思いますので、他の参考書から本記事へアクセスした場合も最後まで読んでいただけると幸いです。
6.2.1 周波数応答関数の大きさを用いる方法（半値幅法）
多自由度系のコンプライアンス（変位／力の周波数応答関数）は式(3.116)です。r次の共振周波数近傍では、このr次の振動モードが支配的な振動成分となります。そのため、この共振峰（共振周波数のピーク値周辺）は近似的に1自由度系とみなすことができます。例えば、下図のようなコンプライアンスの場合、各共振峰は1自由度系とみなすことができます。
$$
G(ω)= x_j / F_i e^{jωt} = (\sum_{r=1}^{N} \frac{ φ_{ri} φ_{rj} }{ -m_r ω^2 + j c_r ω + k_r }) (r=1～N) \tag{3.116}
$$
図1 共振峰が離れている場合のコンプライアンス
なお、下図のような2～3次の共振峰のように、2つの共振周波数が近く、共振峰が近接している場合は、お互いの振動成分が影響しているため1自由度系とはみなせないので注意が必要です。
図2 共振峰が近接する場合のコンプライアンス
式(3.116)が1自由度とみなせるということは、∑や添字rを省略することができるので、式(6.1)で表すことができます。
$$
G(ω) = \frac{ 1/K }{ 1 – β^2 + 2jζβ } , β=\frac{ ω }{Ω} \tag{6.1}
$$
ここで、式(6.1)のコンプライアンスの大きさは式(6.2)となります。
$$
|G(ω)| = \frac{ 1/K }{ \sqrt{ (1 – β^2)^2 + (2ζβ)^2 } } \tag{6.2}
$$
「モード同定」や「モード特性の同定」は測定データ（実験データ）に対して適用します。そのため、この記事のようにn自由度のバネマスモデルで求めた周波数応答関数（コンプライアンスなど）に対してモード同定するのは本来の使い方ではありません。（バネマスモデルであれば、質量行列と剛性行列から共振周波数やモードベクトルなどのモード特性が簡単に求まるので、わざわざ周波数応答関数をモード同定してモード特性を求めるのはナンセンスです。）
ただし、「モード同定」の原理を理解するときにはバネマスモデルを用いた方が、結果の検証などをしやすいのでバネマスモデルを用いて解説をします。
ここからが半値幅法の説明になります。ちなみに半値幅法とは、周波数応答関数からモード減衰比を推定する手法として良く知られていますが、モード減衰比がわかるとモード特性全てを求めることができます。
図1の1次共振峰の拡大図を図3に示します。
図3 1次共振峰の拡大図(m1の自己コンプライアンス)
式(6.2)の最大値|G|maxの1/√2になる2点(最大値から-3dBとなる2点)の周波数f1とf2を読み取り、共振周波数f0からモード減衰比ζを求めます。
$$
ζ= \frac{ f2-f1 }{ 2f_0 } \tag{6.7.1}
$$
なお、共振周波数f0を使わないでも式(6.7.2)で求めることができます。
$$
ζ= \frac{ f2-f1 }{ f2+f1 } \tag{6.7.1}
$$
また、コンプライアンスの最大値は式(6.5)で表すことができます。
$$
|G(ω)|_{max} ≒ \frac{ 1 }{ 2Kζ } \tag{6.5}
$$
ここで、このコンプライアンスG(ω)が自己コンプライアンスの場合を基準として、モードベクトルの振幅を1とするので、式(6.5)のKはモード剛性kと等価として扱えます。自己コンプライアンスとは、例えば、m1を加振した時のm1の応答から求めた伝達関数のことを指します。
モード剛性kが求まると、下式より簡単にモード質量が求めることができます。
$$
Ω^2 m = k
$$
ここまでで、共振周波数（最大ピーク値の周波数なのですぐわかる）、モード減衰比ζ、モード剛性k、モード質量mが判明しています。あとモードベクトルさえ判明すれば全てのモード特性を取得できたことになります。
自己コンプライアンスの場合はモードベクトルの成分が1となりました。
例えば、図3の場合では共振周波数9.432Hzのときに、m1に対応するモードベクトルの成分が1です。
ではm2、m3、….はどのように求めるのでしょうか？
それは相互コンプライアンスGis（加振点がiで応答点がs）を測定して求めることになります。φi=1なので、モードベクトルの成分φsは式(6.11)となります。
$$
K_{is}= \frac{ k }{ φ_s } \tag{6.11}
$$
また、相互コンプライアンスの最大値との関係は式(6.12)となります。
$$
|G_{is}(ω)|_{max} ≒ \frac{ 1 }{ 2K_{is}ζ } \tag{6.12}
$$
したがって、モードベクトルの成分φsは式(6.13)で求めることができます。
$$
φ_s ≒ 2kζ|G_{is}(ω)|_{max} \tag{6.13}
$$
例題
では、図3を例題に解いてみましょう。
答え（バネマスモデルで設定したモード特性）
図3の例題は、4自由度のバネマスモデルでモード特性は下記です。
1次の共振周波数：9.43205759496419Hz
モード減衰比：0.1%
モード質量：1
モード剛性：3512.14651375420
モードベクトル：
v=[-0.345365005476218
-0.487723195715426
-0.554079400516416
-0.579521453681323];
なお、半値幅法ではモードベクトルの加振点が1として正規化されるので、モードベクトル、モード質量およびモード剛性の値はここの設定値である「答え」とは変わってしまいます。
なので、m1を加振して、m1のモードベクトルの成分が1になるようにしておきましょう。
モード減衰比：0.1%
モード質量：8.38384691873869
モード剛性：29445.2987274969
モードベクトル：
v=[1
1.41219633715614
1.60432988788892
1.67799703065524];
半値幅法での求めた値
共振周波数f0：9.4320Hz
最大値が1/√2になる周波数f1：9.4226Hz
最大値が1/√2になる周波数f2：9.4414Hz
式(6.7.1)よりモード減衰比ζは約0.1%ということになりますね。
f0=9.4320; f1=9.4226; f2=9.4414; zeta_hanchi=(f2 - f1)/(f0*2) zeta_hanchi = 9.9661e-04
式(6.5)より、モード剛性kが求まり、モード質量mも求まります。
モード剛性：29545.6919999358
モード質量：8.41253424855546
式(6.13)よりモードベクトルを求めてみましょう。
phys =
1.0000
1.4122
1.6043
1.6780
正しくモード特性を同定できていることがわかります。
実験データで半値幅法を使うときの注意点
この記事ではMATLABで周波数応答関数を求めているので、周波数分解能が細かくできます。しかし、実験データを使う場合は細かくできない場合がありますので、ご注意ください。細かくできないと、半値となる周波数がわからず、半値幅法が適用できません。
実験データで周波数分解能を細かくしたければ、時刻歴データを長めに測定し、FFTの点数（ライン数）を大きくとることで周波数分解能を細かくできます。もし、インパルス加振をしていて、時刻歴データが1s程度しか測定してなかった場合は、疑似的に無音(振幅が0のデータ)を足すことで、周波数分解能を上げることができます。
MATLABプログラム（例題の解答）
実行コード
計算周波数(freq)の周波数分解能やモード減衰比(modal_dampimg)を変更して、ほかのケースでも試してみてください。
clear all;clc;close all
freq=9.43-0.1:0.0001:9.43+0.1; %計算周波数の定義
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
for ii1=1:length(m_vec)
V(:,ii1)=V(:,ii1)/V(1,ii1); % 加振点m1のモードベクトルの成分が1になるように正規化
end
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.1*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,:)+V(ii1,ii2)*V(ii,ii2)./( -mr(ii2).*w.^2 + 1i*cr(ii2)*w +kr(ii2) )*F(ii);
% 横ﾍﾞｸﾄﾙ =横ﾍﾞｸﾄﾙ + ｽｶﾗｰ ×ｽｶﾗｰ ./( ｽｶﾗｰ×横ベクトル + ｽｶﾗｰ×横ベクトル +ｽｶﾗｰ）×ｽｶﾗｰ
end
end
figure
semilogy(freq,abs(Xj(1,:)),'k')
hold on
legend('モード法')
xlabel('周波数 Hz')
ylabel('変位 [m]')
figure
semilogy(freq,abs(Xj(1,:)),'k')
hold on
% legend('モード法')
xlabel('周波数 Hz')
ylabel('変位 [m]')
xlim([freq(1) freq(end)])
semilogy([8.5 10],[1 1]*max(abs(Xj(1,:)))/sqrt(2),'--k')
% % % % 半値幅法
f0=9.4320;
f1=9.4226;
f2=9.4414;
zeta_hanchi=(f2 - f1)/(f0*2);
Gmax=max(abs(Xj(1,:)));
k_hanchi=1/(2*zeta_hanchi*Gmax);
m_hanchi=k_hanchi/(2*pi*f0)^2;
% Gis=max(abs(Xj(:,:)),[],2);
% phys=2*k_hanchi*zeta_hanchi*Gis
[Gis,id]=max( abs( imag(Xj(:,:)) ) ,[],2 );
Gis=Gis .* imag(-Xj(:,id(1)))./abs(Xj(:,id(1))) ;
phys=2*k_hanchi*zeta_hanchi*Gis;
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
コメント