Differential Equation Problem using ODE45

1 次查看（过去 30 天）
NURUL AINA SYAHIRAH
评论： Sam Chak 2024-2-5
In below script, I was using the equation of dF/dz [F=CO2, CO, H2O, H2, DME, CH3OH] which written as A(1) until A(6) to solve the molar flowrate of each component. It was assumed that the feed only contained CO2 and H2 in the inlet flowrate. The mole fraction of CO2 and H2 was assumed to 0.25 and 0.75 with a total flowrate of 0.2667 kmol/s [960 kmol/hr].
In this case, the script was successfully run. However, the mole flowrate obtained for each component was NaN that results to no reaction occured/error. Appreciate if you could help me and give suggestion to improve the script. Thank you in advance.
The script:
function A = test_A(Z,F)
%Fixed parameters
Ptot=50;
T=523;
pc=1900;
%Defining constant
FT=F(1)+F(2)+F(3)+F(4)+F(5)+F(6)+F(7);
RT=523*8.314;
%Feed based on inlet as 0.2667 kmol/s
F(1)=0.25*0.2667;
F(2)=0.75*0.2667;
F(3)=0*0.2667;
F(4)=0*0.2667;
F(5)=0*0.2667;
F(6)=0*0.2667;
F(7)=0*0.2667;
%Partial pressure for each component
pCO2=(F(1)/FT)*Ptot;
pH2=(F(2)/FT)*Ptot;
pCO=(F(3)/FT)*Ptot;
pH2O=(F(4)/FT)*Ptot;
pDME=(F(5)/FT)*Ptot;
pCH3OH=(F(6)/FT)*Ptot
pN2=(F(7)/FT)*Ptot;
%Kinetic parameters
k1=35.45*exp(-1.7069e4/(RT));
k2=7.3976*exp(-2.0436e4/(RT));
k3=8.2894e4*exp(-5.2940e4/(RT));
%Adsorption constant
kH2=0.249*exp(3.4394e4/(RT));
kCO2=1.02*10^-7*exp(6.74*10^4/(RT));
kCO=7.99*10^-7*exp(5.81*10^4/(RT));
%Equilibrium constant
kp1=exp((4213/T)-(5.752*log(T))-(1.707e-3*T)+(2.682e-6*(T^2))-(7.232e-10*(T^3))+17.6);
kp2=exp((2167/T)-(0.5194*log(T))+(1.037e-3*T)-(2.331e-7*(T^2))-1.2777);
kp3=exp((4019/T)+(3.707*log(T))-(2.783e-3*T)+(3.81e-7*(T^2))-(6.561e4/(T^3))-26.64);
%Rate of Reaction
r1=k1*(pCO2*pH2*(1-((pCH3OH*pH2O)/(kp1*pCO2*pH2^3))))/(1+kCO2*pCO2+kCO*pCO+sqrt(kH2*pH2))^3;
r2=k2*((((pCO2*pH2)/(kp2*pCO))-pH2O)/((1+(kCO2*pCO2)+(kCO*pCO)+(sqrt(kH2*pH2)))^3));
r3=k3*((pCH3OH^2/pH2O)-(pDME/kp3)); %MeOH dehydration
%Trans-membrane molar flux
permN2=1.3333e-8;
permH2O=1e-6;
prN2=0*Ptot;
prH2O=0*Ptot;
ppN2=0.79320;
ppH2O=0.20680;
JN2=permN2*(prN2-pN2);
JH2O=permH2O*(prH2O-pH2O);
%Mass balance for each component
vf=0.33;
Dsi=0.0198;
Dmo=14e-3;
%Mass balance for each component
A(1)=(pc*(1-vf)*(r1-r2+r3)*(pi*Dmo^2)-(JH2O*pi*Dmo));
A(2)=(-pc*(1-vf)*(3*r1-r2)*((pi/4)*(Dsi^2-Dmo^2)));
A(3)=(pc*(1-vf)*(r2)*((pi/4)*(Dsi^2-Dmo^2)));
A(4)=(-pc*(1-vf)*(r1-r2)*((pi/4)*(Dsi^2-Dmo^2)));
A(5)=(pc*(1-vf)*(r1-2*r3)*((pi/4)*(Dsi^2-Dmo^2)));
A(6)=(pc*(1-vf)*(r3)*((pi/4)*(Dsi^2-Dmo^2)));
A=A';
  1 个评论
Torsten
Torsten 2024-1-2
编辑：Torsten 2024-1-2
You didn't include the script, but only the function to be used by ode45.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Sam Chak
Sam Chak 2024-1-2
Ran in:
Hi @NURUL AINA SYAHIRAH
First things first, check the values computed by the differential equations. All differential states return either NaN or Inf because the equations contain r2, or r3, or both. In case you don't know, r2 returns Inf due to a division by zero problem (pCO = 0), and r3 returns NaN due to the indeterminate term in pCH3OH^2/pH2O. Fix these two issues and the code should run.
Zspan = [0 1];
F0 = [1 1 1 1 1 1];
[Z, F] = ode45(@test_A, Zspan, F0);
plot(Z, F)
function A = test_A(Z,F)
%Fixed parameters
Ptot=50;
T=523;
pc=1900;
%Feed based on inlet as 0.2667 kmol/s
F(1)=0.25*0.2667;
F(2)=0.75*0.2667;
F(3)=0*0.2667;
F(4)=0*0.2667;
F(5)=0*0.2667;
F(6)=0*0.2667;
F(7)=0*0.2667;
%Defining constant
FT=F(1)+F(2)+F(3)+F(4)+F(5)+F(6)+F(7);
RT=523*8.314;
%Partial pressure for each component
pCO2=(F(1)/FT)*Ptot; % 12.5
pH2=(F(2)/FT)*Ptot; % 37.5
pCO=(F(3)/FT)*Ptot + 0.1; % 0 plus a small non-zero value
pH2O=(F(4)/FT)*Ptot + 0.1; % 0 plus a small non-zero value
pDME=(F(5)/FT)*Ptot; % 0
pCH3OH=(F(6)/FT)*Ptot; % 0
pN2=(F(7)/FT)*Ptot; % 0
%Kinetic parameters
k1=35.45*exp(-1.7069e4/(RT));
k2=7.3976*exp(-2.0436e4/(RT));
k3=8.2894e4*exp(-5.2940e4/(RT));
%Adsorption constant
kH2=0.249*exp(3.4394e4/(RT));
kCO2=1.02*10^-7*exp(6.74*10^4/(RT));
kCO=7.99*10^-7*exp(5.81*10^4/(RT));
%Equilibrium constant
kp1=exp((4213/T)-(5.752*log(T))-(1.707e-3*T)+(2.682e-6*(T^2))-(7.232e-10*(T^3))+17.6);
kp2=exp((2167/T)-(0.5194*log(T))+(1.037e-3*T)-(2.331e-7*(T^2))-1.2777);
kp3=exp((4019/T)+(3.707*log(T))-(2.783e-3*T)+(3.81e-7*(T^2))-(6.561e4/(T^3))-26.64);
%Rate of Reaction
r1=k1*(pCO2*pH2*(1-((pCH3OH*pH2O)/(kp1*pCO2*pH2^3))))/(1+kCO2*pCO2+kCO*pCO+sqrt(kH2*pH2))^3;
r2=k2*( (((pCO2*pH2)/(kp2*pCO)) - pH2O) / ((1 + kCO2*pCO2 + kCO*pCO + sqrt(kH2*pH2))^3) );
r3=k3*((pCH3OH^2/pH2O)-(pDME/kp3)); %MeOH dehydration
%Trans-membrane molar flux
permN2=1.3333e-8;
permH2O=1e-6;
prN2=0*Ptot;
prH2O=0*Ptot;
ppN2=0.79320;
ppH2O=0.20680;
JN2=permN2*(prN2-pN2);
JH2O=permH2O*(prH2O-pH2O);
%Mass balance for each component
vf=0.33;
Dsi=0.0198;
Dmo=14e-3;
%Mass balance for each component
A(1)=( pc*(1-vf)*(r1 - r2 + r3)*(pi*Dmo^2)-(JH2O*pi*Dmo)); % NaN
A(2)=(-pc*(1-vf)*(3*r1 - r2) *((pi/4)*(Dsi^2-Dmo^2))); % Inf
A(3)=( pc*(1-vf)*(r2)*((pi/4) *(Dsi^2-Dmo^2))); % Inf
A(4)=(-pc*(1-vf)*(r1 - r2) *((pi/4)*(Dsi^2-Dmo^2))); % Inf
A(5)=( pc*(1-vf)*(r1 - 2*r3) *((pi/4)*(Dsi^2-Dmo^2))); % NaN
A(6)=( pc*(1-vf)*(r3)*((pi/4) *(Dsi^2-Dmo^2))); % NaN
A=A';
end

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by