Home » Calculations in Chemical Engineering Blog » Plug-flow reactor – modeling and simulation

# 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 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):

In the above equation Ci 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. Caf 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 Fv (m3/s) and the cross-section area of the reactor S (m2):

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

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:

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

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:

### Second-order irreversible reaction

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

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

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

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

# 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