MATLABで学ぶ振動工学減衰系の自由振動 - MATLABパイセンが教える振動・騒音・音響・機械工学

はじめに
1. 関連

はじめに

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

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

参考書解説一覧

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

運動方程式

前回紹介しておりますが、運動方程式は式(2.4)に示します。

$m\ddot{x(t)} + c\dot{x(t)} +kx(t) = 0 \tag{2.4}$

式(2.4)を満たす解は式(2.42)です。

$x =x_0 e^{λt} \tag{2.42}$

式(2.42)を式(2.4)に代入すると式(2.45)が得られます。

$\left\{ \begin{array}{l} x_0&=0\\ mλ^2+cλ+k&=0 \end{array} \right. \tag{2.45}$

x0=0は系が静止していることを意味するので、振動した状態であるためには式(2.45)の下式が成立している必要があります。式(2.45)の下式をλについて解くと、中学校で習う「解の公式」から式(2.46)が得られます。

$\begin{align*} x & = {-c \pm \sqrt{c^2-4mk} \over 2m} \\ & = -{c \over 2m}\pm \sqrt{({c \over 2m})^2-{k \over m}} \\ & = \sqrt{{k \over m}}({-c \over {2 \sqrt{mk}}}) \pm \sqrt{{k \over m}} \sqrt{{c^2 \over 4mk}-1} \\ \tag{2.46} \end{align*}$

振動工学ではおなじみですが、下記に置き換えます。Ω＝√(k/m)は減衰のない角共振振動数ですよね。なので、Ωを不減衰固有角振動数と呼びます。

$\left\{ \begin{array}{l} c_c=2 \sqrt{mk}=2m\sqrt{\frac{k}{m}}=2mΩ\\ ζ=c/c_c \end{array} \right. \tag{2.47}$

すると、式(2.46)は式(2.48)で書き換えられます。

$λ = -Ωζ \pm Ω\sqrt{ζ^2-1} \tag{2.48}$

$\left\{ \begin{array}{l} λ_1 = -Ωζ + Ω\sqrt{ζ^2-1} \\ λ_2 = -Ωζ – Ω\sqrt{ζ^2-1} \end{array} \right. \tag{2.48}$

式(2.48)の一般的な解は式(2.49)となる。ここで、X1とX2は未定係数です。

$x = X_1 e^{λ_1 t} + X_2 e^{λ_2 t} \tag{2.49}$

この式(2.49)は式(2.28)の平方根の中身の正負によって、大きく現象が変化するので、場合分けして考えることになります。

減衰自由振動

ζ<1の場合

ここで、下記のような量を導入します。

$\left\{ \begin{array}{l} σ = Ωζ \\ ω_d = Ω\sqrt{ 1 – ζ^2} \end{array} \right. \tag{2.51}$

すると、式(2.51)を式(2.48)に代入して式(2.52)が得られる。

$λ_1, λ_2 = -σ \pm jω_d \tag{2.52}$

式(2.52)を式(2.49)に代入して、式(2.53)および式(2.54)が得られます。

$x = e^{-σ t}(X_1 e^{ jω_d t} + X_2 e^{ -jω_d t} ) \tag{2.53}$
$x = e^{-σ t}( (X_1+X_2) \cos{ω_d t} + j (X_1-X_2) \sin{ω_d t} ) \tag{2.54}$

ここで、X1とX2を下記のように変形します。

$X_1 = (X_C + j X_S)/2 ,X_2 = (X_C – j X_S)/2 , \tan{φ} = \frac{X_S}{X_C} , x_0 = \sqrt{X_C^2 + X_S^2} \tag{2.56}$

式(2.56)を式(2.54)に代入すると式(2.60)になります。

$x = x_0 e^{-σ t}( \cos{ω_d t} \cos{φ} – \sin{ω_d t} \sin{φ} ) = x_0 e^{-σ t} \cos{( ω_d t +φ )} \tag{2.60}$

さあMATLABでプロットしてみましょう

clear all;clc;close all;
 
 
t=0:0.01:3;
k=1000;
m=1;
 
 
RGB=jet(11);
c=0;
for s=0:0.1:1 %c/cc;
    c=1+c;
Omega=sqrt(k/m);
wd=Omega*sqrt(1-s^2);
sigma=Omega*s;
phy=0;
xo=1;
y=xo*exp(-sigma*t).*cos(wd*t+phy);
plot(t,y,'color',RGB(c,:))
grid on
hold on
xlabel('Time [s]')
ylabel('Ampli. [m]')
end
s=0:0.1:1;
legend(split(num2str(s)));