Chemical Reactor

Simulation of a continuous stirred tank reactor (CSTR) with consecutive exothermic reactions.

You can also find this example as a single file in the GitHub repository.

Chemical Reaction System

The reactor contains two consecutive reactions:

\[A \xrightarrow{k_1} B \xrightarrow{k_2} C\]

Both reactions are:

• Exothermic (release heat)

• Temperature-dependent (Arrhenius kinetics)

• Highly nonlinear

Mathematical Model

The system is described by three ODEs:

Concentration of A:

\[\frac{dC_A}{dt} = \frac{C_{A,in} - C_A}{\tau} - k_1 e^{-E_1/(RT)} C_A\]

Concentration of B:

\[\frac{dC_B}{dt} = \frac{-C_B}{\tau} + k_1 e^{-E_1/(RT)} C_A - k_2 e^{-E_2/(RT)} C_B\]

Reactor Temperature:

\[\frac{dT}{dt} = \frac{T_{in} - T}{\tau} + \frac{-\Delta H_1}{\rho C_p} k_1 e^{-E_1/(RT)} C_A + \frac{-\Delta H_2}{\rho C_p} k_2 e^{-E_2/(RT)} C_B - \frac{UA(T-T_c)}{V\rho C_p}\]

Why is this System Stiff?

The system is stiff because:

1. Exponential temperature dependence creates vastly different timescales

2. Fast reactions (high \(k\) values) vs. slow heat transfer

3. Tight coupling between concentration and temperature

Standard explicit solvers would require extremely small timesteps. We use GEAR52A, an implicit solver designed for stiff systems.

This stiff ODE system uses the ODE block with the GEAR52A solver, designed specifically for stiff differential equations.

import numpy as np
import matplotlib.pyplot as plt

# Apply PathSim docs matplotlib style for consistent, theme-friendly figures
plt.style.use('../pathsim_docs.mplstyle')

from pathsim import Simulation, Connection
from pathsim.blocks import ODE, Source, Scope
from pathsim.solvers import GEAR52A

System Parameters

We define the physical and chemical parameters for the reactor.

# Initial conditions
Ca_0 = 1.0    # Initial concentration of A [mol/L]
Cb_0 = 0.0    # Initial concentration of B [mol/L]
T_0 = 300.0   # Initial temperature [K]

# System parameters
Tc = 280.0    # Coolant temperature [K]
tau = 1.0     # Residence time [s]
k1_0 = 1e4    # Pre-exponential factor 1 [1/s]
k2_0 = 1e3    # Pre-exponential factor 2 [1/s]
E1 = 5e4      # Activation energy 1 [J/mol]
E2 = 5.5e4    # Activation energy 2 [J/mol]
dH1 = -5e4    # Reaction enthalpy 1 [J/mol]
dH2 = -5.2e4  # Reaction enthalpy 2 [J/mol]
rho = 1000.0  # Density [kg/m³]
Cp = 4.184    # Heat capacity [J/(g·K)]
U = 1000.0    # Heat transfer coefficient [W/(m²·K)]
V = 0.1       # Reactor volume [m³]
A = 0.1       # Heat transfer area [m²]
R = 8.314     # Gas constant [J/(mol·K)]

Input Signals

The feed concentration and temperature vary with time to simulate disturbances.

# Time-varying inputs
Src_Ca = Source(lambda t: 2.0 + np.sin(0.5*t))  # Feed concentration A
Src_T = Source(lambda t: 280.0 * (1 - 0.8 * np.exp(-0.6*t)))  # Feed temperature

ODE Function

We define the right-hand side of the ODEs, implementing the reactor dynamics with Arrhenius kinetics.

def reaction_rates(x, u, t):
    """CSTR dynamics with consecutive reactions

    Parameters
    ----------
    x : array [Ca, Cb, T]
        State variables
    u : array [Ca_in, T_in]
        Input variables
    t : float
        Time

    Returns
    -------
    dx_dt : array
        Time derivatives
    """
    # Unpack states
    Ca, Cb, T = x

    # Unpack inputs
    Ca_in, T_in = u

    # Reaction rate constants (Arrhenius)
    k1 = k1_0 * np.exp(-E1/(R*T))
    k2 = k2_0 * np.exp(-E2/(R*T))

    # Concentration dynamics
    dCa_dt = (Ca_in - Ca)/tau - k1*Ca
    dCb_dt = -Cb/tau + k1*Ca - k2*Cb

    # Temperature dynamics
    Q_reaction1 = (-dH1/(rho*Cp)) * k1 * Ca
    Q_reaction2 = (-dH2/(rho*Cp)) * k2 * Cb
    Q_cooling = U*A*(T - Tc)/(V*rho*Cp)

    dT_dt = (T_in - T)/tau + Q_reaction1 + Q_reaction2 - Q_cooling

    return np.array([dCa_dt, dCb_dt, dT_dt])

System Setup

We create the ODE block and scope for recording the state variables.

# Define system blocks
CSTR = ODE(reaction_rates, np.array([Ca_0, Cb_0, T_0]))
Sco = Scope(labels=['Ca', 'Cb', 'T'])

blocks = [CSTR, Src_Ca, Src_T, Sco]

connections = [
    Connection(CSTR, Sco),         # Ca output
    Connection(CSTR[1], Sco[1]),   # Cb output
    Connection(CSTR[2], Sco[2]),   # T output
    Connection(Src_Ca, CSTR[0]),   # Ca_in input
    Connection(Src_T, CSTR[1])     # T_in input
]

Simulation with Stiff Solver

We use GEAR52A, an adaptive order adaptive timestep implicit GEAR method optimized for stiff ODEs.

# Initialize simulation with stiff solver
Sim = Simulation(
    blocks,
    connections,
    dt=0.001,
    log=True,
    Solver=GEAR52A,
    tolerance_lte_abs=1e-6,
    tolerance_lte_rel=1e-4
)

# Run simulation for 20 seconds
Sim.run(20)

11:39:18 - INFO - LOGGING (log: True)
11:39:18 - INFO - BLOCKS (total: 4, dynamic: 1, static: 3, eventful: 0)
11:39:18 - INFO - GRAPH (nodes: 4, edges: 5, alg. depth: 1, loop depth: 0, runtime: 0.035ms)
11:39:18 - INFO - STARTING -> TRANSIENT (Duration: 20.00s)
11:39:18 - INFO - --------------------   1% | 0.0s<1.1s | 337.5 it/s
11:39:18 - INFO - ####----------------  20% | 0.0s<0.0s | 1154.6 it/s
11:39:18 - INFO - ########------------  41% | 0.0s<0.0s | 1134.5 it/s
11:39:18 - INFO - ############--------  61% | 0.1s<0.0s | 1132.7 it/s
11:39:18 - INFO - ################----  81% | 0.1s<0.0s | 1117.3 it/s
11:39:18 - INFO - #################### 100% | 0.1s<--:-- | 1239.9 it/s
11:39:18 - INFO - FINISHED -> TRANSIENT (total steps: 59, successful: 52, runtime: 65.00 ms)

{'total_steps': 59, 'successful_steps': 52, 'runtime_ms': 65.0030010001501}

Results: Reactor Dynamics

The plots show:

• Ca (blue): Concentration of reactant A

• Cb (orange): Concentration of intermediate B

• T (green): Reactor temperature

Notice the complex dynamics with multiple timescales and the strong coupling between concentrations and temperature.

# Plot results
Sco.plot(".-", lw=1.5)
plt.show()

Analysis: Phase Portrait

Let’s examine the relationship between concentration and temperature in phase space.

# Get simulation data
time, [Ca, Cb, T] = Sco.read()

# Phase portraits
fig, axes = plt.subplots(1, 2, figsize=(9, 4), dpi=200)

# Ca vs T
axes[0].plot(T, Ca, linewidth=1.5)
axes[0].set_xlabel('Temperature [K]')
axes[0].set_ylabel('Concentration Ca [mol/L]')
axes[0].set_title('Phase Portrait: Ca vs T')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Cb vs T
axes[1].plot(T, Cb, linewidth=1.5, color='orange')
axes[1].set_xlabel('Temperature [K]')
axes[1].set_ylabel('Concentration Cb [mol/L]')
axes[1].set_title('Phase Portrait: Cb vs T')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()