Constant pressure adiabatic stirred batch reactor with variable heat capacities

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

CODE


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]:
[<matplotlib.lines.Line2D at 0xb3c5eb8>]

Widgets


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]:
<function __main__.CPressureAdiabaticBatchReactor>

In [ ]: