MATLABで学ぶ振動工学　質量付加による伝達関数の変化

2021.05.03


スポンサーリンク

はじめに

久しぶりに振動工学についての記事を更新します。
本当は「構造最適化」について説明しようと思っていました。
しかし、いきなり構造最適化を扱うよりも基礎的な内容を扱った方が良いと思い改まって、本日は伝達関数合成法について説明します。

質量付加による伝達関数の変化

例題として、構造体の点hに力Fが作用しており、点jの変位振動Xの関係を伝達関数Gjhで表せるとする。そこに図1のように点iに質量mが付加されたとする。このとき質量付加により伝達関数は下式のようになる。ちなみに、理論式は下記を参考にしています。

 

$$
newG_{jh}=G_{jh}-\frac{G_{ji}G_{ih}}{G_{ii}-\frac{1}{ω^2 m}}
$$

 

図1　質量付加による伝達関数の変化

スポンサーリンク

例題による検証

では、2章の理論が正しいかを例題を通して検証したいと思います。

○例題
図2のような4自由度のバネマス系でm2に力Fが作用し、m3に質量mが付加される。このときの質量付加前後の伝達関数G_4,2を求めなさい。ここで、質量はm1=m2=m3=m4=1 [kg]、m=5 [kg]、バネ定数はk1=k2=k3=k4=10^3とする。

図2　例題

さて、MATLABで解いて2章が正しいかを検証してみましょう。

 

スポンサーリンク

 

結果はこのようになります。

図3　検証結果

質量付加後の計算結果（2章の理論式と直接法で解いた結果）が一致しています。2章の結果は正しいことがわかります。
これがわかると何が良いって？
これがわかると、質量付加で構造最適化やロバスト最適化をする方法を検討できるようになります。
まあ詳しくは後日説明したいと思います。

 

今日はこの辺でGood luck!

MATLABプログラム

○実行コード

clear all
close all
% % % % 初期設定
k_vec=[1 1 1 1]*10^3;
[K]=eval_Kmatrix(k_vec);
freq=0.5:0.1:20;
f=[0;1;0;0];
C=zeros(length(k_vec));

% % % % 質量付加前の質量行列
m_vec=[1 1 1 1];
[M]=eval_Mmatrix(m_vec);
before_x=eval_direct_x(M,K,C,f,freq);
% % % % % % 
figure(1)
loglog(freq,abs(before_x(4,:)),'k') %力fの成分が1なので、xを伝達関数Gは同じになるのでxを表示
hold on

% % % % % % 質量付加後の質量行列
m=5;
m_vec2=[1 1 1+m 1];
[M2]=eval_Mmatrix(m_vec2);
after_x=eval_direct_x(M2,K,C,f,freq);
% % % % % % 
figure(1)
loglog(freq,abs(after_x(4,:)),'b')

% % % % % % 2章の理論式（質量付加による伝達関数の変化）
% % % 相反定理を使えば
before_x2=eval_direct_x(M,K,C,[0;0;1;0],freq);
newG=before_x(4,:)-( before_x2(4,:).*before_x(3,:) )./ (before_x2(3,:)-1./((2*pi*freq).^2*m));
figure(1)
loglog(freq,abs(newG),'r--')
xlim([freq(1) freq(end)])
legend('質量付加前','質量付加後（直接法で解いた）','質量付加後（2章の理論式）')
xlabel('周波数 Hz')
ylabel('伝達関数 m/N')

○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

function x=eval_direct_x(M,K,C,f,freq)

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

end

 

コメント