No values for RK method

2 views (last 30 days)
Lee Ek Hao on 6 Mar 2021
Edited: KALYAN ACHARJYA on 6 Mar 2021
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

KALYAN ACHARJYA on 6 Mar 2021
Edited: KALYAN ACHARJYA on 6 Mar 2021
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! :)