clear all
clc
format long
global T0 Rg E1 E2 E_2 Eg k10 k20 k2_0 kgo Deo k1 k2 k_2 Ke Kg Cb Ct Ce a
a=1;
T0=800;
Rg=8.314;
E1=10^5;
E2=10^5;
E_2=10^2;
Eg=10^4;
k10=0.1;
k20=0.1;
k2_0=0.1;
kgo=0.1;
Deo=10^-10;
k1 = k10*exp(-E1/(Rg*T0));
k2 = k20*exp(-E2/(Rg*T0));
k_2 = k2_0*exp(-E_2/(Rg*T0));
Ke=(2.01*10^(-6))*T0^2-(1.43*10^(-3))*T0+0.2544;
Kg=kgo*exp(Eg/(Rg*T0));
Cb =0;
Ct=20;
Ce=(7.67131719*10^(-47))*T0^(14.5144484);
r=linspace(0,100e-4,5);
t=linspace(0,300,100);
n = numel(r);
CC0 = zeros(n,1);
CO0 = ones(n,1);
CD0 = ones(n,1);
CM0 = zeros(n,1);
CA0 = zeros(n,1);
CO0(1)=15;
CD0(1)=5;
y0 = [CC0;CO0;CD0;CM0;CA0];
[T, Y] = ode15s(@(t,y)fun(t,y,r,n),t,y0)
plot(t,Y)
set(gca,'YLim',[0 1e-2],'xLim',[0 200]);
xlabel('time')
ylabel('concentration')
function DyDt=fun(t,y,r,n)
global T0 Rg E1 E2 E_2 Eg k10 k20 k2_0 kgo Deo k1 k2 k_2 Ke Kg Cb Ct Ce a
CC = zeros(n,1);
CO = zeros(n,1);
CD = zeros(n,1);
CM = zeros(n,1);
CA = zeros(n,1);
DCCDt = zeros(n,1);
DCODt = zeros(n,1);
DCDDt = zeros(n,1);
DCMDt = zeros(n,1);
DCADt = zeros(n,1);
DyDt = zeros(5*n,1);
rhalf = zeros(n-1,1);
DCCDr = zeros(n-1,1);
D2CCDr2 = zeros(n-1,1);
CC = y(1:n);
CO = y(n+1:2*n);
CD = y(2*n+1:3*n);
CM = y(3*n+1:4*n);
CA = y(4*n+1:5*n);
rhalf(1:n-1)=(r(1:n-1)+r(2:n))/2;
DCCDr(1)=0;
D2CCDr2(1) = 1.0/(rhalf(1)-r(1))*(CC(2)-CC(1))/(r(2)-r(1));
DCCDt(n)=Cb
DCODt(n) =-k1*(CO(n)^2 -( 1/Ke)*CD(n)^2);
DCDDt(n) =k1*(CO(n)^2 - (1/Ke)*CD(n)^2)-k2*CD(n)+k_2*CA(n)*CM(n)*CC(n);
DCMDt(n) =k2*CD(n)-k_2*CA(n)*CM(n)*CC(n);
DCADt(n) =k2*CD(n)-k_2*CA(n)*CM(n)*CC(n);
for i=2:n-1
DCCDr(i) = ((r(i)-r(i-1))/(r(i+1)-r(i))*(CC(i+1)-CC(i))+(r(i+1)-r(i))/(r(i)-r(i-1))*(r(i)-CC(i-1)))/(r(i+1)-r(i-1));
D2CCDr2(i) = (rhalf(i)*(CC(i+1)-CC(i))/(r(i+1)-r(i))-rhalf(i-1)*(CC(i)-CC(i-1))/(r(i)-r(i-1)))/(rhalf(i)-rhalf(i-1));
end
for i=1:n-1
De=Deo*exp(a*CM(i)/Ct);
DCCDt(i) =2*De*DCCDr(i)+De*r(i)*D2CCDr2(i)+k2*CD(i)-k_2*CA(i)*CM(i)*CC(i);
DCODt(i) =-k1*(CO(i)^2 -( 1/Ke)*CD(i)^2);
DCDDt(i) =k1*(CO(i)^2 - (1/Ke)*CD(i)^2)-k2*CD(i)+k_2*CA(i)*CM(i)*CC(i);
DCMDt(i) =k2*CD(i)-k_2*CA(i)*CM(i)*CC(i);
DCADt(i) =k2*CD(i)-k_2*CA(i)*CM(i)*CC(i);
end
DyDt = [DCCDt;DCODt;DCDDt;DCMDt;DCADt];
end
0 Comments
Sign in to comment.