はじめに
振動工学やモード解析をご存じの方であれば、一度は読んだことがあると言われているバイブル的な書籍「モード解析入門」について再復習しております。今回はモード解析入門の6.2.2～6.2.3項の「周波数応答関数の虚部を用いる方法」P.347から復習します。長松先生は「半値幅法」という言葉は使っていませんが、一般的には6.2.1～6.2.3は「半値幅法」と呼ばれているモード減衰比ζを求めるための手法に分類できると思います。
なお、振動工学の参考書はどの本も似たり寄ったりで内容が重複しています。例えば、「実用モード解析入門」は長松先生の息子さんが著者ということもあり、重複している部分が多いです。
各参考書で重複している部分は本記事で大まかに理解できると思いますので、他の参考書から本記事へアクセスした場合も最後まで読んでいただけると幸いです。
6.2.2 周波数応答関数の虚部を用いる方法（最大虚数部法＋半値幅法）
6.2.1項の復習になりますが、隣り合う共振周波数が十分に離れている場合はその共振峰を1自由度とみなせます。そのため、コンプライアンス（変位／力の周波数応答関数）は式(6.1)で表すことができます。
$$
G(ω) = \frac{ 1/K }{ 1 – β^2 + 2jζβ } , β=\frac{ ω }{Ω} \tag{6.1}
$$
式(6.1)の実部と虚部は、式(6.14)で表すことができます。
$$
G_R(ω) = \frac{ (1 – β^2)/K }{ (1 – β^2)^2 + (2ζβ)^2 } , G_I(ω) = \frac{ -2ζβ/K }{ (1 – β^2)^2 + (2ζβ)^2 } \tag{6.14}
$$
ここからは6.2.1項の半値幅法と同じ流れです。ちなみに半値幅法とは、周波数応答関数からモード減衰比を推定する手法として良く知られていますが、モード減衰比がわかるとモード特性全てを求めることができます。
1次共振峰の拡大図を図1に示します。
図1 1次共振峰の拡大図(m1の自己コンプライアンス)
図1に示すように式(6.16)が成立していますね。
$$
|G_I(ω)|_{max}=|G(ω)|_{max} ≒ \frac{ 1 }{ 2Kζ } \tag{6.16}
$$
式(6.16)の最大値|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}
$$
ここで、このコンプライアンス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_{I is}(ω)|_{max} ≒ \frac{ 1 }{ 2K_{is}ζ } \tag{6.12}
$$
したがって、モードベクトルの成分φsは式(6.13)で求めることができます。
$$
φ_s ≒ 2kζ|G_{I 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];
6.2.2項での算出方法で求めた値
共振周波数f0：9.4321Hz
最大値が1/√2になる周波数f1：9.4260Hz
最大値が1/√2になる周波数f2：9.4381Hz
式(6.7.1)よりモード減衰比ζは約0.1%ということになりますね。ただ、6.2.1項の方が精度が高かったですね。
f0=9.4321; f1=9.4260; f2=9.4381; zeta_hanchi=(f2 - f1)/(f0*2) zeta_hanchi = 6.4143e-04
式(6.5)より、モード剛性kが求まり、モード質量mも求まります。
モード剛性：45907.0580387381
モード質量：13.0708229511585
式(6.13)よりモードベクトルを求めてみましょう。
phys =
1.0000
1.4122
1.6043
1.6780
正しくモード特性を同定できていることがわかります。
実験データで半値幅法を使うときの注意点
この記事ではMATLABで周波数応答関数を求めているので、周波数分解能が細かくできます。しかし、実験データを使う場合は細かくできない場合がありますので、ご注意ください。細かくできないと、半値となる周波数がわからず、半値幅法が適用できません。
実験データで周波数分解能を細かくしたければ、時刻歴データを長めに測定し、FFTの点数（ライン数）を大きくとることで周波数分解能を細かくできます。もし、インパルス加振をしていて、時刻歴データが1s程度しか測定してなかった場合は、疑似的に無音(振幅が0のデータ)を足すことで、周波数分解能を上げることができます。
6.2.3 周波数応答関数の実部と虚部を用いる方法
図2に示すように、実部と虚部が等しくなる周波数において、式(6.18)が成立します。
$$
β=\frac{ ω }{Ω} , Δβ^2≒4ζ^2 \tag{6.18}
$$
図2 実部と虚部が等しくなる周波数
つまり、式(6.18)からモード減衰比ζが求めることができます。あとは、6.2.1項の半値幅法や式(6.16)からの流れと同様です。
6.2.3項での算出方法で求めた値
式(6.18)より求めたモード減衰比ζは約0.1%ということになります。6.2.2項の算出方法よりも精度が高くなりました。
zeta_623 = 9.9660e-04
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
xlabel('周波数 Hz')
ylabel('変位 [m]')
xlim([freq(1) freq(end)])
semilogy(freq,abs(real(Xj(1,:))),'r')
semilogy(freq,abs(imag(Xj(1,:))),'b')
legend('|G|','|G_R|','|G_I|')
semilogy([8.5 10],[1 1]*max( abs(imag(Xj(1,:))) )/sqrt(2),'--k')
% % % % % % % % % % % 6.2.2項 最大虚数部法&半値幅法
f0=9.4321;
f1=9.4260;
f2=9.4381;
zeta_hanchi=(f2 - f1)/(f0*2)
Gmax=max(abs(imag(Xj(1,:))));
k_hanchi=1/(2*zeta_hanchi*Gmax);
m_hanchi=k_hanchi/(2*pi*f0)^2;
% Gis=max( abs( imag(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;
% % % % % % % % % % % 6.2.3項
Omega=2*pi*f0;
w1=2*pi*9.4227;
w2=2*pi*9.4415;
beta1=w1/Omega;
beta2=w2/Omega;
delta_beta=(beta2-beta1);
zeta_623=sqrt( delta_beta^2/4 )
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
コメント