CAChemE.org - License (CC-BY)
In [1]:
# Loading our favorite python libraries
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import fsolve
plt.style.use('bmh')
Let's go...
In [13]:
#set constants, ideal gas law
Rg = 0.082 # atm·L/(mol·K)
alpha = np.array([-1, -1, 2, 1, 0]) # Stoichiometric matrix (the first element is the key component)
# kinetic data
k0 = np.exp(8.2) # kinetic constant (L/(mol·h))
Ea = 1000 * Rg # activation energy /R (K)
#thermodynamic data
Cpm = np.array([10, 15, 15, 12, 20]) # average heat capacity for each component (J/(mol·K)) [1×5]
DH_Tref = -12500 * alpha[2] # enthalpy change of the reaction (J/(mol of component C))
Tref = 293.15 # reference temperature for enthalpy (K)
# Set operating conditions of the reactor
P = 1 # constant pressure (Pa)
#Initial state of the reactor
mt = 500 # total amount of gas mixture fed into the reactor (g)
PM0 = 40 # the average molecular weight of the feed mixture (g/mol)
Nt0 = mt / PM0 # initial total number of moles of the gas mixture (mol)
N0 = np.array([4, 4, 0, 0, Nt0 - 4 * 2]) # inital number of moles for each component (mol) [1×5]
T0 = 473.15
t0 = 0
#Calculate paremeters for the reaction-reactor system
#Calculate the expansion factor for the key component:
'''
$\epsilon_{ik} = \frac{N_{k0}}{N_{total,0}} \frac{1}{(-\alpha_{ik})} \sum_{j=1}^S\alpha_{ij}$
'''
epsilon = N0[0] / Nt0 * (1 / (-alpha[0])) * np.sum([alpha]) # dimensionless
# Define the differential equation system
def dydXk(y,Xk):
t = y[0]
T = y[1]
# Write the Arrehnius law:
'''$k(T)=k_0e^{(-E_a/(RT))}$'''
k = k0 * np.exp(-Ea / (Rg * T))
# Express the volume as a function of the problem variables:
'''$V=\frac{N_{total,0}R_{g}T}{p}(1 + \varepsilon {X_k})$'''
V = Nt0 * Rg * T / P * (1 + epsilon * Xk) #L
# Combine the mole balance for the key component and the kinetic law
'''$$ \frac{{dt}}{{d{X_k}}} =
\frac{V} {{{N_{k0}k(T)}{{\left( {1 - {X_k}} \right)}^2}}} $$ '''
# Assume $Cp_j$ are independent of the temperature:
'''$Cp_j \neq f(T), Cp_j=Cp_j^m$'''
Cp = Cpm #(mol) [1×5]
# Calculate the heat capacity change on reaction,
'''$$ \Delta C_{p}=\sum_{j} \alpha_{ij}Cp_j $$'''
DCp = np.sum(alpha * Cp) # J/(mol·K), [1×5]=[1×5]×[5×1]
# As $Cp_j$ is constant then $\Delta Cp$ is also constant, and the enthalpy of the reaction can be expressed as $\Delta H=\Delta H^{ref}+\Delta Cp(T-T^{ref})$
DH = DH_Tref + DCp * (T - Tref) #J/mol
# Write the Energy Balance an solve for
'''$\;\frac{dT}{dX_k}$
$$ \frac{dT}{dX_k} =
\frac{N_{k0}(-\Delta H(T))}{ \sum \limits_{j=1}^S N_{j0}C{p_j} + N_{k0}\Delta Cp{X_k}}$$
'''
balance_1 = V / (N0[0] * k * (1 - Xk) ** 2)
balance_2 = N0[0] * (-DH)/(np.sum(N0*Cp)+N0[0] * DCp * Xk)
return (balance_1, balance_2)
#Set the range for the independent variable
#Xk0=0
#Xkf=0.9
Xkk = np.linspace(0, 0.9)
# Set initial values for the dependent variables
y0 = [t0, T0] #first dependent variable := time, second dependent variable := temperature
# Solve the ordinary differential equation system
# call the "ode45" MATLAB command
Y = odeint(dydXk, y0, Xkk)
# Plot the solution
t = Y[:,0]
T = Y[:,1]
fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(20,10))
# Time vs conversion
ax1.plot(t, Xkk)
ax1.set_xlabel('time / h',)
ax1.set_ylabel('Conversion')
ax1.set_title('''Constant pressure adiabatic stirred batch reactor with variable heat capacities
A+B -> 2C+D''')
#ax1.set_xlim([0,0.9])
#ax1.set_ylim([0,4])
# Temperature vs conversion
ax2.plot(T, Xkk)
ax1.set_xlabel('Temperaure / K',)
ax1.set_ylabel('Conversion')
Out[13]:
In [54]:
from IPython.html.widgets import interact
from IPython.display import clear_output, display, HTML
In [68]:
def CPressureAdiabaticBatchReactor (T0):
#set constants, ideal gas law
Rg = 0.082 # atm·L/(mol·K)
alpha = np.array([-1, -1, 2, 1, 0]) # Stoichiometric matrix (the first element is the key component)
# kinetic data
k0 = np.exp(8.2) # kinetic constant (L/(mol·h))
Ea = 1000 * Rg # activation energy /R (K)
#thermodynamic data
Cpm = np.array([10, 15, 15, 12, 20]) # average heat capacity for each component (J/(mol·K)) [1×5]
DH_Tref = -12500 * alpha[2] # enthalpy change of the reaction (J/(mol of component C))
Tref = 293.15 # reference temperature for enthalpy (K)
# Set operating conditions of the reactor
P = 1 # constant pressure (Pa)
#Initial state of the reactor
mt = 500 # total amount of gas mixture fed into the reactor (g)
PM0 = 40 # the average molecular weight of the feed mixture (g/mol)
Nt0 = mt / PM0 # initial total number of moles of the gas mixture (mol)
N0 = np.array([4, 4, 0, 0, Nt0 - 4 * 2]) # inital number of moles for each component (mol) [1×5]
#T0 = 473.15
t0 = 0
#Calculate paremeters for the reaction-reactor system
#Calculate the expansion factor for the key component:
'''
$\epsilon_{ik} = \frac{N_{k0}}{N_{total,0}} \frac{1}{(-\alpha_{ik})} \sum_{j=1}^S\alpha_{ij}$
'''
epsilon = N0[0] / Nt0 * (1 / (-alpha[0])) * np.sum([alpha]) # dimensionless
# Define the differential equation system
def dydXk(y,Xk):
t = y[0]
T = y[1]
# Write the Arrehnius law:
'''$k(T)=k_0e^{(-E_a/(RT))}$'''
k = k0 * np.exp(-Ea / (Rg * T))
# Express the volume as a function of the problem variables:
'''$V=\frac{N_{total,0}R_{g}T}{p}(1 + \varepsilon {X_k})$'''
V = Nt0 * Rg * T / P * (1 + epsilon * Xk) #L
# Combine the mole balance for the key component and the kinetic law
'''$$ \frac{{dt}}{{d{X_k}}} =
\frac{V} {{{N_{k0}k(T)}{{\left( {1 - {X_k}} \right)}^2}}} $$ '''
# Assume $Cp_j$ are independent of the temperature:
'''$Cp_j \neq f(T), Cp_j=Cp_j^m$'''
Cp = Cpm #(mol) [1×5]
# Calculate the heat capacity change on reaction,
'''$$ \Delta C_{p}=\sum_{j} \alpha_{ij}Cp_j $$'''
DCp = np.sum(alpha * Cp) # J/(mol·K), [1×5]=[1×5]×[5×1]
# As $Cp_j$ is constant then $\Delta Cp$ is also constant, and the enthalpy of the reaction can be expressed as $\Delta H=\Delta H^{ref}+\Delta Cp(T-T^{ref})$
DH = DH_Tref + DCp * (T - Tref) #J/mol
# Write the Energy Balance an solve for
'''$\;\frac{dT}{dX_k}$
$$ \frac{dT}{dX_k} =
\frac{N_{k0}(-\Delta H(T))}{ \sum \limits_{j=1}^S N_{j0}C{p_j} + N_{k0}\Delta Cp{X_k}}$$
'''
balance_1 = V / (N0[0] * k * (1 - Xk) ** 2)
balance_2 = N0[0] * (-DH)/(np.sum(N0*Cp)+N0[0] * DCp * Xk)
return (balance_1, balance_2)
#Set the range for the independent variable
#Xk0=0
#Xkf=0.9
Xkk = np.linspace(0, 0.9)
# Set initial values for the dependent variables
y0 = [t0, T0] #first dependent variable := time, second dependent variable := temperature
# Solve the ordinary differential equation system
# call the "ode45" MATLAB command
Y = odeint(dydXk, y0, Xkk)
# Plot the solution
t = Y[:,0]
T = Y[:,1]
fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(20,10))
# Time vs conversion
ax1.plot(t, Xkk)
ax1.set_xlabel('time / h',)
ax1.set_ylabel('Conversion')
ax1.set_title('''Constant pressure adiabatic stirred batch reactor with variable heat capacities
A+B -> 2C+D''')
#ax1.set_xlim([0,0.9])
#ax1.set_ylim([0,4])
# Temperature vs conversion
ax2.plot(T, Xkk)
ax1.set_xlabel('Temperaure / K',)
ax1.set_ylabel('Conversion')
In [69]:
#Widgets Interactivos
interact(CPressureAdiabaticBatchReactor,
t = (0, 5, 0.1), T = (1500, 2000, 1))
Out[69]:
In [ ]: