Plug-flow reactor - modeling and simulation

Plug-flow reactor model

Plug-flow reactor (PFR, sometimes called piston-flow reactor or tubular reactor) is an idealized reactor in which all particles in a given cross-section have identical velocity and direction of motion. In PFR there is no mixing or back flow, so the fluid flows like a plug, from inlet to outlet (Figure below).

Plug-flow reactor scheme — Fig. 1. Scheme of a plug flow reactor

Plug flow reactor model is formulated based on mass balance and heat balance in a differential volume of a fluid (Fig. 1). If we assume that the process is isothermal (heat of reaction is negligible), then only mass balance is considered. We will assume steady state conditions, which mean that concentrations of reactant do not change over time. It is a typical way of operation of plug flow reactors.

The mathematical model of PFR can be written as follows (Levenspiel, 1998):

$u\frac{dC_{i}}{dx}=v_{i}}r$

$C_{i}(0)=C_{if}$

$0\le x\le L$

In the above equation C_i is reactant i concentration, u denotes fluid velocity, ν_i is the stoichiometric coefficient, r is the reaction rate and x is the position in the reactor. C_a_f is reactant A concentration at the reactor inlet (in the feed stream) and L is the reactor length.

The fluid velocity u is calculated based on the volumetric flow rate F_v (m³/s) and the cross-section area of the reactor S (m²):

$u=F_v/S$

In an ideal plug flow reactor, all the fluid particles have been inside the reactor for exactly the same amount of time (Fogler, 2010). We call this value a mean residence and calculate as

$\tau=L/u$

The residence time data is commonly used in chemical reactors engineering to make predictions of conversion and exit concentrations.

First-order irreversible reaction

Let’s consider simple decomposition reaction:

$A\to 2B$

When the reaction is irreversible and first-order, then we have:

$u\frac{dC_{a}}{dx}=-kC_{a}$

where k is a kinetic constant. In general, kinetic constant depends on temperature. Commonly an Arrhenius equation is used to describe this relationship. Here, we assume isothermal conditions, so will not use this dependency.

For first-order irreversible reaction, the model can be solved analytically. The solution is following:

$C_{a}=C_{af}exp(-x*k/u)$

Second-order irreversible reaction

As an example of a second-order irreversible reaction, let’s use the following one:

$2A\to B$

When the reaction is irreversible and second-order, then we have:

$u\frac{dC_{a}}{dx}=-2k*(C_{a})^2$

Online simulations in chemical engineering

Get started for FREE

Simulation of plug flow reactor

First-order reaction

First, let’s import necessary libraries: numpy for matrix operations, scipy for solving differential equations and matplotlib for making a chart.

Python

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

Next, we will define reaction function. In this example we use first-order irreversible reaction.

Python

def reaction(ca, k):
  r = k * ca
  return r

The function plug_flow describes the differential equation.

Python

def plug_flow(ca, x, k, ni, u):
  r = reaction(ca, k)
  dcadx = ni * r / u
  return dcadx

Snippet below presents the definition of model parameters.

Python

# kinetic constant
k = 0.1
# reactor length
L = 10
x = np.linspace(0, L, 101)
ca0 = 0.6
ni = -1
u = 0.5

Next, we will use odeint function from scipy library to solve the differential equation.

Python

sol = odeint(plug_flow, ca0, x, args=(k, ni, u))

Finally, we will plot the results together with the analytical solution for the comparison.

Python

plt.plot(x, sol)
plt.plot(x[::4], [ca0 * np.exp(-x*k/u) for x in x[::4]],'o')
plt.xlabel('x')
plt.ylabel('$C_{a}$')
plt.legend(['numerical', 'analytical'])
plt.plot

Concentration distribution of substrate along the plug-flow reactor length (first-order reaction) — Fig. 2. Distribution of A concentration along reactor length for first-order reaction – comparison of numerical and analytical solution

Second-order reaction

Obtaining the solution is very similar. We need to change the reaction function and coefficient ν. The final code looks as follows:

Python

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
def reaction(ca, k, n):
  r = k * ca ** n
  return r

def plug_flow(ca, x, k, ni, u, n):
  r = reaction(ca, k, n)
  dcadx = ni * r / u
  return dcadx

# kinetic constant
k = 0.1
# reactor length
L = 10
x = np.linspace(0, L, 101)
# inlet A concentration
ca0 = 0.6
# reaction order
n = 2
# stoichiometric coefficient
ni = -2
# fluid velocity
u = 0.5

sol = odeint(plug_flow, ca0, x, args=(k, ni, u, n))

plt.plot(x, sol)
plt.xlabel('x, m')
plt.ylabel('$c_a, mol/m^3$')
plt.legend(['numerical'])
plt.plot