Help with odes and two different reactions

7 views (last 30 days)
Tom Goodland
Tom Goodland on 14 Jan 2022
Commented: Torsten on 18 Jan 2022
๐‘‘๐ถ๐ด/๐‘‘๐‘ก = ๐‘Ÿ1๐ด = โˆ’๐‘˜๐ด๐ถ๐ด๐ถ๐ต ( 9 ) The first two have the same rate constant
๐‘‘๐ถ๐ต/๐‘‘๐‘ก = ๐‘Ÿ1๐ต = ๐‘Ÿ1๐ด = โˆ’๐‘˜๐ด๐ถ๐ด๐ถ๐ต ( 10 )
d๐ถ๐‘†/๐‘‘๐‘ก = ๐‘Ÿ2๐‘† = โˆ’๐‘˜๐‘†๐ถ๐‘† ( 11 )
A + B โ†’ C + ยฝ D (gas)
S โ†’ 3 D (gas) + misc liquids and solids
These are the two reactions, do I have to do two seperate odes because of the different reactions (and rate constants)? Or is there a way to work out all the differentials in the same ode with initial concentrations known?
If any of my code is needed to answer the question let me know.
Thanks
  1 Comment
Torsten
Torsten on 14 Jan 2022
Look at the second example (Solve Stiff ODE) to see how you can solve all three ODEs in one call to the integrator.

Sign in to comment.

Answers (1)

HWIK
HWIK on 14 Jan 2022
Just write a function for the ODEs and feed it to the solver:
function out = name_your_function(tspan, inputs)
%Get your inputs and assign them to variables:
%(keep the same order as in your initial conditions!)
ca = inputs(1);
cb = inputs(2);
cs = inputs(3);
%Define your relationships:
%choose your values of k1 & k2:
k1 = ;
k2 = ;
%Reaction rates:
r1 = k1*ca*cb;
r2 = k2*cs;
%Define your ODEs:
dca = -r1;
dcb = -r1;
dcs = -r2;
%Assign your ourput:
out = [dca; dcb; dcs];
end
Solve over a given time span and for whatever initial conditions you have, using whatever solver seems to suit your system best.
  33 Comments
Tom Goodland
Tom Goodland on 17 Jan 2022
I am thinking some type of loop would be best for this challenge but I'm not sure what, does anyone know what loop code would be best for the problem I described in the 2 previous comments on this thread?
Thanks
Torsten
Torsten on 18 Jan 2022
You mean something like this ?
function main
tspan = [0 10]; %just as an example
initial_conditions = [4.3 5.1 3 0 4.4 422]; %its the values of the variables you are solving for at t=0
%there must be as many inputs to the initial conditions as variables solved
%Solve your ODE function (I picked ode45 in this case,could be another):
Tstop_left = 450.0;
Tstop_right = 495.0;
Tstop = (Tstop_left + Tstop_right)/2.0;
while Tstop_right - Tstop_left >= 0.1
options = odeset('RelTol',1e-8,'AbsTol',1e-8,'Events',@(t,y)myEventsFcn(t,y,Tstop));
iflag = 1;
[t1,y1,te,ye,ie] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),iflag), tspan, initial_conditions,options);
iflag = 2;
options = odeset('RelTol',1e-8,'AbsTol',1e-8,'InitialStep',1e-8);
[t2,y2] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),iflag), [te tspan(end)], ye, options);
Tmax = max(y2(:,6))
if Tmax > 500
Tstop_left = Tstop_left;
Tstop_right = Tstop;
Tstop = (Tstop_left + Tstop_right)/2;
else
Tstop_left = Tstop;
Tstop_right = Tstop_right;
Tstop = (Tstop_left + Tstop_right)/2;
end
end
Tstop
t = [t1;te;t2];
y = [y1;ye;y2];
%output t is your time range, and y an array with all your variables where:
ca = y(:,1);
cb = y(:,2);
cs = y(:,3);
cd = y(:,4);
P = y(:,5);
T_2 = y(:,6);
figure(1)
plot(t,[ca,cb,cs,cd])
figure(2)
plot(t,P)
figure(3)
plot(t,T_2)
end
function [value,isterminal,direction] = myEventsFcn(t,y,Tstop)
value = y(6) - Tstop; % The value that we want to be zero
isterminal = 1; % Halt integration
direction = 0; % The zero can be approached from either direction
end
function out = fun(tspan,ca,cb,cs,cd,P,T_2,iflag)
%Get your inputs and assign them to variables:
%(keep the same order as in your initial conditions!)
%Define your relationships:
k_0_1 = 4E14;
k_0_2 = 1E84;
Ea1 = 128000;
Ea2 = 800000;
R_gas = 8.314;
V0 = 4000;% L
VH = 5000;
Cv1 = 3360; % mol / h * atm
Cv2 = 53600; % mol / h * atm
Hr1 = -45400; % J / mol
Hr2 = -3.2E5;
Ta = 373;
if iflag == 1
UA = 0;
elseif iflag == 2
UA = 2.77E6;
end
NCp = 1.26E7; % J / K
%choose your values of k1 & k2:
k1 = k_0_1*exp(-Ea1/(R_gas*T_2));
k2 = k_0_2*exp(-Ea2/(R_gas*T_2));
%Reaction rates:
r1 = k1*ca*cb;
r2 = k2*cs;
%Fd
Fd = (-0.5*r1-3*r2)*V0;
%
if Fd < 11.4 % mol / h
Fvent = Fd;
else
if P < 28.2 % atm
Fvent = (P - 1)*Cv1;
else
Fvent = (P - 1)*(Cv1 + Cv2);
end
end
%Define your ODEs:
dca = -r1;
dcb = -r1;
dcs = -r2;
dcd = r1/2 + 3*r2;
dp = (Fd - Fvent)*((R_gas*T_2)/VH);
dT = (V0*(r1*-Hr1 + r2*-Hr2) - UA*(T_2 - Ta))/NCp;
%Assign your ourput:
out = [dca; dcb; dcs; dcd;dp;dT];
end

Sign in to comment.

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Tags

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!