MATLABで学ぶ振動工学 減衰系の自由振動と強制振動

はじめに

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

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

参考書解説一覧

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

mech-eng-uni.com

2023.01.07

3.2 減衰系の自由振動

3.2.0 減衰の種類

振動工学を考える上で「減衰」は非常に重要な項目です。
しかし、「減衰」とGoogleなどで検索すると多数の種類があることがわかります。

あらかじめ「減衰」についての考え方や種類を頭に入れておくと、3.2.0項以降の説明がわかりやすくなるので、ここで説明します。

一般的に、外力が作用しない場合の振動（自由振動）は、振動している状態から徐々に振幅が減少し、静止状態へと推移します。この振幅が徐々に減少していく現象を「減衰」と呼びます。（もしくは「減衰が作用している」と呼びます。） この「減衰」をモデル化（定式化）する考え方として、比例粘性減衰、構造減衰、モード減衰などがあります。

ただし、実際の減衰は様々な種類が複合的に絡み合って生じるので、単一の減衰（例えば、構造減衰だけ）で現象を表現することには未だに課題があります。 しかし、現状ではこの課題を解決する術がなく、単一の減衰でモデル化されています。

3.2.1 運動方程式

運動方程式の細かい導出方法については、「モード解析入門」を参照ください。
ここでは、3.2.0項に関連する基本的な内容だけを説明します。

多自由度系の自由振動の運動方程式は下記になります。

$$[M](\ddot{x})+[C](\dot{x})+[K](x)=(0) \tag{3.57}$$

ここで、[ ]は行列を意味し、( )は縦ベクトルを意味します。
[C]は減衰行列ですが、[C]をどのようにモデル化するのかというのが本記事での取り上げたい内容です。

ちなみに、下記のように速度に比例した減衰を粘性減衰と呼びます。

$$[C](\dot{x})$$

一般粘性減衰として[C]自体の行列要素をモデル化することもできますが、それだと計算コストがかかります。 上記については「モード解析の歴史」について説明するときにでも解説したいと思います。

3.2.2 比例粘性減衰

比例粘性減衰とは、減衰行列[C]が質量行列[M]と剛性行列[K]に比例する成分で表現できるという仮定の基での減衰のモデル化方法です。 つまり、減衰行列[C]はαとβの係数を用いて式(3.63)のように表現できます。

$$[C]=α[M]+β[K] \tag{3.63}$$

ちなみに、有限要素法などの構造解析の分野ではレイリー減衰と呼ばれます。 なぜか振動工学の参考書ではレイリー減衰ではなく、比例粘性減衰と書かれているケースが多いですが、比例粘性減衰もレイリー減衰も同じです。

そして、β=0の場合を質量比例減衰と呼び、α=0の場合を剛性比例減衰と呼びます。 剛性比例減衰の最も有名なものが構造減衰（ヒステリシス減衰）です。構造減衰は下式でモデル化されます。なお、Gは構造減衰係数です。

$$[C](\dot{x})= iG[K](x)$$

3.2.2.1 モード減衰（モード減衰比）

この項目は参考文献の「モード解析入門」には記載がありませんが、項目を分けたほうがわかりやすいので、独立させて解説します。

モード法やモード解析とは

モード法やモード解析とは下記①～③で振動現象を表現する理論です。

　①振動モード（固有ベクトル）
　②共振周波数（固有値）
　③モード減衰比

振動モードや共振周波数がわからない方は下記や2章を参照ください。

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

長松昭男著「モード解析入門」の3.1節：不減衰系の自由振動についてMATLABのコードの示しつつ解説いたします。

mech-eng-uni.com

2021.12.26

3.1節で説明した「固有ベクトルの直交性」から、r次の固有ベクトルφrを質量行列Mや剛性行列Kに左右から乗じると、式(3.38)のようにmrやkrのようにある値になります。一方、左右に異なる時数の固有ベクトルを質量行列Mや剛性行列Kに乗じると0になります。

$$(φ_r)^T [M] (φ_r)= m_r,　(φ_r)^T [K] (φ_r)= k_r,　(r=1～N)  \tag{3.38}$$

このmrをモード質量、krをモード剛性と呼びます。このモード質量mrやモード剛性krはスカラーです。スカラーでは扱いにくいので、下記のように扱う参考書もあります。

$$[Φ]^T [M] [Φ]= [m_m],　[Φ]^T [K] [Φ]= [k_m],　\tag{3.38.2}$$

ここで、[Φ]はモードベクトル行列（もしくはモード行列）と呼び、モードベクトルを下記のように配列させたものを指します。

$$[Φ]= [(φ_1), (φ_2), …, (φ_N)],　(r=1～N) 　\tag{3.38.3}$$

上記の専門用語を理解した上でモード法を厳密に説明すると、「モード法はモード質量、モード剛性、固有値（共振周波数）、モード減衰比を使って計算をする方法」を意味します。

モード減衰比

モード減衰の考え方は、モード質量やモード剛性と同様に、r次の固有ベクトルφrを減衰行列Cに左右から乗じると、式(3.88)のように”モード次元の減衰を表す値”が求まるだろう！というものです。
この”モード次元の減衰を表す値”をモード減衰係数crとなると呼びます。（モード減衰係数crを略して「モード減衰」と呼ぶ人もいます。）

$$(φ_r)^T [C] (φ_r)= c_r　(r=1～N)  \tag{3.88}$$

（減衰行列[C]自体をモデル化することは非常に難しいので、モード減衰係数crを使って減衰行列[C]を使わないで解こう！という考え方とも言いかえれます）

では、このモード減衰係数crはどのように求めるのか？

まず、r次のモード質量mrとモード剛性krを用いて、式(3.89.1)のccrを導入します。

$$c_{cr}=2 \sqrt{m_r k_r}　(r=1～N)  \tag{3.89.1}$$

そして、モード減衰比ζrという値を導入して、式(3.89.2)のように定義します。

$$c_r=ζ_r c_{cr}　(r=1～N)  \tag{3.89.2}$$

モード質量やモード剛性は、質量行列や剛性行列を構成すれば求めることができます。
つまり、多自由度系でも理論モデルから導出することができます。また、有限要素法でも質量行列や剛性行列を構成することができるので、モード質量やモード剛性は求まります。

しかし、モード減衰比は理論的に求めることができません。そのため、モード減衰比は実験から求めるか、過去の経験やエイヤで決めることになります。

モード減衰比さえ理解できれば、いわゆる「モード減衰」は理解できたと言えます。ただ、補足情報としてモード減衰率を紹介しておきます。

r次のモード減衰率は下記のように表現できます。

$$σ_r=Ω_r ζ_r,　Ω_r=\sqrt{k_r/m_r}　(r=1～N)  \tag{3.93.1}$$

減衰が作用すると、共振周波数（固有振動数もしくは固有角振動数）は少し小さくなります。専門的にはこれを式(3.93.2)の減衰角振動数で表します。

$$ω_{dr}=Ω_r \sqrt{1-{ζ_r}^2}　(r=1～N)  \tag{3.93.2}$$

3.3 強制振動

詳しく述べませんでしたが、「減衰系の自由振動」の共振周波数と固有ベクトル（振動モード）は「不減衰系の自由振動」の結果と一致します。減衰が関係する「減衰角振動数」はモード減衰比さえ判明すれば求まります。そして、モード減衰比は実験値を使うか、過去の経験上で決めるか、エイヤで決めて計算することになります。そのため、「減衰自由振動系」では解くべき問題がありません。

一方、外力が作用する強制振動では、モード法で解くことのメリットが大きいです。

モード法でN自由度モデルを解く前に、重要な項目について説明します。

モード解析やモード法は「振動モードの重ね合わせで、全ての振動を表現できる」との仮定しています。ただし、（例えば）3自由度系のモデルでは3次までの振動モードが現れますが、この1～3次の振動モードがすべて均等に励起されるわけではありません。
（1次が33.3%、2次が33.3%、3次が33.3%の影響というわけではない）

つまり、振動変位や振動加速度を振動モードで表現する場合、式(3.45)のように振動モードの影響度を示す係数が必要となります。

$$(x)=ξ_1(φ_1)+ξ_2(φ_2)+ξ_3(φ_3)+….+ξ_r(φ_r)+…+ξ_N(φ_N)=[Φ](ξ)　 \tag{3.45}$$

上記の係数ξrは変位ベクトル(x)の中にr次の固有モード成分がどの程度含まれているかを示します。

式(3.45)を用いた周波数応答関数の定式を確認しましょう。
モード座標系の運動方程式は式(3.114)となります。ここで(f)は外力ベクトルを意味します。

$$m_r \ddot{ξ_r} + c_r \dot{ξ_r} + k_r ξ_r =(φ_{ri})F_i e^{jωt}　(r=1～N)  \tag{3.114}$$

ここで、φriはr次の固有ベクトル(φr)のi行目の要素です。上記から式(3.115)が得られます。

$$ξ_r = \frac{ (φ_{ri})F_i  }{ -m_r ω^2 + j c_r ω + k_r }e^{jωt}　(r=1～N)  \tag{3.115}$$

式(3.45)に式(3.115)を代入すると

$$(x) = \sum_{\substack{r=1}} \frac{ φ_{ri} F_i  }{ -m_r ω^2 + j c_r ω + k_r }(φ_r)e^{jωt}　(r=1～N)  \tag{3.116}$$

上記は点iに加振力が作用したときの変位ベクトル(x)を意味するが、(x)には全点の変位が含まれています。そこで、点jの応答だけを下式のように考えます。

$$x_j = (\sum_{\substack{r=1}} \frac{ φ_{ri} φ_{rj}  }{ -m_r ω^2 + j c_r ω + k_r })F_i e^{jωt}　(r=1～N)  \tag{3.116}$$

周波数応答関数G(ω)は下式となります。

$$G(ω)= x_j / F_i e^{jωt} = (\sum_{\substack{r=1}} \frac{ φ_{ri} φ_{rj}  }{ -m_r ω^2 + j c_r ω + k_r })　(r=1～N)  \tag{3.116}$$

加振力を1として変位ベクトル(x)を求めたいと思います。
ちなみに、加振力を1にすると周数応答関数G(ω)と変位の結果は一致します。もっと深く勉強したい方は、この記事のMATLABプログラムを参考に、周数応答関数を求めるプログラムへと修正してみてください。

MATLABプログラム(3.4節 数値例)

MATLABプログラムの結果

まずは直接法とモード法の比較結果を示します。

モード法も直接法も結果が一致していますね。ちなみに質量行列のM(1,1)を加振したときのM(1,1)の応答です。
直接法は計算（コード）が非常にシンプルなので、直接法の結果と一致していれば、モード法のコードも正しいことが推察できます。

では、ほかの応答点の結果も見てみましょう。

加振点や応答点がモードの節ではない場合、共振周波数で応答が卓越します。その様子が再現できていますね。

今回はこのへんでGood luck

MATLABコード

実行ファイル

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
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
clear all;clc;close all
 
freq=0.1:0.1:100; %計算周波数の定義
w=2*pi*freq; %角振動数
 
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); %固有振動数
 
mr=diag(V.'*M*V); %モード質量行列
kr=diag(V.'*K*V); %モード剛性行列
c_cr=2*sqrt(mr.*kr);
modal_dampimg=0.1*0.01; %モード減衰比（とりあえず0.1%と設定）
cr=c_cr*modal_dampimg; %モード減衰係数
 
F=zeros(length(m_vec),1); F(1)=1; %外力ベクトル(全周波数帯を1で加振する)
Xj=zeros(length(m_vec),length(freq)); %変位ベクトル（計算前にベクトルを定義してメモリを確保）
 
% % % % モード法
ii=1; %加振点 （加振点が複数ある場合はfor ii=1:kength(F)のループを追加すればよい ）
for ii1=1:length(m_vec) %応答点
    for ii2=1:length(wn)
        Xj(ii1,:)=Xj(ii1,:)+V(ii1,ii2)*V(ii,ii2)./( -mr(ii2).*w.^2 + 1i*cr(ii2)*w +kr(ii2) )*F(ii);
%       横ﾍﾞｸﾄﾙ　=横ﾍﾞｸﾄﾙ  +  ｽｶﾗｰ　　×ｽｶﾗｰ　　　./(　ｽｶﾗｰ×横ベクトル + ｽｶﾗｰ×横ベクトル +ｽｶﾗｰ）×ｽｶﾗｰ
    end
end
 
% % % % 直接法
X_direct=zeros(length(m_vec),length(freq));
C=(V.')\ diag(cr) /(V);
for ii1=1:length(freq)
    X_direct(:,ii1)=(-w(ii1)^2*M+1i*w(ii1)*C+K)\F;
end
 
 
figure
semilogy(freq,abs(Xj(1,:)),'k')
hold on
semilogy(freq,abs(X_direct(1,:)),'r--')
legend('モード法','直接法')
xlabel('周波数 Hz')
ylabel('変位 [m]')
xlim([0 25])
 
figure
semilogy(freq,abs(Xj(1,:)),'k')
hold on
semilogy(freq,abs(Xj(50,:)),'b')
semilogy(freq,abs(Xj(100,:)),'r')
% semilogy([fn fn].',[ones(100,1)*10^0 ones(100,1)*10^-10].','k--')
grid on
legend('x1','x50','x100')
xlabel('周波数 Hz')
ylabel('変位 [m]')
xlim([0 25])

functionファイル

1
2
3
4
5
function [M]=eval_Mmatrix(m_vec)
 
M=diag(m_vec);
 
end

1
2
3
4
5
6
7
8
9
10
11
12
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