Solving ODE 45 of chemical reaction mole balance dF/dV

Im trying to solve an ODE of dF/dV derived from a chemical reaction as shown below but the figure plotted does not seem to fit the answer of change in flow rate over the volume of the reactor as expected

close all
clear 
clc

%Set Conditions
T=318;
P=15;
R=8.205745e-5;
C=P/(R*T);

%Rate Constant Functions
K=(0.076*T^2)/(exp(189*((1/T)-(1/90))));
k_1f=0.06*(exp(40*((1/502)-(1/T))));
k_2=0.00202*(exp(4030*((1/441)-(1/T))));
k_1b=k_1f/K;

%Set Parameter
para.k1f=k_1f;
para.k1b=k_1b;
para.k2=k_2;

para.T=T;
para.P=P;
para.R=R;
para.C=C;

%Initial Condition
F0=[1.5 3 0 0 2.1];
FT=sum(F0);
para.FT=FT;

%Volume of Reaction
Vspan=[0 10];

%Solve ODE
options = odeset('NonNegative',1:length(F0));
[V,F] = ode45(@myOde,Vspan,F0,options,para);

% Plot results
figure
plot(V,F,'LineWidth',1.5)
xlabel('Volume')
ylabel('Molar Flow Rate')
legend('C2H4','H2O','C2H5OH','C2H52O','N2')

function dFdV = myOde(V,F,para)

% Initialization
n=length(F);
dFdV=zeros(n,1); % initialize dydt as a column vector

% Assign state vectors to each column
C2H4=F(1);
H2O=F(2);
C2H5OH=F(3);
C2H52O=F(4);
N2=F(5);

% Parameters
k1f=para.k1f;
k1b=para.k1b;
k2=para.k2;

t=para.T;
p=para.P;
r=para.R;
c=para.C;

ft=para.FT;

% ODE-related expressions
r1f=(k1f*c^2*C2H4*H2O)/ft^2;
r1b=(k1b*c*C2H5OH)/ft;
r2=(k2*c*C2H5OH)/ft;

% ODEs
dFdV = zeros(5,1);
    dFdV(1)=-r1f+r1b;
    dFdV(2)=-r1f+r1b+((r2)/2);
    dFdV(3)=r1f-r1b-r2;
    dFdV(4)=((r2)/2);
    dFdV(5)=0;
end

[Sam Chak]

Hi @Proone

What are the expected solutions for the five changes in flow rates over the volume of the reactor? Please show the closed-form solutions using equations or formulas so that we can test them for validating numerical solutions.

Answer 1

[Alan Stevens]

You had the units of R in m^3.atm/(mol.K), but should probably be in L.atm/(mol.K) , since other volumes are in L.

Does this look more like what you expect to see?

%Set Conditions
T=318;                                      % K
P=15;                                       % atm
R=8.205745e-2;                              % L.atm/(mol.K)
C=P/(R*T);                                  % mol/L

%Rate Constant Functions
K=(0.076*T^2)/(exp(189*((1/T)-(1/90))));    % L/mol
k_1f=0.06*(exp(40*((1/502)-(1/T))));        % L/(mol.s)
k_2=0.00202*(exp(4030*((1/441)-(1/T))));    % 1/s
k_1b=k_1f/K;                                % 1/s

%Set Parameter
para.k1f=k_1f;                              % L/(mol.s) 
para.k1b=k_1b;                              % 1/s
para.k2=k_2;                                % 1/s

para.T=T;                                   % K
para.P=P;                                   % atm
para.R=R;                                   % L.atm/(mol.K)
para.C=C;                                   % mol/L

%Initial Condition
F0=[1.5 3 0 0 2.1];                         % mol/s
FT=sum(F0);                                 % mol/s
para.FT=FT;                                 % mol/s

%Volume of Reaction
Vspan=[0 10000];                               % L

%Solve ODE
options = odeset('NonNegative',1:length(F0));
[V,F] = ode45(@myOde,Vspan,F0,options,para);

C2H4=F(:,1);
H2O=F(:,2);
C2H5OH=F(:,3);
C2H52O=F(:,4);
N2=F(:,5);

% Plot results
figure
subplot(2,2,1)
plot(V,C2H4,'LineWidth',1.5)
xlabel('Volume [L]')
ylabel('C2H4 Molar Flow Rate [mol/s]')
subplot(2,2,2)
plot(V,H2O,'LineWidth',1.5)
xlabel('Volume [L]')
ylabel('H2O Molar Flow Rate [mol/s]')
subplot(2,2,3)
plot(V,C2H5OH,'LineWidth',1.5)
xlabel('Volume [L]')
ylabel('C2H5OH Molar Flow Rate [mol/s]')
subplot(2,2,4)
plot(V,C2H52O,'LineWidth',1.5)
xlabel('Volume [L]')
ylabel('C2H52O Molar Flow Rate [mol/s]')

%legend('C2H4','H2O','C2H5OH','C2H52O','N2')

function dFdV = myOde(~,F,para)

% Initialization
% n=length(F);
% dFdV=zeros(n,1); % initialize dydt as a column vector

% Assign state vectors to each column         mol/s
C2H4=F(1);
H2O=F(2);
C2H5OH=F(3);
C2H52O=F(4);
N2=F(5);

% Parameters
k1f=para.k1f;                                   % L/(mol.s)
k1b=para.k1b;                                   % 1/s
k2=para.k2;                                     % 1/s

t=para.T;                                       % K
p=para.P;                                       % atm
r=para.R;                                       % L.atm/(mol.K)
c=para.C;                                       % mol/L

ft=para.FT;                                     % mol/s

% ODE-related expressions                       [c/ft] = s/L
r1f=(k1f*c^2*C2H4*H2O)/ft^2;                    % (mol/s)/L
r1b=(k1b*c*C2H5OH)/ft;                          % (mol/s)/L                       
r2=(k2*c*C2H5OH)/ft;                            % (mol/s)/L

% ODEs
dFdV = zeros(5,1);
    dFdV(1)=-r1f+r1b;                           % (mol/s)/L
    dFdV(2)=-r1f+r1b+((r2)/2);                  % (mol/s)/L
    dFdV(3)=r1f-r1b-r2;                         % (mol/s)/L
    dFdV(4)=((r2)/2);                           % (m0l/s)/L
    dFdV(5)=0;                                  % (mol/s)/L
end