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

5 views (last 30 days)
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);
Error using mupadengine/feval_internal
System does not seem to contain a first-order differential equation for the field 'D(Ch2o)(t)'.

Error in odeToVectorField>mupadOdeToVectorField (line 171)
T = feval_internal(symengine,'symobj::odeToVectorField',sys,x,stringInput);

Error in odeToVectorField (line 119)
sol = mupadOdeToVectorField(varargin);
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)'.''

Answers (1)

Sam Chak
Sam Chak on 20 Dec 2022
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)
ans = 1×10
1.0e-09 * 0.0528 0 0.1585 0 0 0 0.0028 0 0 0
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)
V = 
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);
Warning: Failure at t=9.538715e-07. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.388132e-21) at time t.
plot(z, y), grid on, xlabel('z')

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!