Trying to plot variables as a function of time with a coupled ode

I have to use 5 ODEs (Ca, Cb, Cs, P, T) all coupled togeher to Plot the reactor temperature, the head space pressure, and the concentrations of A, B, S, D as function of time. I am not familiar with odes and the recommended ode to use is ode15s because of stiffness. There are two reactions for this process: A + B → C + ½ D (gas) S → 3 D (gas) + misc liquids and solids

I am trying to use these equations to make this coupled ode

r1𝐴 = −𝑘𝐴𝐶𝐴𝐶B

𝑟2𝑆 = −𝑘𝑆𝐶s

Assuming that the volume of the batch does not change, the mole balances for the reaction liquid yield: 𝑑𝐶𝐴/𝑑𝑡 = 𝑟1𝐴 = −𝑘𝐴𝐶𝐴𝐶𝐵 , d𝐶𝐵/𝑑𝑡 = 𝑟1𝐵 = 𝑟1𝐴 = −𝑘𝐴𝐶𝐴𝐶𝐵 , 𝑑𝐶𝑆/𝑑𝑡 = 𝑟2𝑆 = −𝑘𝑆𝐶S

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

I believe to calculate temperature as a function of time I can use one of these equations:

k(𝑇) = 𝑘0𝑒𝑥𝑝 (− 𝐸𝑎 /RT)

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

I am struggling to be able to couple all these odes together to be able to calculate all these different variables as a function of time. I would be very appreciative if someone could help me understand how to do this and help with the code. If to help you need any extra data just ask and i'll give it if available.

Thanks

Answer 1

[Abhishek Chakram]

Hi Tom Goodland,

It appears to me that you are facing difficulty in plotting variables as a function of time with a coupled ode. One of the ways of achieving this will be using "ode15s" and passing a custom function to it. Here is a sample code for the same:

y0 = [10; 10; 10; 1; 300]; % Initial values of CA, CB, CS, P, T
tspan = [0 100]; % Time span from 0 to 100 seconds

options = odeset('RelTol',1e-6,'AbsTol',1e-9); % Set the solver options
[t,y] = ode15s(@myODEs,tspan,y0,options); % Solve the ODEs
plot(t,y(:,1),'r',t,y(:,2),'g',t,y(:,3),'b') % Plot CA, CB, and CS vs time
xlabel('Time (s)')
ylabel('Concentration (mol/L)')
legend('CA','CB','CS')

figure % Create a new figure
plot(t,y(:,4),'k') % Plot P vs time
xlabel('Time (s)')
ylabel('Pressure (Pa)')

figure % Create a new figure
plot(t,y(:,5),'m') % Plot T vs time
xlabel('Time (s)')
ylabel('Temperature (K)')

function dydt = myODEs(t,y)
% Define the parameters and constants
kA = 0.005;
kS = 0.1;
S = 100;
P = 0;
E = 10;
C = 0;
R = 8.314;
V0 = 1;
VH = 1; % Assume some value for the head space volume
FD = 0.01; % Assume some value for the molar flow rate of D
Fvent = 0.02; % Assume some value for the venting flow rate
U = 0.5; % Assume some value for the overall heat transfer coefficient
A = 1; % Assume some value for the heat transfer area
Ta = 300; % Assume some value for the ambient temperature
Cp = 1; % Assume some value for the heat capacity of the liquid
dHrxn1 = -100; % Assume some value for the enthalpy change of reaction 1
dHrxn2 = -200; % Assume some value for the enthalpy change of reaction 2
k0 = 1; % Assume some value for the pre-exponential factor
Ea = 10000; % Assume some value for the activation energy

% Define the reaction rates
r1A = -kA*y(1)*y(2); % y(1) is CA, y(2) is CB
r2S = -kS*y(3); % y(3) is CS

% Define the ODEs
dydt = zeros(5,1); % Initialize the output vector
dydt(1) = r1A; % dCA/dt
dydt(2) = r1A; % dCB/dt
dydt(3) = r2S; % dCS/dt
dydt(4) = (FD - Fvent)*R*y(5)/VH; % dP/dt, y(5) is T
dydt(5) = (V0*(r1A*dHrxn1 + r2S*dHrxn2) - U*A*(y(5) - Ta))/(sum(y(1:3))*Cp); % dT/dt
end

Kindly refer to the following documentation to know more about the functions used:

• ode15s: https://www.mathworks.com/help/matlab/ref/ode15s.html
• function: https://www.mathworks.com/help/matlab/ref/function.html

Best Regards,

Abhishek Chakram