import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# Constants
= 8.314 # J/mol K
R = 273 # K
TR
# A -> B + C
# 0: A; 1: B; 2: C
= np.array ([-70.0, -50.0, -40.])*1000 # J/mol
HF = np.array ([40.0, 25.0, 15.0]) # J/mol K
CP
= 31.4 * 1e3 # J/mol
E = lambda t: 0.133 * np.exp( (E/R) * (1/450.0 - 1/t) ) # dm3/ kgcat s
k
def pfr (w, x, *args):
= x[0]
X = args
(ca0, fa0, T0, epsilon, delta_hr_tr, delta_cp, theta)
# Calculate T using energy balance
= T0 + X * (-delta_hr_tr)/( np.sum(theta * CP) + X * delta_cp)
T = ca0 * (( 1 - X ) / ( 1 + epsilon * X )) * (T0/T)
ca
= k(T) * ca # -r_A
rate = rate/fa0
dxdw
return [dxdw]
# Data
= np.array ([-1.0, 1.0, 1.0]) # stoichiometric coefficients
nu
= 20 # dm3/s
V_0 = 10 # atm
P0 = 450 # K
T0
= P0 * V_0/ (R * T0)
fa0 = fa0/V_0
ca0
# inlet mole fraction
= np.array ([1.0, 0.0, 0.0])
y0 = np.array([1.0, 0.0, 0.0])
theta
# Heat of reaction at 298K
= np.sum(HF * nu)
delta_hr_tr = np.sum(CP * nu)
delta_cp
= y0[0]
ya0 = ya0 * np.sum(nu)
epsilon
= (ca0, fa0, T0, epsilon, delta_hr_tr, delta_cp, theta)
args
# initial condition
= np.array([0])
initial_conditions = 50 # kg
w_final
= solve_ivp(pfr,
sol 0, w_final],
[
initial_conditions, =args,
args=True)
dense_output
= np.linspace(0,w_final, 1000)
w = sol.sol(w)[0]
x = T0 + x * (-delta_hr_tr)/( np.sum(theta * CP) + x * delta_cp) T
Solutions to workshop 07: Non-isothermal reactor design
Lecture notes for chemical reaction engineering
Try following problems from Fogler 5e P 11-5, P 11-6, P 12-6, P 12-21
We will go through some of these problems in the workshop.
P 11-5
The elementary, irreversible gas-phase reaction
is carried out adiabatically in a PFR packed with a catalyst. Pure A enters the reactor at a volumetric flow rate of 20 dm3/s, at a pressure of 10 atm, and a temperature of 450 K.
Additional information:
; ;
; ;
All heats of formation are referenced to 273 K.
Plot and then analyze the conversion and temperature down the plug-flow reactor until an 80% conversion (if possible) is reached. (The maximum catalyst weight that can be packed into the PFR is 50 kg.) Assume that .
Vary the inlet temperature and describe what you find.
Plot the heat that must be removed along the reactor ( Q vs. V) to maintain isothermal operation.
Now take the pressure drop into account in the PBR with . The reactor can be packed with one of two particle sizes. Choose one.
- Plot and then analyze the temperature, conversion, and pressure along the length of the reactor. Vary the parameters and to learn the ranges of values in which they dramatically affect the conversion.
Parameter | Value |
---|---|
273.00 K | |
10.00 bar | |
-20000.00 | |
0.00 | |
1.00 | |
2.673e-03 | |
0.053 mol/s |
Energy balance:
T = 450.00 + 500.00 X
Catalyst weight required to achieve 80% conversion: 43.14 kg
Temperature at 80% conversion: 850.18 K
Exit (W = 50.00 kg) conversion: 0.91
Exit (W = 50.00 kg) temperature: 903.19 K
="Conversion")
plt.plot(w,x, label
0,50)
plt.xlim(0,1)
plt.ylim(
plt.grid()
plt.legend()
'Catalyst weight ($kg$)')
plt.xlabel('Conversion')
plt.ylabel(
plt.show()
='Temperature')
plt.plot(w,T, label
0,50)
plt.xlim(
= (np.max(T) - np.min(T)) * 0.05
head_margin = np.min(T)
ylim_lower = np.max(T) + head_margin
ylim_upper
plt.ylim(ylim_lower, ylim_upper)
plt.grid()
plt.legend()
'Catalyst weight ($kg$)')
plt.xlabel('Temperature (K)')
plt.ylabel(
plt.show()
Effect of temperature
= np.arange(400, 501, 1)
T0_range = np.zeros(len(T0_range))
X_final = np.zeros(len(T0_range))
T_final = np.zeros(len(T0_range))
W_p8
for i in range(len(T0_range)):
= T0_range[i]
T0 = (ca0, fa0, T0, epsilon, delta_hr_tr, delta_cp, theta)
args
# initial condition
= np.array([0])
initial_conditions = 50 # kg
w_final
= solve_ivp(pfr,
sol 0, w_final],
[
initial_conditions, =args,
args=True)
dense_output
= np.linspace(0,w_final, 1000)
w = sol.sol(w)[0]
x = T0 + x * (-delta_hr_tr)/( np.sum(theta * CP) + x * delta_cp)
T
# add values to array for reporting
= x[-1]
X_final[i] = T[-1]
T_final[i]
# see if you reach 80% conversion
if x[-1] > 0.8:
= w[x>0.8][0] W_p8[i]
='Conversion')
plt.plot(T0_range,X_final, label
0],T0_range[-1])
plt.xlim(T0_range[0,1)
plt.ylim(
plt.grid()
plt.legend()
'Inlet temperature (K)')
plt.xlabel('Conversion')
plt.ylabel(
plt.show()
='Outlet temperature')
plt.plot(T0_range,T_final, label
0],T0_range[-1])
plt.xlim(T0_range[
= (np.max(T_final) - np.min(T_final)) * 0.05
head_margin = np.min(T_final) - head_margin
ylim_lower = np.max(T_final) + head_margin
ylim_upper print (np.max(T_final) + head_margin)
plt.ylim(ylim_lower, ylim_upper)
plt.grid()
plt.legend()
'Inlet temperature (K)')
plt.xlabel('Temperature (K)')
plt.ylabel(
plt.show()
1022.0467173119888
> 0], W_p8[W_p8 > 0], label='Catalyst weight (kg)')
plt.plot(T0_range[W_p8
0],T0_range[-1])
plt.xlim(T0_range[min(W_p8), np.max(W_p8))
plt.ylim(np.
plt.legend()
'Inlet temperature (K)')
plt.xlabel('Catalyst weight (kg)')
plt.ylabel(
0, 1,
plt.fill_between(T0_range, =W_p8==0,
where='gray',
color=0.3,
alpha=plt.gca().get_xaxis_transform())
transform
= np.min(T0_range) + 0.5 * (T0_range[W_p8 > 0][0] - np.min(T0_range))
x_pos = np.min(W_p8) + 0.5 * (np.max(W_p8) - np.min(W_p8))
y_pos
plt.text(x_pos, y_pos, 'Conversion < 0.8',
='center',
horizontalalignment='center',
verticalalignment='black',
color=10)
fontsize
plt.show()
Minimum inlet temperature required to achieve 80% conversion: 439 K
Heat that must be removed along the reactor to maintain isothermal operation:
= 450
T0 = (ca0, fa0, T0, epsilon, delta_hr_tr, delta_cp, theta)
args
# initial condition
= np.array([0])
initial_conditions = 50 # kg
w_final
= solve_ivp(pfr,
sol 0, w_final],
[
initial_conditions, =args,
args=True)
dense_output
= np.linspace(0,w_final, 1000)
w = sol.sol(w)[0]
x = T0 + x * (-delta_hr_tr)/( np.sum(theta * CP) + x * delta_cp)
T = delta_hr_tr + delta_cp * (T - TR)
delta_h = ca0 * (( 1 - x ) / ( 1 + epsilon * x )) * (T0/T)
ca = k(T) * ca # -r_A
rate = -rate*delta_h heat_removed
='Heat removed')
plt.plot(w, heat_removed, label
min(w), np.max(w))
plt.xlim(np.
= (np.max(heat_removed) - np.min(heat_removed)) * 0.05
head_margin = np.max(heat_removed) + head_margin
ylim_upper min(heat_removed), ylim_upper)
plt.ylim(np.
plt.legend()
'Catalyst weight (kg)')
plt.xlabel('Heat removed (J/s)')
plt.ylabel(
plt.show()
Pressure drop:
def pfr_with_pressure_drop (w, x, *args):
= x
X,p = args
(ca0, fa0, T0, epsilon, delta_hr_tr, delta_cp, theta, alpha)
# Calculate T using energy balance
= T0 + X * (-delta_hr_tr)/( np.sum(theta * CP) + X * delta_cp)
T = ca0 * p * (( 1 - X ) / ( 1 + epsilon * X )) * (T0/T)
ca
= k(T) * ca # -r_A
rate = rate/fa0
dxdw = -(alpha/2*p) * (T/T0) * (1 + epsilon * X)
dpdw
return [dxdw, dpdw]
= 450
T0
# initial condition
= np.array([0.0, 1.0])
initial_conditions = 50 # kg
w_final
= 0.0075 # 1/kg
alpha1 = (ca0, fa0, T0, epsilon, delta_hr_tr, delta_cp, theta, alpha1)
args
= solve_ivp(pfr_with_pressure_drop,
sol1 0, w_final],
[
initial_conditions, =args,
args=True)
dense_output
= np.linspace(0,w_final, 1000)
w = sol1.sol(w)[0]
x1 = T0 + x1 * (-delta_hr_tr)/( np.sum(theta * CP) + x1 * delta_cp)
T1 = sol1.sol(w)[1]
p1
= 0.019 # 1/kg
alpha2 = (ca0, fa0, T0, epsilon, delta_hr_tr, delta_cp, theta, alpha2)
args
= solve_ivp(pfr_with_pressure_drop,
sol2 0, w_final],
[
initial_conditions, =args,
args=True)
dense_output
= np.linspace(0,w_final, 1000)
w = sol2.sol(w)[0]
x2 = T0 + x2 * (-delta_hr_tr)/( np.sum(theta * CP) + x2 * delta_cp)
T2 = sol2.sol(w)[1] p2
=f'$\\alpha$ = {alpha1}')
plt.plot(w,x1, label=f'$\\alpha$ = {alpha2}')
plt.plot(w,x2, label
0,50)
plt.xlim(0,1)
plt.ylim(
plt.grid()
plt.legend()
'Catalyst weight ($kg$)')
plt.xlabel('Conversion')
plt.ylabel(
plt.show()
=f'$\\alpha$ = {alpha1}')
plt.plot(w,T1, label=f'$\\alpha$ = {alpha2}')
plt.plot(w,T2, label
0,50)
plt.xlim(
= np.min(np.concatenate((T1, T2)))
min_t = np.max(np.concatenate((T1, T2)))
max_t
= (max_t - min_t) * 0.05
head_margin = min_t + head_margin
ylim_lower = max_t + head_margin
ylim_upper
plt.ylim(ylim_lower, ylim_upper)
plt.grid()
plt.legend()
'Catalyst weight ($kg$)')
plt.xlabel('Temperature (K)')
plt.ylabel(
plt.show()
P 11-6
The irreversible endothermic vapor-phase reaction follows an elementary rate law
and is carried out adiabatically in a 500-dm3 PFR. Species A is fed to the reactor at a rate of 10 mol/min and a pressure of 2 atm. An inert stream is also fed to the reactor at 2 atm, as shown in Figure P11-6 B . The entrance temperature of both streams is 1100 K.
Additional information:
(T in degrees Kelvin);
; ; ;
First derive an expression for as a function of and .
Sketch the conversion and temperature profiles for the case when no inerts are present. Using a dashed line, sketch the profiles when a moderate amount of inerts are added. Using a dotted line, sketch the profiles when a large amount of inerts are added. Qualitative sketches are fine. Describe the similarities and differences between the curves.
Sketch or plot and then analyze the exit conversion as a function of . Is there a ratio of the entering molar flow rates of inerts (I) to A (i.e., at which the conversion is at a maximum? Explain why there “is” or “is not” a maximum.
What would change in parts (b) and (c) if reactions were exothermic and reversible with and at 1100 K?
Sketch or plot FB for parts (c) and (d), and describe what you find.
Plot the heat that must be removed along the reactor ( Q vs. V) to maintain isothermal operation for pure A fed and an exothermic reaction.
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.optimize import minimize
= lambda t: np.exp(34.34 - 34222/t)
k116 def pfr116 (V, x, *args):
= x[0]
X = args
(ca0, fa0, T0, epsilon, delta_hr_tr, sumCp)
# Calculate T using energy balance
= (- delta_hr_tr * X + sumCp* T0)/(sumCp)
T = ca0 * (( 1 - X ) / ( 1 + epsilon * X )) * (T0/T)
ca
= k116(T) * ca # -r_A
rate = rate/fa0
dxdw
return [dxdw]
# data
= 0.082 # dm^3 atm / mol K
R # A -> B + C
= 170 #J/mol K
cpA = 90 #J/mol K
cpB = 80 #J/mol K
cpC = 200 #J/mol K
cpI = 80000 # J/mol
delta_hr_tr
= 500 # dm^3
volume = 10 # mol/min
fa0 = 2 # atm
P0 = 1100 # K
T0
= {}
results = [0, 10, 100]
thetaIs
for thetaI in thetaIs:
= fa0 * thetaI # mol/min
fi0 = fa0 + fi0
ft
= P0/(R * T0)
ca0 = P0/(R * T0)
ci0
= (ca0 + ci0)/(thetaI + 1)
ca01 = 1/(1 + thetaI)
epsilon
= cpA + cpI * thetaI
sumCp
= (ca01, fa0, T0, epsilon, delta_hr_tr, sumCp)
args = np.array([0.0])
initial_conditions
= solve_ivp(pfr116,
sol116 0, volume],
[
initial_conditions, =args,
args=True)
dense_output
= np.linspace(0,volume, 1000)
v = sol116.sol(v)[0]
x = (- delta_hr_tr * x + sumCp* T0)/(sumCp)
T
= {'v': v, 'x': x, 'T': T} results[thetaI]
for thetaI, data in results.items():
'v'],
plt.plot(data['x'],
data[=f'$\\Theta_I$ = {thetaI}')
label
0,volume)
plt.xlim(0,1)
plt.ylim(
plt.grid()
plt.legend()
'Volume ($dm^3$)')
plt.xlabel('Conversion, $X$')
plt.ylabel(
plt.show()
= []
min_t = []
max_t
for thetaI, data in results.items():
= data['v']
v = data['T']
T
plt.plot(v, T, =f'$\\Theta_I$ = {thetaI}')
label
= np.min(T)
m
min_t.append(m)
= np.max(T)
m
max_t.append(m)
0,volume)
plt.xlim(
= np.min(min_t)
min_temp = np.max(max_t)
max_temp = (max_temp - min_temp) * 0.05
head_margin = min_temp + head_margin
ylim_lower = max_temp + head_margin
ylim_upper
plt.ylim(ylim_lower,ylim_upper)
plt.grid()
plt.legend()
'Volume ($dm^3$)')
plt.xlabel('Temperature (K)')
plt.ylabel(
plt.show()
Optimal :
= np.arange(0,15 + 0.5,0.5)
thetaIs = []
conv for thetaI in thetaIs:
= fa0 * thetaI # mol/min
fi0 = fa0 + fi0
ft
= P0/(R * T0)
ca0 = P0/(R * T0)
ci0
= (ca0 + ci0)/(thetaI + 1)
ca01 = 1/(1 + thetaI)
epsilon
= cpA + cpI * thetaI
sumCp
= (ca01, fa0, T0, epsilon, delta_hr_tr, sumCp)
args = np.array([0.0])
initial_conditions
= solve_ivp(pfr116,
sol116 0, volume],
[
initial_conditions, =args,
args=True)
dense_output
= np.linspace(0,volume, 1000)
v = sol116.sol(v)[0]
x = (- delta_hr_tr * x + sumCp* T0)/(sumCp)
T
-1])
conv.append(x[
plt.plot(thetaIs,conv)
min(thetaIs), np.max(thetaIs))
plt.xlim(np.
= np.min(conv)
min_x = np.max(conv)
max_x = (max_x - min_x) * 0.05
head_margin = min_x + head_margin
ylim_lower = max_x + head_margin
ylim_upper
plt.grid()
'$\\Theta_I$')
plt.xlabel('Conversion, $X$')
plt.ylabel(
plt.show()
def objective(thetaI):
= P0 / (R * T0)
ca0 = P0 / (R * T0)
ci0 = (ca0 + ci0) / (thetaI + 1)
ca01 = 1 / (1 + thetaI)
epsilon = cpA + cpI * thetaI
sumCp = (ca01, fa0, T0, epsilon, delta_hr_tr, sumCp)
args = np.array([0.0])
initial_conditions
= solve_ivp(pfr116, [0, volume], initial_conditions, args=args, dense_output=True)
sol = sol.y[0, -1] # Get the final conversion
final_x return -final_x # Minimize the negative of the final conversion to maximize it
# Constants
= 0.082 # dm^3 atm / mol K
R = 170 # J/mol K
cpA = 200 # J/mol K
cpI = 80000 # J/mol
delta_hr_tr = 500 # dm^3
volume = 10 # mol/min
fa0 = 2 # atm
P0 = 1100 # K
T0
# Optimization setup
= (0, 100) # Define bounds for thetaI as a tuple (min, max)
thetaI_bounds = minimize(objective, x0=[10], bounds=[thetaI_bounds], method='L-BFGS-B')
result
# Display results
= result.x[0] optimal_thetaI
The optimal is 8.187. Maximum conversion = 0.915
P 12-6
The endothermic liquid-phase elementary reaction
proceeds, substantially, to completion in a single steam-jacketed, continuous-stirred reactor (Table P12-6 B ). From the following data, calculate the steady-state reactor temperature:
Reactor volume: 125 gal;
Steam jacket area: 10 ft2
Jacket steam: 150 psig (365.9 F saturation temperature)
Overall heat-transfer coefficient of jacket, U: 150
Agitator shaft horsepower: 25 hp
Heat of reaction, Btu/lb-mol of A (independent of temperature)
P 12-21
The irreversible liquid-phase reactions
Reaction 1:
Reaction 2:
are carried out in a PFR with heat exchange. The following temperature profiles were obtained for the reactor and the coolant stream:
The concentrations of A, B, C, and D were measured at the point down the reactor where the liquid temperature, T, reached a maximum, and they were found to be CA = 0.1, CB = 0.2, CC= 0.5, and CD= 1.5, all in mol/dm3. The product of the overall heat-transfer coefficient and the heat-exchanger area per unit volume, Ua, is 10 . The entering molar flow rate of A is 10 mol/s.
Additional information
,
; at 400 K
;
- What is the activation energy for Reaction (1)?
Citation
@online{utikar2024,
author = {Utikar, Ranjeet},
title = {Solutions to Workshop 07: {Non-isothermal} Reactor Design},
date = {2024-03-24},
url = {https://cre.smilelab.dev/content/workshops/07-non-isothermal-reactor-design/solutions.html},
langid = {en}
}