Help with odes and two different reactions

𝑑𝐶𝐴/𝑑𝑡 = 𝑟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

[Torsten]

Look at the second example (Solve Stiff ODE) to see how you can solve all three ODEs in one call to the integrator.

https://de.mathworks.com/help/matlab/ref/ode15s.html

Answer 1

[HWIK]

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.

[Tom Goodland]

Thanks, for the question I have to plot A,B,S,D concentrations, aswell as temperature and head space pressure as a function of time (at the moment i'm focusing on the concentrations). What would I have to alter in the code you supplied to be able to plot cA etc. as a function of time?

[HWIK]

What I wrote is just the function, to solve it you just define what time range you want to plot over and define your initial conditions, which you then feed to the ode solver:

tspan = [0 200]; %just as an example
initial_conditions = [3 2 1]; %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):
[t , y] = ode45(@name_your_function, tspan initial_conditions);
%output t is your time range, and y an array with all your variables where:
ca = y(:,1);
cb = y(:,2); %and so on, you can then plot(t, ca)

If you want to solve for pressure and temperature you will have to add those in the code of the function I added before in a similar way as I did with the rate laws

[Torsten]

tspan = [0 200]; %just as an example
initial_conditions = [4.3 5.1 3]; %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):
[t,y] = ode15s(@(t,y)name_your_function(t,y(1),y(2),y(3)), tspan, initial_conditions);
%output t is your time range, and y an array with all your variables where:
ca = y(:,1);
cb = y(:,2);
cs = y(:,3);
plot(t,[ca,cb,cs]);

function out = name_your_function(tspan,ca,cb,cs);
%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;
%choose your values of k1 & k2:
k1 = k_0_1*exp(-Ea1/(R_gas*422));
k2 = k_0_2*exp(-Ea2/(R_gas*422));
%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

[HWIK]

When assigning variables in your ode function dont assign them to a constant, use the format I followed in the code I gave, assign them to the input of the function, if not the ode solver will be solving for the same thing all the time.

When you are calling the function in the solver use @ before

[Tom Goodland]

thorsten I ran your code and thanks it gave a good curve for A and B, however for S it is a straight horizontal line, I know S shouldn't be like this due to it turning into D (reaction 2). Do you know why it turns up on the graph without any change in conc. over time? Is there something extra I have to specify?

Thanks

[Torsten]

The reaction rate r2 is only about 2.82e-15 mol/(l*s). Thus a decrease in concentration for substance S is not visible in the graphics.

Maybe your reaction constants are wrong.

[Tom Goodland]

I have to plot the concentration of D as a function of time aswell, as cD doesn't show up in any of the rate constants, to be able to plot this do I have to code cD/dt in terms of the other components considering the mole ratios?

[HWIK]

As long as what you put in your question are elementary reactions, dcD/dt = r1/2 + 3*r2

[Tom Goodland]

Thanks I only need to plot temperature and head space pressure over time now. Below is the derivative for head space pressure:

dp /d𝑡 = (𝐹𝐷 − 𝐹𝑣𝑒𝑛𝑡) 𝑅𝑇/𝑉H

𝐹𝐷 = (−0.5𝑟1𝐴 − 3𝑟2𝑆 )V0

𝐹𝑣𝑒𝑛𝑡 = 𝐹𝐷 when FD < 11,400 mol/h

𝐹𝑣𝑒𝑛𝑡 = ∆𝑃𝐶𝑉1 = (𝑃 − 1 atm)𝐶𝑣1 when P < 28.2 atm

𝐹𝑣𝑒𝑛𝑡 = ∆𝑃(𝐶𝑣1 + 𝐶𝑣2 ) = (𝑃 − 1 atm)(𝐶𝑣1 + 𝐶𝑣2 ) when P > 28.2 atm

Below is the derivative for temperature:

d𝑇 /d𝑡 = (V0(𝑟1𝐴∆𝐻𝑟𝑥𝑛,1 + 𝑟2𝑆∆𝐻𝑟𝑥𝑛,2) − 𝑈𝐴(𝑇 − 𝑇𝑎 ) )/∑ 𝑁𝑖𝐶𝑃i

Do you know how I could code these derivatives into odes and then plot the graphs? Some of these values in the equation are known, if you need to know what values are known let me know. I would appreciate any help.

Thanks

[HWIK]

follow the same structure as the current one, just add temperature and pressure to the input and output variables, and define their relationships. For the conditionals in the pressure use an if statement.

[Tom Goodland]

should I add two more ys to the same ode 15 or make a new ode15 for these two?

[HWIK]

Just add the expressions for dP and dT in the function you are adding to the ode solver, you only need one function, and then just update the output: out = [dca ; dcb;... ; dP; dT].

Similarly add the inputs in the same order: P = inputs(5); T = inputs(6)

[Torsten]

If it worked with three, it will also work with five :-)

[Tom Goodland]

I've tried to add the head space pressure and temperature to the ode so I've added stuff to both the function and script. however I've been getting errors:

tspan = [0 200]; %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):

[t,y] = ode15s(@(t,y)thorsten(t,y(1),y(2),y(3),y(4),y(5),y(6)), tspan, initial_conditions);

%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);

plot(t,[ca,cb,cs,cd,P,T_2]);

Here's the script:

function out = thorsten(tspan,ca,cb,cs,cd,P,T_2)

%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;

T_2 = 422; % K

Ta = 373;

UA = 0;

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;

% Fvent

Fvent = Fd;

if Fd < 11,400; % mol / h

Fvent = (P - 1)*Cv1;

if P < 28.2 % atm

end

Fvent = (P - 1)*(Cv1 + Cv2);

if P > 28.2

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

I've named the temperature T_2 because I've used T previously in different work. dT is representing the change in T_2.

The errors im getting are:

Index exceeds the number of array elements. Index must not exceed 1.

Error in thorsten (line 43)

dT = (V0*(r1*Hr1 + r2*Hr2) - UA(T_2 - Ta))/NCp;

Error in untitled2>@(t,y)thorsten(t,y(1),y(2),y(3),y(4),y(5),y(6)) (line 5)

[t,y] = ode15s(@(t,y)thorsten(t,y(1),y(2),y(3),y(4),y(5),y(6)), tspan, initial_conditions);

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode15s (line 152)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in untitled2 (line 5)

[t,y] = ode15s(@(t,y)thorsten(t,y(1),y(2),y(3),y(4),y(5),y(6)), tspan, initial_conditions);

Do you know why im getting these errors, any help would be appreciated.

Thanks

[Torsten]

function main
tspan = [0 200]; %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):
[t,y] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6)), tspan, initial_conditions);
%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 out = fun(tspan,ca,cb,cs,cd,P,T_2)
%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;
UA = 0;
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

Check whether Fvent is calculated correctly.

[HWIK]

Use semicolons in your out variable in the function, it must be a column vector.

In the solver you have one input for your initial conditions, therefore your thorsten function must too, and extract those variables at the begining of the function from the array you input it (just like I put in the initial answer...).

You are missing a * in the dT variable after the UA.

Also, do not define the variables you are trying to solve as constants, like you did for T_2, or else it will remain as a constant.

[Tom Goodland]

Thorsten is this what you mean by 'Get your inputs and assign them to variables':

function out = func(tspan,ca,cb,cs,cd,P,T_2)

%Get your inputs and assign them to variables:

ca = r1/-k1*cb;

cb = r1/-k1*ca;

cs = r2/-k2;

[Tom Goodland]

I'm still getting a bunch of errors when I run that code. Do you know why this is?

Thanks

[Torsten]

I'm still getting a bunch of errors when I run that code. Do you know why this is?

No, on my PC the code runs without error. There is nothing to be changed in the code.

[Tom Goodland]

I'm meant to be plotting these variables over time for a reactor explosion, plotting until the reactor explodes when pressure reaches 4.4 atm or temperature reaches 600 K. Do you know why the temperature graph decreases, could it be due to the very small r2s value. When the srcipt is ran does the r2s value change or stay constant?

Thanks

[Torsten]

Ua = 0, Hr1, Hr2 are both negative, r1 and r2 are positive, so dT/dt < 0.

r2s is positive, but becomes smaller with time: from 2.8e-15 to 5.3e-27 (so all appr. 0).

But you should learn to debug by yourself. Just remove the semicolon after the variable you want to see, run the code and look at the values printed to the command window.

[Tom Goodland]

For another part I have to do these plots all again for: where the cooling water is turned on (UA = 2.77E6 J/(h K)) once the temperature in the reactor reaches 455 K. How would I calculate the initial conditions for this, is there something on matlab to choose 455K and find the CA, CB, CS, and CD values corresponding on the graphs?

Thanks

[Tom Goodland]

for the first part you've helped me with matlab says

Warning: Failure at t=3.426260e+00. Unable to meet integration tolerances without reducing the step size below the

smallest value allowed (7.105427e-15) at time t.

> In ode15s (line 655)

In Twob (line 6)

I don't see this as much of a problem as the temperature has already shot up past 500K (reactor exploded) by then, however when I run the code the pressure time plot is a straight horizontal line and it doesn't shoot up or increase when the temperature increases.

dp = (Fd - Fvent)*((R_gas*T_2)/VH) in the code so would I be right in saying that as temperature increases the pressure should aswell due to T_2 being multuplied by R_gas in the equation.

Thanks

[Torsten]

For another part I have to do these plots all again for: where the cooling water is turned on (UA = 2.77E6 J/(h K)) once the temperature in the reactor reaches 455 K. How would I calculate the initial conditions for this, is there something on matlab to choose 455K and find the CA, CB, CS, and CD values corresponding on the graphs?

Define an event for the integrator if T reaches 455 K.

https://de.mathworks.com/help/matlab/math/ode-event-location.html

dp = (Fd - Fvent)*((R_gas*T_2)/VH) in the code so would I be right in saying that as temperature increases the pressure should aswell due to T_2 being multuplied by R_gas in the equation.

Only if Fd > Fvent

[Tom Goodland]

[t,y,te,ye,ie] = odeXY(odefun,tspan,y0,options)

The three additional outputs returned by the solver correspond to the detected events:

• te is a column vector of the times at which events occurred.
• ye contains the solution value at each of the event times in te.
• ie contains indices into the vector returned by the event function. The values indicate which event the solver detected.

What would ye be in this case would it be the new equation for dT/dt with the new UA value?

[Torsten]

No, it's the solution vector [ca,cb,cs,...] at the time when the event occured.

You will have to call odeXY with these values as start values for a new period of integration and with the equation for dT/dt changed. Hand a flag to odefun which contains the number of the integration period and which serves you to set up the dT/dt equation correctly:

function dydt = odefun(t,ca,cb,cs,cd,P,T_2,flag)
...
    if flag==1
      dT = ...;
    elseif flag ==2
      dT = ...;
    end

[Tom Goodland]

what are the flag 1 and 2 parts, if flag = 1 is the UA = 0 part, and if flag = 2 UA = 2.77E6?

[Tom Goodland]

Also what code could I use to find out the time when the temperature reaches 455K to be able to find the [ca,cb,...] new initial conditions?

[Torsten]

This is a rough guideline.

If errors pop up, first make an attempt to correct them on your own, please.

function main
  tspan = [0 200]; %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
  options = odeset('Events',@myEventsFcn);
  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;
  [t2,y2] = ode15s(@(t,y)fun(t,y(1),y(2),y(3),y(4),y(5),y(6),iflag), [te tspan(end)], ye);
  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)
  value = y(6) - 455; % 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

[Tom Goodland]

Thanks so much it worked nicely. I now have to find out the maximum time the cooling water can be turned off before being reinstated to avoid the reactor exploding (happens at T = 500 K or P = 45 atm). I want to use the backbone of the code and have a function keep on picking temperatures to turn on the cooling water above 455 K. The objective for me is to find out how late I can turn the cooling water on and the reactor not explode. Do you know what code or function I could use that runs the ode with different event details for when the UA changes and gives me a graph for each different temperature when the cooling water is turned on?

I'd appreciate any help

[Tom Goodland]

This is the temp graph for when the cooling water is turned on at 455 K, I need to find out the highest temperature you can turn the water on (or time if that makes it easier for matlab) to left of the peak to avoid explosion.

[Tom Goodland]

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]

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