I've tried to simulate a catalytic reactor by solving 5 second order differential non-linear equations. I've used 'odetovectorfield' which I'm not pretty sure about that
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
clc, clear, close all
syms Cch4(z) Ch2o(z) Cco(z) Ch2(z) Cco2(z)
P = 5*10^5; %[Pa]
T = 1123.15; %[K]
Rg=8.31447; %[J / mol*K]
rhocat = 1900; % [Kg_cat / m3_reactor]
e=1-(rhocat/2900);
Dp=0.02;
a=6*(1-e)/Dp;
sym={'-g' '--r' '-.b'};
cases={'isoT','adb','cooling'};
% R1 R2 R3
nu = [-1 , 0 , -1 %CH4
-1 , -1 , -2 %H2O
1 , -1 , 0 %CO
3 , 1 , 4 %H2
0 , 1 , 1]; %CO2
a_b=[1.932, 1.7888;
1.718 , 1.7884;
2.454 ,1.7370;
11.275 ,1.6936;
1.248, 1.8033];
%Kinetic constant parameters: pre-exp factor (1st column) and activation
%energy (2 column)
Par=[4.225e15 , 1.955e6 , 1.020e15 , 8.23e-5 , 6.12e-9 , 6.65e-4 , 1.77e5;
240.1 , 67.13 , 243.9 , -70.65 , -82.90 , -38.28 , 88.68 ]; %[kJ/mol]
y_in=[50, 150, 0, 2.63, 0, T]; %[mol/s], [K], vector of initial values
C0=y_in(1:5)./((sum(y_in(1:5)))*Rg*T*P); %[mol/m^3] initial concentration
dC=[0 0 0 0 0];
BC=[C0 dC];
C=[Cch4(z),Ch2o(z),Cco(z),Ch2(z),Cco2(z),T];
Di=(a_b(:,1).*(10^-10)).*(C(6).^(a_b(:,2)));
p = (Rg*C(6)).*C(1:5);
%specific heats
cp= 1.386*C(6) + 1648 ;
%Entalpy of reaction for each reaction [kJ/mol]
H=[-4.47e13*C(6).^(-4.459)+226.9;
-271.4.*C(6).^(-0.2977);
99.52.*C(6).^(0.0937)];
% KIN, EQ and absorptionEQ constants
Kp1 = exp(30.481-27.187e3/C(6));
Kp2 = exp(-3.924+4.291e3/C(6));
Kp3 = exp(26.891-23.258e3/C(6));
Kp=[Kp1,Kp2,Kp3];
for i=1:3
k(i)=Par(1,i)*exp(-Par(2,i)/(Rg*C(6)));
end
for i=1:4
Keq(i)=Par(1,i+3)*exp(-Par(2,i+3)/(Rg*C(6)));
end
%Rate of rxns
DEN = 1+Keq(1)*p(3)+Keq(2)*p(4)+Keq(3)*p(1)+Keq(4)*p(2)/p(4);
Ri = 1000/3600.*[k(1)*(p(1)*p(2)-(p(4))^3*p(3)/Kp(1))./((p(4)^2.5).*DEN^2);
k(2)*(p(3)*p(2)-p(4)*p(5)/Kp(2))./(p(4).*DEN^2);
k(3)*(p(1)*(p(2))^2-(p(4))^4*p(5)/Kp(3))./((p(4)^3.5).*DEN^2)].*(rhocat/a); %[mol/s/m^3]
r = nu*Ri;
eq1=diff(C(1),2)==(4*r(1))/(Di(1)*0.125);
eq2=diff(C(2),2)==(4*r(2))/(Di(2)*0.125);
eq3=diff(C(3),2)==(4*r(3))/(Di(3)*0.125);
eq4=diff(C(4),2)==(4*r(4))/(Di(4)*0.125);
eq5=diff(C(5),2)==(4*r(5))/(Di(5)*0.125);
Vars=[Cch4(z);Ch2o(z);Cco(z);Ch2(z);Cco2(z)];
V=odeToVectorField(eq1,eq2,eq3,eq4,eq5);
M=matlabFunction(V,'Vars',{'z','Y'});
[z,y]=ode45(M,[0,8.1487],BC);
plot(z,y);
unfortunately ive encounter the error below in the odetovectorfield line and i couldnt figure out the problem so if someone has any idea that would be really helpfull for me to do this assignment.
''System does not seem to contain a first-order differential equation for the field 'D(Ch2o)(t)'.''
0 Kommentare
Antworten (1)
Sam Chak
am 20 Dez. 2022
Hi @Mahdi
The system can be run now. However, the simulation shows that the system blows up, which may indicate that the system contains a singularity somewhere. On this part, you need to check your equations.
syms Cch4(z) Ch2o(z) Cco(z) Ch2(z) Cco2(z) z
P = 5*10^5; % [Pa]
T = 1123.15; % [K]
Rg = 8.31447; % [J / mol*K]
rhocat = 1900; % [Kg_cat / m3_reactor]
e = 1-(rhocat/2900);
Dp = 0.02;
a = 6*(1-e)/Dp;
sym = {'-g' '--r' '-.b'};
cases = {'isoT','adb','cooling'};
% R1 R2 R3
nu = [-1 , 0 , -1 %CH4
-1 , -1 , -2 %H2O
1 , -1 , 0 %CO
3 , 1 , 4 %H2
0 , 1 , 1]; %CO2
a_b = [1.932, 1.7888;
1.718, 1.7884;
2.454, 1.7370;
11.275, 1.6936;
1.248, 1.8033];
% Kinetic constant parameters: pre-exp factor (1st column) and activation
% energy (2 column)
Par = [4.225e15, 1.955e6, 1.020e15, 8.23e-5, 6.12e-9, 6.65e-4, 1.77e5;
240.1, 67.13, 243.9, -70.65, -82.90, -38.28, 88.68]; %[kJ/mol]
y_in = [50, 150, 0, 2.63, 0, T]; %[mol/s], [K], vector of initial values
C0 = y_in(1:5)./((sum(y_in(1:5)))*Rg*T*P); %[mol/m^3] initial concentration
dC = [0 0 0 0 0];
BC = [C0(1) dC(1) C0(2) dC(2) C0(3) dC(3) C0(4) dC(4) C0(5) dC(5)];
double(BC)
C = [Cch4(z), Ch2o(z), Cco(z), Ch2(z), Cco2(z), T];
Di = (a_b(:,1).*(10^-10)).*(C(6).^(a_b(:,2)));
p = (Rg*C(6)).*C(1:5);
% specific heats
cp = 1.386*C(6) + 1648;
% Entalpy of reaction for each reaction [kJ/mol]
H = [-4.47e13*C(6).^(-4.459)+226.9;
-271.40.*C(6).^(-0.2977);
99.520.*C(6).^(0.0937)];
% KIN, EQ and absorptionEQ constants
Kp1 = exp(30.481 - 27.187e3/C(6));
Kp2 = exp(-3.924 + 4.2910e3/C(6));
Kp3 = exp(26.891 - 23.258e3/C(6));
Kp = [Kp1, Kp2, Kp3];
for i = 1:3
k(i) = Par(1,i)*exp(-Par(2,i)/(Rg*C(6)));
end
for i = 1:4
Keq(i) = Par(1,i+3)*exp(-Par(2,i+3)/(Rg*C(6)));
end
% Rate of rxns
DEN = 1+Keq(1)*p(3)+Keq(2)*p(4)+Keq(3)*p(1)+Keq(4)*p(2)/p(4);
Ri = 1000/3600.*[k(1)*(p(1)*p(2)-(p(4))^3*p(3)/Kp(1))./((p(4)^2.5).*DEN^2);
k(2)*(p(3)*p(2)-p(4)*p(5)/Kp(2))./(p(4).*DEN^2);
k(3)*(p(1)*(p(2))^2-(p(4))^4*p(5)/Kp(3))./((p(4)^3.5).*DEN^2)].*(rhocat/a); %[mol/s/m^3]
r = nu*Ri;
% Differential equations
eq1 = diff(Cch4,z,2)==(4*r(1))/(Di(1)*0.125);
eq2 = diff(Ch2o,z,2)==(4*r(2))/(Di(2)*0.125);
eq3 = diff(Cco,z,2)==(4*r(3))/(Di(3)*0.125);
eq4 = diff(Ch2,z,2)==(4*r(4))/(Di(4)*0.125);
eq5 = diff(Cco2,z,2)==(4*r(5))/(Di(5)*0.125);
V = odeToVectorField(eq1, eq2, eq3, eq4, eq5)
Vars = [Cch4(z); Ch2o(z); Cco(z); Ch2(z); Cco2(z)];
M = matlabFunction(V, 'Vars', {'z', 'Y'});
[z, y] = ode45(M, [0, 8.1487], BC);
plot(z, y), grid on, xlabel('z')
0 Kommentare
Siehe auch
Kategorien
Mehr zu Equation Solving finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!