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 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponses (1)
Sam Chak
le 20 Déc 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 commentaires
Voir également
Catégories
En savoir plus sur Equation Solving dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!