1. Reaction in series

Written on December, 2024 by Bolton Tran

Consider a gas-phase reaction in series: $A \overset{k_1}{\rightarrow} B \overset{k_2}{\rightarrow} C$ taking place in an isothermal batch reactor, we want to solve for the concentrations of $A$, $B$, and $C$ versus time.

We write the rate of change of concentration of each species with respect to time as follows:

\begin{align} \frac{dC_A}{dt} &= -k_1 C_A\\ \frac{dC_B}{dt} &= k_1 C_A - k_2 C_B\\ \frac{dC_C}{dt} &= k_2 C_B \end{align}

We can solve the above ordinary differential equations (ODEs) once values for rate constants $k_1$ and $k_2$ are known. We can use an approximate form of transition state theory (TST) to calculate $k_1$ and $k_2$:

\begin{align} k_1 &= \frac{k_BT}{h}\exp{\left(-\frac{E^\dagger_1}{k_BT}\right)}\\ k_2 &= \frac{k_BT}{h}\exp{\left(-\frac{E^\dagger_2}{k_BT}\right)} \end{align}

The activation energies $E^\dagger_1$ and $E^\dagger_2$ and the temperature $T$ control reaction rates. From a simple collision theory perspective, larger activation energies slow reaction down since the probability of particles’ collision having enough energy for chemical reaction is lower; higher reaction temperature speeds reaction up since the average kinetic energy of particles and the chance of them colliding increase. Equivalent TST interpretation involves invoking the shift in a quasi-equilibrium between the transition state and the initial state.

$T$ (K): 400

$E^\dagger_1$ (kJ/mol): 140

$E^\dagger_2$ (kJ/mol): 150

Solving ODEs is a must-have skill for chemical engineers. While many ODE solver in different coding languages (MATLAB, Wolfram, Scipy) are available, knowing how to set up an ODE problem in one code allows for easy transfer to another code. Demonstration using Python/Scipy is shown below:

import numpy as np
from scipy.integrate import odeint #module for ODE solver

#Constants
kB=8.3144/1000 #kJ/mol.K
h=4.135667696E-15*96.49/3600 #kJ/mol.h

#Solve kinetic ODE at a given temperature and EAs
def solveODE(EA1, EA2, T):
    #Rate constants
    k1=(kB*T/h)*np.exp(-EA1/kB/T)
    k2=(kB*T/h)*np.exp(-EA2/kB/T)  
    #Define ODEs
    def dcdt(c,t):
        #c is array for concentrations
        cA=c[0]; cB=c[1]; cC=c[2]
        dcAdt=-k1*cA
        dcBdt=k1*cA - k2*cB
        dcCdt=k2*cB
        return [dcAdt, dcBdt, dcCdt]
    #Initial parameters
    c0=[2,0,0] #mol/L
    #solve ODE
    cs=odeint(dcdt, c0, time)
    return cs

#Define a time range
t0=0 
t1=50
time=np.linspace(t0,t1,100) #hr

#Invoke ODEsolve function for a given set of parameters
para=[150,160,400] #EA1, EA2, T
result=solveODE(*para) #solve
print(result)

For anyone who is curious: the interactive plot was created by hosting a python-Flask (a Web Server Gateway Interface application) script on https://www.pythonanywhere.com/. The script contains python functions that perform mathematical operations. The sliders send values of ($E^\dagger_1$, $E^\dagger_2$, $T$) to the PythonAnywhere app, which performs the above Scipy-ODE solver and outputs results as json arrays. The json arrays are then parsed back to Plotly written in JavaScript, which created the plot interactively.

This interactive exercise is heavily inspired by the LearnChemE group at the University of Colorado Boulder. Check out their library for other cool interaction simulations for chemical engineering students.