Differential Equation Problem using ODE45
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
NURUL AINA SYAHIRAH
le 2 Jan 2024
Commenté : Sam Chak
le 5 Fév 2024
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 commentaire
Réponse acceptée
Sam Chak
le 2 Jan 2024
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
2 commentaires
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!