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

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)'.''

Answer 1

[Sam Chak]

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)

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')