Zenn
🌊

MATLABで学ぶ振動工学　N自由度モデル（モード解析）2

2021/03/14に公開
#
matlab
tech

はじめに

前回の記事を読んでない方はこちらを見てから、本記事を読んでいただけると幸いです。
https://zenn.dev/mech_eng_uni/articles/e37fc1fb154ce9

前回の記事ではMK行列（質量行列と剛性行列）のプログラムを紹介しました。
今回は100自由度モデルを直接法で解く方法と、モード法で解く方法について紹介します。

理論

直接法

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

(-w^2 [M] + [K]) (x) = (F)
(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')

結果

どちらの結果も同じ結果が得られます。

「同じ結果が得られるなら直接法で良くない？」って思う方もいると思います。
そうですね。
このプログラムだけだとモード法のメリットが伝わりにくいですよね。
次回以降にその辺について解説する予定なので、少しお待ちください。

今回はこのへんでGood luck

私の個人ブログもよろしくお願いします。
https://mech-eng-uni.com/

Discussion