ODE System not Functioning Correctly

1 次查看（过去 30 天）

显示更早的评论

Can Cakiroglu 2019-5-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/464093-ode-system-not-functioning-correctly

关闭： MATLAB Answer Bot 2021-8-20

在 MATLAB Online 中打开

Hello Matlab Community,

I'm a chemical engineering student, and I have to solve an ODE system. However, while one of my values change, the other values seem to stay constant, although they are deffinetly functions that are not constant. Could there be any problems in the calculation of Matlab? And if yes, how could it be solved?

The function codes are as following:

Function File

function func = fchem(W, Variables)
X = Variables(1);
T = Variables(2);
P = Variables(3);
Fao = 1487.351; % kmol/h
Tr = 273; % K
X0 = 0;
P0 = 150; %bar
T0 = 673; %Kelvin
a = (3 * 1.5) / ((3 * 1.5) + 1);
b1 = (0.5 + (0.5 / 1.5)) / 1;
b2 = (-0.5+(1.5*1.5))/1;
R = 8.314 / 1000; % kJ/mol*K
K30 = 0.0307;
k20 = 2.12 * (10 ^ 13);  % should be 2.12e13 I suspect
E2 = 72100.82; % kJ/mol
E3 = -80989.98; % kJ/mol
omega = 1.564;
alpha = 0.640;
 
Bo = 0.254626; %atm/m
rho_c = 1600; %kg/m^3
Phi = 0.5;
D = 0.05; %Diameter of tube, m
Ac = pi * (D^2)/4; %m^2
alpha_pressure = (2*Bo)/Ac*rho_c*(1-Phi)*P0;
c1= 2.691122;
c2= 5.519265;
c3= 1.848863;
c4= 2001.6;
c5= 2.6899;
k2o = k20*exp(-E2/(R*T));
z = 0.5*X-1;
gamma = 1.5;
a_NH3 = P*z;
a_N2 = P*(a/3*gamma)*(1-b2*z);
a_H2 = P*a*(1-b1*z);
Ka = (1/T^c1)*10^(-c2*10^(-5)*T + c3 * 10^(-7)* (T^2) + (c4/T)  + c5);
K3 = K30*exp(-E3/(R*T));
ra = (k2o*(a_N2*(Ka^2)-(a_NH3^2)/(a_H2^3)))/(1+K3*a_NH3/(a_H2^omega))^(2*alpha);
Cp_NH3 = 6.70 + 0.00630*T;
Cp_H2 = 6.62 + 0.00081*T;
Cp_N2 = 6.50 + 0.00100*T;
delta_H_Tr = -91630;
delta_Cp = 2*Cp_NH3-3*Cp_H2-Cp_N2;
delta_H = delta_H_Tr + (-12.96*(T-Tr)+0.00917*((T-Tr)^2)/2);
dXdW = -ra/Fao;
dTdW = (ra*delta_H)/(Fao*((Cp_N2+3*Cp_H2+2*Cp_NH3)+delta_Cp*X));
dPdW = ((-alpha_pressure/2)/P)*(T/T0)*(1-0.5*X);
func = [dXdW; dTdW; dPdW];
end

And this is the calling script

clc
 Wspan = [0 500]; % Range for Catalyst Weight
 y0 = [0.; 673.; 150]; % Initial values X,T,and P respectively
[W y]=ode45(@fchem,Wspan,y0);
z=size(y);
plot(W,y);