## 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 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 *C _{i}* is reactant

*i*concentration,

*u*denotes fluid velocity,

*ν*is the stoichiometric coefficient,

_{i}*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

^{3}/s) and the cross section area of the reactor S (m

^{2}):

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.

```
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.

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

The function *plug_flow *describes the differential equation.

```
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.

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

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

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

```
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('')
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:

```
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('')
plt.legend(['numerical'])
plt.plot
```