MATLAB: Solving simultaneous 3 ODE : I can’t use the results in the equations themselves.

chemical enginnering reaction kineticspacked bed reactorssimultaneous 3 ode solution

Hi everyone.
I want to solve 3 ODE equations simultaneously and I need to use the values calculated in the ODE itself. My variable is z, and dependent values are X, T, P. In the function file, I want to use the resulting X, T, P values but I obtained NaN values only. Also, if I write dep = dep' at the end of function file in order to convert it into column vector, I received an error that I should convert it to column vector.
Here is my main program :
clc , clear all
global nu a b c d e f h A Ea R eps rocat D Mw dcat mugas f0
%%data
nu = [-1 -0.5 0 1]; % stoichiometric matrix
R = 8.314; % [J.mol-1.K]
T0 = 703.15; %feed temperature [K]
F0 = 1.2676e3; % [mol.s-1] calculated by hand 
P0 = 101e3; % feed pressure [Pa]
y0 = [0.095 0.115 0.79 0 ]; % initial feed composition
f0 = F0.*y0; % inlet feed concentration of components [mol.s-1]
robed = 5e2; %bed/bulk density [kg.m-3]
eps = 0.40; %porosity
rocat = robed/(1-eps); %catalyst density [kg.m-3]
dcat = 8e-3; % catalyst diameter [m]
mugas = 36.6*10^-6; %[Pa.s] 
Mw = [64.066 16.000 14.007 80.006]*10^-3; % molecular weight [kg.mol-1]
X0 = 0;
a = [21.43029 30.03235 19.50583 24.02503];
b = [74.35094 8.772972 19.88705 119.4607];
c = [-57.75217 -3.988133 -8.598535 -94.38686];
d = [16.35534 0.788313 1.369784 26.96237];
e = [0.086731 -0.741599 0.527601 -0.117517];
f = [-305.7688 -11.32468 -4.935202 -407.8526];
h = [-296.8422 0 0 -395.7654];
a = a';
b = b';
c = c';
d = d';
e = e';
f = f';
h = h';
A =  5.3428;
Ea = 7.9670e+04; %[J.mol-1.K-1]
%%initial condition    
dep0(1) = X0;
dep0(2) = T0;
dep0(3) = P0;
%%ODE solver
L = 5; % [m]

D = 0.5; % [m]
wcat = (3.14/4.*D.^2).*L./robed; % catalyst weight [kg]
zrange = [0 L];
[z, dep] = ode15s(@newodesolv, zrange, dep0);
X = dep(1);
T = dep(2);
P = dep(3);
fi = f0 - nu.*f0(1)*X;
fprintf('%12.4e %12.4e %12.4e %12.4e\n', wcat, X, T, P)
fprintf('%12.4e %12.4e %12.4e %12.4e\n', fi)
————————————————————————————————
And this is my function file.
————————————————————————————————
function dep = newodesolv(z,dep)
%dep(1) = X
%dep(2) = T
%dep(3) = P
global nu a b c d e f h A Ea R eps rocat D Mw dcat mugas f0
X = dep(1);
T = dep(2);
P = dep(3);
fi = f0 - nu.*f0(1)*dep(1);
y = fi/(fi(1)+fi(2)+fi(3)+fi(4));
Mwmix = (Mw*y'); % molecular weight of mixture
rogas = (P * Mwmix)/(R*T); %[kg.m-3]
t = T/1000; %[K]
cp = a+b.*t+c.*(t^2)+d.*(t^3)+e./(t^2); % [J.mol-1.K-1] 
Formation_ent = (a.*t+(b.*t^2)/2+(c.*t^3)/3+(d.*t^4)/4-(e./t)+ f - h)*1000; % [J.mol-1]

delta_H_rxn = (nu*Formation_ent)*10^3; % [J.mol-1]
k = A*exp(-Ea/(R*T))*1000/101325; %[mol.kg,cat-1.Pa-1]
Kp = 7.97*10^-5*exp(12100./T); % [MPa^-0.5]
ficp = fi*cp; %[J.K-1.s-1]
fiMw = fi*Mw'; % [kg.s-1] 
r = k.*((P*y(2)*y(1)^(0.5))/(y(4)^(0.5))-((y(4)^(1.5))/((y(1)^(1.5))*(Kp^2)*10^6))); % [mol.kg,cat-1.s-1]
dep(1) = r.*(1-eps).*rocat*3.14/4*(D.^2)./f0(1);
dep(2) = -(r.*delta_H_rxn.*(1-eps).*rocat.*(3.14/4)*(D.^2))/(ficp);
dep(3) = -(fiMw)./((rogas).*(3.14/4).*(D.^2).*dcat)*...
    ((1-eps)/(eps^3)).*((150*(1-eps).*mugas./dcat) + (1.75.*(fiMw)./(3.14/4*(D.^2))));
dep;
Could you please help me? Thank you.

Best Answer

First, remove the semicolons from your dep(1)...dep(3) statements. It could be that one or more of them are returning matrices. If so, you need to figure out why and how to fix them so they will return elements of column vectors.
The NaNs are usually the result of (0/0) calculations. One NaN in the wrong place can ruin all your calculations. You need to find out what is generating them and fix it.
I don’t know what you’re doing, or what you can or can’t change to get the result you want. Only you can determine effectively what’s wrong with your code and how to fix it.
I suggest you set zrange = [0 0.1 0.2]; until you get your code running successfully.
Related Question