ode45 not solving past the initial conditions (reactor dynamics)
3 views (last 30 days)
Show older comments
Charles Morton
on 31 Dec 2021
Answered: Sulaymon Eshkabilov
on 31 Dec 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
Accepted Answer
Sulaymon Eshkabilov
on 31 Dec 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 Comments
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!