Azzera filtri
Azzera filtri

Use of ODE45 for concentration plot help

27 visualizzazioni (ultimi 30 giorni)
Daniel Henry
Daniel Henry il 14 Feb 2023
Commentato: Daniel Henry il 21 Mar 2023
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

Risposte (1)

Bjorn Gustavsson
Bjorn Gustavsson il 14 Feb 2023
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
  5 Commenti
Bjorn Gustavsson
Bjorn Gustavsson il 21 Mar 2023
@Daniel Henry, did our suggestions solve your problem?
Daniel Henry
Daniel Henry il 21 Mar 2023
Hi there! Yes with a bit of rejigging I got the code to work. Thanks very much for your help.

Accedi per commentare.

Prodotti


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by