Multiple reactions - 2013

Class date(s):	06 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]

Textbook references

• F2011: Chapter 8
• F2006: Chapter 6

Suggested problems

Will be posted soon

Class materials

06 March 2013 (08B-2)

• Audio and video recording of the class.

07 March 2013

• Audio and 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.

# 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

• Audio and video recording of the class.

13 March 2013

• Audio and video recording of the class.

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

• Audio and 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.

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