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

はじめに

振動工学やモード解析をご存じの方であれば、一度は読んだことがあると言われているバイブル的な書籍「モード解析入門」について再復習しております。

[商品価格に関しましては、リンクが作成された時点と現時点で情報が変更されている場合がございます。]

【中古】 モード解析入門 長松昭男 【古本・古書】 価格:9400円(税込、送料別) (2021/10/13時点)

  再復習で学んだことをMATLABコードを示しながら解説したいと思います。 今回は3章の多自由度系、3.1節の「不減衰系の自由振動」のP.79から復習します。まだ過去の記事を読んでいない方、この章や節以外を復習したい方は下記記事を参照ください。下記記事が目次となっており、各章や節にアクセスできるようにしています。 

長松昭男著「モード解析入門」のMATLABコード付き解説
長松昭男先生の「モード解析入門」をMATLABコード付きで解説しています。

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は暫定的(適当)に定義しました。

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自由度モデルを解いてみます。

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は下記のようになります。

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

コメント