🌊
MATLABで学ぶ振動工学 N自由度モデル(モード解析)2
はじめに
前回の記事を読んでない方はこちらを見てから、本記事を読んでいただけると幸いです。
前回の記事ではMK行列(質量行列と剛性行列)のプログラムを紹介しました。
今回は100自由度モデルを直接法で解く方法と、モード法で解く方法について紹介します。
理論
直接法
直接法は運動方程式から逆行列計算で直接的に変位xベクトルを求める方法です。
モード法
モード法の理論についても詳しく説明したいのですが、KaTeXでの式の書き方に慣れていないため、下式だけ示します。なお、全ての振動は振動モードの重ね合わせで表現できる、というのがモード法の理論で、モードの重ね合わせで変位xベクトルを求めることができます。
(モード法を用いることで、連成系が非連成系になるなどのメリットがあるのですが、その辺の議論は次回以降で記載したいと思います。)
プログラム
直接法
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
私の個人ブログもよろしくお願いします。
Discussion