Use of ODE45 for concentration plot help
11 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
0 commentaires
Réponses (1)
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
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!