ode45 not solving past the initial conditions (reactor dynamics)

3 views (last 30 days)
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
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);

More Answers (0)

Tags

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!