ode45 not solving past the initial conditions (reactor dynamics)

2 次查看（过去 30 天）
Charles Morton
Charles Morton 2021-12-31
回答： Sulaymon Eshkabilov 2021-12-31
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

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

采纳的回答

Sulaymon Eshkabilov
Sulaymon Eshkabilov 2021-12-31
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 个）

类别

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

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by