% Integrate the ODE
% -----------------
% The independent variable always requires an initial and final value:
indep_start = 0.0; % m^3
indep_final = 5.0; % m^3
% Set initial condition(s): for integrating variables (dependent variables)
X_depnt_zero = 0.0; % i.e. X(V=0) = 0.0
% Other parameters. The "param" is just a variable in MATLAB.
% It is called a structured variable, or just a "struct"
% It can have an arbitrary number of sub-variables attached to it.
% In this example we have two of them.
param.T_0 = 330; % feed temperature [K]
param.FT0 = 163000/3600; % mol/s (was kmol/hour originally)
% Integrate the ODE(s):
[indep, depnt] = ode45(@pfr_example, [indep_start, indep_final], [X_depnt_zero], odeset(), param);
% Deal with the integrated output to show interesting plots
X = depnt(:,1);
T = 43.3.*X + param.T_0; % what was the temperature profile?
rA_over_FA0 = zeros(numel(X), 1); % what was the rate profile?
C_A = zeros(numel(X), 1); % what was the concentration profile?
for i = 1:numel(X)
[rA_over_FA0(i), C_A(i)] = pfr_example([], X(i), param);
end % the above can be done more efficiently in a single line, but the code would be too confusing
% Plot the results
f=figure;
set(f, 'Color', [1,1,1])
subplot(2, 2, 1)
plot(indep, X); grid
xlabel('Volume, V [kg]', 'FontWeight', 'bold')
ylabel('Conversion, X [-]', 'FontWeight', 'bold')
subplot(2, 2, 2)
plot(indep, T); grid
xlabel('Volume, V [kg]', 'FontWeight', 'bold')
ylabel('Temperature profile [K]', 'FontWeight', 'bold')
subplot(2, 2, 3)
plot(indep, C_A); grid
xlabel('Volume, V [kg]', 'FontWeight', 'bold')
ylabel('Concentration C_A profile [K]', 'FontWeight', 'bold')
subplot(2, 2, 4)
plot(indep, rA_over_FA0); grid
xlabel('Volume, V [kg]', 'FontWeight', 'bold')
ylabel('(Reaction rate/FA0) profile [1/m^3]', 'FontWeight', 'bold')
% Now plot one of the most important figures we saw earlier in the course:
% F_A0 / (-rA) on the y-axis, against conversion X on the x-axis. This plot
% is used to size various reactors.
% The material leaves the reactor at equilibrium; let's not plot
% that far out, because it distorts the scale. So plot to 95% of
% equilibrium
f = figure; set(f, 'Color', [1,1,1])
index = find(X>0.95 * max(X), 1);
plot(X(1:index), 1./rA_over_FA0(1:index)); grid
xlabel('Conversion, X [-]', 'FontWeight', 'bold')
ylabel('FA0/(-r_A) profile [m^3]', 'FontWeight', 'bold')
% Updated code to find the optimum inlet temperature:
temperatures = 300:3:400;
conv_at_exit = zeros(size(temperatures));
i = 1;
for T =
param.T_0 = T; % feed temperature [K]
[indep, depnt] = ode45(@pfr_example, [indep_start, indep_final], ...
[X_depnt_zero], odeset(), param);
conv_at_exit(i) = depnt(end, 1);
i = i + 1;
end
plot(temperatures, conv_at_exit); grid;
xlabel('Inlet temperature')
ylabel('Conversion at reactor exit')