Professor Mohamed Hafez

Professor Mohamed Hafez

TA: Edward Tavernetti

COSMOS – Cluster 3

1. Semenov

%Semenov Law

clear; close all; clc;

dt = 1/20;

h = 1;

alpha = 20;

Beta = 4;

u_infty = 1;

u(1) = 0; %u0

k = 0;

while k < 300

k = k+1;

u(k+1) = -h*(u(k)-u_infty)*dt+alpha*exp(-Beta/u(k))*dt+u(k);

end

t = linspace(0,dt*(k+1),k+1);

plot(t,u)

xlabel('Time'),ylabel('Temperature'),

title(['T_0 = ',num2str(u(1)), ', T_\infty = 1']),

steady = alpha*exp(-Beta/u(k))-h*(u(k)-u_infty)

2. Frank-Kamenetskii

alpha = 20;

Beta = 4;

n = 21;

j = 1;

u(1,j) = 1; %Left boundary

for i = 2:n-1

u(i,j) = 1; %initial condition.

end

u(n,j) = 1; %Right boundary

dx = 1/(n-1);

K = 86*dx^2/2;

dt = 1/100 ;

while j < 1600

u(1,j+1) = 1; %Left boundary

for i = 2:n-1

u(i,j+1) = dt*K*(u(i+1,j)-2*u(i,j)+u(i-1,j))/dx^2+alpha*exp(-Beta/u(i,j))*dt+u(i,j);

end

u(n,j+1) = 1; %Right boundary

j = j+1;

end

t = linspace(0,dt*(j),j);

subplot(1,2,1), plot(t,u(10,:)),

xlabel('Time'),ylabel('Temperature'),

x = linspace(0,1,n);

subplot(1,2,2),surf(t(1:round(j/50):j),x,u(:,1:round(j/50):j,:))

xlabel('time'),ylabel('space')

alpha = 10;

Beta = 4;

n = 21;

j = 1;

u(1,j) = 1; %Left boundary

for i = 2:n-1

u(i,j) = 1; %initial condition.

end

u(n,j) = 1; %Right boundary

dx = 1/(n-1);

K = 86*dx^2/2;

dt = 1/100 ;

while j < 1000

u(1,j+1) = dt*K*(2*u(i+1,j)-2*u(i,j))/dx^2 ...

+dt*alpha*exp(-Beta/u(i,j))+u(i,j);

for i = 2:n-1

u(i,j+1) = dt*K*(u(i+1,j)-2*u(i,j)+u(i-1,j))/dx^2 ...

+dt*alpha*exp(-Beta/u(i,j))+u(i,j);

end

u(n,j+1) = dt*K*(-2*u(i,j)+2*u(i-1,j))/dx^2 ...

+dt*alpha*exp(-Beta/u(i,j))+u(i,j);

j = j+1;

end

t = linspace(0,dt*(j),j);

subplot(1,2,1),plot(t,u(10,:)),

xlabel('Time'),ylabel('Temperature'),

x = linspace(0,1,n);

subplot(1,2,2),surf(t(1:round(j/50):j),x,u(:,1:round(j/50):j,:))

xlabel('time'),ylabel('space')

u_infty = 1;

h = 1;

alpha = 10;

Beta = 4;

n = 21;

j = 1;

u(1,j) = 1; %Left boundary

for i = 2:n-1

u(i,j) = 1; %initial condition.

end

u(n,j) = 1; %Right boundary

dx = 1/(n-1);

K = 86*dx^2/2;

dt = 1/100 ;

while j < 2000

UD = u(i+1,j)+2*dx*h/K*(u_infty-u(i,j));

u(1,j+1) = dt*K*(u(i+1,j)-2*u(i,j)+UD)/dx^2 ...

+dt*alpha*exp(-Beta/u(i,j))+u(i,j);

for i = 2:n-1

u(i,j+1) = dt*K*(u(i+1,j)-2*u(i,j)+u(i-1,j))/dx^2 ...

+dt*alpha*exp(-Beta/u(i,j))+u(i,j);

end

UD = u(i-1,j)+2*dx*h/K*(u_infty-u(i,j));

u(n,j+1) = dt*K*(UD-2*u(i,j)+u(i-1,j))/dx^2 ...

+dt*alpha*exp(-Beta/u(i,j))+u(i,j);

j = j+1;

end

t = linspace(0,dt*(j),j);

subplot(1,2,1),plot(t,u(10,:)),

xlabel('Time'),ylabel('Temperature'),

x = linspace(0,1,n);

subplot(1,2,2),surf(t(1:round(j/50):j),x,u(:,1:round(j/50):j,:))

xlabel('time'),ylabel('space')

%Chemical Reaction

%Semenov Law 2

clear; close all; clc;

dt = 1/100;

h = 0;

alpha0 = 20;

Beta = 5;

u_infty = 1;

u(1) = 1; %u0

k = 0;

while k < 6000

k = k+1;

alpha = alpha0/(cosh((dt*k/20)^4))^2;

A(k) = alpha;

u(k+1) = -h*(u(k)-u_infty)*dt+alpha*exp(-Beta/u(k))*dt+u(k);

end

t = linspace(0,dt*(k+1),k+1);

plot(1:length(A),A)

plot(t,u)

xlabel('Time'),ylabel('Temperature')