Multiple reactions - 2013

From Introduction to Reactor Design: 3K4
Jump to: navigation, search
Class date(s): 06 March to 14 March
Download video: Link (plays in Google Chrome) [142 M]

Download video: Link (plays in Google Chrome) [367 M]

Download video: Link (plays in Google Chrome) [211 M]

Download video: Link (plays in Google Chrome) [331 M]

Download video: Link (plays in Google Chrome) [243 M]

Textbook references

  • F2011: Chapter 8
  • F2006: Chapter 6

Suggested problems

F2011 F2006
Problem 8-12 Problem 6-12
Problem 8-14 (covered in class) Problem 6-15 (covered in class)
Problem 8-18 (set up equations) Problem 6-21 (set up equations)

Class materials

06 March 2013 (08B-2)

07 March 2013 (08C)

Polymath code for example in class. Make sure you plot the instantaneous selectivity, overall selectivity and yield over time. Compare these 3 plots during the batch to understand what each of these 3 variables mean.

# ODEs
d(CA) / d(t) = -k1*CA 
d(CB) / d(t) = k1*CA - k2*CB
d(CC) / d(t) = k2*CB
 
# Initial conditions
CA(0) = 2 # mol/L
CB(0) = 0 # mol/L
CC(0) = 0 # mol/L
 
# Algebraic equations
k1 = 0.5 # 1/hr
k2 = 0.2 # 1/hr
 
# The 3 important algebraic variables: plot these 3 against time and interpret them.
S_DU = if (t>0.001) then (k1*CA - k2*CB) / (k2*CB)  else 0
Overall_SDU = if (t>0.001) then CB/CC else 0
Yield = if (t>0.001) then CB / (2 - CA) else 0
 
# Independent variable details
t(0) = 0
t(f) = 3.1  # hours

11 March 2013 (09A)

13 March 2013 (09B)

Code for the CSTR example:

tau = 0:0.05:10;
 
CA0 = 2;  % mol/L
k1 = 0.5; % 1/hr
k2 = 0.2; % 1/hr
 
CA = CA0 ./ (1 + k1 .* tau);
CB = tau .* k1 .* CA ./ (1 + k2 .* tau);
CC = tau .* k2 .* CB;
 
instant_selectivity = (k1.*CA - k2.*CB) ./ (k2.*CB);
overall_selectivity = CB ./ CC;
overall_yield = CB ./ (CA0 - CA);
conversion = (CA0 - CA)./CA0;
 
plot(tau, CA, tau, CB, tau, CC)
grid on
xlabel('\tau')
ylabel('Concentrations [mol/L]')
 
figure 
plot(tau, overall_selectivity)
xlabel('\tau')
ylabel('Overall Selectivity')
grid on
 
figure 
plot(tau, overall_yield)
xlabel('\tau')
ylabel('Overall Yield')
grid on
 
figure
plot(tau, conversion)
xlabel('\tau')
ylabel('Conversion')
hold on
grid on

14 March 2013 (09C)

Despite the fact that Polymath code is shorter to write, I still recommend you use MATLAB or Python. For example, comparing two simulations and generating plots is so much easier in MATLAB than Polymath.

MATLAB Polymath

pfr.m

function d_depnt__d_indep = pfr(indep, depnt)
% Dynamic balance for the reactor
% 
%    indep: the independent ODE variable, such as time or length
%    depnt: a vector of dependent variables
% 
%    Returns d(depnt)/d(indep) = a vector of ODEs
 
% Assign some variables for convenience of notation
FA = depnt(1);
FB = depnt(2);
FC = depnt(3);
FD = depnt(4);
FE = depnt(5);
FG = depnt(6);
FW = depnt(7);
y  = depnt(8);
 
% Constant(s)
k1 = 0.014;    % L^{0.5} / mol^{0.5} / s
k2 = 0.007;    % L/(mol.s)
k3 = 0.14;     % 1/s
k4 = 0.45;     % L/(mol.s)
alpha = 0.002; % 1/L
CT0 = 1.0;     % mol/L
FA0 = 10;      % mol/s
FB0 = 5.0;     % mol/s
FT0 = FA0 + FB0;
 
FT = FA + FB + FC + FD + FE + FW + FG;
 
CA = CT0 * FA/FT * y;
CB = CT0 * FB/FT * y;
CC = CT0 * FC/FT * y;
CD = CT0 * FD/FT * y;
CE = CT0 * FE/FT * y;
CG = CT0 * FG/FT * y;
CW = CT0 * FW/FT * y;
 
% Reaction 1: A + 0.5B -> C
r1A = -k1*(CA)*(CB)^(0.5);
r1B = 0.5*r1A;
r1C = -r1A;
 
%# Reaction 2: 2A -> D
r2A = -k2*(CA)^2;
r2D = -r2A/2;
 
% Reaction 3: C -> E + W
r3C = -k3*(CC);
r3E = -r3C;
r3W = -r3C;
 
% Reaction 4: D + W -> G + C
r4D = -k4*(CD)*(CW);
r4W = r4D;
r4G = -r4D;
r4C = -r4D;
 
% Output from this ODE function must be a COLUMN vector, with n rows
n = numel(depnt);
d_depnt__d_indep = zeros(n,1);
d_depnt__d_indep(1) = r1A + r2A;
d_depnt__d_indep(2) = r1B;
d_depnt__d_indep(3) = r1C + r3C + r4C;
d_depnt__d_indep(4) = r2D + r4D;
d_depnt__d_indep(5) = r3E;
d_depnt__d_indep(6) = r3W + r4W;
d_depnt__d_indep(7) = r4G;
d_depnt__d_indep(8) = -alpha/(2*y) * (FT / FT0);

ODE_driver.m

% Integrate the ODE
% -----------------
 
% The independent variable: requires an initial and final value:
indep_start = 0.0;   % kg
indep_final = 500.0; % kg
 
% Set initial condition(s) for dependent variables
FA_depnt_zero = 10.0;   % i.e. FA(W=0) = 10.0
FB_depnt_zero = 5.0;    % i.e. FB(W=0) = 10.0
FC_depnt_zero = 0.0;    % i.e. FC(W=0) = 10.0
FD_depnt_zero = 0.0;    % etc
FE_depnt_zero = 0.0;
FG_depnt_zero = 0.0;
FW_depnt_zero = 0.0;
y_depnt_zero = 1.0;     % i.e. y(W=0) = 1.0
 
IC = [FA_depnt_zero, FB_depnt_zero, FC_depnt_zero, FD_depnt_zero, ...
      FE_depnt_zero FG_depnt_zero, FW_depnt_zero, y_depnt_zero];
 
% Integrate the ODE(s):
[indep, depnt] = ode45(@pfr, [indep_start, indep_final], IC);
 
% Calculate Yields and Selectivities
FA = depnt(:,1);
FC = depnt(:,3);
FD = depnt(:,4);
FE = depnt(:,5);
 
Yield_C = FC ./ (FA_depnt_zero - FA);
S_CE = FC./FE;
S_CD = FC./FD;
 
% Plot the results:
clf;
plot(indep, depnt(:,1), indep, depnt(:,2), indep, depnt(:,3), ...
     indep, depnt(:,4), indep, depnt(:,5), indep, depnt(:,6), ...
     indep, depnt(:,7), indep, depnt(:,8))
grid('on')
hold('on')
plot(indep, depnt(:,2), 'g')
xlabel('Catalyst weight, W [kg]')
ylabel('Concentrations and pressure drop')
legend('FA', 'FB', 'FC', 'FD', 'FE', 'FG', 'FW', 'y')
k1 = 0.014    # L^{0.5} / mol^{0.5} / s
k2 = 0.007    # L/(mol.s)
k3 = 0.14     # 1/s
k4 = 0.45     # L/(mol.s)
alpha = 0.002 # 1/L
CT0 = 1.0     # mol/L
FA0 = 10      # mol/s
FB0 = 5.0     # mol/s
FT0 = FA0 + FB0
 
# Concentration functions (isothermal conditions)
CA = CT0 * FA/FT * y
CB = CT0 * FB/FT * y
CC = CT0 * FC/FT * y
CD = CT0 * FD/FT * y
CE = CT0 * FE/FT * y
CG = CT0 * FG/FT * y
CW = CT0 * FW/FT * y
 
FT = FA + FB + FC + FD + FE + FW + FG
 
# Reaction 1: A + 0.5B -> C
r1A = -k1*(CA)*(CB)^(0.5)
r1B = 0.5*r1A
r1C = -r1A
 
# Reaction 2: 2A -> D
r2A = -k2*(CA)^2
r2D = -r2A/2
 
# Reaction 3: C -> E + W
r3C = -k3*(CC)
r3E = -r3C
r3W = -r3C
 
# Reaction 4: D + W -> G + C
r4D = -k4*(CD)*(CW)
r4W = r4D
r4G = -r4D
r4C = -r4D
 
# ODE's
d(FA) / d(W) = r1A + r2A
FA(0) = 10
 
d(FB) / d(W) = r1B
FB(0) = 5
 
d(FC) / d(W) = r1C + r3C + r4C
FC(0) = 0
 
d(FD) / d(W) = r2D + r4D
FD(0) = 0
 
d(FE) / d(W) = r3E
FE(0) = 0
 
d(FW) / d(W) = r3W + r4W
FW(0) = 0
 
d(FG) / d(W) = r4G
FG(0) = 0
 
W(0) = 0    # kg
W(f) = 500 # kg
 
Yield_C = if (W>0.001) then (FC / (FA0 - FA)) else (0)
S_CE = if (W>0.001) then (FC/FE) else (0)
S_CD = if (W>0.001) then (FC/FD) else (0)
 
d(y) / d(W) = -alpha/(2*y) * (FT / FT0)
y(0) = 1.0