Untitled

avatar
unknown
plain_text
2 years ago
9.3 kB
30
Indexable
%% Initiate the script
clc; clear; close all

%% Define the parameters
L = 0.2;       % m
m = 4;       % kg
k = 1600;       % N/m
c = 6;       % Ns/m  
mud = 1.8;
mus = 2.1;
r = 0;
vd= 0.0001;

%% Q1 Plot the friction model
Fd=mud*m*10;
xdot=linspace(-10e-3,10e-3,100);
figure(1); clf
hold on
plot(xdot,ffriction(xdot),'linewidth',2) % we apply the friction funtion 
rectx=[-vd,vd,vd,-vd];
recty=[-84,-84,84,84];
patch(rectx,recty,'red');
    xlabel('xdot in Seconds')
    ylabel('Ffriction in N')
    set(gca,'fontsize',11)
    set(gcf,'color','w','units','centimeters','outerposition',[5 5 15 10])
hold off

exportgraphics(gcf,'Figure_1.png','Resolution',300);    % Save a png image of the first graph with a resolution of 300 dpi

%% Q3
%Initial conditions
IC = [0; 0];
tt = linspace(0,10,2e3);  % define the time vector
opts = odeset('RelTol',1e-4,'AbsTol',1e-6);

[~,z_0] = ode45(@(t,z)  EX_1_ode_1(t,z,m,k,c,mud,mus), tt, IC, opts);   %we use ode45 to solve the differential equation for each case

%% Q4 Simulation
% Harmonic force parameters
A = 200;               % N
     
% Initial conditions
w = [5,20,35];      % rad/s

% Simulate the three cases for each frequency
[~,z_1] = ode45(@(t,z)  EX_1_ode_2(t,z,m,k,c,mud,mus,200,w(1)), tt, IC, opts); 
[~,z_2] = ode45(@(t,z)  EX_1_ode_2(t,z,m,k,c,mud,mus,200,w(2)), tt, IC, opts);
[~,z_3] = ode45(@(t,z)  EX_1_ode_2(t,z,m,k,c,mud,mus,200,w(3)), tt, IC, opts);

%% Q5
w = 5:0.1:35;
z_min = zeros(size(w));

for i = 1:1:length(w)
    [~,z_w]= ode45(@(t,z) EX_1_ode_2(t,z,m,k,c,mud,mus,A,w(i)), tt, IC, opts);
    z_min(i)=min(z_w(:,1));
%     clear z_w
end

idx = find(z_min<=-L);  %if the mass falls from the table
w_min = w(idx(1));
idx_max=length(idx);
w_max = w(idx(idx_max));

[~,z_w_min]= ode45(@(t,z) EX_1_ode_2(t,z,m,k,c,mud,mus,A,w_min), tt, IC, opts);
[~,z_w_max]= ode45(@(t,z) EX_1_ode_2(t,z,m,k,c,mud,mus,A,w_max), tt, IC, opts);


%% Q6
k_a = 700;
k_b = -700;
w = 5:0.1:35;
z_min_2 = zeros(size(w));
z_min_3 = zeros(size(w));

for i = 1:1:length(w)
    [~,z_w_2]= ode45(@(t,z) EX_1_ode_3(t,z,m,k,k_a,c,mud,mus,200,w(i)), tt, IC, opts);
    [~,z_w_3]= ode45(@(t,z) EX_1_ode_3(t,z,m,k,k_b,c,mud,mus,200,w(i)), tt, IC, opts);
    z_min_2(i)=min(z_w_2(:,1));
    z_min_3(i)=min(z_w_3(:,1));
%     clear z_w_2
%     clear z_w_3
end

idx_a = find(z_min_2<=-L);
w_min_2 = w(idx_a(1));
idx_a_max=length(idx_a);
w_max_2 = w(idx_a(idx_a_max));

idx_b = find(z_min_3<=-L);  
w_min_3 = w(idx_b(1));
idx_b_max=length(idx_b);
w_max_3 = w(idx_b(idx_b_max));

[~,z_w_2_min]= ode45(@(t,z) EX_1_ode_3(t,z,m,k,k_a,c,mud,mus,200,w_min_2), tt, IC, opts);
[~,z_w_2_max]= ode45(@(t,z) EX_1_ode_3(t,z,m,k,k_a,c,mud,mus,200,w_max_2), tt, IC, opts);

[~,z_w_3_min]= ode45(@(t,z) EX_1_ode_3(t,z,m,k,k_b,c,mud,mus,200,w_min_3), tt, IC, opts);
[~,z_w_3_max]= ode45(@(t,z) EX_1_ode_3(t,z,m,k,k_b,c,mud,mus,200,w_max_3), tt, IC, opts);

%% Data visulaization 
hold on

figure(2);      %case where F=-250 N ans we only have 1 spring
yyaxis left
    plot(tt,z_0(:,1),'linewidth',2);
    xlabel('$t$(s)','interpreter','Latex')
    ylabel('$x$(m)','interpreter','Latex')
    set(gca,'FontSize',11)
yyaxis right
    plot(tt,z_0(:,2),'linewidth',2);
    ylabel('$\dot{x}$(m/s)','interpreter','Latex')
    set(gca,'FontSize',11)
set(gcf,'color','w','units','centimeters','outerposition',[5 5 15 10])

figure(3);  % plot of the velocity with added inital conditions for the frequency
yyaxis left
    plot(tt,z_1(:,1),'linewidth',2);
    xlabel('$t$(s)','interpreter','Latex')
    ylabel('$x$(m)','interpreter','Latex')
    set(gca,'FontSize',11)
yyaxis right
    plot(tt,z_1(:,2),'linewidth',2);
    ylabel('$\dot{x}$(m/s)','interpreter','Latex')
    set(gca,'FontSize',11)
set(gcf,'color','w','units','centimeters','outerposition',[5 5 15 10])
hold off

figure(4);
yyaxis left
    plot(tt,z_2(:,1),'linewidth',2);
    xlabel('$t$(s)','interpreter','Latex')
    ylabel('$x$(m)','interpreter','Latex')
    set(gca,'FontSize',11)
yyaxis right
    plot(tt,z_2(:,2),'linewidth',2);
    ylabel('$\dot{x}$(m/s)','interpreter','Latex')
    set(gca,'FontSize',11)
set(gcf,'color','w','units','centimeters','outerposition',[5 5 15 10])

figure(5);
yyaxis left
    plot(tt,z_3(:,1),'linewidth',2);
    xlabel('$t$(s)','interpreter','Latex')
    ylabel('$x$(m)','interpreter','Latex')
    set(gca,'FontSize',11)
yyaxis right
    plot(tt,z_3(:,2),'linewidth',2);
    ylabel('$\dot{x}$(m/s)','interpreter','Latex')
    set(gca,'FontSize',11)
set(gcf,'color','w','units','centimeters','outerposition',[5 5 15 10])

figure(6);      % plot for w min and max to make sure the mass won't fall from the table
yyaxis left
    plot(tt,z_w_min(:,1),'linewidth',2);
    xlabel('$t$(s)','interpreter','Latex')
    ylabel('$x$(m)','interpreter','Latex')
    set(gca,'FontSize',11)
yyaxis right
    plot(tt,z_w_min(:,2),'linewidth',2);
    ylabel('$\dot{x}$(m/s)','interpreter','Latex')
    set(gca,'FontSize',11)
set(gcf,'color','w','units','centimeters','outerposition',[5 5 15 10])

figure(7);
yyaxis left
    plot(tt,z_w_max(:,1),'linewidth',2);
    xlabel('$t$(s)','interpreter','Latex')
    ylabel('$x$(m)','interpreter','Latex')
    set(gca,'FontSize',11)
yyaxis right
    plot(tt,z_w_max(:,2),'linewidth',2);
    ylabel('$\dot{x}$(m/s)','interpreter','Latex')
    set(gca,'FontSize',11)
set(gcf,'color','w','units','centimeters','outerposition',[5 5 15 10])

figure(8);      % plot of the velocity with the second spring for w min and max 
yyaxis left
    plot(tt,z_w_2_min(:,1),'linewidth',2);
    xlabel('$t$(s)','interpreter','Latex')
    ylabel('$x$(m)','interpreter','Latex')
    set(gca,'FontSize',11)
yyaxis right
    plot(tt,z_w_2_min(:,2),'linewidth',2);
    ylabel('$\dot{x}$(m/s)','interpreter','Latex')
    set(gca,'FontSize',11)
set(gcf,'color','w','units','centimeters','outerposition',[5 5 15 10])

figure(9);
yyaxis left
    plot(tt,z_w_2_max(:,1),'linewidth',2);
    xlabel('$t$(s)','interpreter','Latex')
    ylabel('$x$(m)','interpreter','Latex')
    set(gca,'FontSize',11)
yyaxis right
    plot(tt,z_w_2_max(:,2),'linewidth',2);
    ylabel('$\dot{x}$(m/s)','interpreter','Latex')
    set(gca,'FontSize',11)
set(gcf,'color','w','units','centimeters','outerposition',[5 5 15 10])

figure(10);
yyaxis left
    plot(tt,z_w_3_min(:,1),'linewidth',2);
    xlabel('$t$(s)','interpreter','Latex')
    ylabel('$x$(m)','interpreter','Latex')
    set(gca,'FontSize',11)
yyaxis right
    plot(tt,z_w_3_min(:,2),'linewidth',2);
    ylabel('$\dot{x}$(m/s)','interpreter','Latex')
    set(gca,'FontSize',11)
set(gcf,'color','w','units','centimeters','outerposition',[5 5 15 10])

figure(11);
yyaxis left
    plot(tt,z_w_3_max(:,1),'linewidth',2);
    xlabel('$t$(s)','interpreter','Latex')
    ylabel('$x$(m)','interpreter','Latex')
    set(gca,'FontSize',11)
yyaxis right
    plot(tt,z_w_3_max(:,2),'linewidth',2);
    ylabel('$\dot{x}$(m/s)','interpreter','Latex')
    set(gca,'FontSize',11)
set(gcf,'color','w','units','centimeters','outerposition',[5 5 15 10])

figure(12);
plot(w,z_min,'linewidth',2);
xline(w_min,'--b','linewidth',1)
xline(w_max,'--b','linewidth',1)
yline(-L,'--r','linewidth',1);
xlabel('$w$(rad/s)','interpreter','Latex')
ylabel('$x$(m)','interpreter','Latex')
title('minimum postition $x$ in fonction of frequencies $w$','interpreter','Latex')
exportgraphics(gcf,'Figure_6.png','Resolution',300);

figure(13);     % next 3 plots are the plots where we determine w min and max for each case 
plot(w,z_min_2,'linewidth',2);
xline(w_min_2,'--b','linewidth',1)
xline(w_max_2,'--b','linewidth',1)
yline(-L,'--r','linewidth',1);
xlabel('$w$(rad/s)','interpreter','Latex')
ylabel('$x$(m)','interpreter','Latex')
title('minimum postition $x$ in fonction of frequencies $w$','interpreter','Latex')
exportgraphics(gcf,'Figure_6.png','Resolution',300);

figure(14);
plot(w,z_min_3,'linewidth',2);
xline(w_min_3,'--b','linewidth',1)
xline(w_max_3,'--b','linewidth',1)
yline(-L,'--r','linewidth',1);
xlabel('$w$(rad/s)','interpreter','Latex')
ylabel('$x$(m)','interpreter','Latex')
title('minimum postition $x$ in fonction of frequencies $w$','interpreter','Latex')
exportgraphics(gcf,'Figure_6.png','Resolution',300);

function y=ffriction(v) % Karnopp friction model without Fext for Q1
Fd=72;
vd=0.0001;                                                       
if(abs(v)>vd)
    y=-Fd*sign(v);
end
end

function dzdt2 = EX_1_ode_1(~,z,m,k,c,mud,mus)% order reduction Q3
vd=10e-4;
g=10;
dzdt2(1,1) = z(2);
if (abs(z(2))<=vd)
    dzdt2(2,1) = (-min(max(-mus*m*g,-250-k*z(1)-c*z(2)),mus*m*g)-k*z(1)-c*z(2)-250)/m;
elseif (abs(z(2))>vd)
    dzdt2(2,1)=(-mud*m*g*sign(z(2))-k*z(1)-c*z(2)-250)/m;
end
end

function dzdt = EX_1_ode_2(t,z,m,k,c,mud,mus,A,w) % order reduction Q4 and Q5
vd=10e-4;
g=10;
dzdt(1,1) = z(2);
if (abs(z(2))<=vd)
    dzdt(2,1) = (-min(max(-mus*m*g,A*sin(w*t)-k*z(1)-c*z(2)),mus*m*g)-k*z(1)-c*z(2)+A*sin(w*t))/m;
elseif (abs(z(2))>vd)
    dzdt(2,1)=(-mud*m*g*sign(z(2))-k*z(1)-c*z(2)+A*sin(w*t))/m;
end
end

function dzdt3 = EX_1_ode_3(t,z,m,k,k2,c,mud,mus,A,w)% order reduction Q6. I have put K2 in this function.
vd=10e-4;
g=10;
dzdt3(1,1) = z(2);
if (abs(z(2))<=vd)
    dzdt3(2,1) = (-min(max(-mus*m*g,A*sin(w*t)-k*z(1)-c*z(2)-k2*(z(1))^2),mus*m*g)-k*z(1)-c*z(2)-k2*(z(1))^2+A*sin(w*t))/m;
elseif (abs(z(2))>vd)
    dzdt3(2,1)=(-mud*m*g*sign(z(2))-k*z(1)-c*z(2)-k2*(z(1))^2+A*sin(w*t))/m;
end
end
Editor is loading...
Leave a Comment