# Use of ODE45 for concentration plot help

조회 수: 48 (최근 30일)
Daniel Henry 2023년 2월 14일
댓글: Daniel Henry 2023년 3월 21일
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

댓글을 달려면 로그인하십시오.

### 답변 (1개)

Bjorn Gustavsson 2023년 2월 14일
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이전 댓글 3개 표시이전 댓글 3개 숨기기
Bjorn Gustavsson 2023년 3월 21일
@Daniel Henry, did our suggestions solve your problem?
Daniel Henry 2023년 3월 21일
Hi there! Yes with a bit of rejigging I got the code to work. Thanks very much for your help.

댓글을 달려면 로그인하십시오.

### 카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

R2021b

### Community Treasure Hunt

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

Start Hunting!

Translated by