Non physically fearsible results from a 10 chemical reaction system (8 ODEs)

Question

I'm trying to script a 10 chemical reaction system (image link below) with 8 components (so 8 ODEs) in order to get the reactor volume, but I'm getting negative values for most of the components. It may be something to do with the molar balances or the script itself so if anyone with some chemical background can look into it, I would apreciate it!

Data:

2,00e-06    48000   2,50e-13    -175000 0   0   0   0,76
23,20       182000  8,30e-14    -186000 2,30e-12    -124000 1   0,40
5,20e-07    68000   3,60e-14    -187000 0   0   0   0,85
0,00011     104000  4,00e-13    -168000 0   0   1   0,55
0,17        157000  4,50e-13    -166000 0   0   0   0,37
0,06        166000  1,60e-13    -211000 0   0   1   0,96
12000000    226000  0   0   0   0   0   0
9300        300000  0   0   0   0   0   0
0,00019     173000  0   0   0   0   1   1
0,026       220000  0   0   0   0   1   1

Function (OCML.m):

function dFdV = OCML( V , F0 , par , T , Pt)

%% Initial conditions

F = sum(F0);
y = F0./F;
pp = y * Pt;

p_CH4  = pp(1)   ;
p_O2   = pp(2)   ;
p_CO2  = pp(3)   ;
p_H2O  = pp(4)   ;
p_C2H6 = pp(5)   ;
p_CO   = pp(6)   ;
p_C2H4 = pp(7)   ;
p_H2   = pp(8)   ;

%% Parameters

ko = par(:,1)       ;
Ea = par(:,2)       ;
KCO2 = par(:,3)     ;
DHadCO2 = par(:,4)  ;
KO2 = par(:,5)      ;
DHadO2 = par(:,6)   ;
m = par(:,7)        ;
n = par(:,8)        ;
R = 8.314462        ;
ro_cat = 3600000    ; % g/m3

%% Stoichiometric coefficients matrix

N = [ -1   -2  1  2  0  0  0  0 ;
      -2 -0.5  0  1  1  0  0  0 ;
      -1   -1  0  1  0  1  0  1 ;
       0 -0.5  1  0  0  1  0  0 ;
       0 -0.5  0  1 -1  0  1  0 ;
       0   -2  0  2  0  2 -1  0 ;
       0    0  0  0 -1  0  1  1 ;
       0    0  0 -2  0  2 -1  4 ;
       0    0  1 -1  0 -1  0  1 ;
       0    0 -1  1  0  1  0 -1 ; ];

%% Reaction Rates

% Reaction 1

exp1 = exp(-Ea(1)/(R*T))                                ;
exp_ad1 = exp(-DHadCO2(1)/(R*T))                        ;
denom1 = 1 + KCO2(1) * exp_ad1 * p_CO2                  ;
r1 = ( ko(1) * exp1 * p_CH4^m(1) * p_O2^n(1) ) / denom1 ;

% Reaction 2

exp2 = exp(-Ea(2)/(R*T))                                                        ;
exp_ad2a = exp(-DHadO2(2)/(R*T))                                                ;
exp_ad2b = exp(-DHadCO2(2)/(R*T))                                               ;
denom2 = 1 + KO2(2) * exp_ad2a * p_O2 + KCO2(2) * exp_ad2b * p_CO2              ;
r2 = ( ko(2) * exp2 * ( KO2(2) * exp_ad2a * p_O2 )^n(2) * p_CH4 ) / denom2^2    ;

% Reaction 3

exp3 = exp(-Ea(3)/(R*T))                                ;
exp_ad3 = exp(-DHadCO2(3)/(R*T))                        ;
denom3 = 1 + KCO2(3) * exp_ad3 * p_CO2                  ;
r3 = ( ko(3) * exp3 * p_CH4^m(3) * p_O2^n(3) ) / denom3 ;

% Reaction 4

exp4 = exp(-Ea(4)/(R*T))                                ;
exp_ad4 = exp(-DHadCO2(4)/(R*T))                        ;
denom4 = 1 + KCO2(4) * exp_ad4 * p_CO2                  ;
r4 = ( ko(4) * exp4 * p_CO^m(4) * p_O2^n(4) ) / denom4  ;

% Reaction 5

exp5 = exp(-Ea(5)/(R*T))                                    ;
exp_ad5 = exp(-DHadCO2(5)/(R*T))                            ;
denom5 = 1 + KCO2(5) * exp_ad5 * p_CO2                      ;
r5 = ( ko(5) * exp5 * p_C2H6^m(5) * p_O2^n(5) ) / denom5    ;

% Reaction 6

exp6 = exp(-Ea(6)/(R*T))                                        ;
exp_ad6 = exp(-DHadCO2(6)/(R*T))                                ;
denom6 = 1 + KCO2(6) * exp_ad6 * p_CO2                          ;
r6 = ( ko(6) * exp6 * p_C2H4^m(6) * p_O2^n(6) ) / denom6    ;

% Reaction 7

exp7 = exp(-Ea(7)/(R*T))                ;
r7 = ko(7)* exp7 * p_C2H6  * ro_cat     ;

% Reaction 8

exp8 = exp(-Ea(8)/(R*T))                    ;
r8 = ko(8)* exp8 * p_C2H4^m(8) * p_H2O^n(8) ;

% Reaction 9

exp9 = exp(-Ea(9)/(R*T))                    ;
r9 = ko(9)* exp9 * p_CO^m(9) * p_H2O^n(9)   ;

% Reaction 10

exp10 = exp(-Ea(10)/(R*T))                     ;
r10 = ko(10)* exp10 * p_CO2^m(10) * p_H2^n(10) ;

r_vector = [ r1 r2 r3 r4 r5 r6 r7 r8 r9 r10 ]'  ;

%% Molar balances (dFx/dV)

dF_CH4dV  = sum(N(:,1).*r_vector)*ro_cat     ;
dF_O2dV   = sum(N(:,2).*r_vector)*ro_cat     ;
dF_CO2dV  = sum(N(:,3).*r_vector)*ro_cat     ;
dF_H2OdV  = sum(N(:,4).*r_vector)*ro_cat     ;
dF_C2H6dV = sum(N(:,5).*r_vector)*ro_cat     ;
dF_COdV   = sum(N(:,6).*r_vector)*ro_cat     ;
dF_C2H4dV = sum(N(:,7).*r_vector)*ro_cat     ;
dF_H2dV   = sum(N(:,8).*r_vector)*ro_cat     ;

dFdV = [ dF_CH4dV dF_O2dV dF_CO2dV dF_H2OdV dF_C2H6dV dF_COdV dF_C2H4dV dF_H2dV ]' ;
end

Script:

clear; clc; close all ; format long

%% Data

par=importdata('data.mat');

T = 273.15 + 850;
Pt = 101325;

% Methane and O2 feed
F0 = [ 27.77777778 9.259259259 0 0 0 0 0 0 ];

% Volume vector

V = 0.0:0.1:1 ;

% ODE solver

OCML_model = @(V,F0)[OCML(V,F0,par,T,Pt)]     ;
[ V, F_vector ] = ode15s(OCML_model,V,F0)      ;

Answer 1

if you have not mass-balance equations which involve fractions, like Michaëlis-Menten enzymatic kinetics, they usually serve as some quasi-steady state approximations, which might fail under low concentrations of the reagents. This could be one reason of negative values.

As a putty, you can use options = odeset('NonNegative',1);. Which might give you still corrupted results.

But first check the equations, and start first with plain mass-action, and grow your system gradually. Then you'll see after which addition it fails. This is the main merit of minimal working example concept.