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

はじめに

振動工学やモード解析をご存じの方であれば、一度は読んだことがあると言われているバイブル的な書籍「モード解析入門」について再復習しております。

[商品価格に関しましては、リンクが作成された時点と現時点で情報が変更されている場合がございます。]

【中古】 モード解析入門 長松昭男 【古本・古書】
価格:9400円(税込、送料別) (2021/10/13時点)

 

再復習で学んだことをMATLABコードを示しながら解説したいと思います。

今回は2.6節の「周波数応答関数」のP.67から復習します。

まだ過去の記事を読んでいない方、この章や節以外を復習したい方は下記記事を参照ください。下記記事が目次となっており、各章や節にアクセスできるようにしています。

長松昭男著「モード解析入門」のMATLABコード付き解説
長松昭男先生の「モード解析入門」をMATLABコード付きで解説しています。

周波数応答関数とは

入力と出力の比を伝達関数と呼びます。
この伝達関数は時間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

コメント