# Difference between revisions of "Isothermal reactor design - 2013"

Class date(s): 04 February to 27 February

• F2011: Chapter 5 and 6
• F2006: Chapter 4

### 07 February 2013 (05C)

to see the effect on pressure drop in the packed bed.

### 25 February 2013 (07A)

The example covered in class is based on example 4-8 in F2006 and example 6-2 in F2011.

### MATLAB

In a file called membrane.m:

function d_depnt__d_indep = membranem(indep, depnt)

% Dynamic balance for the packed bed reactor (PBR); demo problem class 05C
%
%    indep: the independent ODE variable, such as time or length
%    depnt: a vector of dependent variables
%
%    X = depnt(1) = the conversion
%    y = depnt(2) = the pressure ratio = P/P_0 = y
%
%    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);

% Constants
kDiff = 0.005;  % s^{-1}
k = 0.01;       % s^{-1}
KC = 50;        % mol.m^{-3}
P0 = 830600;    % Pa
T0 = 500;       % K
R  = 8.314;     % J/(mol.K)

% Algebraic equations
FT = FA + FB + FC;
CT0 = P0 / (R * T0);
CA = CT0 * FA / FT;
CB = CT0 * FB / FT;
CC = CT0 * FC / FT;
RB = kDiff * CB;
rA = -k * (CA - CB * CC / KC);
rB = -rA;
rC = -rA;

% 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) = rA;
d_depnt__d_indep(2) = rB - RB;
d_depnt__d_indep(3) = rC;


In a separate file (any name), for example: ode_driver.m, which will "drive" the ODE solver:

% Integrate the ODE
% -----------------

% The independent variable always requires an initial and final value:
indep_start = 0.0;  % m^3
indep_final = 0.4;  % m^3

% Set initial condition(s): for integrating variables (dependent variables)
FA_depnt_zero = 0.25;   % i.e. FA(V=0) = 15 mol/min = 0.25 mol/s
FB_depnt_zero = 0.0;    % i.e. FB(V=0) = 0 mol/s
FC_depnt_zero = 0.0;    % i.e. FC(V=0) = 0 mol/s

% Integrate the ODE(s):
[V, depnt] = ode45(@membranem, [indep_start, indep_final], [FA_depnt_zero, FB_depnt_zero, FC_depnt_zero]);

% Plot the results:
clf;
plot(V, depnt(:,1), 'b')
grid('on')
hold('on')
plot(V, depnt(:,2), 'g')
plot(V, depnt(:,3), 'r')
xlabel('Reactor volume, V [m^3]')
ylabel('F_A, F_B and F_C')
legend('F_A', 'F_B', 'F_C')


### Python

import numpy as np
from scipy import integrate
from matplotlib.pylab import (plot, grid, xlabel, ylabel, show, legend)

def membrane(indep, depnt):
"""
Dynamic balance for the membrane reactor; class 07A

indep: the independent ODE variable, such as time or length
depnt: a vector of dependent variables

X = depnt[0] = the conversion
y = depnt[1] = the pressure ratio = P/P_0 = y

Returns d(depnt)/d(indep) = a vector of ODEs
"""

# Assign some variables for convenience of notation
FA = depnt[0]
FB = depnt[1]
FC = depnt[2]

# Constants
kDiff = 0.005  # s^{-1}
k = 0.01       # s^{-1}
KC = 50        # mol.m^{-3}
P0 = 830600    # Pa
T0 = 500       # K
R  = 8.314     # J/(mol.K)

# Algebraic equations
FT = FA + FB + FC
CT0 = P0 / (R * T0)
CA = CT0 * FA / FT
CB = CT0 * FB / FT
CC = CT0 * FC / FT
R_B = kDiff * CB
rA = -k * (CA - CB * CC / KC)
rB = -rA
rC = -rA

# Output from this ODE function must be a COLUMN vector, with n rows
n = len(depnt)
d_depnt__d_indep = np.zeros((n,1))
d_depnt__d_indep[0] = rA
d_depnt__d_indep[1] = rB - R_B
d_depnt__d_indep[2] = rC
return d_depnt__d_indep

# The "driver" that will integrate the ODE(s):
# -----------
# Start by specifying the integrator:
# use vode with "backward differentiation formula"
r = integrate.ode(membrane).set_integrator('vode', method='bdf')

# Set the independent variable's range
indep_start = 0.0  # m^3
indep_final = 0.4  # m^3
delta = 0.01       # the results will be shown only at these delta points
num_steps = np.floor((indep_final - indep_start)/delta) + 1 # Number of steps: 1 extra for initial condition

# Set initial condition(s): for integrating variables (dependent variables)
FA_depnt_zero = 0.25   # i.e. FA(V=0) = 15 mol/min = 0.25 mol/s
FB_depnt_zero = 0.0    # i.e. FB(V=0) = 0 mol/s
FC_depnt_zero = 0.0    # i.e. FC(V=0) = 0 mol/s

r.set_initial_value([FA_depnt_zero, FB_depnt_zero, FC_depnt_zero], indep_start)

# Create vectors to store trajectories
V = np.zeros((num_steps, 1))
FA = np.zeros((num_steps, 1))
FB = np.zeros((num_steps, 1))
FC = np.zeros((num_steps, 1))
V[0] = indep_start
FA[0] = FA_depnt_zero
FB[0] = FB_depnt_zero
FC[0] = FC_depnt_zero

# Integrate the ODE(s) across each delta
k = 1
while r.successful() and k < num_steps:
r.integrate(r.t + delta)

# Store the results to plot later
V[k] = r.t
FA[k] = r.y[0]
FB[k] = r.y[1]
FC[k] = r.y[2]
k += 1

# All done!  Plot the trajectories:
plot(V, FA)
grid('on')
plot(V, FB)
plot(V, FC)
xlabel('Volume, $V$ [$m^3$]')
ylabel('$F_A, F_B, F_C$ [mol/s]')
legend(['$F_A$', '$F_B$', '$F_C$'])
show()


### Polymath

d(FA)/d(V) = rA
d(FB)/d(V) = rB - RB
d(FC)/d(V) = rC

FA(0) = 0.25  # mol/s
FB(0) = 0.0   # mol/s
FC(0) = 0.0   # mol/s

# Independent variable
V(0) = 0
V(f) = 0.4  #m^3

# Constants
kDiff = 0.005  # s^{-1}
k = 0.01       # s^{-1}
KC = 50        # mol.m^{-3}
P0 = 830600    # Pa
T0 = 500       # K
R  = 8.314     # J/(mol.K)

# Algebraic equations
FT = FA + FB + FC
CT0 = P0 / (R * T0)
CA = CT0 * FA / FT
CB = CT0 * FB / FT
CC = CT0 * FC / FT
RB = kDiff * CB
rA = -k * (CA - CB * CC / KC)
rB = -rA
rC = -rA