% Catalytic Reactor Design
reg_no = 409083;
dob_day = 29 ;
% Extract values
n = mod(reg_no, 100); % Last two digits of registration number
a = dob_day;
b = a - 1;
% Determine rounding based on even/odd reg_no
if mod(reg_no, 2) == 0
round_figs = 4;
else
round_figs = 5;
end
% GIVEN CONSTANTS
Fao = 5.0; % mol/min
To = 450; % K
Po = 10; % atm
Cpa = 35.0; % J/mol.K
Cao = Po / (0.0821 * To); % mol/dm^3
Cbo = 0.55555; % mol/dm^3
theta_b = Cbo / Cao;
R = 8.314; % J/mol.K
Ea = 41800; % J/mol
dHrx = -45000; % J/mol
alpha = 0.015; % kg^-1
% Calculation of k1
k1 = 0.5 * sqrt((n + a) / (n + b));
k1 = round(k1, round_figs);
disp(['Calculated Rate Constant k1: ', num2str(k1)]);
% --- Rate constant as function of T (Arrhenius equation) ---
k_func = @(T) k1 .* exp((Ea / R) .* (1 / To - 1 ./ T));
% Rate law Function
rate_func = @(X, T, yP) ...
k_func(T) .* ...
sqrt(max(Cao .* (1 - X) .* yP .* (To ./ T), 0)) .* ...
(max(Cao .* (theta_b - X) .* yP .* (To ./ T), 0)).^1.5;
% ODE System
% y(1) = X, y(2) = T, y(3) = yP
odefun = @(W, y) [
rate_func(y(1), y(2), y(3)) / Fao;
(-rate_func(y(1), y(2), y(3)) * dHrx) / (Fao * Cpa);
-alpha/2 * (1 - y(1)) * y(3) * (y(2) / To)
];
% Solve the ODEs
Wspan = [0 20];
y0 = [0 To 1];
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-9, 'NonNegative', [1 2 3]);
[W, Y] = ode45(odefun, Wspan, y0, options);
% Extract results
X = min(max(Y(:,1), 0), 1); % Clamp X between 0 and 1
T = Y(:,2);
yP = max(Y(:,3), 0); % Avoid negative yP
% Concentrations
Ca = max(Cao .* (1 - X) .* yP .* (To ./ T), 0);
Cb = max(Cao .* (theta_b - X) .* yP .* (To ./ T), 0);
Cc = max(2 * Cao .* X .* yP .* (To ./ T), 0);
Cd = max(Cao .* X .* yP .* (To ./ T), 0);
% Reaction rate and heat
ra = rate_func(X, T, yP);
Q = ra .* dHrx;
% Pressure and Temperature arrays
P = yP * Po;
T0_array = To * ones(size(W));
P0_array = Po * ones(size(W));
k_arr = k_func(T);
% Results Table
ExtendedTable = table(W, X, T, T0_array, P, P0_array, yP, Q, ra, ...
Ca, Cb, Cc, Cd, k_arr, ...
'VariableNames', {'W_kg', 'Conversion_X', 'T_K', 'T0_K', ...
'P_atm', 'P0_atm', 'y', 'Q_J_per_kg', ...
'Rate_of_Reaction_mol_kg_min', 'Conc_A', 'Conc_B', 'Conc_C', 'Conc_D',
'Rate_Constant_k'});
disp('--- Extended Reactor Conditions and Concentrations ---');
disp(ExtendedTable(1:min(1000, height(ExtendedTable)), :));
% PLOTTING SECTION
% Plot 1: Conversion and Temperature
figure;
subplot(2,1,1);
yyaxis left
plot(W, X, 'b-', 'LineWidth', 2);
ylabel('Conversion (X)');
ylim([0 1.1]);
yyaxis right
plot(W, T, 'r-', 'LineWidth', 2);
ylabel('Temperature (K)');
xlabel('Catalyst Weight (kg)');
title('Conversion and Temperature Profiles');
legend('Conversion', 'Temperature', 'Location', 'best');
grid on;
% Plot 2: Pressure and Heat
subplot(2,1,2);
yyaxis left
plot(W, P, 'g-', 'LineWidth', 2);
ylabel('Pressure (atm)');
yyaxis right
plot(W, Q, 'm-', 'LineWidth', 2);
ylabel('Heat Generated (J/kg)');
xlabel('Catalyst Weight (kg)');
title('Pressure and Heat Generation Profiles');
legend('Pressure', 'Heat Generated', 'Location', 'best');
grid on;
% Plot 3: Concentrations
figure;
plot(W, Ca, 'b-', W, Cb, 'g--', W, Cc, 'r-.', W, Cd, 'k:', 'LineWidth', 2);
xlabel('Catalyst Weight (kg)');
ylabel('Concentration (mol/dm^3)');
legend('[A]', '[B]', '[C]', '[D]', 'Location', 'best');
title('Concentration Profiles along Catalyst Weight');
grid on;
filename = '[Link] [Link]'; % Name of the Excel file
writetable(ExtendedTable, filename);
disp(['Table successfully saved to ', filename]);