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
