MATLABで学ぶ振動工学 周波数応答関数

2021.12.19

はじめに

振動工学やモード解析をご存じの方であれば、一度は読んだことがあると言われているバイブル的な書籍「モード解析入門」について再復習しております。今回はモード解析入門の2.6節の「周波数応答関数」のP.67から復習します。なお、振動工学の参考書はどの本も似たり寄ったりで内容が重複しています。例えば、「実用モード解析入門」は長松先生の息子さんが著者ということもあり、重複している部分が多いです。

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

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

周波数応答関数とは

入力と出力の比を伝達関数と呼びます。
この伝達関数は時間tの関数であることや角周波数ωの関数である場合など、様々な形態があります。
周波数応答関数とは、角周波数ωや周波数freqを変数とした伝達関数を指します。

周波数応答関数の種類

振動や音響で扱う周波数応答関数には種類があります。
以下に一覧表を示します。なお、振動工学では変位／力（コンプライアンス）を基準として理論展開することが多く、表2.3の関係式の項目でもコンプライアンスGを基準としたときの比較をしています。

表2.3 周波数応答関数の種類

種類 名称 関係式 単位(SI)
変位/力 コンプライアンス
レセプタンス
アドミッタンス
動柔性		 G m/N
速度/力 モビリティ jωG m/(Ns)
加速度/力 アクセレランス
イナータンス		 ーω2G m/(Ns2)
力/変位 動剛性 1/G N/m
力/速度 機械インピーダンス ーj/(ωG) Ns/m
力/加速度 動質量 ーj/(ω2G) Ns2/m

先ほど、振動工学の理論展開ではコンプライアンスを用いられることが多いと述べました。
しかし、実測では入力をインパクトハンマー（力）として、応答を加速度センサーで測定することが多く、実用上ではアクセレランスを使うことが多いです。

将来的に説明するモード同定で混乱を生じることがあります。モード同定は周波数応答関数から共振周波数、モード減衰、モードベクトルを特定します。しかし、理論ではコンプライアンスを用いられていることが多く、実測の周波数応答関数をコンプライアンスに変形する必要があります。
コンプライアンスへの変形は表2.3からわかるかと思いますが、アクセレランスからコンプライアンスに変形する場合は、アクセレランス÷(－ω2)で求まります。

ちなみに、コンプライアンスGは式(2.125)です。

$$
G(ω) = frac{X(ω)}{F} =frac{1/k}{1-β^2+2jζβ} \tag{2.125}
$$

図示・プロットの種類

周波数応答関数は大きく分けて下記の3種の図示方法があります。

1つ目：ボード線図
2つ目：コクアド線図
3つ目：ナイキスト線図

一般的にはボード線図が使われることが多いです。
この理由は、現代では実測の測定周波数が広いため、測定データに多数の共振周波数が含まれることが理由かと思います。
もしナイキスト線図で表すと、図がぐちゃぐちゃになります。
コクアド線図だと対数目盛が扱えないので、値の大きい周波数帯域しか比較できなくなってしまいます。

まあ、上記のような説明をしてもわからないかと思いますので、MATLABで図示してみましょうか。

図　ボード線図

図　コクアド線図

 

図　ナイキスト線図

 

どうですか？
見やすさは個人の印象によるかと思いますが、ナイキスト線図は見にくそうですよね。

 

今日はこのへんでGood luck

MATLABプログラム

メインコード（実行ファイル）

 

clear all;close all;clc

freq=1:1000;
w=2*pi*freq;
m_vec=logspace(1,4,100);
% k_vec=logspace(5,13,100);
k_vec=ones(1,100)*10^8;

[M]=eval_Mmatrix(m_vec);
[K]=eval_Kmatrix(k_vec);
[C]=eval_Kmatrix(k_vec)*10^-5;
F=zeros(length(m_vec),1); F(1)=1;

x=zeros(length(m_vec),length(freq));
for ii1=1:length(freq)
    x(:,ii1)=inv(-w(ii1)^2*M+j*w(ii1)*C+K)*F;
end



% % % ボード線図
figure
subplot(211)
semilogy(freq,abs(x(1,:)),'k')
hold on
grid on
title('ボード線図')
legend('x1/F')
ylabel('応答倍率 [m/N]')
xlabel('周波数 [Hz]')

subplot(212)
plot(freq,angle(x(1,:)),'k')
hold on
grid on
legend('x1/F')
ylabel('位相 [rad]')
xlabel('周波数 [Hz]')


% % % コクアド線図
figure
subplot(211)
plot(freq,real(x(1,:)),'k')
hold on
grid on
title('コクアド線図')
legend('x1/F')
ylabel('応答倍率の実部 [m/N]')
xlabel('周波数 [Hz]')

subplot(212)
plot(freq,imag(x(1,:)),'k')
hold on
grid on
legend('x1/F')
ylabel('応答倍率の虚部 [m/N]')
xlabel('周波数 [Hz]')

% % % ナイキスト線図
figure
plot(real(x(1,:)),imag(x(1,:)),'k')
hold on
grid on
title('ナイキスト線図')
legend('x1/F')
ylabel('応答倍率の虚部 [m/N]')
xlabel('応答倍率の実部 [Hz]')

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

コメント