1自由度振動の数値シミュレーションについて書こうかな

この記事は、WMMC Advent Calendar 2022 - Adventar12日目の記事です。

こんばんは。kuroです。

皆様いかがお過ごしでしょうか？

昨日の記事は、pizzagatakasugiさんのデータサイエンスコンペティションでした。

データコンペ反省会 - pizzagatakasugi’s blog (hatenablog.com)

早稲田が、主催しているデータサイエンスコンペティションがあるとは、私も知らなかったです。機械工学科にいると、データサイエンス系のお話は、あまり耳に入ってこないので、新鮮な気持ちで、記事を読ませていただきました。

さて、皆さんが、モノづくり系の記事をバンバン書いていらっしゃいますが、その流れに逆行し、理論系？？の記事を書こうと思います。具体的には、1自由度系振動の数値シミュレーションについて書いていきます。よって、マウス全く関係ないです。あと、記事が堅苦しいので、苦手な人は、ブラウザバックを。

なーぜ振動学なんじゃという方がいらっしゃると思うのですが、去年のように狂った？記事を書く気力はもう残っていないので、老害として、多少まともな記事を残しとこうかなと思ったからです（後、私の備忘録も兼ねてます....）。

機械工学科の方ならば、振動学は多少目にするでしょう。ただ、講義の振動学は、数式ゴリゴリにいじって、いまいち何やってるかわからない人が多いかと思います（機械工学科じゃない人は、こんな現象もあるんだなという感じで見て頂けると幸いです。）

なので、数値計算を通して、力学的な振動現象のイメージを持っていただけたらなと思います。なお、簡単のため、厳密に正しくない表現等を用いることがあるので、そこだけご留意を。

1自由度系の減衰ありの自由振動の運動方程式の立式および解
数値シミュレーションについて
1自由度減衰あり自由振動のシミュレーション
1自由度強制振動系のシミュレーション
最後に

1自由度系の減衰ありの自由振動の運動方程式の立式および解

　　　図1のような力学モデルを扱います。

　　　　　　　　　図1 減衰ありの1自由度振動系の力学モデル

図1を見ると、ばねとおもりとダンパーがついていることが分かります。

「自由度」とは、簡単に言うと、いくつの方向に動けるかみたいなものです。厳密には正しくないですが、おもりの数と思ってください。図1は重りが1個（厳密には、1次元的な方向のみ運動する）なので1自由度振動系と呼ばれます。

　機械科じゃない方のために、ばねとダンパの機能を簡単に説明します。ばねは、与えられた変位（位置の変化）に比例した復元力と呼ばれる力を変位の方向と逆におもりに与えます。また、ダンパは、振動を抑える役割を担い、おもりの速度に比例した粘性力をおもりに与え、振動を抑制させます。

　ダンパが無いモデルでしたら、高校物理で見かける単振動モデルとなりますが、それじゃつまらないので、ダンパをつけてみました。

　さて、高校物理でおなじみ、運動の解析をするときは、運動方程式を立てます。おもりの質量を $M$ , ばね定数を $k$ , 減衰係数を $c$ と仮定します。減衰係数とは、ばね定数のダンパ版だと思ってください。

運動方程式は、質量と力が必要なので、図1のおもりにかかっている力を抜き出します。すると、おもりには図2のような力がかかります。

　　　　　　　　　　　　　図2 おもりにかかる力の様子　　

　　図2より，運動方程式は、以下のように表されます。

$M\dfrac {d^2x}{dt^2}+c\dfrac {dx}{dt}+kx=0 \tag{1}$

ここで、解析のため、式(1)の両辺を $M$ で割り、式を無次元化します。 $\zeta=\dfrac{c}{\sqrt{2Mk}}$ , $\omega=\sqrt{\dfrac{k}{M}}$ と仮定すると、式(1)は、

　　　　　 $\dfrac {d^2x}{dt^2}+2\zeta\omega\dfrac {dx}{dt}+\omega^2x=0 \tag{2}$

　　ここで、あまり数学的なことに首を突っ込んでも、発狂必至なので、式(2)の形の方程式を解く解き方は、微分方程式等の教科書を参照してください。簡単に説明すると、式(2)の解を $x=e^{{\lambda}t}$ と置いて、式(2)に代入し、指数関数部を消去後、

$\lambda$ の二次方程式を、 $\zeta$ について場合分けして、 $\lambda$ について解きます。この二次方程式のことを特性方程式と呼びます。なお、式(2)のような形の方程式の場合、特定方程式の判別式を解けばわかるはずですが、 $0\lt\zeta\lt1$ , $\zeta=1$ , $\zeta\gt1$ の3つの場合で、式(2)の一般解が変化します。今回は、数値計算の対象とする $0\lt\zeta\lt1$ の解のみを式(3)に示しておきます。ただし、 $C_1$ および $C_2$ は任意の定数とします。

　　　　 $x(t)=exp(-\zeta\omega t)(C_1cos(\omega\sqrt{1-\zeta^2}t+C_2 sin(\omega\sqrt{1-\zeta^2}t))\tag{3}$

式(3)を見ると $t$ を無限大に飛ばすと、変位が0になっていく様子が分かります。

数値シミュレーションについて

　数値シミュレーションとは、微分方程式を「数値的」に解くことです。コンピュータは、解析的な積分なんぞ出来ません。よって、コンピュータは微分方程式を解くときは、被積分関数を離散的に積分する「数値積分」を用いて解きます。また、数値計算ソフトとして、MATLABを今回は使用します。

　数値積分の手法は、様々な種類がありますが、有名なものとして、「オイラー法」と「4次のルンゲクッタ法」等があります。「4次のルンゲクッタ法」は、精度の高い数値積分を実現しますが、少し説明が難しいので、今回は、簡便な手法である「オイラー法」を用いて、振動現象を解析していきたいと思います。

　オイラー法は、簡単に説明すると、テーラ展開を利用して、ある時刻 $t$ から、 $\Delta t$ ずれた時刻の $x(t+\Delta t)$ を計算することを繰り返し、ドミノ倒しのように積分値を求める方法です。

　具体例を通して説明します。式(4)のような微分方程式（初期条件 $x(0)=1$ ）をオイラー法で解くことを考えます。

　　　　　　　　　　　　 $\dfrac{dx}{dt}=x\tag{4}$

MATLABのソースコードですが、見にくいかもしれませんがご容赦を。

t_end=10;      %シミュレーション時間
    delta_t=0.0001;    %数値計算の時間刻みΔt(解析的な微分ではΔt→0)
   


    n=t_end/delta_t;    %2次元配列の行数を指定

    x=zeros(n,1);       %数値解における変位を格納する配列
   
    x_anal=zeros(n,1);  %解析解における変位を格納する配列
    
    
    t=zeros(n,1);       %時刻を格納する配列

    x(1,1)=1;           %初期変位
    

    for i=1:n-1     %オイラー法による微分方程式の数値計算
        %下の行にて1次の項までのテーラー展開にて、x(t+Δt)を計算
        x(i+1,1)=x(i,1)+x(i,1)*delta_t;
       
        t(i+1,1)=t(i,1)+delta_t;        

    end

    for i=1:n       %解析的な解を格納する

        x_anal(i,1)=exp(t(i,1));

    end

    plot(t(1:n,1),x(1:n,1));        %数値解のプロット
    hold on                         %グラフの重ね合わせコマンド
    plot(t(1:n,1),x_anal(1:n,1));        %解析解のプロット
    legend("数値解","解析解");      %凡例の追加
    xlabel("時刻(s)");
    ylabel("変位(m)")

このコードを実行し、数値解と解析解を図3のように重ね合わせます。

図3を見ると、式(4)の解析解と数値解がほとんど同じであることが見て取れます。本当は、数値計算の妥当性の検証と称して、Δtを細かくして、数値解の収束性を見るのですが、研究ではないので省略します。

1自由度減衰あり自由振動のシミュレーション

　さて、前置きはこれくらいにして、式(2)の数値シミュレーションに移りましょう。図1のモデルにおいて、初期変位 $x(0)=1$ , $v(0)=0$ として、オイラー法を用いて計算を行います。コードを以下のように作成します。なお、オイラー法では、2階の微分方程式を直接計算することはできないので、式(2)を以下のように、1階の連立微分方程式に変形します。

$\dfrac{dx}{dt}=v\\ \dfrac{dv}{dt}=-2\zeta\omega v-\omega^2 x \tag{5}$

式(5)の連立一階微分方程式を解くコードは以下の通りです。

t_end=100;   %シミュレーション時間
 delta_t=0.0001; %数値計算の時間刻みΔt(解析的な微分ではΔt→0)
 M=1;     %おもりの質量
 k=1;    %ばね定数
 c=0.2;   %減衰係数。減衰比が1以下となるように調整

 teta=c/(2*sqrt(M*k));   %減衰比
 omega=sqrt(k/M);     %系の固有角振動数


 n=t_end/delta_t; %2次元配列の行数を指定

 x=zeros(n,1);    %数値解における変位を格納する配列
 v=zeros(n,1);    %数値解における速度を格納する変数
 x_anal=zeros(n,1);  %解析解における変位を格納する配列
 
 
 t=zeros(n,1);    %時刻を格納する配列

 x(1,1)=1;     %初期変位
 v(1,1)=0;           %初期速度

    for i=1:n-1     %オイラー法による微分方程式の数値計算
        %下の行にて1次の項までのテーラー展開にて、x(t+Δt)を計算
        x(i+1,1)=x(i,1)+v(i,1)*delta_t;
        v(i+1,1)=v(i,1)+(-2*teta*omega*v(i,1)-(omega^2)*x(i,1))*delta_t;
        t(i+1,1)=t(i,1)+delta_t;        

    end

    for i=1:n       %解析的な解を格納する

        x_anal(i,1)=exp(-teta*omega*t(i,1))*(cos(omega*sqrt(1-teta^2)*t(i,1))+(teta/(sqrt(1-teta^2)))*sin(omega*sqrt(1-teta^2)*t(i,1)));

    end

    plot(t(1:n,1),x(1:n,1));        %数値解のプロット
    hold on                         %グラフの重ね合わせコマンド
    plot(t(1:n,1),x_anal(1:n,1));        %解析解のプロット
    legend("数値解","解析解");      %凡例の追加
    xlabel("時刻(s)");
    ylabel("変位(m)")

このコードで、シミュレーションを回した結果は図4の通りです。

図4を見ると、質点が振動しながら減衰している様子がわかります。これが、俗にいう不足減衰です。この振動の周期ですが、 $\zeta$ の値を変えると、少し変わるはずです。気になる人は、やってみてください（見た目では、わからないレベルかもしれませんが）。

1自由度強制振動系のシミュレーション

ここで、図1の質点に正弦外力を加えてみます。

すると式(2)のようになっていた運動方程式が外力が加わったことにより式(6)のように変化します。

　 $\dfrac {d^2x}{dt^2}+2\zeta\omega\dfrac {dx}{dt}+\omega^2x=\dfrac{f_0}{M} sin(\omega _R t) \tag{6}$

式(6)を解くのは、本当に煩雑なので、ここでは省略します。このプロセスは、微分方程式や振動工学の参考書に乗ってるはずです。

数値計算コードを乗せておきます。なお、自由振動のシミュレーションの時と、ばねマスダンパのパラメータは変えずに、 $f_0=1$ N, $\omega_R=2$ rad/sとして、シミュレーションを行いました。なお、簡単のため、 $x(0)=0$ , $v(0)=0$ としてシミュレーションを行います。

   t_end=100;   %シミュレーション時間
 delta_t=0.0001; %数値計算の時間刻みΔt(解析的な微分ではΔt→0)
 M=1;     %おもりの質量
 k=1;    %ばね定数
 c=0.2;   %減衰比が1以下となるように調整
 f_0=0.1;    %外力振幅
 omega_R=2; %外力角周波数

 teta=c/(2*sqrt(M*k));   %減衰比
 omega=sqrt(k/M);     %系の固有角振動数

 %%下4行は解析解を構成するために必要な変数
 A=-(f_0/M)*(2*teta*omega_R*omega)/(4*(teta^2*omega_R^2*omega^2)+(omega^2-omega_R^2)^2);
 B=(f_0/M)*((omega^2-omega_R^2))/(4*(teta^2*omega_R^2*omega^2)+(omega^2-omega_R^2)^2);
 C_1=-A;
 C_2=(-teta*omega*A-B*omega_R)/(omega*sqrt(1-teta^2));




 n=t_end/delta_t; %2次元配列の行数を指定

 x=zeros(n,1);    %数値解における変位を格納する配列
 v=zeros(n,1);    %数値解における速度を格納する変数
 x_anal=zeros(n,1);  %解析解における変位を格納する配列
 
 
 t=zeros(n,1);    %時刻を格納する配列

 x(1,1)=0;     %初期変位
 v(1,1)=0;           %初期速度

    for i=1:n-1     %オイラー法による微分方程式の数値計算
        %下の行にて1次の項までのテーラー展開にて、x(t+Δt)を計算
        x(i+1,1)=x(i,1)+v(i,1)*delta_t;
        v(i+1,1)=v(i,1)+(-2*teta*omega*v(i,1)-(omega^2)*x(i,1)+(f_0/M)*sin(omega_R*t(i,1)))*delta_t;
        t(i+1,1)=t(i,1)+delta_t;        

    end

    for i=1:n       %解析的な解を格納する

        x_anal(i,1)=exp(-teta*omega*t(i,1))*(C_1*cos(omega*sqrt(1-teta^2)*t(i,1))+C_2*sin(omega*sqrt(1-teta^2)*t(i,1)))+A*cos(omega_R*t(i,1))+B*sin(omega_R*t(i,1));

    end

    plot(t(1:n,1),x_anal(1:n,1));        %解析解のプロット
    hold on                         %グラフの重ね合わせコマンド
    plot(t(1:n,1),x(1:n,1));        %数値解のプロット
    legend("解析解","数値解");      %凡例の追加
    xlabel("時刻(s)");
    ylabel("変位(m)");

このコードを実行し、解析解と数値解を重ね合わせます。

　ここで、図5を眺めてみます。私の受けた振動学の授業では、式(6)の解の形として、正弦波の形を仮定していました。それは、ダンパが付いており、十分時間が経てば正しいのです。しかし、図5を見ると最初から、正弦波になるわけではないことが分かります。むしろ、最初の数十秒の間は、波形が乱れていることが分かります。これは、式(6)の初期条件により決まる解（過渡解）の影響といえます。過渡解とは式(6)の右辺を0にした時の解です（今回は、式(3)と同じ形をしています）。

　このことは、微分方程式の教科書に載っていることですが、式(6)のような非斉次微分方程式の解は、非斉次微分方程式の右辺を0にした時の一般解と、非斉次方程式そのものの特殊解の重ね合わせにより一般解が定まります。

　そして、今回の条件の場合、過渡解が式(3)の形をしており、 $t$ を無限大に飛ばすと0ヘ収束します。このことから、図5において、十分時間が経過した定常状態においては、波形は、正弦波に限りなく近くなります。しかしながら、強制振動において、最初から、波形が正弦波になるとは、限らないのです。

　私は、振動学を学んでいるときに、過渡解の存在を忘れて一回沼りました。私の学科の機械工学科の方は、授業の時、あまり過渡解を意識することがありませんでしたので、授業を理解するという観点で、過渡解の意識はあまりしなくて良いと思います。

ただ、研究で振動学を扱う際には、過渡解の存在を意識すると、振動の世界がだいぶ違って見えると思います。言われれば、当たり前のことなのですが、解析の時に存在を忘れがちなので、覚えておくと、振動現象の本質的な部分が見えやすくなると思います。

最後に

　ここまで、長々と書いて参りましたが、お付き合いいただきありがとうございました。面白みもない記事だったと思いますが、機械工学科の方は、振動学に対するイメージ、それ以外の方は、教養の力学で扱った振動現象が実際には、どのような現象なのかということを、わかっていただけたら幸いです。

　次回は、太郎さんの、「HM-StarterKitについて(初学者向け)」です。

例のマウスの学習キットについてですかね？私は、標準マウスを作成するときに、はんだ付けが下手すぎて、動かなかったので、そういうキットを使うのは、初学者の方ならありかもしれませんね。

それでは、また。