Simple Reactors

Batch reactors

Cases

TP: constant T, P
TV: constant T, V
HP: constant H, P (adiabatic, constant pressure)
UV: constant U, V (adiabatic, constant volume)

Species equation

Simple mass balance, where $m_i$ is the mass of species $i$, $V$ is the reactor volume (possibly changing in time), and $\dot{m}_i'''$ is the net species mass production rate per unit volume.

$$\frac{dm_i}{dt} = \dot{m}_i'''V.$$

Now, $m_i= my_i$, where $y_i$ is mass fraction and $m$ is the total mass of the reactor, which is constant. Note, if the reactor volume is constant, then so is $\rho$. This gives

$$\frac{dmy_i}{dt} = m\frac{dy_i}{dt} = \dot{m}_i'''V,$$

$$\frac{dy_i}{dt} = \frac{\dot{m}_i'''}{\rho}.$$

This equation holds for any of the four cases listed.

We can write $\dot{m}_i''' = \dot{\omega}_iM_i$, where $\omega_i$ is a species molar production rate with units of kmol/m$^3$*s. Cantera provides $\omega_i$ via the property net_production_rates.

Energy equations

An energy equation is not needed for cases TP, or TV. For cases HP and UV, in a code like Cantera, we can set the gas state as either gas.HPY = H, P, Y, or gas.UVY = U, V, Y, where $V$ is $1/\rho$.
- Temperature is then available simply as T = gas.T, where inverting $H(T)=H$ to solve for $T$ is done internally by Cantera.
Otherwise, we can simply write $$\frac{dh}{dt} = 0,$$ or $$\frac{du}{dt} = 0.$$

Temperature equation

Consider $dh/dt=0$.

$$h = h(T, y_i),$$

$$dh = \underbrace{\frac{\partial h}{\partial T}}_{c_p}dT = \sum_i\underbrace{\frac{\partial h}{\partial y_i}}_{h_i}dy_i.$$

Then divide through by $dt$ and solve for $dT/dt$ with $dh/dt = 0$.

$$\frac{dT}{dt} = -\frac{1}{\rho c_p}\sum_ih_i\dot{m}_i'''.$$

For $du/dt=0$, we have

$$\frac{dT}{dt} = -\frac{1}{\rho c_v}\sum_iu_i\dot{m}_i'''.$$

CSTR

Here we have an unsteady reactor with an inlet and outlet stream. The reactor is well mixed so the outlet composition is the same as in the reactor. ### Species equation.
The mass balance in this case gives: $$\frac{dm_i}{dt} = \dot{m}_i^{in} - \dot{m}_i^{out} + \dot{m}_i'''V.$$
Assume the mass in the reactor is constant in time (this is not true in general). Also assume that $\dot{m}^{in}=\dot{m}^{out}$, which may also not be true.

$$\frac{dy_i}{dt} = \frac{y_i^{in}-y_i}{\tau} + \frac{\dot{m}_i'''}{\rho}.$$

Here, $\tau=m/\dot{m}$ is a reactor residence time.
When solved to steady state, the resulting composition is independent of the two assumptions listed above.
If we set $\frac{dy_i}{dt}=0$, we can solve for species compositions using a nonlinear algebraic solver, rather than time advance the ODEs given using an ODE solver. The ODE solution will be more robust to convergence.

Energy equations

An enthalpy equation can be written as $$\frac{dh}{dt} = \frac{h^{in}-h}{\tau}.$$
A similar equation for internal energy can be formulated, as well as a temperature equation using the approach shown above.
At steady state, we simply have $h-h^{in}$.

CSTR versus batch

Note, the CSTR equations for species and enthalpy are very similar to the respective equations for batch reactor, but with the additional mixing term involving $\tau$.
Hence, a single solver can be written to solve both types of reactors.
It is common to solve a CSTR to steady state, in which case, some measure of the approach to steady state is needed.

Example code

Solve a Batch reactor HP problem

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.integrate import odeint
import cantera as ct

#-------------------- setup

P    = 101235                      # Pa
T0   = 1500                        # K
x0   = "CH4:1, O2:2, N2:7.52"      # -
trun = 0.005                       # run time (s)
nt   = 1000                        # number of times to record

#-------------------- define rate function to integrate

def rhsf(y, t):
    gas.HPY = h0, P, y
    return gas.net_production_rates * gas.molecular_weights / gas.density

#-------------------- initial conditions and constant enthalpy

gas = ct.Solution('gri30.yaml')
gas.TPX = T0, P, x0
y0 = gas.Y
h0 = gas.h

#-------------------- solve problem

t = np.linspace(0, trun, nt)

y = odeint(rhsf, y0, t)

#-------------------- recover the temperature profile corresponding to y(t)

T = np.empty(nt)
for i in range(nt):
    gas.HPY = h0, P, y[i,:]
    T[i] = gas.T
    
#-------------------- plot results

plt.rc('font', size=14)

plt.figure(figsize=(8,5))
plt.plot(t, T)
plt.xlabel('t (s)')
plt.ylabel('T (K)');

plt.figure(figsize=(8,5))
plt.plot(t, y[:,gas.species_index("CH4")])
plt.plot(t, y[:,gas.species_index("O2")])
plt.plot(t, y[:,gas.species_index("CO2")])
plt.plot(t, y[:,gas.species_index("CO")])
plt.legend(['CH4', 'O2', 'CO2', 'CO'], frameon=False)
plt.xlabel('t (s)')
plt.ylabel('y')
plt.xlim([0,trun]);