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

動吸振器の問題を解く(1)の続き

前回はルンゲ＝クッタ法で2自由度の2階の微分方程式（連成振動）の数値計算が行えるところまで確かめた。
今回は、周波数特性のグラフ化まで行う。
主系に対する従系の質量比、ダンパ減衰比、バネ剛性比を変えた状態で周波数特性を比較することで、目的に応じた動吸振器の設計が可能になるだろう。

#ループの打ち切り

　　　　　　 %% Store extreme values
            if i>2 && ym(i,j,iRk)<ym(i-1,j,iRk) && ym(i-1,j,iRk)>ym(i-2,j,iRk)
                yex(s,j,iRk)=ym(i-1,j,iRk);
                s=s+1;
            end

            %% Check the decrese of extreme values 
            if s>10 && yex(s-1,j,iRk)<(yex(s-2,j,iRk)+yex(s-3,j,iRk)+yex(s-4,j,iRk)+yex(s-5,j,iRk))/4.8
                fin(j,iRk)=s;
                break
            else 
            [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

一番下のサンプルコードの、この部分でループの打ち切り判定を行う
おおまかな意味は

①極大値格納
②if 極大値が（手前４つの極大値の和/4.8）より小さい
　→減衰しているのでループ打ち切り（打ち切り点を記録）
　else 計算続行

である。
この部分を時間ループ、周波数ループ、パラメータ比のループ（Rm,Rc,Rk）の5つの入れ子でループさせる。
打ち切りなし

打ち切りあり

これでかなり計算量が減るので計算時間が格段に減った。

#パラメータ比ごとの計算
パラメータ比の宣言がこちらである。

Rm = 1;
Rc = 1;
Rk = [0.1,1,10];

ベクトルになっているので、好きに値を変更・追加できる。

パラメータ比ごとに計算させるために、変位、速度の変数が

ym = zeros(Tn+1,Fn,RmN,RcN,RkN);
ya = zeros(Tn+1,Fn,RmN,RcN,RkN);
vm = zeros(Tn+1,Fn,RmN,RcN,RkN);
va = zeros(Tn+1,Fn,RmN,RcN,RkN);

このような形の5次元配列になっている。
各要素ごとにforループを増やせば良いというわけだ。

計算数はTn,Fnで宣言している。ここを増やすことで時間、および周波数の計算数が増加し、精度の高い結果が得られる。

#周波数特性
各周波数ごとに変位の最大値を格納していき、周波数に対応させてプロットさせる。
ここの計算で、無次元化した変位をもとの変位に戻している。（$y = x/(F/k_m)$より)

xm(j,iRm,iRc,iRk) = F/km*max(ym(:,j,iRm,iRc,iRk));
xa(j,iRm,iRc,iRk) = F/km*max(ya(:,j,iRm,iRc,iRk));

周波数は$f = R_s/T_m$ なので、それに対応させたプロットをすれば、無次元化ではなくもとの微分方程式の値が算出される。

#サンプルコード

clear
close
hold off

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

Tn = 5000;
Fn = 500;

Rm = [1];
Rc = [1.2];
Rk = [0.1,1,10];

mm = 10;
cm = 0.001;
km = 30;
F = 1; % Force

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


%% Loop preparation

% Declaration
ym = zeros(Tn+1,Fn,RmN,RcN,RkN);
ya = zeros(Tn+1,Fn,RmN,RcN,RkN);
vm = zeros(Tn+1,Fn,RmN,RcN,RkN);
va = zeros(Tn+1,Fn,RmN,RcN,RkN);


% Loop variable
RmN = length(Rm);
RcN = length(Rc);
RkN = length(Rk);
dtau = 0.1 *sm;
Rs = RsLow:(RsHigh-RsLow)/Fn:RsHigh;

% For plot
freq = Rs/Tm;
N_freq = sm/(2*pi);

% max variable
xm=zeros(Fn+1,RmN,RcN,RkN);
xa=zeros(Fn+1,RmN,RcN,RkN);

% Loop cutoff variable
yex= zeros(Tn,Fn,RmN,RcN,RkN);
maxI = 0;
fin=zeros(Fn,RmN,RcN,RkN);

%% Calculation
for iRm = 1:RmN
for iRc = 1:RcN
for iRk = 1:RkN
    for j=1:Fn % frequency loop
        s=1;
        for i=1:Tn % time loop
            %% Store extreme values
            if i>2 && ym(i,j,iRk)<ym(i-1,j,iRk) && ym(i-1,j,iRk)>ym(i-2,j,iRk)
                yex(s,j,iRk)=ym(i-1,j,iRk);
                s=s+1;
            end
            %% Check the decrese of extreme values 
            if s>10 && yex(s-1,j,iRk)<(yex(s-2,j,iRk)+yex(s-3,j,iRk)+yex(s-4,j,iRk)+yex(s-5,j,iRk))/4.8
                fin(j,iRk)=s;
                break
            else 
            [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
        xm(j,iRm,iRc,iRk) = F/km*max(ym(:,j,iRm,iRc,iRk));
        xa(j,iRm,iRc,iRk) = F/km*max(ya(:,j,iRm,iRc,iRk));
    end
end
end
end


%% Plot

%%%%% Separate type %%%%%
for iRk = 1:RkN
    subplot(RkN,1,iRk)
    plot(freq,xm(:,1,1,iRk))
    xline(N_freq,'--','DisplayName','Natural Frequency')
    title(['Rk=',num2str(Rk(iRk))])
end

%%%%% Same area %%%%%
% hold on
% for iRk = 1:RkN
%     plot(freq,xm(:,1,1,iRk),'DisplayName',strcat('Rk =',num2str(Rk(iRk))))
% end
% xline(N_freq,'--','DisplayName','Natural Frequency')
% legend()
% hold off

%% Functions
 
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通りのプロットを用意したので順番に見てみよう。
今回はRkだけを3通りの値で変化させて周波数特性を出力した。
Rk = [ 0.1 , 1 , 10 ]

Rk = 0.1およびRk = 10にはピークが1つしか無いように見える。
しかし、よく見ると極小のピークが立っている。
また、破線が主系のみの固有振動数である。つまり、従系（動吸振器）を追加することで共振周波数がシフトすることが確認できた。

どちらのグラフも縦軸が変位、横軸が周波数でプロットしている。
このときの値は無次元化したものではなく、もとの微分方程式の値である。
すなわち、実際に設計に反映できる値である。

#まとめ

今回の記事では動吸振器をモデルにした2自由度振動系の数値計算をmatlabで行った。
個人的にはmatlabはpythonよりも使いやすいと思う。
pythonと比べたとき、
・標準ライブラリが非常に豊富なので、追加インストールなしで使える。
・計算が速い（特に配列の計算）
・変数の値を表示してくれる。(↓画像)
・複数のタスクを開いておける（計算とplotを分割出来る）
この4点が非常に優秀だと感じる。

matlabの計算性能は非常に高いので、膨大なデータの計算をゴリゴリ回したい方は是非挑戦してみてほしい。

（完）