ode45 not solving past the initial conditions (reactor dynamics)
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Charles Morton
le 31 Déc 2021
Réponse apportée : Sulaymon Eshkabilov
le 31 Déc 2021
I first attepted to plot y(1) through y(6) but no data showed up on the figure. Then I checked the array [z, y] of the ode45 output and only the initial conditions showed up in row 1 with the rest of the array filled in with NaN. I am pretty new to ode45 but I got ode45 to solve a similar less complicated problem when I was first trying it out. Hopefully there is just some typo but I cant find it and I think this is my last option. Thanks!
zspan=[0 1];
z0=[3.1135; 0.76830; 0; 0; 0; 0.1617; 0];
[z, y]=ode45(@diffy, zspan, z0)
%y(1)=Fh2
%y(2)=Fco2
%y(3)=Fm
%y(4)=Fw
%y(5)=Fdme
%y(6)=Fco
%y(7)=H
%T=y(7)/(y(1)*cph+y(2)*cpco2+y(3)*cpm+y(4)*cpw+y(5)*cpdme+y(6)*cpco)
function dy=diffy(~,y)
%thermo kJ/mol
dh1=-56.33;
dh2=-21.255;
dh3=-40.9;
R=8.3145; %kJ/kmol/K
%reactor
eps=0.33;
rho=1900;
vdot=2.37;
Re=0.019;
Ri=0.01;
U=20;
Tv=475;
Permw=3*10^-7;
Permh=2.2*10^-7;
%Temp and cp kJ/kmol/K
cph=14.5*2;
cpco2=43.99;
cpm=87.49;
cpw=2.25*18;
cpdme=125000;
cpco=27.995*1.185;
T=y(7)/(y(1)*cph+y(2)*cpco2+y(3)*cpm+y(4)*cpw+y(5)*cpdme+y(6)*cpco);
%catalyst kinetics
k1=35.45*exp(-1.7069*10^4/R./T);
k2=8.2894*10^4*exp(-5.2940*10^4/R./T);
k3=7.3976*exp(-2.0436*10^4/R./T);
Kh=0.249*exp(3.4394*10^4/R./T);
Kco2=1.02*10^-7*exp(6.74*10^4/R./T);
Kco=7.99*10^-7*exp(5.81*10^4/R./T);
K1=exp(4213./T-5.752*log(T)-1.707*10^-3*T+2.682*10^-6*T.^2-7.232*10^-10.*T.^3+17.6);
K2=exp(4019./T+3.707*log(T)-2.783*10^-3*T+3.8*10^-7*T.^2-6.561*10^4./T.^3-26.64);
K3=10^(2167./T-0.5194*log10(T)+1.037*10^-3.*T-2.331*10^-7.*T.^2-1.2777);
%pressure conversion
Ph=y(1)./(vdot*R*T);
Pco2=y(2)./(vdot*R*T);
Pm=y(3)./(vdot*R*T);
Pw=y(4)./(vdot*R*T);
Pdme=y(5)./(vdot*R*T);
Pco=y(6)./(vdot*R*T);
%reactions
r1=k1*(Pco2*Ph*(1-Pm*Pw/K1/Pco2/Ph^3))/(1+Kco2*Pco2+Kco*Pco+sqrt(Kh*Ph))^3;
r2=k2*(Pm^2/Pw-Pdme/K2);
r3=k3*(Pw-(Pco2*Ph/K3/Pco))/(1+Kco2*Pco2+Kco*Pco+sqrt(Kh*Ph));
%perm
pw=0.05*Pw;
ph=0.05*Ph;
Jw=Permw*(Pw-pw);
Jh=Permh*(Ph-ph);
%deez
d=[0 0 0 0 0 0 0];
d(1)=-rho*(1-eps)*pi*(Re^2-Ri^2)*(3*r1-r3)-Jh*2*pi*Ri;
d(2)=-rho*(1-eps)*pi*(Re^2-Ri^2)*(r1-r3);
d(3)=rho*(1-eps)*pi*(Re^2-Ri^2)*(r1-2*r2);
d(4)=rho*(1-eps)*pi*(Re^2-Ri^2)*(r1-r3+r2)-Jw*2*pi*Ri;
d(5)=rho*(1-eps)*pi*(Re^2-Ri^2)*r2;
d(6)=-rho*(1-eps)*pi*(Re^2-Ri^2)*r3;
d(7)=(-dh1*r1-dh2*r2-dh3*r3)*rho*(1-eps)*pi*(Re^2-Ri^2)-U*(T-Tv)*2*pi*Re;
dy=d';
end
1 commentaire
Réponse acceptée
Sulaymon Eshkabilov
le 31 Déc 2021
The problem is not only with T = 0 and K2 = 0, there must be some other pitfalls with the formulations. Please check the formulations and also, if necessary, you may need to adjust error tolerances with odeset options, e.g.:
OPTs = odeset('reltol', 1e-10, 'abstol', 1e-16);
[z, y]=ode45(@diffy, zspan, z0, OPTs);
0 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!