MATLABで学ぶモード同定法 モード円適合/カーブフィット

はじめに

振動工学やモード解析をご存じの方であれば、一度は読んだことがあると言われているバイブル的な書籍「モード解析入門」について再復習しております。今回はモード解析入門の6.2.4項の「モード円適合」P.349から復習します。モード円適合はカーブフィットとも呼ばれます。モード解析入門ではかなりアッサリと説明されているのですが、この記事では内容を補足して説明したいと思います。

なお、振動工学の参考書はどの本も似たり寄ったりで内容が重複しています。例えば、「実用モード解析入門」は長松先生の息子さんが著者ということもあり、重複している部分が多いです。

各参考書で重複している部分は本記事で大まかに理解できると思いますので、他の参考書から本記事へアクセスした場合も最後まで読んでいただけると幸いです。

参考書解説一覧

MATLABとNVH(Noise, Vibration and Harshness)をこよなく愛する人生の先輩、それが私（MATLABパイセン）です。NVHを極めるにあたり、周辺の知識も網羅的に勉強することになり、その知識を共有すべく、本HPを運営しています。日本の製造業を応援すべく、機械エンジニアの「車輪の再開発」の防止と業務効率化の手助けをモットーに活動。専門である振動騒音工学や音響工学について、プログラムを示しながら解説。大学の授業よりもわかりやすく説明することを目指す。質問頂ければ、回答したいと思います。

mech-eng-uni.com

2023.01.07

6.2.4 モード円適合（カーブフィット）

まず、下図のようなFRFを例として考えます。

図1 FRF

図1に示した赤枠の周波数帯域のナイキスト線図を図2に示します。

図2 ナイキスト線図（5～15Hz）

図2に示すように、FRFのナイキスト線図は赤線となります。周波数分解能によりますが、ナイキスト線図ではカクカクの円となります。（教科書などでは黒点線のように綺麗な円で説明されていますが…）

モード円適合では、図2に示すようなナイキスト線からモード特性（共振周波数およびモード減衰、モードベクトル）を同定すること意味します。なお、図1および図2は1次共振周波数の共振峰が他の共振周波数の影響を受けていないと考えられます。このような場合にのみモード円適合は適用可能です。

他の振動モード（共振周波数）の影響を受けない共振峰のナイキスト線図は図3のような円となります。

図3 コンプライアンスのナイキスト線図（他の振動モードの影響が無視できる場合）

$$
中心 = (R, I – \frac{ 1 }{4Kζ}) \tag{6.21}
$$

また、これは下記のような関係があります。

$$
a=2R, b=2I – \frac{ 1 }{2Kζ}, c=\frac{ I }{2Kζ}-R^2-I^2 \tag{6.25}
$$

$$
円の半径=\sqrt{c+\frac{ a^2+b^2}{4}}, 中心=(\frac{ a }{2}, \frac{ b }{2}) \tag{6.27}
$$

ここで、コンプライアンスの実部をx、虚部をyとして、係数a,b,cを式(6.30)の最小二乗法の関係で表すことができます。なお、l（小文字のL）はデータ数を意味します。

\[
\left[ \begin{array}{ccc}
∑x_i^2 & ∑x_iy_i & ∑x_i \\
∑x_iy_i & ∑y_i^2 & ∑y_i \\
∑x_i & ∑y_i & l \\
\end{array} \right]
\left[ \begin{array}{c} a  \\ b  \\ c \end{array} \right]
=\left[ \begin{array}{c} ∑(x_i^3+x_iy_i^2)  \\ ∑(x_i^2y_i+y_i^3)  \\ ∑(x_i^2+y_i^2) \end{array} \right] \tag{6.30}
\]

式(6.30)により、a,b,cが求まるので、図2のような黒点線の円が描けるようになる。
しかし、綺麗な円が描けるだけではモード特性（共振周波数、モード減衰比、モードベクトル）が求まるわけではありません。

ではここから、モード特性を求める手順を示していきたいと思います。

図4 コンプライアンスのナイキスト線図 part2

まず、図4に示す共振周波数fnを求める必要があります。
一般的は、9.4Hzと9.5Hzの中間なので、（9.4Hz＋9.5Hz）/2＝9.45Hzと求めることが多いです。ただし、これは長松先生がモード解析入門を1990年代の初頭の話であって、現在のように数値演算が簡単に高精度に行えるのであれば、もう少し賢く求めたいと思います。

では、どのように共振周波数fnを求めるのか？
式(6.27)より円の中心が判明しているので、共振周波数fnがとる円上の位置はわかっています。ただし、周波数との対応が不明なので、円上の位置と周波数の対応をとる必要があります。
図4のように、共振周波数fnを跨ぐ二つの周波数（9.4Hzと9.5Hz）が判明しているので、下式のように線形補完して求めることができます。

$$
fn=\frac{θ_d}{θ_d+θ_e}df + fd
$$

なお、角度を線形補完しているので、誤差が生じますが、カーブフィットをするときの周波数分解能dfは大きくても2Hz程度なので、この誤差は許容できるかと思います。

共振周波数fnを求めるとモード減衰比が求まります。
モード減衰比の算出方法を下記に示します。

$$
βd=\frac{fd}{fn}, βe=\frac{fe}{fn} \tag{B9.2}
$$

$$
ζ=\frac{βe^2-βd^2}{2(βd\tan{\frac{θd}{2}}+βe\tan{\frac{θe}{2}})} \tag{B9.5}
$$

後の流れは半値幅と同様です。念のために再度説明したいと思います。

コンプライアンスG(ω)が自己コンプライアンスの場合を基準として、モードベクトルの振幅を1とするので、式(6.5)のKはモード剛性kと等価として扱えます。自己コンプライアンスとは、例えば、m1を加振した時のm1の応答から求めた伝達関数のことを指します。

モード剛性kが求まると、下式より簡単にモード質量が求めることができます。

$$
Ω^2 m = k
$$

ここまでで、共振周波数fn、モード減衰比ζ、モード剛性k、モード質量mが判明しています。あとモードベクトルさえ判明すれば全てのモード特性を取得できたことになります。

自己コンプライアンスの場合はモードベクトルの成分が1となりました。
一方、相互コンプライアンスGis（加振点がiで応答点がs）に関しては、下記手順で求めることになります。φi=1なので、モードベクトルの成分φsは式(6.11)となります。

$$
K_{is}= \frac{ k  }{ φ_s }   \tag{6.11}
$$

また、相互コンプライアンスの共振周波数の値との関係は式(6.12)となります。

$$
|G_{is}(ωn)| ≒ \frac{ 1  }{ 2K_{is}ζ }   \tag{6.12}
$$

したがって、モードベクトルの成分φsは式(6.13)で求めることができます。
$$
φ_s ≒ 2kζ|G_{is}(ωn)|  \tag{6.13}
$$

プログラムの検証

4dofモデルでプログラムの精度を検証したいと思います。表1に理論値を示します。この理論値を「モード円適合」の理論で求めることができるのか検証していきます。なお、周波数分解能dfが小さいと精度が高くなることが知られているので、周波数分解能との関係も検証したいと思います。

理論値

同定結果（周波数分解能が0.1Hzのとき）

同定結果（周波数分解能が1Hzのとき）

同定結果（周波数分解能が2Hzのとき）

精度比較

まず、共振周波数の推定精度を表1で比較します。表1に示すように、理論値を基準（正解値）として誤差を算出しています。やはり、周波数分解能が小さいほど誤差が小さいことがわかりますね。個人的には、周波数分解能が2Hzのときの誤差が最大8.2%なのが驚きでした。もっと低いと思っていたので、周波数分解能の設定は重要ということがわかりますね。

表1　共振周波数の推定精度

次に、モード減衰比の推定精度を表2で比較します。表2に示すように、周波数分解能が大きくなると誤差が大きくなり、推定精度が低下していることがわかります。個人的な感覚としては、モード減衰比は誤差60～70%程度であれば許容可能かなーという感じです。（モード減衰比は最悪のケースで倍ぐらい違ってもしょうがないかなーという感覚でいます。）

なので、周波数分解能2Hzの4次共振のモード減衰比の推定精度は許容できないレベルで誤差が大きいですね。

表2　モード減衰比の推定精度

最後にモードベクトルの精度を検証します。モードベクトルは図5に示すようにクロスMACを求めて精度を検証しました。図中の文字が小さくなってしまいましたが、左から周波数分解能0.1Hz、中央が1Hz、右が2Hzです。モードベクトルは周波数分解能1Hzであれば十分な精度で求められていることがわかります。周波数分解能2Hzのときは、(理論の4次)×（モード円適合の2次）が0.5と類似度が高く、モードベクトルの同定がうまくいっていないことがわかります。

図5 クロスMACによるモードベクトルの精度検証結果

MATLABプログラム

実行プログラム

clear all;clc;close all


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); %固有振動数

for ii10=1:4
freq=fn(ii10)-7:2:fn(ii10)+7; %計算周波数の定義
w=2*pi*freq; %角振動数



mr=diag(V.'*M*V); %モード質量行列
kr=diag(V.'*K*V); %モード剛性行列
c_cr=2*sqrt(mr.*kr);
modal_dampimg=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))


% % % % % % % % ここからが「モード円適合」のプログラム
[zeta1,R1,I1,k1,m1,modal_vector1,fn_curve1] = modal_circle_fit(Xj(1,:),freq,[1 1]);
[zeta2,R2,I2,k2,m2,modal_vector2,fn_curve2] = modal_circle_fit(Xj(2,:),freq,[2 1],k1,m1);
[zeta3,R3,I3,k3,m3,modal_vector3,fn_curve3] = modal_circle_fit(Xj(3,:),freq,[3 1],k1,m1);
[zeta4,R4,I4,k4,m4,modal_vector4,fn_curve4] = modal_circle_fit(Xj(4,:),freq,[4 1],k1,m1);

ii10

modal_vector=[
modal_vector1
modal_vector2
modal_vector3
modal_vector4
]

zeta_temp=[
zeta1
zeta2
zeta3
zeta4
];

fn_curve_temp=[fn_curve1
    fn_curve2
    fn_curve3
    fn_curve4];

zeta=mean(zeta_temp);
fn_curve=mean(fn_curve_temp)
zeta*100

end

functionファイル

function [zeta,R,I,k,m,modal_vector,fn_curve]=modal_circle_fit(Xj,freq,out_in,modalk,modalm)
% % Xj コンプライアンス
% % out_in 入力と応答のid、例えば1点応答-2番入力の場合はout_in=[1,2]
df=freq(2) - freq(1);

xi = real(Xj(1,:));
yi = imag(Xj(1,:));

if nargin < 4
    modalk=0;
    modalm=0;
end


A = [sum(xi.^2) sum(xi.*yi) sum(xi)
    sum(xi.*yi) sum(yi.^2) sum(yi)
    sum(xi) sum(yi) length(xi)
    ]; %式(6.30)

B = [sum(xi.^3+xi.*yi.^2)
    sum(xi.^2.*yi+yi.^3)
    sum(xi.^2+yi.^2)
    ]; %式(6.30)

abc=A\B; %式(6.30)

figure
plot(real(Xj(1,:)),imag(Xj(1,:)))
hold on
grid on
axis equal
title([num2str(freq(1)) '～' num2str(freq(end)) 'Hzまで'])


plot(abc(1)/2,abc(2)/2,'ko') %式(6.27)

r=sqrt(abc(3) + (abc(1)^2+abc(2)^2)/4); %式(6.27)
for ii1=0:0.05:2*pi
    plot(abc(1)/2+r*cos(ii1),abc(2)/2+r*sin(ii1),'k.') %式(6.27)
end

% % % % モード円適合で算出された共振周波数の「円上の角度」
fn_theta=atan( abc(2)/abc(1) ); %「円上の角度」≒　共振周波数
if ( abc(2)/2+r*sin(fn_theta) )^2 < ( abc(2)/2+r*sin(-fn_theta) )^2
    fn_theta=-fn_theta;
end
plot(abc(1)/2+r*cos(fn_theta),abc(2)/2+r*sin(fn_theta),'ko')
Gmax_fn=( abc(1)/2+r*cos(fn_theta) ) + i*( abc(2)/2+r*sin(fn_theta) ); %　「円上の角度」のコンプライアンスの値（共振周波数におけるFRFの値）

% % % % 入力データ|Xj|maxの角度　　|G|maxの角度
[value,id_max]=max(abs(Xj(1,:)));
Gmax_theta = atan( imag(Xj(1,id_max))/real(Xj(1,id_max)) );

% % % % 「円上の角度」と隣り合うデータの角度を求める（片方は|G|maxの角度）
if Gmax_theta >= fn_theta
    Gmax_theta_next = atan( imag(Xj(1,id_max+1))/real(Xj(1,id_max+1)) );
    id_next=id_max+1;
else
    Gmax_theta_next = atan( imag(Xj(1,id_max-1))/real(Xj(1,id_max-1)) );
    id_next=id_max-1;
end

% % % % 「円上の角度」と隣り合う角度を求める
dtheta1 = abs( Gmax_theta - fn_theta );
dtheta2 = abs( Gmax_theta - Gmax_theta_next );
% % % % 「円上の角度」から求めた共振周波数（線形補完しているので、共振周波数には誤差があるが、一般的には許容範囲内で収まるはず）
fn_curve=(dtheta1/(dtheta2+dtheta1))*df + freq(id_max);

% % % % atanで角度を求めているので、座標上の第三象限では正負の符号がおかしくなるので、調整する必要がある
if Gmax_theta>0; Gmax_theta = -Gmax_theta; end
if fn_theta>0; fn_theta = -fn_theta; end
if Gmax_theta_next>0; Gmax_theta_next = -Gmax_theta_next; end

theta_d = abs( Gmax_theta - fn_theta );
theta_e = abs( Gmax_theta_next - fn_theta );


beta_d = freq(id_max) / fn_curve; %式(B9.2)
beta_e = freq(id_next) / fn_curve; %式(B9.2)

zeta = abs(beta_e^2 - beta_d^2) / ( 2*(beta_d*tan(theta_d)+beta_e*tan(theta_e)) ); %式(B9.5)

% % % 式(6.25)に代入
R=abc(1)/2;
K=1/(2*zeta)*sqrt( 1/( 4*(abc(3)+R^2) + abc(2)^2) ); % FRFが自己コンプライアンスの場合はk=K;
I=1/2*(abc(2)+1/(2*K*zeta));

% out_in,modalk)
if out_in(1)==out_in(2)
    modal_vector=1;
    k=K;
    m=k/(2*pi*fn_curve)^2;
else
    Gmax_fn
%     Gmax_fn=Gmax_fn .* imag(Gmax_fn)./abs(Gmax_fn) 
    modal_vector=2*modalk*zeta*imag(-Gmax_fn);
    k=modalk;
    m=modalm;
end

% [Gis,id]=max(  abs( imag(Xj(:,:)) ) ,[],2  );
% Gis=Gis .* imag(-Xj(:,id(1)))./abs(Xj(:,id(1))) ;
% phys=2*k_hanchi*zeta_hanchi*Gis;