No difference in return values when using ode45

Hi,

I've been trying to model the reaction of a Water Gas Shift reaction using some kinetics I found by using ode45. I have used it before but I must be doing something wrong because the values of x(1)-x(5) are not changing. My code looks like this:

% Script calling Moe Kinetics Function
%
% Ptot=8.85;   % (bar) ptot
% xco =.0499;    % Mole Fractions
% xco2=.0844;
% xh20=.2913;
% xh2 =.4532;
% xch4=.1212;
%
% Pco =.0501*8.85 = x(1) % Partial Pressures (bar)
% Pco2=.0830*8.85 = x(2)
% Ph2o=.2824*8.85 = x(3)
% Ph2 =.4485*8.85 = x(4)
% Pch4=.1360*8.85 = x(5)
% T = 628.35 K = x(6)
R = 8.3144598*(10^-5); % (m^3*bar)/(K*mol)
% Initial condition Partial Pressures (Pa) and Initial temperature (K)
ic=[.1212*8.85,.4532*8.85,.2913*8.85,.0844*8.85,.0499*8.85,628.35];
Vspan=[0:.001:10];
[V, x] = ode45('project_wgs',Vspan,ic);
% Concentrations of species
x1=x(:,1);
x2=x(:,2);
x3=x(:,3);
x4=x(:,4);
x5=x(:,5);
T=x(:,6);  % Temperature (K)
% Converting concentration to partial pressure (bar)
x11 = x1.*R.*T;
x22 = x2.*R.*T;
x33 = x3.*R.*T;
x44 = x4.*R.*T;
x55 = x5.*R.*T;

P = x11+x22+x33+x44+x55;  % Total Pressure

function [dCdV] = project_wgs (V,x)
R = 8.3144598*(10^-5); % (m^3*bar)/(K*mol)
q = 5.888*(10^5)/60;    % (m^3/min)
% Pco =.0501*8.85 = x(1) % Partial Pressures (bar)
% Pco2=.0830*8.85 = x(2)
% Ph2o=.2824*8.85 = x(3)
% Ph2 =.4485*8.85 = x(4)
% Pch4=.1360*8.85 = x(5)
% T = 628.35 K    = x(6)

% Rate equation in (mol/(m^3*min))
r = ((7.26*exp(-15429/R*x(6))*x(1)*x(3)) - (551.6*exp(-53482/R*x(6))*x(2)*x(4)))*1000*7633.65;

dCdV(1,1)  = r/q;
dCdV(2,1)  = r/q;
dCdV(3,1)  = r/q;
dCdV(4,1)  = r/q;
dCdV(5,1)  = r/q;
dCdV(6,1)  = x(6);

Answer 1

[Steven Lord]

One part of your calculation for r is exp(-15429/R*x(6)). With R around 1e-5 and x(6) initially around 600, this term in your expression for r is roughly exp(-1e11). That underflows to 0. Similarly, the other part of your calculation for r, exp(-53482/R*x(6)), also underflows. So r is equal to 0 and so is r/q.

[Jrevs]

Thank you, I fixed that issue but also found a few other problems with my math and fixed those. I am now using conversion (X) as my single variable along with temperature that I want ode45 to solve for. However, it returns as zero every time. I have tried adjust Fa0 incase that was the issue but I am stumped. I attached the code below

Any suggestions?

if true
  clear all, close all
% Script calling Moe Kinetics Function
%
% Initial total pressure, Ptot=8.85 (bar)
% Initial mole fractions
% xco =.0499 Mole Fractions
% xco2=.0844
% xh20=.2913
% xh2 =.4532
% xch4=.1212Pco =.0501*8.85;
Pco =.0501*8.85;
Pco2=.0830*8.85;
Ph2o=.2824*8.85;
Ph2 =.4485*8.85;
Pch4=.1360*8.85;
% Initial condition Partial Pressures (Pa) and Initial temperature (K)
ic=[0,628.35];
Vspan=[0:.1:100];
[V, x] = ode45('project_wgs',Vspan,ic);
% Partial Pressures
X=x(:,1);
T=x(:,2);
% x22=x(:,2);
% x33=x(:,3);
% x44=x(:,4);
% x55=x(:,5);
% T=x(:,6);  % Temperature (K)
P_CO = Pco*(1-X);
P_H2O = Ph2o-(Pco*X);
P_H2 = Ph2+(Pco*X);
P_CO2 = Pco2+(Pco*X);
P_CH4 = Pch4;
% Converting from partial pressure to concentration
% x11 = x(1)/(R*x(6));
% x22 = x(2)/(R*x(6));
% x33 = x(3)/(R*x(6));
% x44 = x(4)/(R*x(6));
% x55 = x(5)/(R*x(6));
% Total Pressure
P = P_CO+P_H2O+P_H2+P_CO2+P_CH4;
% W=Weigh of catalyst (Kg)
W = 7633.65*X;

end

[Jrevs]

if true
  function [dXdV] = project_wgs (V,x)
% Initial Partial Pressures (bar)
Pco =.0501*8.85;
Pco2=.0830*8.85;
Ph2o=.2824*8.85;
Ph2 =.4485*8.85;
Pch4=.1360*8.85;
Fa0=2.498*(10^4)*1000/60; % Molar Flow (gmol/min)
% Specific Heats (kj/kg*K)
Cpco=1.1;
Cpco2=1.9;
Cph2o=2.109;
Cph2=14.56;
Cpch4=3.34;
% Temperature (K)
T0 = 628.35;
% Gas Constant (m^3*Pa)/(K*mol)
R = 8.3144598;
X=x(1);
r = ((7.26*exp(-15429/R*x(2))*((Ph2o-(Pco*X))*Pco*(1-X))) - (551.6*exp(-53482/R*x(2))*(Ph2+(Pco*X))*(Pco2+(Pco*X))))*1000*7633.65;
x(2) = T0 + (((X*41.1)+(T0*((Cpco+((Ph2o/Pco)*Cph2o)+((Pco2/Pco)*Cpco2)+((Ph2/Pco)*Cph2)+((Pch4/Pco)*Cpch4)))))/(Cpco+((Ph2o/Pco)*Cph2o)+((Pco2/Pco)*Cpco2)+((Ph2/Pco)*Cph2)+((Pch4/Pco)*Cpch4)));
dXdV(1,1)  = r/Fa0;
dXdV(2,1)  = x(2); % Temperature (K)
end