r_damped_oscillator.m 4.05 KB
function r_damped_oscillator

x0 = -0.191865;
u0 = 0.0;

dt = 0.01;
t  = 0.0:dt:200;
N  = length(t);

x1(1) = x0;
u1(1) = u0;

x2(1) = x0;
u2(1) = u0;

x3(1) = x0;
u3(1) = u0;

x4(1) = x0;
u4(1) = u0;

x5(1) = x0;
u5(1) = u0;

D    = 0.99;                    % m
h    = 0.033;                   % m
A    = h.*D;                    % m^2
S    = pi.*(D./2).^2;           % m^2

T    = 12.5;                    % s

rhod = 550;                     % kg m^-3
rhow = 1000;                    % kg m^-3
m    = S.*h.*rhod;              % kg
d    = m./(rhow.*S);            % m
As   = d.*D;                    % Submerged cross section
cd   = 0.75;                     % -
ca   = 0.07;                    % -
ma   = ca.*rhow.*S.*d;          % kg
k    = (m + ma).*(2.*pi./T).^2; % N m^-1 (kg s^-2)
%k    = 30.*31.0;
L0   = 2.98;                    % m

disp(['D  = ',num2str(D)]);
disp(['h  = ',num2str(h)]);
disp(['d  = ',num2str(d)]);
disp(['k  = ',num2str(k)]);
disp(['m  = ',num2str(m)]);
disp(['ma = ',num2str(ma)]);
disp(['A  = ',num2str(A)]);
disp(['As = ',num2str(As)]);
disp(['S  = ',num2str(S)]);

%T = sqrt(k./m); disp(T);

%% Version 1 - First-order accurate (Euler method)
% Compute x(n+1) based on n(n) and u(n) and u(n) based x(n)
for n = 1:N-1;
    a_spring = -1.*k.*x1(n)./m;
    u1(n+1) = a_spring.*dt + u1(n);
    x1(n+1) = x1(n) + u1(n).*dt + a_spring.*dt.^2./2;
end

%% Version 2
% Compute x(n+1) based on x(n) and u(n+1) and u(n+1) based on x(n)
for n = 1:N-1;
    a_spring = -1.*k.*x2(n)./m;
    u2(n+1) = a_spring.*dt + u2(n);
    x2(n+1) = x2(n) + u2(n+1).*dt + a_spring.*dt.^2./2;
end

%% Version 3 - Second-order accurate
% Compute x(n+1) based on x(n) and u(n) and u(n+1) based on x(n)
for n = 1:N-1;
    astmp1  = -1.*k.*x3(n)./m;
    xtmp    = x3(n) + u3(n).*dt + 0.5.*astmp1.*dt.^2;
    astmp2  = -1.*k.*xtmp./m;
    as      = 0.5.*(astmp1 + astmp2);
    u3(n+1) = as.*dt + u3(n);
    x3(n+1) = x3(n) + u3(n).*dt + 0.5.*as.*dt.^2; 
end

%% Version 4 - Second-order accurate with quadratic damping
for n = 1:N-1;
    adtmp1  = -1.*0.5.*rhow.*cd.*A.*sign(u4(n)).*u4(n).^2./(m + ma);
    astmp1  = -1.*k.*x4(n)./(m + ma);
    xtmp    = x4(n) + u4(n).*dt + 0.5.*(astmp1 + adtmp1).*dt.^2;
    astmp2  = -1.*k.*xtmp./(m + ma);
    as      = 0.5.*(astmp1 + astmp2);
    
    u4(n+1) = (as + adtmp1).*dt + u4(n);
    adtmp2  = -1.*0.5.*rhow.*cd.*A.*sign(u4(n+1)).*u4(n+1).^2./(m + ma);
    ad      = 0.5.*(adtmp1 + adtmp2);
    x4(n+1) = x4(n) + u4(n).*dt + 0.5.*(as + ad).*dt.^2;
end

%% Version 5 - Second-order accurate with quadratic damping
for n = 1:N-1;
    adtmp1  = -1.*0.5.*rhow.*cd.*A.*sign(u5(n)).*u5(n).^2./(m + ma);
    theta1  = atan(x5(n)./L0);
    astmp1  = -1.*k.*x5(n)./(m + ma);
    %astmp1  = -1.*k.*(x5(n) - L0.*sin(theta1))./(m + ma);
    %astmp1  = -1.*7.5.*x5(n).^2./(m + ma);
    xtmp    = x5(n) + u5(n).*dt + 0.5.*(astmp1 + adtmp1).*dt.^2;
    theta2  = atan(xtmp./L0);
    astmp2  = -1.*k.*xtmp./(m + ma);
    %astmp2  = -1.*k.*(xtmp - L0.*sin(theta2))./(m + ma);
    %astmp2  = -1.*7.5.*xtmp.^2./(m + ma);
    as      = 0.5.*(astmp1 + astmp2);
    
    u5(n+1) = (as + adtmp1).*dt + u5(n);
    adtmp2  = -1.*0.5.*rhow.*cd.*A.*sign(u5(n+1)).*u5(n+1).^2./(m + ma);
    ad      = 0.5.*(adtmp1 + adtmp2);
    x5(n+1) = x5(n) + u5(n).*dt + 0.5.*(as + adtmp1).*dt.^2;
end

% figure(1); clf
% subplot(3,1,1)
% hold on
% plot(t,x3,'-k','LineWidth',1);
% plot(t,x0.*ones(size(t)),'--k');
% plot(t,-x0.*ones(size(t)),'--k');
% hold off
% ylim(1.2.*[-1.*abs(x0) abs(x0)])
% 
% subplot(3,1,2)
% hold on
% plot(t,x1,'-b','LineWidth',1);
% plot(t,x0.*ones(size(t)),'--k');
% plot(t,-x0.*ones(size(t)),'--k');
% hold off
% ylim(1.2.*[-1.*abs(x0) abs(x0)])
% 
% subplot(3,1,3)
% hold on
% plot(t,x2,'-r','LineWidth',1);
% plot(t,x0.*ones(size(t)),'--k');
% plot(t,-x0.*ones(size(t)),'--k');
% hold off
% ylim(1.2.*[-1.*abs(x0) abs(x0)])



%load decay_001.dat

figure(2); clf
hold on
plot(t,x4,'-r');
plot(t,x5,'-b');
%plot((decay_001(:,1)-342.16)./10,decay_001(:,2)./100,'-k');
hold off

% figure(3); clf
% hold on
% x = 0.001:0.1:2;
% theta = atan(x./L0);
% plot(x,(x + L0.*sin(theta)),'-k');
% plot(x,x,'--k');
% hold off