No values for RK method

27 次查看（过去 30 天）
Lee Ek Hao
Lee Ek Hao 2021-3-6
编辑： KALYAN ACHARJYA 2021-3-6
I have tried to solve a set of coupled first order differential using RK method. However, there is nothing being plotted on my graph. After checking all my codes and values, I still cannot understand where have I go wrong. Would really appreciate if someone is able to look into it.
clear all, clc
%Constants
pg=2.41;
Fgo=1.247;
Fso=4.8;
vgo=0.739;
Cgo=24.446;
Fsup=0.4/0.6*Fso;
Tgo=723;
d=0.5;
A=pi*(d^2)/4;
pb=2200;
ps=(pb)*0.4/160;
ko=1.31*(10^8);
Ea=219900;
R=8.314;
n=0.681;
b=4;
dp=75*(10^-6);
Ap=pi*(dp^2);
Vp=(pi/6)*(dp^3); %Vol. of particle m3
pp=4000; %Particle density kg/m3
Sm=(pb*Ap)/(pp*Vp); %Heat transfer per unit vol. m2/m3
vbar=0.0003; %Superficial gas velocity m/s
%Heat of reaction J/mol
%HOR=(5034082686347630884407*Ts)/22517998136852480000 - 766695369022715765625/(140737488355328*Ts) + (64148441052234189*Ts^2)/1125899906842624000 - (82760967131826871*Ts^3)/6755399441055744000000 + (134900196956432657181*Ts^4)/112589990684262400000000000000 + 31193507971733846684869577643645007001/362427180012640665600000000000000;
%Differential Eq
%dXgdz=(A*ps*k*(Cg^n))/(b*Fgo);
% =(A*ps*(ko*exp(-Ea/(R*Ts)))*(((Cgo*(1-Xg)*Tgo)/((1+(2*Xg))*Tg))^n))/(b*Fgo);
%dXsdz=-A*ps*k*(Cg^n)/Fso;
% =-A*ps*(ko*exp(-Ea/(R*Ts)))*(((Cgo*(1-Xg)*Tgo)/((1+(2*Xg))*Tg))^n))/Fso;
%dTgdz=(A*Sm*h*(Ts-Tg))/(Fgo*(CpM+Xg*(2*CpW+CpC+CpM)));
% =(A*Sm*(((-0.007215+(8.015*(10^-5))*Tg+(5.477*10^-9)*Tg^2+(-1.053*10^-11)*Tg^3)/dp)*0.33*((dp*vbar*pg/((0.002545+(4.549*(10^-5))*Tg+(-8.649*(10^-9))*(Tg^2))))^(1/3)))*(Ts-Tg))/(Fgo*(CpM+Xg*(2*CpW+CpC+CpM)));
%dTsdz=((A*Sm*h*(Ts-Tg))+(Fgo*dXgdz*HOR))/(Fgo*(CpFE2+Xs*(2*CpFE-CpFE2))+Fsup*Cpsup);
% =((A*Sm*(((-0.007215+(8.015*(10^-5))*Tg+(5.477*10^-9)*Tg^2+(-1.053*10^-11)*Tg^3)/dp)*0.33*((dp*vbar*pg/((0.002545+(4.549*(10^-5))*Tg+(-8.649*(10^-9))*(Tg^2))))^(1/3)))*(Ts-Tg))+(Fgo*fXg*HOR))/(Fgo*(CpFE2+Xs*(2*CpFE-CpFE2))+Fsup*Cpsup);
%RK method
fXg=@(z,Ts,Xg,Tg) (A*ps*(ko*exp(-Ea/(R*Ts)))*(((Cgo*(1-Xg)*Tgo)/((1+(2*Xg))*Tg))^n))/(b*Fgo);
fXs=@(z,Ts,Xg,Tg) -(A*ps*(ko*exp(-Ea/(R*Ts)))*(((Cgo*(1-Xg)*Tgo)/((1+(2*Xg))*Tg))^n))/Fso;
fTg=@(z,Ts,Xg,Tg) (A*Sm*(((-0.007215+(8.015*(10^-5))*Tg+(5.477*10^-9)*Tg^2+(-1.053*10^-11)*Tg^3)/dp)*0.33*((dp*vbar*pg/((0.002545+(4.549*(10^-5))*Tg+(-8.649*(10^-9))*(Tg^2))/1000))^(1/3)))*(Ts-Tg))/(Fgo*((36.154+(-0.05111)*Tg+(2.214*(10^-4))*(Tg^2)+(-1.8243*(10^-7))*(Tg^3)+(4.898*(10^-11))*(Tg^4))+Xg*(2*(33.763+(-5.945*(10^-3))*Tg+(2.235*(10^-5))*(Tg^2)+(-9.962*(10^-9))*(Tg^3)+(1.097*(10^-12))*(Tg^4))+(29.268+(-0.02236)*Tg+(2.652*(10^-4))*(Tg^2)+(-4.153*(10^-7))*(Tg^3)+(2.005*(10^-10))*(Tg^4))+(36.154+(-0.05111)*Tg+(2.214*(10^-4))*(Tg^2)+(-1.8243*(10^-7))*(Tg^3)+(4.898*(10^-11))*(Tg^4)))));
fTs=@(z,Ts,Tg,Xs,Xg) ((A*Sm*(((-0.007215+(8.015*(10^-5))*Tg+(5.477*10^-9)*Tg^2+(-1.053*10^-11)*Tg^3)/dp)*0.33*((dp*vbar*pg/((0.002545+(4.549*(10^-5))*Tg+(-8.649*(10^-9))*(Tg^2))/1000))^(1/3)))*(Ts-Tg))+(Fgo*(A*ps*(ko*exp(-Ea/(R*Ts)))*(((Cgo*(1-Xg)*Tgo)/((1+(2*Xg))*Tg))^n))/(b*Fgo)*((5034082686347630884407*Ts)/22517998136852480000 - 766695369022715765625/(140737488355328*Ts) + (64148441052234189*Ts^2)/1125899906842624000 - (82760967131826871*Ts^3)/6755399441055744000000 + (134900196956432657181*Ts^4)/112589990684262400000000000000 + 31193507971733846684869577643645007001/362427180012640665600000000000000)))/(Fgo*((110.9362+(32.04714*(Ts/1000))+(-9.192333*(Ts/1000)^2)+(0.901506*(Ts/1000)^3)+(5.433677/(Ts/1000)^2))+Xs*(2*(45.7512+(18.78553*(Ts/1000))+(-5.952201*(Ts/1000)^2)+(0.852779*(Ts/1000)^3)+(-0.081265/(Ts/1000)^2))-(110.9362+(32.04714*(Ts/1000))+(-9.192333*(Ts/1000)^2)+(0.901506*(Ts/1000)^3)+(5.433677/(Ts/1000)^2))))+Fsup*(146.5551+(35.91295*(Ts/1000))+(-0.183978*(Ts/1000)^2)+(0.031409*(Ts/1000)^3)+(-3.659941/(Ts/1000)^2)));
%Initial conditions
z(1)=0;
Xg(1)=0;
Xs(1)=1;
Tg(1)=723;
Ts(1)=1173;
%Step size
ss=0.0001;
zfinal=1;
N=ceil(zfinal/ss);
%Update loop
for i=1:N
%Update length
z(i+1)=z(i)+ss;
%Update variables
k1Xg=fXg(z(i) ,Ts(i) ,Xg(i) ,Tg(i));
k1Xs=fXs(z(i) ,Ts(i) ,Xg(i) ,Tg(i));
k1Tg=fTg(z(i) ,Ts(i) ,Xg(i) ,Tg(i));
k1Ts=fTs(z(i) ,Ts(i) ,Xg(i) ,Tg(i) ,Xs(i)) ;
k2Xg=fXg(z(i)+ss/2,Ts(i)+ss/2*k1Ts,Xg(i)+ss/2*k1Xg,Tg(i)+ss/2*k1Tg) ;
k2Xs=fXs(z(i)+ss/2,Ts(i)+ss/2*k1Ts,Xg(i)+ss/2*k1Xg,Tg(i)+ss/2*k1Tg) ;
k2Tg=fTg(z(i)+ss/2,Ts(i)+ss/2*k1Ts,Xg(i)+ss/2*k1Xg,Tg(i)+ss/2*k1Tg) ;
k2Ts=fTs(z(i)+ss/2,Ts(i)+ss/2*k1Ts,Xg(i)+ss/2*k1Xg,Tg(i)+ss/2*k1Tg,Xs(i)+ss/2*k1Xs);
k3Xg=fXg(z(i)+ss/2,Ts(i)+ss/2*k2Ts,Xg(i)+ss/2*k2Xg,Tg(i)+ss/2*k2Tg) ;
k3Xs=fXs(z(i)+ss/2,Ts(i)+ss/2*k2Ts,Xg(i)+ss/2*k2Xg,Tg(i)+ss/2*k2Tg) ;
k3Tg=fTg(z(i)+ss/2,Ts(i)+ss/2*k2Ts,Xg(i)+ss/2*k2Xg,Tg(i)+ss/2*k2Tg) ;
k3Ts=fTs(z(i)+ss/2,Ts(i)+ss/2*k2Ts,Xg(i)+ss/2*k2Xg,Tg(i)+ss/2*k2Tg,Xs(i)+ss/2*k2Xs);
k4Xg=fXg(z(i)+ss ,Ts(i)+ss *k3Ts,Xg(i)+ss *k3Xg,Tg(i)+ss *k3Tg) ;
k4Xs=fXs(z(i)+ss ,Ts(i)+ss *k3Ts,Xg(i)+ss *k3Xg,Tg(i)+ss *k3Tg) ;
k4Tg=fTg(z(i)+ss ,Ts(i)+ss *k3Ts,Xg(i)+ss *k3Xg,Tg(i)+ss *k3Tg) ;
k4Ts=fTs(z(i)+ss ,Ts(i)+ss *k3Ts,Xg(i)+ss *k3Xg,Tg(i)+ss *k3Tg,Xs(i)+ss *k3Xs);
Xg(i+1)=Xg(i)+ss/6*(k1Xg + 2*k2Xg + k2Xg*k3Xg + k4Xg);
Xs(i+1)=Xs(i)+ss/6*(k1Xs + 2*k2Xs + k2Xs*k3Xs + k4Xs);
Tg(i+1)=Tg(i)+ss/6*(k1Tg + 2*k2Tg + k2Tg*k3Xg + k4Tg);
Ts(i+1)=Ts(i)+ss/6*(k1Ts + 2*k2Ts + k2Ts*k3Ts + k4Ts);
end
%Plot the solution
figure (1); clf(1)
hold on
plot(z,Xg,'-r','displayname','Xg')
plot(z,Xs,'-b','displayname','Xs')
legend;
ylabel('conversion');
xlabel('Length/m');
hold off
figure (2); clf(2)
hold on
plot(z,Tg,'-r','displayname','Tg')
plot(z,Ts,'-b','displayname','Ts')
legend;
ylabel('Temperature');
xlabel('Length/m');
hold off

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

回答（1 个）

KALYAN ACHARJYA
KALYAN ACHARJYA 2021-3-6
编辑：KALYAN ACHARJYA 2021-3-6
There is no numerical data within Xg and Xs (all are NaN), same for 2nd figure also, hence there is no plots. Your next goal should be to find out why Xg, Xs or others have NaN data, may be because wrong implementation expression. Check carefully, this is the best way to learn.
Good Luck! :)

类别

在 Help Center 和 File Exchange 中查找有关 Line Plots 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by