Multiple reactions  2013
Class date(s):  06 March  
 
 
 
 
Textbook references
 F2011: Chapter 8
 F2006: Chapter 6
Suggested problems
Will be posted soon
Class materials
06 March 2013 (08B2)
07 March 2013
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
13 March 2013
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
Despite the fact that Polymath code is shorter to write, I still recommend you use MATLAB or Python. For example, comparing simulations 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
%
% X = depnt(1) = the conversion
% y = depnt(2) = the pressure ratio = P/P_0 = y
%
% 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 always requires an initial and final value:
indep_start = 0.0; % kg
indep_final = 500.0; % kg
% Set initial condition(s): for integrating variables (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
