Use of ODE45 for concentration plot help

Currently writing a code for the modelling of a chemical reactor for teh reaction of NH3+HCl -> NH4Cl by solving a first order ODE to gain a concentration plot.

I have tried using ODE45 but I cannot get the code to work to consider all three components in the system.

My boundary conditions are initial reactant concentrations are both 0.0109 mol/m3 and my constant k = 2.02*10^10 m3/s. Final reactant concentration = 0 assuming 100% converion.

The kinetic equation is just -ra = k*Cnh3*Chcl

Any help is appreciated, my code is below

%% Matlab code for Irreversible first order ammonium chloride synthesis reaction
%% Using ode45

function xdot = f(t,x)
% Rate constants 
k1 = 2.02;
% Reactants 
cA = x(1);
cB = x(2);
% Reaction rates 
r1 = k1*cA*cB;

% Differential equations 
xdot(1) = -r1;
xdot(2) = r1;
xdot = [xdot(1); xdot(2)];

%% Initial conditions 
cA0 = 0.0109;
cB0 = 0.0109;
x0 = [cA0; cB0];

%% Time span
tspan = [0 6];

%% Solving ODE
[t,x] = ode45(@f,tspan,x0);

%% Plotting
plot(t,x(:,1),'r-',t,x(:,2),'b-');
title('Ammonium Chloride Synthesis (Irreversible)');
xlabel('time');

end

Answer 1

[Bjorn Gustavsson]

To the best of my understanding you should include all three species in the ode-function. Perhaps something like this:

function xdot = f(t,x)
% Rate constants 
k1 = 2.02;
% Reactants 
cA = x(1);
cB = x(2);
cC = x(3);
% Reaction rates 
r1 = k1*cA*cB;
% Differential equations 
xdot = zeros(3,1);
xdot(1) = -r1;
xdot(2) = -r1;
xdot(3) = r1;

Then you specify the initial conditions for all 3 species:

%% Initial conditions 
cA0 = 0.0109;
cB0 = 0.0109;
cC0 = 0;
x0 = [cA0; cB0; cC0];
%% Time span
tspan = [0 160]; % It seems necessary to increase time-span - slow reaction?
%% Solving ODE
[t,x] = ode45(@f,tspan,x0);
%% Plotting
ph = plot(t,x);
set(ph(1),'color','r');
set(ph(2),'color','b')
set(ph(3),'color','g','linewidth',2)
title('Ammonium Chloride Synthesis (Irreversible)');
xlabel('time');

Then I get curves that seem familiar to me - but I'm not an expert in this type of chemistry so I don't know if reaction-response should be on these time-scales or if they are completely off.

HTH

[Daniel Henry]

Hi, thanks for your reply, I get an error message when I try to run this code? I tried changing the vartiable names but it still doesnt work

[Torsten]

Then take the code @Bjorn Gustavsson supplied. It works.

The function part must follow the script part (see below).

%% Initial conditions 
cA0 = 0.0109;
cB0 = 0.0109;
cC0 = 0;
x0 = [cA0; cB0; cC0];
%% Time span
tspan = [0 160]; % It seems necessary to increase time-span - slow reaction?
%% Solving ODE
[t,x] = ode45(@f,tspan,x0);
%% Plotting
ph = plot(t,x);
set(ph(1),'color','r');
set(ph(2),'color','b')
set(ph(3),'color','g','linewidth',2)
title('Ammonium Chloride Synthesis (Irreversible)');
xlabel('time');

function xdot = f(t,x)
  % Rate constants 
  k1 = 2.02;
  % Reactants 
  cA = x(1);
  cB = x(2);
  cC = x(3);
  % Reaction rates 
  r1 = k1*cA*cB;
  % Differential equations 
  xdot = zeros(3,1);
  xdot(1) = -r1;
  xdot(2) = -r1;
  xdot(3) = r1;
end

[Bjorn Gustavsson]

@Daniel Henry, my habit is to put all types of ODE-functions in separate files (in my relevant directory/toolbox), and then have separate scripts for each work-task. That way I find it easier to use, reuse, and reuse the same functions over and over again. But for that to apply in this specific case you would better use another name for the ode-function than a simple f. Perhaps something like: ode_ammonium_chloride_synth, or ode_NH4Cl_synth - to make it easier to remember and search for.

[Bjorn Gustavsson]

@Daniel Henry, did our suggestions solve your problem?

[Daniel Henry]

Hi there! Yes with a bit of rejigging I got the code to work. Thanks very much for your help.