Loading [MathJax]/extensions/tex2jax.js

MATLABで学ぶ振動工学 固有モードの直交性

2021.12.26

はじめに

振動工学やモード解析をご存じの方であれば、一度は読んだことがあると言われているバイブル的な書籍「モード解析入門」について再復習しております。今回はモード解析入門の3章の多自由度系、3.1節の「不減衰系の自由振動」のP.79から復習します。なお、振動工学の参考書はどの本も似たり寄ったりで内容が重複しています。例えば、「実用モード解析入門」は長松先生の息子さんが著者ということもあり、重複している部分が多いです。

各参考書で重複している部分は本記事で大まかに理解できると思いますので、他の参考書から本記事へアクセスした場合も最後まで読んでいただけると幸いです。

参考書解説一覧
MATLABとNVH(Noise, Vibration and Harshness)をこよなく愛する人生の先輩、それが私（MATLABパイセン）です。NVHを極めるにあたり、周辺の知識も網羅的に勉強することになり、その知識を共有すべく、本HPを運営しています。日本の製造業を応援すべく、機械エンジニアの「車輪の再開発」の防止と業務効率化の手助けをモットーに活動。専門である振動騒音工学や音響工学について、プログラムを示しながら解説。大学の授業よりもわかりやすく説明することを目指す。質問頂ければ、回答したいと思います。
mech-eng-uni.com
2023.01.07

3.1 不減衰系の自由振動

「モード解析入門」の運動方程式はかなり力を入れて解説されています。 長松先生の運動方程式や機械系物理学の思考の深さについての熱の入れ具合は、下記の文献を読めばよく理解できるかと思います。

  ただ、「モード解析入門」の運動方程式は丁寧なんですけど、 読み進めるには重たいんですよね。 なので、この記事では手短に解説したいと思います。

3.1.3 多自由度系

下記のような100自由度モデルを例として考えます。

上記モデルのような多自由度系での質量行列と剛性行列のMATLABプログラムを紹介します。 下記のm_vecおよびk_vecは下式を意味します。m_vecおよびk_vecを入力として、質量行列Mおよび剛性行列Kを出力します。

質量行列

質量行列Mは下のようになります。

MATLABプログラム(functionファイル)
function [M]=eval_Mmatrix(m_vec)

M=diag(m_vec);

end

剛性行列

剛性行列Kは下のようになります。

MATLABプログラム(functionファイル)
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

3.1.4 固有振動数と固有モード

固有振動数と固有モードは固有値解析すればわかります。固有値解析は下記のようにeigを使います。（eigsを使う場合もあります） ここで、m_vecとk_vecは暫定的（適当）に定義しました。

1
2
3
4
5
6
7
8
9
10
clear all
 
m_vec=ones(1,100);
k_vec=(1:100)*10^3;
[M]=eval_Mmatrix(m_vec);
[K]=eval_Kmatrix(k_vec);
 
[V,D]=eig(K,M); %固有値D、固有ベクトルV
wn=sqrt(diag(D)); %固有角振動数
fn=wn/(2*pi); %固有振動数

3.1.5 固有モード（固有ベクトル）の直交性

固有ベクトルは直交性があります。固有ベクトルは固有モードや振動モードなどとも呼ばれます。 そのため、（n次の固有ベクトル）×（n次の固有ベクトル）＝1となりますが、 （n次の固有ベクトル）×（i次の固有ベクトル）＝0になります。（n≠i） 教科書ではこの辺が詳しく説明されているのですが、イメージをつかむために、まずはMATLAB上で計算してみるのがよいでしょう。 例として下記の4自由度モデルを解いてみます。

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
clear all
 
m_vec=ones(1,4);
k_vec=(1:4)*10^3;
[M]=eval_Mmatrix(m_vec);
[K]=eval_Kmatrix(k_vec);
 
[V,D]=eig(K,M); %固有値D、固有ベクトルV
wn=sqrt(diag(D)); %固有角振動数
fn=wn/(2*pi); %固有振動数
 
 
V1=V(:,1); %1次の固有ベクトル
V2=V(:,2); %2次の固有ベクトル
V3=V(:,3); %3次の固有ベクトル
V4=V(:,4); %4次の固有ベクトル
 
MAC_1_1=V1.'*V1;
MAC_1_2=V1.'*V2;
MAC_1_3=V1.'*V3;
MAC_1_4=V1.'*V4;
MAC_2_1=V2.'*V1;
MAC_2_2=V2.'*V2;
 
MAC=[MAC_1_1
MAC_1_2
MAC_1_3
MAC_1_4
MAC_2_1
MAC_2_2
    ];

一部だけですが検証した結果をMACというベクトルに保存したので、MACの中身を見てみましょう。 MACは下記のようになります。

1
2
3
4
5
6
7
8
MAC =
 
    1.0000
    0.0000
    0.0000
    0.0000
    0.0000
    1.0000

（※実際には0ではなく、10^-17のような非常に小さい値です） 上記結果から、 （1次の固有ベクトル）×（1次の固有ベクトル）＝1、 （2次の固有ベクトル）×（2次の固有ベクトル）＝1、 それ以外は0となっていることが確認できます。 上記結果で固有ベクトルの直交性についてのイメージがつかめたかと思います。 「モード解析入門」の3.1節では、モード質量やモード剛性、モード正規固有モード、モード座標についての説明がありますが、これはモード法を説明するときに解説したいと思います。

今回はこのへんでGood luck

コメント