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

スポンサーリンク

はじめに

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

前回の記事では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)}$$

スポンサーリンク

プログラム

　　直接法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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')

　　モード法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
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