1. Description

This example is intended as an introduction to the nonlinear dynamics of an exothermic continuous stirred-tank reactor. The example has been studied by countless researchers and students since the pioneering work of Amundson and Aris in the 1950's. The particular formulation and parameter values described below are taken from example 2.5 from Seborg, Edgar, Mellichamp and Doyle (SEMD).

Exothermic Reactor

(Diagram By Daniele Pugliesi - Own work, CC BY-SA 3.0, Link)

2. Reaction Kinetics

We assume the kinetics are dominated by a single first order reaction

$$A \xrightarrow{kc_A}{} \text{Products}$$

The reaction rate per unit volume is modeled as the product $kc_A$ where $c_A$ is the concentration of $A$. The rate constant $k(T)$ is a increases with temperature following the Arrehenius law

$$k(t) = k_0 e^{-\frac{E_a}{RT}}$$

$E_a$ is the activation energy, $R$ is the gas constant, $T$ is absolute temperature, and $k_0$ is the pre-exponential factor.

We can see the strong temperature dependence by plotting $k(T)$ versus temperature over typical operating conditions.

In [4]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

Ea  = 72750     # activation energy J/gmol
R   = 8.314     # gas constant J/gmol/K
k0  = 7.2e10    # Arrhenius rate constant 1/min

T = np.linspace(290,400)

# Arrhenius rate expression
def k(T):
    return k0*np.exp(-Ea/R/T)

plt.xlabel('Absolute Temperature [K]')
plt.ylabel('Rate Constant [1/min]')
plt.title('Arrenhius Rate Constant')

This graph shows the reaction rate changes by three orders of magnitude over the range of possible operating temperatures. Because an exothermic reaction releases heat faster at higher temperatures, there is a positive feedback that can potentially result in unstable process behavior.

3. Model Equations and Parameter Values

The model consists of mole and energy balances on the contents of the well-mixed reactor.

\begin{align*} V\frac{dc_A}{dt} & = q(c_{Ai}-c_A)-Vkc_A \\ V\rho C_p\frac{dT}{dt} & = wC_p(T_i-T) + (-\Delta H_R)Vkc_A + UA(T_c-T) \end{align*}

Normalizing the equations to isolate the time rates of change of $c_A$ and $T$ give

\begin{align*} \frac{dc_A}{dt} & = \frac{q}{V}(c_{Ai} - c_A)- kc_A \\ \frac{dT}{dt} & = \frac{q}{V}(T_i - T) + \frac{-\Delta H_R}{\rho C_p}kc_A + \frac{UA}{V\rho C_p}(T_c - T) \end{align*}

which are the equations that will be integrated below.

Quantity Symbol Value Units Comments
Activation Energy $E_a$ 72,750 J/gmol
Arrehnius pre-exponential $k_0$ 7.2 x 1010 1/min
Gas Constant $R$ 8.314 J/gmol/K
Reactor Volume $V$ 100 liters
Density $\rho$ 1000 g/liter
Heat Capacity $C_p$ 0.239 J/g/K
Enthalpy of Reaction $\Delta H_r$ -50,000 J/gmol
Heat Transfer Coefficient $UA$ 50,000 J/min/K
Feed flowrate $q$ 100 liters/min
Feed concentration $c_{A,f}$ 1.0 gmol/liter
Feed temperature $T_f$ 350 K
Initial concentration $c_{A,0}$ 0.5 gmol/liter
Initial temperature $T_0$ 350 K
Coolant temperature $T_c$ 300 K Primary Manipulated Variable
In [5]:
Ea  = 72750     # activation energy J/gmol
R   = 8.314     # gas constant J/gmol/K
k0  = 7.2e10    # Arrhenius rate constant 1/min
V   = 100.0     # Volume [L]
rho = 1000.0    # Density [g/L]
Cp  = 0.239     # Heat capacity [J/g/K]
dHr = -5.0e4    # Enthalpy of reaction [J/mol]
UA  = 5.0e4     # Heat transfer [J/min/K]
q = 100.0       # Flowrate [L/min]
cAi = 1.0       # Inlet feed concentration [mol/L]
Ti  = 350.0     # Inlet feed temperature [K]
cA0 = 0.5;      # Initial concentration [mol/L]
T0  = 350.0;    # Initial temperature [K]
Tc  = 300.0     # Coolant temperature [K]

def deriv(y,t):
    cA,T = y
    dcA = (q/V)*(cAi - cA) - k(T)*cA
    dT = (q/V)*(Ti - T) + (-dHr/rho/Cp)*k(T)*cA + (UA/V/rho/Cp)*(Tc-T)
    return [dcA,dT]

We first define a visualization function that will be reused in later simulations.

In [6]:
# visualization
def plotReactor(t,y):
    plt.xlabel('Time [min]')
    plt.ylabel('Concentration [gmol/liter]')

    plt.xlabel('Time [min]')
    plt.ylabel('Temperature [K]');

4. Transient Behavior

The first simulation assumes a cooling water temperature of 300 K.

In [7]:
from scipy.integrate import odeint

Tc = 300
IC = [cA0,T0]
t = np.linspace(0,10.0,2000)
y = odeint(deriv,IC,t)


5. Effect of Cooling Temperature

The primary means of controlling the reactoris through temperature of the cooling water jacket. The next calculations explore the effect of plus or minus change of 5 K in cooling water temperature on reactor behavior. These simulations reproduce the behavior shown in Example 2.5 SEMD.

In [8]:
Tlist = [295,300,305]

for Tc in Tlist:
    y = odeint(deriv,IC,t)

/Users/jeff/anaconda3/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:106: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  warnings.warn(message, mplDeprecation, stacklevel=1)

6. Interactive Simulation

Executing the following cell provides an interactive tool for exploring the relationship of cooling temperature with reactor behavior. Use it to observe a thermal runaway, sustained osciallations, and low and high conversion steady states.

In [6]:
from ipywidgets import interact

def sim(Tcooling):
    global Tc
    Tc = Tcooling
    y = odeint(deriv,IC,t)

interact(sim,Tcooling = (290.0,310.0));

7. Nullclines

The nullclines of two first-order differential equations are points in the phase plane for which one or the other of the two derivatives are zero.

\begin{align*} V\frac{dc_A}{dt} & = 0 = q(c_{Ai}-c_A)-Vkc_A \\ V\rho C_p\frac{dT}{dt} & = 0 = wC_p(T_i-T) + (-\Delta H_R)Vkc_A + UA(T_c-T) \end{align*}

The intersection of the nullclines correspond to steady states. The relative positions of the nullclines provide valuable insight into the dynamics of a nonlinear system.

In [16]:
# plot nullclines

def plotNullclines():
    T = np.linspace(300.0,460.0,1000)
    plt.plot((q/V)*cAi/((q/V) + k(T)),T,'b:')
    plt.plot(((q/V)*(Ti-T) + (UA/V/rho/Cp)*(Tc-T))/((dHr/rho/Cp)*k(T)),T,'r:')
    plt.legend(['dC/dt = 0','dT/dt = 0'])

8. Phase Plane Analysis

The final analysis is display the simulation in both time and phase plane coordinates. The following cell produces a

In [17]:
def plotPhase(t,y):

def phase(Tcooling=300,cinitial=0.5,Tinitial=350):
    # put slider value into the global Tc
    global Tc
    Tc = Tcooling
    y = odeint(deriv,[cinitial,Tinitial],t)


interact(phase,Tcooling = (290.0,310.0), cinitial=(0,1,.01), Tinitial=(300,400,1))

Suggested Exercises

1. Using the interactive simulation, adjust the cooling temperature over the range of values from 290K to 310K, identify the cooling water temperatures at which you observe a qualitative difference in system behavior. Try to identify the following behaviors:

  • Stable steady state
  • Oscillatory steady state
  • The onset of a thermal runaway condition

2. Repeat the exercise using the phase plane simulation. Describe the underlying behaviors in light of the nullclines.

3. Assume the reactor is start up with a tank of reactant feed with a concentration $c_A = 1$. What is the maximum initial temperature that would avoid thermal runaway during reactor startup?