Difference between revisions of "Multiple reactions - 2013"

From Introduction to Reactor Design: 3K4
Jump to navigation Jump to search
 
(23 intermediate revisions by the same user not shown)
Line 1: Line 1:
{{ClassSidebar
{{ClassSidebarYouTube
| date = 06 March
| date = 06 March to 14 March
| dates_alt_text =  
| dates_alt_text =  
| vimeoID1 = 61228970
| vimeoID1 = yPH7JvEi1Jw
| vimeoID2 = 61319707
| vimeoID2 = SpLniff-kcw
| vimeoID3 =
| vimeoID3 = fX8v4189Riw
| vimeoID4 =
| vimeoID4 = 0W48zitUYVw
| vimeoID5 =
| vimeoID5 = MvRWk6R0oSc
| course_notes_PDF =  
| course_notes_PDF =  
| course_notes_alt = Course notes
| course_notes_alt = Course notes
Line 19: Line 19:
| video_notes2 =
| video_notes2 =
| video_download_link3_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-09A.mp4
| video_download_link3_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-09A.mp4
| video_download_link3_MP4_size = M
| video_download_link3_MP4_size = 211 M
| video_notes3 =
| video_notes3 =
| video_download_link4_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-09B.mp4
| video_download_link4_MP4_size = 331 M
| video_notes4 =
| video_download_link5_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-09C.mp4
| video_download_link5_MP4_size = 243 M
| video_notes5 =
}}__NOTOC__
}}__NOTOC__


Line 29: Line 35:


== Suggested problems ==
== Suggested problems ==
''Will be posted soon''
 
<!--
{| class="wikitable"
{| class="wikitable"
|-
|-
Line 36: Line 41:
! F2006
! F2006
|-
|-
| Problem 7-7 (a)
| Problem 8-12
| Problem 5-6 (a)
| Problem 6-12
|-
|-
| Problem 7-8 (a)
| Problem 8-14 (covered in class)
| Problem 5-7 (a)
| Problem 6-15 (covered in class)
|-
|-
| Problem 7-15
| Problem 8-18 (set up equations)
| Not in this edition
| Problem 6-21 (set up equations)
|} -->
|}


==Class materials ==
==Class materials ==
Line 50: Line 55:
=== 06 March 2013 (08B-2) ===
=== 06 March 2013 (08B-2) ===


* [http://learnche.mcmaster.ca/media/3K4-2013-Class-08B-2.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-08B-2.mp4 video] recording of the class
* [http://learnche.mcmaster.ca/media/3K4-2013-Class-08B-2.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-08B-2.mp4 video] recording of the class.


=== 07 March 2013 ===
=== 07 March 2013 (08C) ===


* [http://learnche.mcmaster.ca/media/3K4-2013-Class-08C.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-08C.mp4 video] recording of the class
* [http://learnche.mcmaster.ca/media/3K4-2013-Class-08C.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-08C.mp4 video] recording of the class.


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.
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.
Line 71: Line 76:
k1 = 0.5 # 1/hr
k1 = 0.5 # 1/hr
k2 = 0.2 # 1/hr
k2 = 0.2 # 1/hr
S_DU = (k1*CA - k2*CB) / (k2*CB+1E-15)
 
Overall_SDU = CB/(CC+1E-15)
# The 3 important algebraic variables: plot these 3 against time and interpret them.
Yield = CB/(2 - CA+1E-15)
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
# Independent variable details
Line 79: Line 86:
t(f) = 3.1  # hours
t(f) = 3.1  # hours
</syntaxhighlight>
</syntaxhighlight>
=== 11 March 2013 (09A) ===
* [http://learnche.mcmaster.ca/media/3K4-2013-Class-09A.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-09A.mp4 video] recording of the class.
=== 13 March 2013 (09B) ===
* [http://learnche.mcmaster.ca/media/3K4-2013-Class-09B.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-09B.mp4 video] recording of the class.
Code for the CSTR example:
<syntaxhighlight lang="matlab">
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
</syntaxhighlight>
=== 14 March 2013 (09C) ===
* [http://learnche.mcmaster.ca/media/3K4-2013-Class-09C.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-09C.mp4 video] recording of the class.
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.
{| class="wikitable"
|-
! MATLAB
! Polymath
|-
|
'''<tt>pfr.m</tt>'''
<syntaxhighlight lang="matlab">
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);
</syntaxhighlight>
'''<tt>ODE_driver.m</tt>'''
<syntaxhighlight lang="matlab">
% 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')
</syntaxhighlight>
| valign="top"|
<syntaxhighlight lang="text">
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
</syntaxhighlight>
|}

Latest revision as of 17:28, 25 January 2017

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