Main

## Temperature Control of a Stirred Reactor

The objective of this case study is to automatically maintain a temperature in stirred tank reactor with a highly exothermic reaction. This is an example of a highly nonlinear process that is prone to exponential run-away when the temperature rises too quickly.

A reactor is used to convert a hazardous chemical A to an acceptable chemical B in waste stream before entering a nearby lake. This particular reactor is dynamically modeled as a Continuously Stirred Tank Reactor (CSTR) with a simplified kinetic mechanism that describes the conversion of reactant A to product B with an irreversible and exothermic reaction. Because the analyzer for product B is not fast enough for real-time control, it is desired to maintain the temperature at a constant set point that maximizes the consumption of A (highest possible temperature).

The exercise involves creating a dynamic first-order model, obtaining tuning parameters, and adjusting the PID tuning to achieve acceptable performance.

#### Doublet Test for FOPDT Fit

Perform the necessary open loop tests to determine a first order plus dead time (FOPDT) model that describes the relationship between cooling jacket temperature (MV=Tc) and reactor temperature (CV=T). Use the file data.txt to extract the information for the FOPDT model. Obtain an FOPDT model and report the values of the three parameters K_p, \tau_p, \theta_p.

Below is a step test example in Python. Modify this to start at 300K for the cooling jacket temperature, step up to 303K, down to 297K, and back to 300K in a doublet test. Fit K_p, \tau_p, \theta_p of an FOPDT model using optimization to achieve the best match to the simulated CSTR data.

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

# define CSTR model
def cstr(x,t,u,Tf,Caf):
# Inputs (3):
# Temperature of cooling jacket (K)
Tc = u
# Tf = Feed Temperature (K)
# Caf = Feed Concentration (mol/m^3)

# States (2):
# Concentration of A in CSTR (mol/m^3)
Ca = x[0]
# Temperature in CSTR (K)
T = x[1]

# Parameters:
# Volumetric Flowrate (m^3/sec)
q = 100
# Volume of CSTR (m^3)
V = 100
# Density of A-B Mixture (kg/m^3)
rho = 1000
# Heat capacity of A-B Mixture (J/kg-K)
Cp = 0.239
# Heat of reaction for A->B (J/mol)
mdelH = 5e4
# E - Activation energy in the Arrhenius Equation (J/mol)
# R - Universal Gas Constant = 8.31451 J/mol-K
EoverR = 8750
# Pre-exponential factor (1/sec)
k0 = 7.2e10
# U - Overall Heat Transfer Coefficient (W/m^2-K)
# A - Area - this value is specific for the U calculation (m^2)
UA = 5e4
# reaction rate
rA = k0*np.exp(-EoverR/T)*Ca

# Calculate concentration derivative
dCadt = q/V*(Caf - Ca) - rA
# Calculate temperature derivative
dTdt = q/V*(Tf - T) \
+ mdelH/(rho*Cp)*rA \
+ UA/V/rho/Cp*(Tc-T)

# Return xdot:
xdot = np.zeros(2)
xdot[1] = dTdt
return xdot

# Steady State Initial Conditions for the States
Ca_ss = 0.87725294608097
T_ss = 324.475443431599
x0 = np.empty(2)
x0[0] = Ca_ss
x0[1] = T_ss

u_ss = 300.0
# Feed Temperature (K)
Tf = 350
# Feed Concentration (mol/m^3)
Caf = 1

# Time Interval (min)
t = np.linspace(0,10,100)

# Store results for plotting
Ca = np.ones(len(t)) * Ca_ss
T = np.ones(len(t)) * T_ss
u = np.ones(len(t)) * u_ss

# Step cooling temperature to 295
u[10:] = 295.0

# Simulate CSTR
for i in range(len(t)-1):
ts = [t[i],t[i+1]]
y = odeint(cstr,x0,ts,args=(u[i+1],Tf,Caf))
Ca[i+1] = y[-1][0]
T[i+1] = y[-1][1]
x0[0] = Ca[i+1]
x0[1] = T[i+1]

# Construct results and save data file
# Column 1 = time
# Column 2 = cooling temperature
# Column 3 = reactor temperature
data = np.vstack((t,u,T)) # vertical stack
data = data.T             # transpose data
np.savetxt('data.txt',data,delimiter=',')

# Plot the results
plt.figure()
plt.subplot(3,1,1)
plt.plot(t,u,'b--',linewidth=3)
plt.ylabel('Cooling T (K)')
plt.legend(['Jacket Temperature'],loc='best')

plt.subplot(3,1,2)
plt.plot(t,Ca,'r-',linewidth=3)
plt.ylabel('Ca (mol/L)')
plt.legend(['Reactor Concentration'],loc='best')

plt.subplot(3,1,3)
plt.plot(t,T,'k.-',linewidth=3)
plt.ylabel('T (K)')
plt.xlabel('Time (min)')
plt.legend(['Reactor Temperature'],loc='best')

plt.show()

#### Design PI Controller

Use the K_p, \tau_p, \theta_p from the FOPDT model in the Internal Model Control (IMC) PI tuning correlation to obtain starting values of K_c and \tau_I for a PI controller. Use an aggressive value for the closed loop time constant, \tau_c = 0.1 \tau_p.

#### Implement PI Controller

Using K_c and \tau_I from the IMC tuning correlation, implement a PI controller with anti-reset windup. Test the set point tracking capability of this controller by plotting the response of the process to steps in set point of the reactor temperature (not cooling jacket temperature) from 300K up to 320K and then down to 280K. Limit the cooling jacket temperature between 250K and 350K. Comment on how the nonlinear behavior of this process impacts your observed set point response performance.

#### Tune PI Controller

Improve the tuning by adjusting K_c and \tau_I by trial and error until the controller displays a 10% to 15% overshoot in response to a reactor temperature set point step from 300K to 320K. Plot this set point step response.

#### Challenge: Lowest Concentration

Step up the set point to achieve a maximum temperature in the reactor without exceeding the maximum allowable temperature of 400K (don’t cause a reactor run-away). What is the lowest concentration that can be achieved without exceeding the maximum allowable temperature? This is a very difficult problem due to the potential for reactor temperature run-away. Derivative action may be necessary to achieve the best performance.