Constant pressure adiabatic stirred batch reactor with variable heat capacities



In [1]:

%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...

CODE



In [11]:

#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)
ax2.set_xlabel('Temperaure / K',)
ax2.set_ylabel('Conversion')




Out[11]:

<matplotlib.text.Text at 0x855d780>



Widgets



In [12]:

from IPython.html.widgets import interact
from IPython.display import clear_output, display, HTML




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)
ax2.set_xlabel('Temperaure / K',)
ax2.set_ylabel('Conversion')




In [10]:

#Widgets Interactivos
T0 = (298., 1000., 1.))







In [ ]: