Packed bed reactor problem

I'm trying to plot varying flowrates against volume but am not sure how cause the flowrates to vary.

My matlab code is below

clc
clear all

global Ea K T R

% Activation energy of the reaction
Ea = 23000; % kcal/mol

% Pre-exponential factor of reaction
K = 426; % kmol/kg s

% Ideal gas constant
R = 1.9872; % kcal/ K kmol

% PBR dimensions
d = 1.07; % m
L = 2.98; % m

% Operating conditions
T = 450 + 273; % K
P = 100; % bar


% Setting inital values of flowrates and temperature in the reactor as Y0
Y0(1) = 816.993; % kmol/hr of nitrogen
Y0(2) = 2450.979; % kmol/hr of hydrogen
Y0(3) = 0; % kmol/hr of ammonia

% Volume span for ode45 and graphing
vol_span = [0,2.67];

% Ode45 solving the system of odes
[Z,Y] = ode45(@HBfunc, vol_span, Y0);

% Plotting a graph of flowrate against volume
figure(1)
plot(Z,Y(:,1:3))
xlabel("Volume(m^3)")
ylabel("Flowrate(kmol/hr)")

% Creating flowrate variables
FN = Y(:,1); % kmol/hr
FH = Y(:,2); % kmol/hr
FAm = Y(:,3); % kmol/hr










function dYdZ = HBfunc(Z,Y)

global Ea K T R 

% Creating flowrate variables that change with volume
FN = Y(1); % kmol/hr
FH = Y(2); % kmol/hr
FAm = Y(3); % kmol/hr


% Adding all flowrates together
Ft = FN + FH + FAm; % kmol/hr

% Reaction rate
r = K*(exp(-Ea/R*T)); % kmol/kg s

% Creating ODE's for all components
dFNdZ = -(r); % kmol/kg s
dFHdZ = -(3*r); % kmol/kg s
dFAmdZ = 2*r; % kmol/kg s

% Creating ODE for temperature
%fill here

% Putting ODE's into dYdZ
dYdZ(1) = dFNdZ;
dYdZ(2) = dFHdZ;
dYdZ(3) = dFAmdZ;

dYdZ = dYdZ';

end

[Torsten]

r = K*exp(-Ea/(R*T)); % kmol/kg s

instead of

r = K*(exp(-Ea/R*T)); % kmol/kg s

And check your units:

Ea/(R*T) has unit kcal/mol / ( kcal/ (K kmol) * K) = 10^3 

Answer 1

[Ishu]

Hi Raheel Razaq,

I understand that you are trying to plot varying flowrates against volume for a packed bed reactor problem.

As in general, reaction flow rate also depends upon "time" and "temperature". And in your function "HBfunc()" you have not used both of these because of why flowrate is not changing properly. To vary flow rates you can use "Ramp", "sin" functionalities. So in "FN", "FH" and "FAm" you have to add one of these functionalities according to your problem statement to vary the flowrate. I have shown an example below.

FN = Y(1) + sin(pi*Z);  
FH = Y(2) + 0.5*sin(pi*Z); 
FAm = Y(3) + 0.2*sin(pi*Z);

And further you have to use these flowrates in the function to calculate the exact results. Like use these while calculating "dFNdZ", "dFHdZ" and "dFAmdZ". You can try these changes at your end and then try to plot the graph. Also accordingly add other dependencies on which your flow rate depends.

For more information you can refer these documentations:

ode45: https://www.mathworks.com/help/matlab/ref/ode45.html

sin: https://www.mathworks.com/help/matlab/ref/sin.html?searchHighlight=sin&s_tid=srchtitle_support_results_1_sin

Hope it helps.