MATLABで学ぶ振動工学 N自由度モデル（モード解析）2 - MATLABパイセンが教える振動・騒音・音響・機械工学

はじめに
理論
プログラム
結果
1. 関連

はじめに

前回の記事を読んでない方はこちらを見てから、本記事を読んでいただけると幸いです。

https://mech-eng-uni.com/matlab%e3%81%a7%e5%ad%a6%e3%81%b6%e6%8c%af%e5%8b%95%e5%b7%a5%e5%ad%a6%e3%80%80n%e8%87%aa%e7%94%b1%e5%ba%a6%e3%83%a2%e3%83%87%e3%83%ab%ef%bc%88%e3%83%a2%e3%83%bc%e3%83%89%e8%a7%a3%e6%9e%90%ef%bc%891/

前回の記事ではMK行列（質量行列と剛性行列）のプログラムを紹介しました。
今回は100自由度モデルを直接法で解く方法と、モード法で解く方法について紹介します。ちなみに、私の愛読している長松先生著の”モード解析入門”を参考にしています。

理論

　　直接法

直接法は運動方程式から逆行列計算で直接的に変位xベクトルを求める方法です。

$$ (x) = (F) (-w^2 [M] + [K])^{-1} $$

　　モード法

モード法の理論についても詳しく説明したいのですが、KaTeXでの式の書き方に慣れていないため、下式だけ示します。なお、全ての振動は振動モードの重ね合わせで表現できる、というのがモード法の理論で、モードの重ね合わせで変位xベクトルを求めることができます。
（モード法を用いることで、連成系が非連成系になるなどのメリットがあるのですが、その辺の議論は次回以降で記載したいと思います。）

$$
(x) =\displaystyle\sum_{r=1}^n \frac{F \phi_{ri}\phi_{rl}}{m_r (w_r^2 – w^2)}
$$

プログラム

　　直接法

clear all;close all;clc

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

[M]=eval_Mmatrix(m_vec);
[K]=eval_Kmatrix(k_vec);
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+K)*F;
end

figure
semilogy(freq,abs(x(1,:)),'k')
hold on
semilogy(freq,abs(x(50,:)),'b')
semilogy(freq,abs(x(100,:)),'r')
grid on
legend('x1/F','x50/F','x100/F')

　　モード法

clear all;close all;clc

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

[M]=eval_Mmatrix(m_vec);
[K]=eval_Kmatrix(k_vec);
F=zeros(length(m_vec),1); F(1)=1;


[V,D]=eig(K,M);
wn=sqrt(diag(D));
fn=wn/(2*pi);

mr=diag(V.'*M*V);
kr=diag(V.'*K*V);

x2=zeros(length(m_vec),length(freq));
for ii1=1:length(wn)
    for ii2=1:length(m_vec)
        x2(ii2,:)=x2(ii2,:)+V(1,ii1)*V(ii2,ii1)./( mr(ii1).*(wn(ii1)^2-w.^2) );
    end
end

figure
semilogy(freq,abs(x2(1,:)),'k')
hold on
semilogy(freq,abs(x2(50,:)),'b')
semilogy(freq,abs(x2(100,:)),'r')
grid on
legend('x1/F','x50/F','x100/F')