Effacer les filtres
Effacer les filtres

Use of ODE45 for concentration plot help

20 vues (au cours des 30 derniers jours)
Daniel Henry
Daniel Henry le 14 Fév 2023
Commenté : Daniel Henry le 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

Réponses (1)

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

Connectez-vous pour commenter.

Produits


Version

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by