動吸振器の問題を解く（前編）

#動吸振器とは

動吸振器というものをご存知だろうか

普通に生きてればなかなか目にしない。私も大学の制御工学の授業で始めてその名を聞いた。

こちらに目を通せば概要はつかめるはず。→動吸振器(ウィキペディア)

要は「振動させたくないものに追加でバネますダンパ系を取り付けて、振動を吸わせよう」
というコンセプトである。

目を通してくれた方は気づいたかもしれない.

今回解くの2自由度系の2階線形微分方程式である。（つまり連成振動）

この記事は動吸振器だけではなく、2自由度の振動問題を扱うのであればある程度参考になると思われる。

#解析
試しに、主系に周期外力がかかったときの、従系のふるまいを考えてみる。

こちらが動吸振器の力学的な計算モデルになる。

今回は主系に周期外力がかかったときの、従系のふるまいを考える。（2自由度強制振動）

###方程式の無次元化

※微分方程式の数学的な処理に興味ねーよ！という方は読み飛ばして頂いて構わない。
　式を扱いやすくするための処理である。

振動方程式を解く際によく用いられるのが「無次元化」という手法だ。

ざっくり説明すると、変数変換を使って、両辺に存在する変数の単位（ここでは、力（N）が基準）を消してやる。
また、振動してしまう主系に、どのような動吸振器を付けてやればいいのかを考えるので、主系のパラメータを基準に変数変換を行う。
「主系に対する従系の」質量・バネ剛性・減衰係数の比をRで書き換えた。

主系を基準にとった無次元化の流れと結果

どうだろうか？格段に見やすくなったと思われる。
余談だがこの式変形は実に気持ち良いので試してみてほしい。（特に合成関数の微分でｙの係数が1になるのはクセになる。）

#数値解析

今回使用するのは4次のルンゲ＝クッタ法。
微分方程式の数値解析の一般的な手法。調べると色々出てくるが、個人的にわかりやすかったサイトはこちら（シキノート様）

微分方程式の数値解析でよくある手法が「2階の微分方程式は1階の連立微分方程式にする」
というものだ。
今回は2階の連立微分方程式なので、計4本の式になる。

ルンゲ=クッタ法の計算部分については割愛。調べるといっぱい出てくる。

###データ構造
さて、今回調べたいのは**「質量比・バネ剛性比・減衰係数比を変化させた際の、周波数特性」**といったところか。

周波数特性を調べたい。ということで、2次元配列を準備する。

ある周波数に対して、時間をゼロから順にルンゲクッタの処理をしてやれば、時間応答がわかる。
振動が収束したところで次の周波数に移動する。
を繰り返すことで周波数特性がわかる（はず）なので、実際に実装してみよう。

#実装(Matlab)

それでは実装をするために、システムの入力と出力を整理する

この条件から、上で示した微分方程式の周波数特性がわかるはずだ。無次元化してあるので、変位x_mやx_aなどは最後に計算し直すとする。

##時間ループ

時間方向のループはどうすれば良いかを考える。
無次元化によって$t=τ/ω_m$の変換がされており、$τ$時間軸内での振動を考えなければいけない。
$t$が$0$~$2π$で変化する時、主系の固有周期$T$($=2π/ω_m$)の１周期分の時間が経過する。
今回は$t$を$0$~$2π$（１波長分）から10点選ぶ形で時間幅を決めてやることにする。
すなわち$dτ=T/10$である。

冗長に書いたが、サンプルコード中のdtau=0.1smの係数（0.1の部分）を大きくすれば時間幅が伸び、小さくすれば縮むだけの話である。

この$dτ$をTmax回forでループさせれば良い。

##周波数ループ

周波数ループも同様に考える。
今回は周期外力の項にのみ周波数の変数が現れている。
$sin{(ω/ω_m)τ}$の$ω$を変えることによって周期外力の周波数を変え、変位を観測する。
**$Rs = ω/ωm$**とすると、$Rs=1$の時に主系固有周波数となる。$Rs=2$は主系固有周波数の2倍である。

つまり、$Rs=0$~$2$として、fmaxで分割すれば周波数を$0$~$2ω_m$までシフトできるというわけだ。

##サンプルコード（パラメータ調整・ループ打ち切りなし）

とりあえず理論はわかったので、一旦コードを書いてみる。
ここでは微分方程式の数値計算が出来る事と、振動のモデリングが出来ることを確認する。

clear

RsLow =0; %"1" is Out of Phase , "0" show both
RsHigh =5; % Freq end

Tmax = 5000;
Fmax = 500;

Rm = 1;
Rc = 1;
Rk = [3,4,5];

mm = 10;
cm = 0.001;
km = 30;

zm = cm / (2*sqrt(mm*km));
sm = sqrt(km/mm);
Tm = 2*pi/ sm; % Period


RmN = length(Rm);
RcN = length(Rc);
RkN = length(Rk);

ym = zeros(Tmax,Fmax,RmN,RcN,RkN);
ya = zeros(Tmax,Fmax,RmN,RcN,RkN);
vm = zeros(Tmax,Fmax,RmN,RcN,RkN);
va = zeros(Tmax,Fmax,RmN,RcN,RkN);

% Loop preparation
dtau = 0.1 *sm;
Rs = RsLow:(RsHigh-RsLow)/Fmax:RsHigh;



%% Calculation

iRm = 1;
iRc = 1;
iRk = 1;

for j=1:Fmax
    for i=1:Tmax-1
        
        [ym(i+1,j,iRm,iRc,iRk),ya(i+1,j,iRm,iRc,iRk),vm(i+1,j,iRm,iRc,iRk),va(i+1,j,iRm,iRc,iRk)] ...
        = RungeKutta(ym(i,j,iRm,iRc,iRk),ya(i,j,iRm,iRc,iRk),vm(i,j,iRm,iRc,iRk),va(i,j,iRm,iRc,iRk),Rm(iRm),Rc(iRc),Rk(iRk),Rs(j),zm,dtau*i,dtau);
    end
end

%% Function

function [x1,x2,u1,u2]=RungeKutta(x1,x2,u1,u2,Rm,Rc,Rk,Rs,zm,tau,dtau)

    [Cx1_1,Cx2_1,Cu1_1,Cu2_1] = DifEqu(x1,x2,u1,u2,Rm,Rc,Rk,Rs,zm,tau);
    [Cx1_2,Cx2_2,Cu1_2,Cu2_2] = DifEqu(x1+Cx1_1*0.5*dtau,x2+Cx2_1*0.5*dtau,u1+Cu1_1*0.5*dtau,u2+Cu2_1*0.5*dtau,Rm,Rc,Rk,Rs,zm,tau+0.5*dtau);
    [Cx1_3,Cx2_3,Cu1_3,Cu2_3] = DifEqu(x1+Cx1_2*0.5*dtau,x2+Cx2_2*0.5*dtau,u1+Cu1_2*0.5*dtau,u2+Cu2_2*0.5*dtau,Rm,Rc,Rk,Rs,zm,tau+0.5*dtau);
    [Cx1_4,Cx2_4,Cu1_4,Cu2_4] = DifEqu(x1+Cx1_3*dtau,x2+Cx2_3*dtau,u1+Cu1_3*dtau,u2+Cu2_3*dtau,Rm,Rc,Rk,Rs,zm,tau+dtau);

    x1 = x1+dtau*(Cx1_1 + 2*Cx1_2 + 2*Cx1_3 + Cx1_4)/6;
    x2 = x2+dtau*(Cx2_1 + 2*Cx2_2 + 2*Cx2_3 + Cx2_4)/6;
    u1 = u1+dtau*(Cu1_1 + 2*Cu1_2 + 2*Cu1_3 + Cu1_4)/6;
    u2 = u2+dtau*(Cu2_1 + 2*Cu2_2 + 2*Cu2_3 + Cu2_4)/6;

end

function [dym,dya,dvm,dva] = DifEqu(ym,ya,vm,va,Rm,Rc,Rk,Rs,zm,tau)
    dym = vm; 
    dya = va;
    dvm = sin(Rs*tau)- 2*zm*vm-ym-Rc*2*zm*(vm-va)-Rk*(ym-ya);
    dva = -1/Rm*(Rc*2*zm*(va-vm)+Rk*(ya-ym));
end

関数2種類でルンゲクッタ法を実装している。
関数RungeKutta内で、4本の微分方程式について、4種類の重みを求め足し合わせている。
そして、ここで参照しているDifEquに解きたい微分方程式を記述すれば良い。

##結果

主系の振幅$y_m$を見てみよう。

縦軸に時間、横軸に周波数。
出力にはmatlabのimage関数を使用した。「黄色くなるほど値が大きい」ので、2つの周波数付近で振幅が大きくなっていることがわかる。

共振と非共振の場合の数値の違いを見る。
下の画像の①②③のラインに沿って$y_m$の変化を確認する。

下の画像が$y_m$の時間変化である。強制振動の変位の波形はこのような形になるのが可視化されて面白い。

非共振の場合は振幅が減少する。共振状態であれば、理想的には振幅は減少せず、増加し続ける。（減衰があればやがて収束する。）

前編では、ルンゲ＝クッタ法で微分方程式の数値解析が行えるところまでを確認した。
宣言した微分方程式の関数DifEquの中身を変えることによって、様々な形に応用が効くので試してみてほしい。

後編では
①変位$y_m$のピークを取得し、周波数特性を出力すること
②変位$y_m$の減衰から、無駄な計算を打ち切ること
③パラーメーター比（Rm,Rc,Rm）を変化させること
の3つを実装する。