Untitled

avatar
unknown
plain_text
a year ago
8.0 kB
32
Indexable
%% Initiate the script
clc; clear; close all

rho=1;
cp=1;
he=0.025;
Ae=0.25*pi;
V=0.05*Ae;
Te=298;
ht=10;
At=16*pi*10^-6;
cq=2*10^-8;
ct=0.06;
j=10^-3;
M=[0.25 0.5 1];
cT=2.2*10^-8;
c=10^-2;

tt = linspace(0,150,2e4);

IC1 = [298; 298 ;0;0];

opts = odeset('RelTol',1e-4,'AbsTol',1e-6);
[~,z_1] = ode45(@(t,z)  EX_3_ode_Students_2022(t,z,V,rho,cp,he,Ae,Te,ht,At,cq,ct,j,M(1),cT,c,1000), tt, IC1, opts); 
[~,z_2] = ode45(@(t,z)  EX_3_ode_Students_2022(t,z,V,rho,cp,he,Ae,Te,ht,At,cq,ct,j,M(2),cT,c,1000), tt, IC1, opts); 
[~,z_3] = ode45(@(t,z)  EX_3_ode_Students_2022(t,z,V,rho,cp,he,Ae,Te,ht,At,cq,ct,j,M(3),cT,c,1000), tt, IC1, opts); 

%% Q3 a)
figure(1)
hold on
plot(tt,z_1(:,1),'linewidth',2);
plot(tt,z_2(:,1),'linewidth',2);
plot(tt,z_3(:,1),'linewidth',2);
xline(16,'--b','linewidth',1)
set(gca,'FontSize',11)
xlabel('$time$(s)','interpreter','Latex');
ylabel('$T$(K)','interpreter','Latex');
title('Temperature as a function of time','interpreter','Latex')
legend("M=0.25 N.m","M=0.5 N.m","M=1 N.m","steady state",'location','best');
legend show

exportgraphics(gcf,'Figure_d.png','Resolution',300);

hold off

%% Q3 b)
T_norm1=(z_1(:,1)-min(z_1(:,1)))./(max(z_1(:,1))-min(z_1(:,1)));
T_norm2=(z_2(:,1)-min(z_2(:,1)))./(max(z_2(:,1))-min(z_2(:,1)));
T_norm3=(z_3(:,1)-min(z_3(:,1)))./(max(z_3(:,1))-min(z_3(:,1)));

figure (2)
hold on

plot (tt, T_norm1,'linewidth',2);
plot (tt, T_norm2,'linewidth',2);
plot (tt, T_norm3,'linewidth',2);
xline(16,'--b','linewidth',1)
xlabel('$time$(s)','interpreter','Latex');
ylabel('$T$(normalized)','interpreter','Latex');
legend("M=0.25 N.m","M=0.5 N.m","M=1 N.m","steady state",'location','best');
legend show
title('Normalized temperature as a function of time','interpreter','Latex')

exportgraphics(gcf,'Figure_d2.png','Resolution',300);

hold off

%% Q3 c)
figure(3)
hold on
plot(tt,z_1(:,4),'linewidth',2);
plot(tt,z_2(:,4),'linewidth',2);
plot(tt,z_3(:,4),'linewidth',2);
set(gca,'FontSize',11)
xlabel('$time$(s)','interpreter','Latex');
ylabel('$Angular Velocity$(rad/sec)','interpreter','Latex');
legend("M=0.25 N.m","M=0.5 N.m","M=1 N.m",'location','best');
legend show
title('Angular Velocity as a function of time','interpreter','Latex')

exportgraphics(gcf,'Figure_3.png','Resolution',300);

hold off


V_norm1=(z_1(:,4)-min(z_1(:,4)))./(max(z_1(:,4))-min(z_1(:,4)));
V_norm2=(z_2(:,4)-min(z_2(:,4)))./(max(z_2(:,4))-min(z_2(:,4)));
V_norm3=(z_3(:,4)-min(z_3(:,4)))./(max(z_3(:,4))-min(z_3(:,4)));

figure (4)
hold on
plot (tt, V_norm1,'linewidth',2);
plot (tt, V_norm2,'linewidth',2);
plot (tt, V_norm3,'linewidth',2);
xlabel('$time$(s)','interpreter','Latex');
ylabel('$Angular Velocity $(normalized)','interpreter','Latex');
legend("M=0.25 N.m","M=0.5 N.m","M=1 N.m",'location','best');
legend show
title('Normalized angular velocity as a function of time','interpreter','Latex')
exportgraphics(gcf,'Figure_4.png','Resolution',300);
hold off

%% Q3 e)
tt2=linspace(0,100,2e3);
[~,z_4] = ode45(@(t,z)  EX_3_ode_Students_2022(t,z,V,rho,cp,he,Ae,Te,ht,At,cq,ct,j,M(1),cT,c,16), tt2, IC1, opts); 
[~,z_5] = ode45(@(t,z)  EX_3_ode_Students_2022(t,z,V,rho,cp,he,Ae,Te,ht,At,cq,ct,j,M(2),cT,c,16), tt2, IC1, opts); 
[~,z_6] = ode45(@(t,z)  EX_3_ode_Students_2022(t,z,V,rho,cp,he,Ae,Te,ht,At,cq,ct,j,M(3),cT,c,16), tt2, IC1, opts);


figure(5)
n=3;
hold on

plot(tt2,z_4(:,1),'linewidth',2);
plot(tt2,z_5(:,1),'linewidth',2);
plot(tt2,z_6(:,1),'linewidth',2);
xline(16,'--b','linewidth',1)
xline(32,'--r','linewidth',1)
set(gca,'FontSize',11)
xlabel('$time$(s)','interpreter','Latex');
ylabel('$T$(K)','interpreter','Latex');
legend("M=0.25 N.m","M=0.5 N.m","M=1 N.m","first steady state","second steady state",'location','best');

legend show
title('Temperature as a function of time','interpreter','Latex')
exportgraphics(gcf,'Figure_5.png','Resolution',300);
hold off

figure(6)
hold on
plot(tt2,z_4(:,4),'linewidth',2);
plot(tt2,z_5(:,4),'linewidth',2);
plot(tt2,z_6(:,4),'linewidth',2);
xline(16,'--b','linewidth',1)
xline(32,'--r','linewidth',1)
set(gca,'FontSize',11)
xlabel('$time$(s)','interpreter','Latex');
ylabel('$Angular Velocity$(rad/sec)','interpreter','Latex');
legend("M=0.25 N.m","M=0.5 N.m","M=1 N.m",'location','best');
legend show
title('Angular velocity as a function of time','interpreter','Latex')

exportgraphics(gcf,'Figure_6.png','Resolution',300);

hold off

%% Q4

 RE_1=100.*(z_4(:,2)-z_4(:,1))./z_4(:,1);
 RE_2=100.*(z_5(:,2)-z_5(:,1))./z_5(:,1);
 RE_3=100.*(z_6(:,2)-z_6(:,1))./z_6(:,1);

 figure (7)
 hold on
 plot(tt2, RE_1,'linewidth',2);
 plot(tt2, RE_2,'linewidth',2);
 plot(tt2, RE_3,'linewidth',2);
 xline(16,'--b','linewidth',1)
 xline(32,'--r','linewidth',1)
 xlabel('$time$(s)','interpreter','Latex');
 ylabel('Relative Error %');
 legend("M=0.25 N.m","M=0.5 N.m","M=1 N.m","first steady state","second steady state",'location','best');
 
 legend show
 title('Error between the temperature in the tank and the one measured by the thermometer ','interpreter','Latex')
 exportgraphics(gcf,'Figure_7.png','Resolution',300);
 hold off

%% Q5
tt3=linspace(0,1000,4e3);
M=linspace(0,2,1e2);
for i=1:length(M)
[~,z_7] = ode45(@(t,z)  EX_3_ode_Students_2022(t,z,V,rho,cp,he,Ae,Te,ht,At,cq,ct,j,M(i),cT,c,10000), tt3, IC1, opts);
SteadyT(i)= z_7(end,1);
SteadyV(i)=z_7(end,4);
end
Mss=interp1(SteadyV, M,100)
Tss= interp1(M,SteadyT,Mss)
figure(8)
yyaxis left
hold on
    yline(Tss,'--b','linewidth',1)
    plot(Mss,Tss,'ro','LineWidth',1.5);
    plot(M,SteadyT,'-','linewidth',2);
    xlabel('$M$(N.m)','interpreter','Latex')
    ylabel('$T$(K)','interpreter','Latex')
    set(gca,'FontSize',11)
yyaxis right
    plot(Mss,100,'go','LineWidth',1.5);
    plot(M,SteadyV,'-','linewidth',2);
    yline(100,'--r','linewidth',1);
    ylabel('$V$(rad/sec)','interpreter','Latex')
    set(gca,'FontSize',11)
set(gcf,'color','w','units','centimeters','outerposition',[5 5 15 10])
title('Final temperature and angular velocity as a function of the input Moment','interpreter','Latex')
exportgraphics(gcf,'Figure_8.png','Resolution',300);

%% Q6

tt = linspace(0,150,2e4);
IC1 = [298; 298 ;0;0];
M = [-0.5,0.5];

[~,z_1] = ode45(@(t,z)  EX_3_2_ode_Students_2022(t,z,V,rho,cp,he,Ae,Te,ht,At,cq,ct,j,M(1),cT,c), tt, IC1, opts); 
[~,z_2] = ode45(@(t,z)  EX_3_2_ode_Students_2022(t,z,V,rho,cp,he,Ae,Te,ht,At,cq,ct,j,M(2),cT,c), tt, IC1, opts);


figure(9)
hold on
plot(tt,z_1(:,1),'linewidth',2);
plot(tt,z_2(:,1),'linewidth',2);
set(gca,'FontSize',11)
xlabel('$time$(s)','interpreter','Latex');
ylabel('$T$(K)','interpreter','Latex');
title('Temperature as a function of time','interpreter','Latex')
legend("M=-0.5 N.m","M=0.5 N.m",'location','best');
legend show
exportgraphics(gcf,'Figure_9.png','Resolution',300);
hold off

figure(10)
hold on
plot(tt,z_1(:,4),'linewidth',2);
plot(tt,z_2(:,4),'linewidth',2);
set(gca,'FontSize',11)
xlabel('$time$(s)','interpreter','Latex');
ylabel('$Angular Velocity$(rad/sec)','interpreter','Latex');
legend("M=-0.5 N.m","M=0.5 N.m",'location','best');
legend show
title('Angular Velocity as a function of time','interpreter','Latex')
exportgraphics(gcf,'Figure_10.png','Resolution',300);
hold off



%% Function
function dzdt = EX_3_ode_Students_2022(t,z,V,rho,cp,he,Ae,Te,ht,At,cq,ct,j,M,cT,c,ts)
    if (t<ts)
        dzdt(1,1) = (he*Ae*(Te-z(1))+ht*At*(z(2)-z(1))+cq*z(4)^4)/(V*rho*cp); %dT/dt
        dzdt(2,1) = (ht*At*(z(1)-z(2)))/ct; %dTt/dt
        dzdt(3,1) = z(4); %omega
        dzdt(4,1)= (M-cT*z(4)^3-c*z(4))/j; %omega dot
    else
        dzdt(1,1) = (he*Ae*(Te-z(1))+ht*At*(z(2)-z(1))+cq*z(4)^4)/(V*rho*cp); %dT/dt
        dzdt(2,1) = (ht*At*(z(1)-z(2)))/ct; %dTt/dt
        dzdt(3,1) = z(4); %omega
        dzdt(4,1)= (-cT*z(4)^3-c*z(4))/j; %omega dot
    end
end

function dzdt = EX_3_2_ode_Students_2022(~,z,V,rho,cp,he,Ae,Te,ht,At,cq,ct,j,M,cT,c)
    dzdt(1,1) = (he*Ae*(Te-z(1))+ht*At*(z(2)-z(1))+cq*z(4)^5)/(V*rho*cp); %dT/dt
    dzdt(2,1) = (ht*At*(z(1)-z(2)))/ct; %dTt/dt
    dzdt(3,1) = z(4); %omega
    dzdt(4,1)= (M-cT*z(4)^2-c*z(4))/j; %omega dot  
end
Editor is loading...
Leave a Comment