Difference between revisions of "Steady-state nonisothermal reactor design - 2013"

From Introduction to Reactor Design: 3K4
Jump to navigation Jump to search
m
 
(30 intermediate revisions by the same user not shown)
Line 1: Line 1:
{{ClassSidebar
{{ClassSidebarYouTube
| date = 18 March
| date = 18 March
| dates_alt_text =  
| dates_alt_text =  
| vimeoID1 = 62121732
| vimeoID1 = FnsyFP1QoP0
| vimeoID2 = 62301489
| vimeoID2 = 8YVhVL5dVIU
| vimeoID3 =
| vimeoID3 = GJxYdIZ-nCI
| vimeoID4 =
| vimeoID4 = -kMr6kk54Pk
| vimeoID5 =
| vimeoID5 = yLXZSyJK05E
| vimeoID6 = bAPUtqIWAsI
| vimeoID7 = 6JIZAhz7Y6o
| vimeoID8 = xy6N_-Ob4GA
| vimeoID9 = gaEvAyo29Yw
| course_notes_PDF =  
| course_notes_PDF =  
| course_notes_alt = Course notes
| course_notes_alt = Course notes
Line 19: Line 23:
| video_notes2 =
| video_notes2 =
| video_download_link3_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-10C.mp4
| video_download_link3_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-10C.mp4
| video_download_link3_MP4_size = M
| video_download_link3_MP4_size = 392 M
| video_notes3 =
| video_notes3 =
| video_download_link4_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-11A.mp4
| video_download_link4_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-11A.mp4
| video_download_link4_MP4_size = M
| video_download_link4_MP4_size = 342 M
| video_notes4 =
| video_notes4 =
| video_download_link4_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-11B.mp4
| video_download_link5_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-11B.mp4
| video_download_link5_MP4_size = M
| video_download_link5_MP4_size = 355 M
| video_notes5 =
| video_notes5 =
}}__NOTOC__
| video_download_link6_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-11C.mp4
| video_download_link6_MP4_size = 392 M
| video_notes6 =
| video_download_link7_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-12A.mp4
| video_download_link7_MP4_size = 331 M
| video_notes7 =
| video_download_link8_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-12B.mp4
| video_download_link8_MP4_size = 556 M
| video_notes8 =
| video_download_link9_MP4 = http://learnche.mcmaster.ca/media/3K4-2013-Class-12C-part1.mp4
| video_download_link9_MP4_size = 180 M
| video_notes9 =}}
__NOTOC__


== Textbook references ==
== Textbook references ==
Line 35: Line 51:


== Suggested problems ==
== Suggested problems ==
''Will be posted soon''
 
<!--
{| class="wikitable"
{| class="wikitable"
|-
|-
Line 42: Line 57:
! F2006
! F2006
|-
|-
| Problem 7-7 (a)
| Problem 12-6
| Problem 5-6 (a)
| Problem 8-5
|-
| Problem 12-15 (a)
| Problem 8-16 (a)
|-
|-
| Problem 7-8 (a)
| Problem 12-16 (a)
| Problem 5-7 (a)
| Problem 8-18 (a)
|-
|-
| Problem 7-15
| Problem 12-24 (set up equations)
| Not in this edition
| Problem 8-26 (set up equations)
|} -->
|}


==Class materials ==
==Class materials ==
Line 56: Line 74:
* '''18 March 2013 (10A)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-10A.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-10A.mp4 video] recording of the class.
* '''18 March 2013 (10A)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-10A.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-10A.mp4 video] recording of the class.
* '''20 March 2013 (10B)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-10B.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-10B.mp4 video] recording of the class; and the [[Media:3K4-2013-Class-10B.pdf|handout used in class]].
* '''20 March 2013 (10B)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-10B.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-10B.mp4 video] recording of the class; and the [[Media:3K4-2013-Class-10B.pdf|handout used in class]].
* '''21 March 2013 (10C)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-10C.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-10C.mp4 video] recording of the class; and the [[Media:3K4-2013-Class-10C.pdf|revised handout]]
* '''25 March 2013 (11A)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-11A.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-11A.mp4 video] recording of the class.
* '''27 March 2013 (11B)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-11B.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-11B.mp4 video] recording of the class.
* '''28 March 2013 (11C)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-11C.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-11C.mp4 video] recording of the class.
* '''01 April 2013 (12A)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-12A.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-12A.mp4 video] recording of the class.
* '''03 April 2013 (12B)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-12B.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-12B.mp4 video] recording of the class.
* '''04 April 2013 (12C-part1)''': [http://learnche.mcmaster.ca/media/3K4-2013-Class-12C-part1.mp3 Audio] and [http://learnche.mcmaster.ca/media/3K4-2013-Class-12C-part1.mp4 video] recording of the class.


==Source code ==
==Source code ==
Line 75: Line 100:


% Constants. Make sure to use SI for consistency
% Constants. Make sure to use SI for consistency
FT0 = 163000/3600;% mol/s (was kmol/hour originally)
FT0 = param.FT0; % mol/s.  Note how we use the "struct" variable to access the total flow
FA0 = 0.9 * FT0;  % mol/s
FA0 = 0.9 * FT0;  % mol/s
T1 = 360;        % K
T1 = 360;        % K
Line 82: Line 107:
R = 8.314;        % J/(mol.K)
R = 8.314;        % J/(mol.K)
HR = -6900;      % J/(mol of n-butane)
HR = -6900;      % J/(mol of n-butane)
CA0 = 9300;      % kmol/m^3
CA0 = 9300;      % mol/m^3
k_1 = 31.1/3600;  % 1/s (was 1/hr originally)
k_1 = 31.1/3600;  % 1/s (was 1/hr originally)
K_Cbase = 3.03;  % [-]
K_Cbase = 3.03;  % [-]
Line 104: Line 129:
</syntaxhighlight>
</syntaxhighlight>


'''<tt>driver.m</tt>'''
'''<tt>driver_ode.m</tt>'''
<syntaxhighlight lang="matlab">
<syntaxhighlight lang="matlab">
% Integrate the ODE
% Integrate the ODE
Line 116: Line 141:
X_depnt_zero = 0.0;  % i.e. X(V=0) = 0.0
X_depnt_zero = 0.0;  % i.e. X(V=0) = 0.0


% Other parameters
% Other parameters. The "param" is just a variable in MATLAB.
param.T_0 = 330;     % feed temperature [K]
% It is called a structured variable, or just a "struct"
% It can have an arbitrary number of sub-variables attached to it.
% In this example we have two of them.
param.T_0 = 330;         % feed temperature [K]
param.FT0 = 163000/3600; % mol/s (was kmol/hour originally)


% Integrate the ODE(s):
% Integrate the ODE(s):
[indep, depnt] = ode45(@pfr_example, [indep_start, indep_final], [X_depnt_zero], optimset(), param);
[indep, depnt] = ode45(@pfr_example, [indep_start, indep_final], [X_depnt_zero], odeset(), param);


% Deal with the integrated output to show interesting plots
% Deal with the integrated output to show interesting plots
Line 129: Line 158:
for i = 1:numel(X)
for i = 1:numel(X)
     [rA_over_FA0(i), C_A(i)] = pfr_example([], X(i), param);
     [rA_over_FA0(i), C_A(i)] = pfr_example([], X(i), param);
end
end % the above can be done more efficiently in a single line, but the code would be too confusing


% Plot the results
% Plot the results
figure
f=figure;
set(f, 'Color', [1,1,1])
subplot(2, 2, 1)
subplot(2, 2, 1)
plot(indep, X); grid
plot(indep, X); grid
xlabel('Volume, V [kg]')
xlabel('Volume, V [kg]', 'FontWeight', 'bold')
ylabel('Conversion, X [-]')
ylabel('Conversion, X [-]', 'FontWeight', 'bold')


subplot(2, 2, 2)
subplot(2, 2, 2)
plot(indep, T); grid
plot(indep, T); grid
xlabel('Volume, V [kg]')
xlabel('Volume, V [kg]', 'FontWeight', 'bold')
ylabel('Temperature profile [K]')
ylabel('Temperature profile [K]', 'FontWeight', 'bold')


subplot(2, 2, 3)
subplot(2, 2, 3)
plot(indep, C_A); grid
plot(indep, C_A); grid
xlabel('Volume, V [kg]')
xlabel('Volume, V [kg]', 'FontWeight', 'bold')
ylabel('Concentration C_A profile [K]')
ylabel('Concentration C_A profile [K]', 'FontWeight', 'bold')


subplot(2, 2, 4)
subplot(2, 2, 4)
plot(indep, rA_over_FA0); grid
plot(indep, rA_over_FA0); grid
xlabel('Volume, V [kg]')
xlabel('Volume, V [kg]', 'FontWeight', 'bold')
ylabel('(Reaction rate/FA0) profile [1/m^3]')
ylabel('(Reaction rate/FA0) profile [1/m^3]', 'FontWeight', 'bold')
 
% Now plot one of the most important figures we saw earlier in the course:
% F_A0 / (-rA) on the y-axis, against conversion X on the x-axis. This plot
% is used to size various reactors.
 
% The material leaves the reactor at equilibrium; let's not plot
% that far out, because it distorts the scale. So plot to 95% of
% equilibrium
f = figure; set(f, 'Color', [1,1,1])
index = find(X>0.95 * max(X), 1);
plot(X(1:index), 1./rA_over_FA0(1:index)); grid
xlabel('Conversion, X [-]', 'FontWeight', 'bold')
ylabel('FA0/(-r_A) profile [m^3]', 'FontWeight', 'bold')
 
% Updated code to find the optimum inlet temperature:
temperatures = 300:3:400;
conv_at_exit = zeros(size(temperatures));
i = 1;
for T =
  param.T_0 = T;        % feed temperature [K]
  [indep, depnt] = ode45(@pfr_example, [indep_start, indep_final], ...
                      [X_depnt_zero], odeset(), param);
  conv_at_exit(i) = depnt(end, 1);
  i = i + 1;
end
plot(temperatures, conv_at_exit); grid;
xlabel('Inlet temperature')
ylabel('Conversion at reactor exit')
</syntaxhighlight>
</syntaxhighlight>

Latest revision as of 19:35, 6 January 2017

Class date(s): 18 March
Download video: Link (plays in Google Chrome) [399 M]

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

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

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

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

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

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

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

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


Textbook references

  • F2011: Chapter 11 and 12
  • F2006: Chapter 8

Suggested problems

F2011 F2006
Problem 12-6 Problem 8-5
Problem 12-15 (a) Problem 8-16 (a)
Problem 12-16 (a) Problem 8-18 (a)
Problem 12-24 (set up equations) Problem 8-26 (set up equations)

Class materials

Source code

Example on 21 March (class 10C)

pfr_example.m

function [d_depnt__d_indep, CA] = pfr_example(indep, depnt, param)
 
% 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
X = depnt(1);

% Constants. Make sure to use SI for consistency
FT0 = param.FT0;  % mol/s.  Note how we use the "struct" variable to access the total flow
FA0 = 0.9 * FT0;  % mol/s
T1 = 360;         % K
T2 = 333;         % K
E = 65700;        % J/mol
R = 8.314;        % J/(mol.K)
HR = -6900;       % J/(mol of n-butane)
CA0 = 9300;       % mol/m^3
k_1 = 31.1/3600;  % 1/s (was 1/hr originally)
K_Cbase = 3.03;   % [-]

% Equations
T = 43.3*X + param.T_0;  % derived in class, from the heat balance
k1 = k_1 * exp(E/R*(1/T1 - 1/T)); % temperature dependent rate constant
KC = K_Cbase * exp(HR/R*(1/T2 - 1/T)); % temperature dependent equilibrium constant
k1R = k1 / KC; % reverse reaction rate constant
CA = CA0 * (1 - X); % from the stoichiometry (differs from Fogler, but we get same result)
CB = CA0 * (0 + X);
r1A = -k1 * CA; % rate expressions derived in class 
r1B = -r1A;
r2B = -k1R * CB;
r2A = -r2B;
rA = r1A + r2A; % total reaction rate for species A

n = numel(depnt);
d_depnt__d_indep = zeros(n,1);
d_depnt__d_indep(1) = -rA / FA0;

driver_ode.m

% Integrate the ODE
% -----------------
 
% The independent variable always requires an initial and final value:
indep_start = 0.0;  % m^3
indep_final = 5.0; % m^3
 
% Set initial condition(s): for integrating variables (dependent variables)
X_depnt_zero = 0.0;   % i.e. X(V=0) = 0.0

% Other parameters. The "param" is just a variable in MATLAB.
% It is called a structured variable, or just a "struct"
% It can have an arbitrary number of sub-variables attached to it.
% In this example we have two of them.
param.T_0 = 330;         % feed temperature [K]
param.FT0 = 163000/3600; % mol/s (was kmol/hour originally)

% Integrate the ODE(s):
[indep, depnt] = ode45(@pfr_example, [indep_start, indep_final], [X_depnt_zero], odeset(), param);

% Deal with the integrated output to show interesting plots
X = depnt(:,1);
T = 43.3.*X + param.T_0;          % what was the temperature profile?
rA_over_FA0 = zeros(numel(X), 1); % what was the rate profile?
C_A = zeros(numel(X), 1);         % what was the concentration profile?
for i = 1:numel(X)
    [rA_over_FA0(i), C_A(i)] = pfr_example([], X(i), param);
end % the above can be done more efficiently in a single line, but the code would be too confusing

% Plot the results
f=figure;
set(f, 'Color', [1,1,1])
subplot(2, 2, 1)
plot(indep, X); grid
xlabel('Volume, V [kg]', 'FontWeight', 'bold')
ylabel('Conversion, X [-]', 'FontWeight', 'bold')

subplot(2, 2, 2)
plot(indep, T); grid
xlabel('Volume, V [kg]', 'FontWeight', 'bold')
ylabel('Temperature profile [K]', 'FontWeight', 'bold')

subplot(2, 2, 3)
plot(indep, C_A); grid
xlabel('Volume, V [kg]', 'FontWeight', 'bold')
ylabel('Concentration C_A profile [K]', 'FontWeight', 'bold')

subplot(2, 2, 4)
plot(indep, rA_over_FA0); grid
xlabel('Volume, V [kg]', 'FontWeight', 'bold')
ylabel('(Reaction rate/FA0) profile [1/m^3]', 'FontWeight', 'bold')

% Now plot one of the most important figures we saw earlier in the course: 
% F_A0 / (-rA) on the y-axis, against conversion X on the x-axis. This plot
% is used to size various reactors.

% The material leaves the reactor at equilibrium; let's not plot
% that far out, because it distorts the scale. So plot to 95% of 
% equilibrium
f = figure; set(f, 'Color', [1,1,1])
index = find(X>0.95 * max(X), 1);
plot(X(1:index), 1./rA_over_FA0(1:index)); grid
xlabel('Conversion, X [-]', 'FontWeight', 'bold')
ylabel('FA0/(-r_A) profile [m^3]', 'FontWeight', 'bold')

% Updated code to find the optimum inlet temperature:
temperatures = 300:3:400;
conv_at_exit = zeros(size(temperatures));
i = 1;
for T = 
   param.T_0 = T;         % feed temperature [K]
   [indep, depnt] = ode45(@pfr_example, [indep_start, indep_final], ...
                       [X_depnt_zero], odeset(), param);
   conv_at_exit(i) = depnt(end, 1);
   i = i + 1;
end
plot(temperatures, conv_at_exit); grid;
xlabel('Inlet temperature')
ylabel('Conversion at reactor exit')