Difference between revisions of "Steady-state nonisothermal reactor design - 2013"
Jump to navigation
Jump to search
Kevin Dunn (talk | contribs) m |
Kevin Dunn (talk | contribs) m |
||
Line 55: | Line 55: | ||
* '''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 ( | * '''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]]. | ||
==Source code == | |||
=== Example on 20 March (class 10B) === | |||
<tt>'''pfr_example.m'''</tt> | |||
<syntaxhighlight lang="matlab"> | |||
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 = 163000/3600;% mol/s (was kmol/hour originally) | |||
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; % kmol/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; | |||
</syntaxhighlight> | |||
'''<tt>driver.m</tt>''' | |||
<syntaxhighlight lang="matlab"> | |||
% 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 | |||
param.T_0 = 330; % feed temperature [K] | |||
% Integrate the ODE(s): | |||
[indep, depnt] = ode45(@pfr_example, [indep_start, indep_final], [X_depnt_zero], optimset(), 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 | |||
% Plot the results | |||
figure | |||
subplot(2, 2, 1) | |||
plot(indep, X); grid | |||
xlabel('Volume, V [kg]') | |||
ylabel('Conversion, X [-]') | |||
subplot(2, 2, 2) | |||
plot(indep, T); grid | |||
xlabel('Volume, V [kg]') | |||
ylabel('Temperature profile [K]') | |||
subplot(2, 2, 3) | |||
plot(indep, C_A); grid | |||
xlabel('Volume, V [kg]') | |||
ylabel('Concentration C_A profile [K]') | |||
subplot(2, 2, 4) | |||
plot(indep, rA_over_FA0); grid | |||
xlabel('Volume, V [kg]') | |||
ylabel('(Reaction rate/FA0) profile [1/m^3]') | |||
</syntaxhighlight> |
Revision as of 16:19, 21 March 2013
Class date(s): | 18 March | ||||
| |||||
| |||||
Textbook references
- F2011: Chapter 11 and 12
- F2006: Chapter 8
Suggested problems
Will be posted soon
Class materials
- 18 March 2013 (10A): Audio and video recording of the class.
- 20 March 2013 (10B): Audio and video recording of the class; and the handout used in class.
Source code
Example on 20 March (class 10B)
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 = 163000/3600;% mol/s (was kmol/hour originally)
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; % kmol/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.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
param.T_0 = 330; % feed temperature [K]
% Integrate the ODE(s):
[indep, depnt] = ode45(@pfr_example, [indep_start, indep_final], [X_depnt_zero], optimset(), 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
% Plot the results
figure
subplot(2, 2, 1)
plot(indep, X); grid
xlabel('Volume, V [kg]')
ylabel('Conversion, X [-]')
subplot(2, 2, 2)
plot(indep, T); grid
xlabel('Volume, V [kg]')
ylabel('Temperature profile [K]')
subplot(2, 2, 3)
plot(indep, C_A); grid
xlabel('Volume, V [kg]')
ylabel('Concentration C_A profile [K]')
subplot(2, 2, 4)
plot(indep, rA_over_FA0); grid
xlabel('Volume, V [kg]')
ylabel('(Reaction rate/FA0) profile [1/m^3]')