Info

이 질문은 마감되었습니다. 편집하거나 답변을 올리려면 질문을 다시 여십시오.

ODE System not Functioning Correctly

조회 수: 2 (최근 30일)
Can Cakiroglu
Can Cakiroglu 2019년 5월 26일
마감: MATLAB Answer Bot 2021년 8월 20일
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);

답변 (0개)

이 질문은 마감되었습니다.

태그

제품

Community Treasure Hunt

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

Start Hunting!

Translated by