Canonical Reactors I¶

Zero-dimensional reactors.
- No spatial dependence.
- Evolve in time (or space for the PFR).
Batch reactors
- Constant TP
- Constant TV
- Constant HP
- Constant UV
Plug flow reactors (PFR)
Perfectly stirred reactor (PSR), also called a continuous stirred tank reactor (CSTR)
- These are flow reactors.

Approach¶

Couple kinetics with conservation laws for heat and mass balances.
These reactors provide convenient, well-characterized baselines for more complicated flow processes.
All equations below are on a mass basis: $h$ is J/kg, etc.

Batch reactor: TP¶

Batch reactor is a closed system $\to$ constant mass.
Specify an initial state: $T$ , $P$ , $y_{i}$ .
Mass balance: ${\dot{m}}_{a c c u m} = {\dot{m}}_{i n} - {\dot{m}}_{o u t} + {\dot{m}}_{g e n e r a t i o n}$
${\dot{m}}_{i n} = {\dot{m}}_{o u t} = 0$ , so we have ${\dot{m}}_{a c c u m} = {\dot{m}}_{g e n e r a t i o n}$
Now, write in the details, where ${\dot{m}}_{i}^{'''}$ is the net mass generation rate of species $i$ with units of kg/m $^{3}$ *s, and $V$ is volume. $\frac{d m_{i}}{d t} = {\dot{m}}_{i}^{'''} V,$

$\frac{d m y_{i}}{d t} = {\dot{m}}_{i}^{'''} V,$

$m = ρ V$ is constant: $\frac{d y_{i}}{d t} = \frac{{\dot{m}}_{i}^{'''} V}{ρ V},$

\frac{d y_{i}}{d t} = \frac{{\dot{m}}_{i}^{'''}}{ρ} .

Also, we have $ρ = \frac{M P}{R T},$

M = \frac{1}{\sum_{i} y_{i} / M_{i}},

{\dot{m}}_{i}^{'''} = {\dot{m}}_{i}^{'''} (T, P, y_{i}) .

That is, we have $N_{s p}$ ODEs for each species mass fraction. All other terms are written in terms of the $y_{i}$ and the known $T$ , $P$ .
The last blue equation indicates the functional dependence of the species rates, which are given in terms of kinetic expressions for the given combustion mechanism.
At constant pressure, the volume will vary. The volume per unit mass is simply $1 / ρ$ , where $ρ$ at any time is given by the ideal gas law in terms of the composition at any time.

Batch reactor: TV¶

If volume (or more likely, volume per unit mass, so that the 0-D system is intensive) is given instead of pressure, then the system density is constant (and is known from the initial state.)
The previous equations hold, but $ρ$ is known.
The ideal gas law can be used to solve for $P$ , which now varies.

Batch reactor: HP¶

For an adiabatic, constant pressure reactor, the previous equations all hold, but now temperature will vary.
The energy conservation equation becomes

h = constant,

\frac{d h}{d t} = 0.

The system enthalpy is known from the initial state.
At any given time, the temperature may be found by solving $h = h (T, y_{i})$ for $T$ given the known $h$ , and the current $y_{i}$ .
In Cantera, if h0 is our enthalpy, this calculation is done internally, and we would simply write:

gas.HPY = h0, P, y
T = gas.T

Temperature equation¶

Alternatively, we can derive a temperature equation as follows.

$\frac{d h}{d t} = 0,$

$h = h (T, y_{i}),$

$d h = \underset{c_{p}}{\underset{⏟}{\frac{\partial h}{\partial T}}} d T + \sum_{i} \underset{h_{i}}{\underset{⏟}{\frac{\partial h}{\partial y_{i}}}} d y_{i},$

Now divide through by $d t$ :

$\frac{d h}{d t} = c_{p} \frac{d T}{d t} + \sum_{i} h_{i} \frac{d y_{i}}{d t} = 0,$

Use $d y_{i} / d t = {\dot{m}}_{i}^{'''} / ρ$ , and solve for $d T / d t$ :

\frac{d T}{d t} = - \frac{1}{ρ c_{p}} \sum_{i} h_{i} {\dot{m}}_{i}^{'''} .

Also, $c_{p} = \sum_{i} y_{i} c_{p, i},$

h_{i} = h_{f, i} (T_{r e f}) + \int_{T_{r e f}}^{T} c_{p, i} d T .

Batch reactor: UV¶

If we have an adiabatic, constant volume batch reactor, then then energy equation is written in terms of internal energy $u$ instead of enthalpy, and $c_{v}$ instead of $c_{p}$ .
Also, $ρ$ is constant as for the TV system above.

u = constant,

\frac{d u}{d t} = 0.

\frac{d T}{d t} = - \frac{1}{ρ c_{v}} \sum_{i} u_{i} {\dot{m}}_{i}^{'''} .

Also, $c_{v} = \sum_{i} y_{i} c_{v, i},$

u_{i} = u_{f, i} (T_{r e f}) + \int_{T_{r e f}}^{T} c_{v, i} d T .

Plug flow reactors (PFR)¶

Instead of a batch reactor evolving in time, we have a flow reactor evolving in space (flow in a tube, say).
Radial mixing is perfect, and the "plug" of fluid moves down the tube without mixing axially.
Hence, the plug behaves as a batch reactor, reacting in time. But the plug is at a different spatial location at each time.
Convert time $t$ to spatial location $z$ :

$\frac{d}{d t} = \frac{d}{d z} \frac{d z}{d t} = v \frac{d}{d z} .$

That is, we transform from time to space using the local velocity: $v = d z / d t \to d t = d z / v$ .
Now, for a PFR, the mass flux is constant (even though the local velocity may vary):

{\dot{m}}^{''} = constant

${\dot{m}}^{''} = ρ v \to v = {\dot{m}}^{''} / ρ$

Hence,

\frac{d}{d t} \to \frac{{\dot{m}}^{''}}{ρ} \frac{d}{d z}

So, in all the previous equations, wherever you see $d / d t$ , replace it with the relation above.
You will have to specify the mass flux.

Perfectly stirred reactor (PSR)¶

This is a flow reactor with an inlet and an outlet.
The outlet composition is the same as the reactor composition.
This can be constant temperature or adiabatic.
Characterized by the reactor size (volume $V$ ) and flow rate $\dot{m}$ .
- These can be combined as $τ = ρ V / \dot{m}$ , where $τ$ is the PSR residence time.
- The PSR has one key parameter, $τ$ , which is the mixing timescale of the reactor.
We can write steady or unsteady versions of the reactor.

Unsteady PSR¶

Assumptions:
- Mass $m$ in the PSR is constant.
- ${\dot{m}}_{i n} = {\dot{m}}_{o u t}$ .
Species balance equation: ${\dot{m}}_{a c c u m} = {\dot{m}}_{i n} - {\dot{m}}_{o u t} + {\dot{m}}_{g e n e r a t i o n},$

$\frac{d m y_{i}}{d t} = \dot{m} y_{i, i n} - \dot{m} y_{i} + {\dot{m}}_{i}^{'''} V,$

$\frac{d y_{i}}{d t} = \frac{\dot{m}}{m} (y_{i, i n} - y_{i}) + \frac{{\dot{m}}_{i}^{'''} V}{m},$

\frac{d y_{i}}{d t} = \frac{(y_{i, i n} - y_{i})}{τ} + \frac{{\dot{m}}_{i}^{'''}}{ρ} .

Energy balance equation.
- If the reactor is isothermal, we just specify a constant temperature.
- For an adiabatic reactor (with constant $m$ and ${\dot{m}}_{i n} = {\dot{m}}_{o u t}$ , we have

$\frac{d h}{d t} = h_{i n} - h .$

If the initial $h$ is equal to $h_{i n}$ (so make the initial composition equal to the inlet composition for consistency), then we have $\frac{d h}{d t} = 0,$ $h = h_{0} = constant .$
Temperature can be found implicitly by solving $h = h (T, y_{i})$ , as noted above.
Alternatively, we can formulate and solve a temperature equation:

Temperature equation¶

The temperature equation parallels the treatment for a batch reactor, where we also had the $d h / d t = 0$ equation:

$\frac{d h}{d t} = c_{p} \frac{d T}{d t} + \sum_{i} h_{i} \frac{d y_{i}}{d t} = 0,$

Solve for $d T / d t$ and use our PSR equation above for $d y_{i} / d t$ :

$\frac{d T}{d t} = - \frac{1}{c_{p}} \sum_{i} h_{i} \frac{d y_{i}}{d t} .$

\frac{d T}{d t} = - \frac{1}{c_{p}} \sum_{i} h_{i} (\frac{(y_{i, i n} - y_{i})}{τ} + \frac{{\dot{m}}_{i}^{'''}}{ρ}) .

Solution approach¶

Note that the equations for an unsteady PSR are very similiar to those for a batch reactor, and we can use the same code to solve both systems.
Often, an unsteady PSR is solved to steady state by solving the ODE system for long enough.
Alternatively, we can directly solve the steady problem by setting the $d / d t$ terms to zero. We then have a coupled system of nonlinear algebraic equations that we can solve using, e.g., Newton's method.

PSR Note¶

Consider again the unsteady PSR equations for species:

$\frac{d y_{i}}{d t} = \frac{(y_{i, i n} - y_{i})}{τ} + \frac{{\dot{m}}_{i}^{'''}}{ρ} .$

At steady state we have $\frac{(y_{i} - y_{i, i n})}{τ} = \frac{{\dot{m}}_{i}^{'''}}{ρ},$

$(mixing term) = (reaction term) .$

That is, mixing balaces reaction.
$τ$ is the mixing timescale and $1 / τ$ is a mixing rate. Hence, we have (mixing rate)=(reaction rate).
As $τ$ decreases, the mixing rate increases, and eventually the chemical reaction rate cannot match it.
- At that point, the PSR blows out.
- The balance still holds, but we get the trivial solution 0 = 0.
- We can decrease $τ$ by increasing $\dot{m}$ , or by decreasing volume $V$ .

Chemical timescale¶

Combustion chemistry is complex and there is no one reaction rate (or no one reaction timescale).
- However, at the point just before blowout, we can take $τ_{0}$ as the characteristic chemical timescale.
- We can then compare $τ_{0}$ to some other mixing timescale to define a Damkohler number: $D a = τ / τ_{0}$ .

Analogy with diffusive mixing¶

Consider the unsteady diffusion equation:

$\frac{\partial y_{i}}{\partial t} = D \frac{\partial^{2} y_{i}}{\partial x^{2}} .$

If we nondimensionalize this equation, then we need a timescale $τ$ , a lengthscale $L$ , and a scale $y_{i, r e f}$ (which cancels).
The nondimensional form is $\frac{1}{τ} \frac{\partial y_{i}^{*}}{\partial t^{*}} = \frac{D}{L^{2}} \frac{\partial^{2} y_{i}^{*}}{\partial x^{* 2}} .$
Here, starred variables are nondimensional.
If we chose appropriate reference values, then the magnitude of terms $\partial y_{i}^{*} / \partial t^{*}$ and $\partial^{2} y_{i}^{*} / \partial x^{* 2}$ are both $O (1)$ .
This gives $τ = \frac{L^{2}}{D}$ as the characteristic diffusion timescale.
Hence, we can write:

$D \frac{\partial^{2} y_{i}}{\partial x^{2}} \sim D \frac{Δ y_{i}}{L^{2}} \sim \frac{Δ y_{i}}{τ}$

But this is just our PSR mixing term. The point is, that the PSR mixing term can be used as an analogy to more complex diffusive mixing in terms of a mixing timescale, where mixing and reaction processes are in a kind of balance.

HP Batch Reactor code¶

In [1]:

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

In [2]:

gas = ct.Solution("gri30.yaml")

#------------- Set batch reactor equation

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

#------------- set initial condition and enthalpy

gas.TPX = 1400, 101325, "CH4:1, O2:2, N2:7.52"
h0 = gas.enthalpy_mass
y0 = gas.Y

#------------- solve the ODE system

nt = 1000
times = np.linspace(0,0.01,nt)
y = odeint(rhsf, y0, times)

#------------- recover the temperature

T = np.zeros(nt)
for i in range(nt):
    gas.HPY = h0, gas.P, y[i,:]
    T[i] = gas.T

#------------- equilibrium state

gas.equilibrate("HP")
Teq = gas.T
yeq = gas.Y

#------------- plot results

plt.rc("font", size=14)
plt.plot(times*1000,T, label="T")
plt.plot(times[-1]*1000,Teq, 'kd', label="Teq")
plt.xlabel("time (ms)")
plt.ylabel("T (K)")
plt.legend(frameon=False);

plt.figure()
species = ["CH4", "O2", "CO2", "CO", "OH"]
for sp in species:
    plt.plot(times*1000,y[:,gas.species_index(sp)], label=sp)
    plt.plot(times[-1]*1000,yeq[gas.species_index(sp)], 'kd', label="eq" if sp==species[-1] else "")
plt.xlabel("time (ms)")
plt.ylabel("y")
plt.legend(frameon=False);

No description has been provided for this image

In [3]:

plt.plot(times*1000, y[:,gas.species_index("CO2")]/y[-1,gas.species_index("CO2")])
plt.plot(times*1000, y[:,gas.species_index("NO")]/y[-1,gas.species_index("NO")])

Out[3]:

[<matplotlib.lines.Line2D at 0x106d81be0>]

In [ ]: