Chemical reaction Rates with changing temperature

Im attampting to model the kenetics of a chemical reaction calculating K (with respect to temperature), concentration of CB, and recator temp simultaniously, I created a while loop that continious calculates that temperature using Eulers Method.

I keep getting the error:

"Error using odearguments

@(T,CB)COEFF.*CB returns a vector of length 29, but the length of initial conditions vector is 1. The vector returned by @(T,CB)COEFF.*CB and the initial conditions vector must have the same number of elements.

Error in CBR_T_reactor_V2 (line 34)

[t,CB1] = ode23(CB_fun(CC,k1), t, CB0);

The following is my code

clc
close all
clear

% let B = H2O2, R = H2O, T = O2, and C = I- (catalyst)
% overall reaction: 2B+C -> T+2R
% Key reactant is B and Key product is T
% Define reaction kinetics and parameters  
V_rect1= 1; % Reactor Volume in liters - 1st run
V_rect2= 1.5; % Reactor Volume in liters - 2nd run
MmCpm= 0.18*4184; % in j/K  
% Time vector in sec
tspan = (0:5:500); 
tspan = tspan';
dt=5;
Ts= 298.15;
T1=Ts; % Inital value 
T2=Ts; % Inital value
iter=1;
while  iter <= 101
    %Solving for CB - Concentration of CB
    t=[0 1500];
    CB0 = 0.09; 
    CC = 0.24; 
    k1= 5.052e8*exp(-59020./(8.314*T1));
    k2= 5.052e8*exp(-59020./(8.314*T2));
    [t,CB1] = ode23(CB_fun(CC,k1), t, CB0);
    [t,CB2]  = ode23(CB_fun(CC,k2), t, CB0);
    
    % Solve for  V_T - Volumetric O2 rate 
    Ps= 1; % pressure in atm
    V_o2= 22.4; % standard gas Volume in Liters at STP
    R= 0.0821;
    rT1=CC*k1*CB1; % (EQ 5)
    rT2=CC.*k2*CB2; % (EQ 5)
    V_T1=(rT1.*V_o2.*R.*Ts*60)/Ps; %(EQ 9)
    V_T2=(rT2.*V_o2.*R.*Ts*60)/Ps; %(EQ 9)
    
    % Solve for  Temperature  
    [dTdt1] =deltaT_fun(V_rect1, rT1);
    [dTdt2] =deltaT_fun(V_rect2, rT2);
    
    T1=T1+dTdt1.*dt;
    T2=T2+dTdt2.*dt;
    iter= iter+ 1;
end 
% Ploting Data 
figure (1)
plot(time1, T1, 'Red ^'); 
hold on
plot(time2, T2, 'blue ^'); 
yyaxis right
plot(time1, V_T1, 'Red +'); 
plot(time2, V_T2, 'blue +'); 
yyaxis left
title( 'Adiabatic Reactor Temp. and O2 rate vs. Time'); 
xlabel('Time (seconds)'); 
ylabel(' Reactor Temp (Kelvin)');

yyaxis right
ylabel('Volumetric O2 Rate (L/min)')
legend('1 L Reactor Tempure', ...
    '1.5 L Reactor Temperature', ...
    'Predicted Volumetic 02 rate')
%% Concentration of CB Function 
function [dCBdt]= CB_fun(CC,k)
N=-2*CC;
Coeff=N*k;
dCBdt = @(t, CB) Coeff.*CB;% dCB/dt (combination of E5 and E7)
end
%% Delta_Recator_Temp  Function 
function [dTdt] =deltaT_fun(Vol_rector,r_T)
cp_w = 4.184; %Mass-based specific heat of water (j/g*Kelvin)
rho_w = 997; %Water mass density kg/m^3
delta_H_RXN=-196*1000; % heat of RXN for O2 produced in (j/Mol)
dTdt= (Vol_rector.*r_T*-delta_H_RXN)./(Vol_rector*cp_w*rho_w);%dT/dt Reactor Temp rate of change (EQ 10)
end

[Torsten]

Maybe you could explain in your own words what you are doing. It's very hard to understand your code.

Answer 1

[Walter Roberson]

% let B = H2O2, R = H2O, T = O2, and C = I- (catalyst)
% overall reaction: 2B+C -> T+2R
% Key reactant is B and Key product is T
% Define reaction kinetics and parameters  
V_rect1= 1; % Reactor Volume in liters - 1st run
V_rect2= 1.5; % Reactor Volume in liters - 2nd run
MmCpm= 0.18*4184; % in j/K  
% Time vector in sec
tspan = (0:5:500); 
tspan = tspan';
dt=5;
Ts= 298.15;
T1=Ts; % Inital value 
T2=Ts; % Inital value
iter=1;
while  iter <= 101
    %Solving for CB - Concentration of CB
    t=[0 1500];
    CB0 = 0.09; 
    CC = 0.24; 
    k1= 5.052e8*exp(-59020./(8.314*T1));
    k2= 5.052e8*exp(-59020./(8.314*T2));
    [t,CB1] = ode23(CB_fun(CC,k1), t, CB0);
    [t,CB2]  = ode23(CB_fun(CC,k2), t, CB0);
    
    % Solve for  V_T - Volumetric O2 rate 
    Ps= 1; % pressure in atm
    V_o2= 22.4; % standard gas Volume in Liters at STP
    R= 0.0821;
    rT1=CC*k1*CB1; % (EQ 5)
    rT2=CC.*k2*CB2; % (EQ 5)
    V_T1=(rT1.*V_o2.*R.*Ts*60)/Ps; %(EQ 9)
    V_T2=(rT2.*V_o2.*R.*Ts*60)/Ps; %(EQ 9)
    
    % Solve for  Temperature  
    [dTdt1] =deltaT_fun(V_rect1, rT1);
    [dTdt2] =deltaT_fun(V_rect2, rT2);
    
    whos
    T1=T1+dTdt1.*dt;
    T2=T2+dTdt2.*dt;
    iter= iter+ 1;
end 

  Name           Size            Bytes  Class     Attributes

  CB0            1x1                 8  double              
  CB1           29x1               232  double              
  CB2           29x1               232  double              
  CC             1x1                 8  double              
  MmCpm          1x1                 8  double              
  Ps             1x1                 8  double              
  R              1x1                 8  double              
  T1             1x1                 8  double              
  T2             1x1                 8  double              
  Ts             1x1                 8  double              
  V_T1          29x1               232  double              
  V_T2          29x1               232  double              
  V_o2           1x1                 8  double              
  V_rect1        1x1                 8  double              
  V_rect2        1x1                 8  double              
  cmdout         1x33               66  char                
  dTdt1         29x1               232  double              
  dTdt2         29x1               232  double              
  dt             1x1                 8  double              
  iter           1x1                 8  double              
  k1             1x1                 8  double              
  k2             1x1                 8  double              
  rT1           29x1               232  double              
  rT2           29x1               232  double              
  t             29x1               232  double              
  tspan        101x1               808  double              

Error using odearguments
@(T,CB)COEFF.*CB returns a vector of length 29, but the length of initial conditions vector is 1. The vector returned by @(T,CB)COEFF.*CB and the initial conditions vector must have the same number of elements.

Error in ode23 (line 102)
odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

% Ploting Data 
figure (1)
plot(time1, T1, 'Red ^'); 
hold on
plot(time2, T2, 'blue ^'); 
yyaxis right
plot(time1, V_T1, 'Red +'); 
plot(time2, V_T2, 'blue +'); 
yyaxis left
title( 'Adiabatic Reactor Temp. and O2 rate vs. Time'); 
xlabel('Time (seconds)'); 
ylabel(' Reactor Temp (Kelvin)');

yyaxis right
ylabel('Volumetric O2 Rate (L/min)')
legend('1 L Reactor Tempure', ...
    '1.5 L Reactor Temperature', ...
    'Predicted Volumetic 02 rate')
%% Concentration of CB Function 
function [dCBdt]= CB_fun(CC,k)
N=-2*CC;
Coeff=N*k;
dCBdt = @(t, CB) Coeff.*CB;% dCB/dt (combination of E5 and E7)
end
%% Delta_Recator_Temp  Function 
function [dTdt] =deltaT_fun(Vol_rector,r_T)
cp_w = 4.184; %Mass-based specific heat of water (j/g*Kelvin)
rho_w = 997; %Water mass density kg/m^3
delta_H_RXN=-196*1000; % heat of RXN for O2 produced in (j/Mol)
dTdt= (Vol_rector.*r_T*-delta_H_RXN)./(Vol_rector*cp_w*rho_w);%dT/dt Reactor Temp rate of change (EQ 10)
end

[Walter Roberson]

    [t,CB1] = ode23(CB_fun(CC,k1), t, CB0);

The first call to ode23 returns a second output that is number of times by number of states. Your input timespan t was a vector of length 2, so the number of output times will be whatever ode23 feels like.

    [t,CB2]  = ode23(CB_fun(CC,k2), t, CB0); 

Here you are passing in the output times from the first ode23 call to become the tspan for the second call, so the number of output times for the second call will be the same as the number of output times for the first call unless the second call exits early due to non-convergence.

    rT1=CC*k1*CB1; % (EQ 5)
    rT2=CC.*k2*CB2; % (EQ 5)

Vectors times something gets you vector rT1 and rT2 outputs.

    [dTdt1] =deltaT_fun(V_rect1, rT1);
    [dTdt2] =deltaT_fun(V_rect2, rT2);

That are passed to delta_fun

    dTdt= (Vol_rector.*r_T*-delta_H_RXN)./(Vol_rector*cp_w*rho_w);%dT/dt Reactor Temp rate of change (EQ 10)

resulting in an output the same size as the vector.

    T1=T1+dTdt1.*dt;
    T2=T2+dTdt2.*dt;

Those vectors dTdt1 and dTdt2 combine with the previously-scalar T1 and T2 giving vector outputs for T1 and T2

Then, next iteration of the loop

    k1= 5.052e8*exp(-59020./(8.314*T1));
    k2= 5.052e8*exp(-59020./(8.314*T2));

Your k1 and k2 become vectors because T1 and T2 are vectors.

    [t,CB1] = ode23(CB_fun(CC,k1), t, CB0);
    [t,CB2]  = ode23(CB_fun(CC,k2), t, CB0);

and the vector k1 and k2 get passed to CB_fun,

function [dCBdt]= CB_fun(CC,k)
N=-2*CC;
Coeff=N*k;
dCBdt = @(t, CB) Coeff.*CB;% dCB/dt (combination of E5 and E7)

where they become a vector of coefficients, leading to an objective fucntion that takes in the scalar state Cb and returns a vector output.

But when you have a scalar state variable for ode*() you need to return a scalar output state, not a vector output state.

[Walter Roberson]

I do not think your code can be repaired. It looks to me as if you are trying to implement a PDE.

The only hope would be if you can replace

    rT1=CC*k1*CB1; % (EQ 5)
    rT2=CC.*k2*CB2; % (EQ 5)

with

    rT1=CC*k1*CB1(end); % (EQ 5)
    rT2=CC.*k2*CB2(end); % (EQ 5)