동 흡진기 문제 해결 (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개의 극대치의 합/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

결과

앞서 언급했듯이 두 가지 플롯을 준비했으므로 차례로 보자.
이번에는 Rk만을 3가지 값으로 변화시켜 주파수 특성을 출력했다.
Rk = [0.1, 1, 10]

Rk = 0.1 및 Rk = 10에는 피크가 하나밖에 없는 것처럼 보인다.
그러나 잘 보면 극소의 피크가 서 있다.
또한, 파선은 주계만의 고유 진동수이다. 즉, 종계(동 흡진기)를 추가함으로써 공진 주파수가 시프트되는 것을 확인할 수 있었다.




두 그래프 모두 세로축이 변위, 가로축이 주파수로 플롯되어 있다.
이 때의 값은 무차원화된 것이 아니고, 원래의 미분 방정식의 값이다.
즉, 실제로 설계에 반영할 수 있는 값이다.

요약

이번 기사에서는 동흡진기를 모델로 한 2 자유도 진동계의 수치 계산을 matlab에서 실시했다.
개인적으로 matlab은 파이썬보다 사용하기 쉽다고 생각합니다.
파이썬과 비교했을 때,
・표준 라이브러리가 매우 풍부하므로, 추가 인스톨 없이 사용할 수 있다.
· 계산이 빠르다 (특히 배열 계산)
· 변수의 값을 표시합니다. (↓ 이미지)
・복수의 태스크를 열어 둔다(계산과 plot를 분할할 수 있다)
이 4점이 매우 우수하다고 느낀다.



matlab의 계산 성능은 매우 높기 때문에 방대한 데이터의 계산을 골고루 돌리고 싶은 분은 꼭 도전해 보길 바란다.

(완료)