bvp4c unsolvable

Hi everyone, I am trying to solve a system of coupled ODEs using bpv4c. The results I have obtained was not what I was expecting as the values are too weird. Also, an error (singular jacobian) will occur when I changed the initial guess.

Here's a script to obtain my results,

xmesh = linspace(0,2,10);
guess=[0 0.9 700 1000];
solinit=bvpinit(xmesh,guess);

sol=bvp4c(@odeEqn,@bc,solinit);

figure(1);
plot(sol.x,sol.y,'-')

Here's the function file for @odeEqn

Everything is well defined except Xgf,Xsf,Tgf and Tsf

function F = odeEqn(z,S)

Xgf=S(1);
Xsf=S(2);
Tgf=S(3);
Tsf=S(4);

dXgfdz=(A/Fgof)*(ps/4)*kf*(Cgf^n);
dXsfdz=(A/Fsof)*(-ps)*kf*(Cgf^n);
dTgfdz=(A*Sm*h*(Tsf-Tgf))/(Fgof*(CpM+Xgf*(2*CpW+CpC+CpM)));
dTsfdz=(A*((Sm*h*(Tsf-Tgf))+((Fgof/A)*dXgfdz*HOR1)))/(Fgof*(CpFE2+Xsf*(2*CpFE-CpFE2))+Fsup*Cpsup);

F=[dXgfdz; dXsfdz; dTgfdz; dTsfdz];

end

Here's the function file for the boundary conditions

function res = bc (Sa,Sb)

res = [Sa(1) 
       Sa(3)-723
       Sb(1)-1
       Sb(4)-1173];

end

Thank you in advance!

Answer 1

[Star Strider]

‘Everything is well defined except Xgf,Xsf,Tgf and Tsf’

And, unfortunately, ‘A’,‘Sm’, ‘h’, ‘Cgf’, and a bunch of others. They need to be passed as additional parameters unless all the code is inside another function. In that situation, they may inherit them from the outer function workspace, however they are best passed as additional parameters. In any event, they are not provided here.

The code otherwise appears to be correct (no obvious problems that I can see). If the other parameters are provided, I’ll again attempt to run this.

[Lee Ek Hao]

Here is the entire code for the function. Thanks for the help! really appreciate it!

function F = odeEqn(z,S)

Xgf=S(1);
Xsf=S(2);
Tgf=S(3);
Tsf=S(4);


%Explicit Information
Tgof=723; %Initial temperature of fuel(methane)-K
Fsof=3; %Initial Fe2O3 feed rate-mol/s
Fsup=(0.4/0.6)*Fsof; %Molar flow rate of support material-mol/s
R=8.314; %Gas constant-J/(mol.K)
d=0.5; %fuel reactor diameter-m
A=pi*(d^2)/4; %Cross-sectional area-m2
pg=0.2411; %gas density-kg/m3
vch4=0.0828; %Vol. flow rate of methane-m3/s
Fgof=1.247; %Initial molar flow rate of fuel(methane)-mol/s
pb=2200; %Bulk density of reactive material in reactor-kg/m3
ps=2.81; %Molar bulk density of reactive material-mol/m3

%Arrhenius
kof=1.31*(10^8); %Pre-exponential factor 
Ea=219900; %Activation energy-J/mol
kf=kof*exp(-Ea/(R*Tsf));

%Cgfn
Cgof=Fgof/vch4; %Initial concentration of fuel(methane)-mol/m3
Cgf=Cgof*((1-Xgf)/(1+(2*Xgf)))*(Tgof/Tgf);
n=0.681; %Order of reaction




%Energy balance parameter of fuel reactor
pp=4000; %Density of oxygen carrier-kg/m3
dp=75*(10^-6); %mean particle diameter-m
Ap=4*pi*(dp/2)^2; %Surface area of a particle-m2
Vp=(4/3)*pi*(dp/2)^3; %Volume of particle-m3
mp=Vp*pp;
Np=pb/mp;
Sm=Np*Ap;
Vbar=vch4/A; %Superficial gas velocity-m/s
mug=(0.001596+((3.439*10^-5)*Tgof)+((-8.14*10^-9)*Tgf^2))/1000; %gas viscosity-kg/(m.s)
Nre=(dp*Vbar*pg)/mug;
kg=-0.007215+(8.015*(10^-5))*Tgf+(5.477*10^-9)*Tgf^2+(-1.053*10^-11)*Tgf^3; %thermal conductivity of gas-W/(m.K)
h=(0.33*(Nre^(1/3))*kg)/dp;

CpM=36.154+(-0.05111)*Tgf+(2.214*(10^-4))*(Tgf^2)+(-1.8243*(10^-7))*(Tgf^3)+(4.898*(10^-11))*(Tgf^4);            %Cp of methane J/(mol.K)
CpW=33.763+(-5.945*(10^-3))*Tgf+(2.235*(10^-5))*(Tgf^2)+(-9.962*(10^-9))*(Tgf^3)+(1.097*(10^-12))*(Tgf^4);       %Cp of Steam
CpC=29.268+(-0.02236)*Tgf+(2.652*(10^-4))*(Tgf^2)+(-4.153*(10^-7))*(Tgf^3)+(2.005*(10^-10))*(Tgf^4);             %Cp of CO2
CpFE2=110.9362+(32.04714*(Tsf/1000))+(-9.192333*(Tsf/1000)^2)+(0.901506*(Tsf/1000)^3)+(5.433677/(Tsf/1000)^2);   %Cp of Fe2O3
CpFE=45.7512+(18.78553*(Tsf/1000))+(-5.952201*(Tsf/1000)^2)+(0.852779*(Tsf/1000)^3)+(-0.081265/(Tsf/1000)^2);    %Cp of FeO
Cpsup=146.5551+(35.91295*(Tsf/1000))+(-0.183978*(Tsf/1000)^2)+(0.031409*(Tsf/1000)^3)+(-3.659941/(Tsf/1000)^2);  %Cp of MgAl2O4

HOR1=(5034082686347630884407*Tsf)/22517998136852480000 - 766695369022715765625/(140737488355328*Tsf) + (64148441052234189*Tsf^2)/1125899906842624000 - (82760967131826871*Tsf^3)/6755399441055744000000 + (134900196956432657181*Tsf^4)/112589990684262400000000000000 + 31193507971733846684869577643645007001/362427180012640665600000000000000; 
%Heat of reaction in fuel reactor-J/mol


%Differential Equation
dXgfdz=(A/Fgof)*(ps/4)*kf*(Cgf^n);
dXsfdz=(A/Fsof)*(-ps)*kf*(Cgf^n);
dTgfdz=(A*Sm*h*(Tsf-Tgf))/(Fgof*(CpM+Xgf*(2*CpW+CpC+CpM)));
dTsfdz=(A*((Sm*h*(Tsf-Tgf))+((Fgof/A)*dXgfdz*HOR1)))/(Fgof*(CpFE2+Xsf*(2*CpFE-CpFE2))+Fsup*Cpsup);

F=[dXgfdz; dXsfdz; dTgfdz; dTsfdz];

end

[Star Strider]

That runs with:

xmesh = linspace(0,2,2600);

since it needs a larger mesh to avoid a singular Jacobian. It produces a result that approaches infinity at about x=1.85.

It appears to be unstable, so now that it runs without error, you need to correct the instability. I am not certain what you are doing, and what appears to be a model of a chemical reactor is outside my areas of expertise anyway, so I cannot help with that. Otherwise, I can help with the code if that appears to be a problem, however everything appears to work correctly, even if it is not producing the correct result.

[Lee Ek Hao]

Thank you for helping!

Could you help on the code to plot the graphs?

Xgf and Xsf vs x on the same graph

Tgf and Tsf vs x on the same graph

[Star Strider]

My pleasure!

Sure!

Here are all the possibilities:

figure(1);
plot(sol.x,sol.y(1:2,:),'-')
legend('Xgf','Xsf', 'Location','NW')

figure(2);
plot(sol.x,sol.y(3:4,:),'-')
legend('Tgf','Tsf', 'Location','NW')

figure(3)
subplot(2,1,1)
plot(sol.x,sol.y(1:2,:),'-')
legend('Xgf','Xsf', 'Location','NW')
subplot(2,1,2)
plot(sol.x,sol.y(3:4,:),'-')
legend('Tgf','Tsf', 'Location','NW')

figure(4)
subplot(2,1,1)
plot(sol.x,real(sol.y(1:2,:)),'-')
hold on
plot(sol.x,imag(sol.y(1:2,:)),'--')
hold off
legend('\Re{Xgf}','\Re{Xsf}','\Im{Xgf}','\Im{Xsf}', 'Location','NW')
subplot(2,1,2)
plot(sol.x,real(sol.y(3:4,:)),'-')
hold on
plot(sol.x,imag(sol.y(3:4,:)),'--')
hold off
legend('\Re{Tgf}','\Re{Tsf}','\Im{Tgf}','\Im{Tsf}', 'Location','NW')

.

[Lee Ek Hao]

Thank you so much!!!!

[Star Strider]

As always, my pleasure!