Steady-state nonisothermal reactor design - 2013
Jump to navigation Jump to search
|Class date(s):||18 March|
- F2011: Chapter 11 and 12
- F2006: Chapter 8
Will be posted soon
- 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.
- 21 March 2013 (10C): Audio and video recording of the class; and the revised handout
Example on 21 March (class 10C)
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;
% 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')