Use of Modelica + Python in Process Systems Engineering Education

Basic analysis of "Seborg reactor"

Bernt Lie

University College of Southeast Norway

Basic import and definitions

In [24]:
from OMPython import ModelicaSystem
import numpy as np
import numpy.random as nr
%matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
LW1 = 2.5
LW2 = LW1/2
Cb1 = (0.3,0.3,1)
Cb2 = (0.7,0.7,1)
Cg1 = (0,0.6,0)
Cg2 = (0.5,0.8,0.5)
Cr1 = "Red"
Cr2 = (1,0.5,0.5)
LS1 = "solid"
LS2 = "dotted"
LS3 = "dashed"
figpath = "../figs/"

Nonlinear reactor model from Seborg et al.

Process diagram

Original state space model

We consider the original Seborg et al. model, written as:

dcAdt=Vi˙V(cA,icA)arCpdTdt=V˙iVCp,i(TiT)+(ΔrH~)rV+Q˙

Here:

r=kcAak=k0exp(ERT)Q˙=UA(TcT)

In the original Seborg et al. model, a=1, while Δr is constant, V˙e=V˙i, and Cp=Cp,i is constant.

Modelica code for original Seborg et al model "ModSeborgCSTRorg" within package "SeborgCSTR"

mo
 model ModSeborgCSTRorg
    // Model of ORiGinal Seborg CSTR in ode form 
    // author:  Bernt Lie
    //          University of Southeast Norway
    //          November 7, 2017
    //
    // Parameters
    parameter Real V = 100 "Reactor volume, L";
    parameter Real rho = 1e3 "Liquid density, g/L";
    parameter Real a = 1 "Stoichiometric constant, -";
    parameter Real EdR = 8750 "Activation temperature, K";
    parameter Real k0 = exp(EdR/350) "Pre-exponential factor, 1/min";
    parameter Real cph = 0.239 "Specific heat capacity of mixture, J.g-1.K-1";
    parameter Real DrHt = -5e4 "Molar enthalpy of reaction, J/mol";
    parameter Real UA = 5e4 "Heat transfer parameter, J/(min.K)";

    // Initial state parameters
    parameter Real cA0 = 0.5 "Initial concentration of A, mol/L";
    parameter Real T0 = 350 "Initial temperature, K";
    // Declaring variables
    // -- states
    Real cA(start = cA0, fixed = true) "Initializing concentration of A in reactor, mol/L";
    Real T(start = T0, fixed = true) "Initializing temperature in reactor, K";
    // -- auxiliary variables
    Real r "Rate of reaction, mol/(L.s)";
    Real k "Reaction 'constant', ...";
    Real Qd "Heat flow rate, J/min";
    // -- input variables
    input Real Vdi "Volumetric flow rate through reactor, L/min";
    input Real cAi "Influent molar concentration of A, mol/L";
    input Real Ti "Influent temperature, K";
    input Real Tc "Cooling temperature', K";
    // -- output variables
    output Real y_T "Reactor temperature, K";
  // Equations constituting the model
  equation
    // Differential equations
    der(cA) = Vdi*(cAi-cA)/V- a*r;
    der(T) = Vdi*(Ti-T)/V + (-DrHt)*r/(rho*cph) + Qd/(rho*V*cph);
    // Algebraic equations
    r = k*cA^a;
    k = k0*exp(-EdR/T);
    Qd = UA*(Tc-T);
    // Outputs
    y_T = T;
  end ModSeborgCSTRorg;

Nominal model study

Instantiating model

In [25]:
sr_org = ModelicaSystem("SeborgCSTR.mo","SeborgCSTR.ModSeborgCSTRorg")
2018-02-05 16:09:00,687 - OMPython - INFO - OMC Server is up and running at file:///c:/users/bernt_~1/appdata/local/temp/openmodelica.port.a27082498c784186bbda1fbf0a7db93c

Checking quantity names for model

In [3]:
pd.DataFrame(sr_org.getQuantities())
Out[3]:
Changeable Description Name Value Variability alias aliasvariable
0 false Initializing temperature in reactor, K T None continuous noAlias None
1 false Initializing concentration of A in reactor, mol/L cA None continuous noAlias None
2 false der(Initializing temperature in reactor, K) der(T) None continuous noAlias None
3 false der(Initializing concentration of A in reactor... der(cA) None continuous noAlias None
4 false None $cse1 None continuous noAlias None
5 false Heat flow rate, J/min Qd None continuous noAlias None
6 true Cooling temperature', K Tc None continuous noAlias None
7 true Influent temperature, K Ti None continuous noAlias None
8 true Volumetric flow rate through reactor, L/min Vdi None continuous noAlias None
9 true Influent molar concentration of A, mol/L cAi None continuous noAlias None
10 false Reaction 'constant', ... k None continuous noAlias None
11 false Rate of reaction, mol/(L.s) r None continuous noAlias None
12 false Reactor temperature, K y_T None continuous noAlias None
13 true Molar enthalpy of reaction, J/mol DrHt -50000.0 parameter noAlias None
14 true Activation temperature, K EdR 8750.0 parameter noAlias None
15 true Initial temperature, K T0 350.0 parameter noAlias None
16 true Heat transfer parameter, J/(min.K) UA 50000.0 parameter noAlias None
17 true Reactor volume, L V 100.0 parameter noAlias None
18 true Stoichiometric constant, - a 1.0 parameter noAlias None
19 true Initial concentration of A, mol/L cA0 0.5 parameter noAlias None
20 true Specific heat capacity of mixture, J.g-1.K-1 cph 0.239 parameter noAlias None
21 false Pre-exponential factor, 1/min k0 None parameter noAlias None
22 true Liquid density, g/L rho 1000.0 parameter noAlias None

Setting simulation length

In [4]:
sr_org.getSimulationOptions()
Out[4]:
{'solver': 'dassl',
 'startTime': 0.0,
 'stepSize': 0.002,
 'stopTime': 1.0,
 'tolerance': 1e-06}
In [5]:
sr_org.setSimulationOptions(stopTime=15,stepSize=0.05)

Checking input

In [6]:
u = sr_org.getInputs()
In [7]:
u
Out[7]:
{'Tc': None, 'Ti': None, 'Vdi': None, 'cAi': None}
In [8]:
u.keys()
Out[8]:
['Vdi', 'cAi', 'Tc', 'Ti']
In [9]:
u.values()
Out[9]:
[None, None, None, None]

Setting nominal inputs

In [10]:
sr_org.setInputs(Tc=300)
sr_org.setInputs(Ti=350)
sr_org.setInputs(Vdi=100)
sr_org.setInputs(cAi=1)
In [11]:
u = sr_org.getInputs()
u.keys()
Out[11]:
['Vdi', 'cAi', 'Tc', 'Ti']
In [12]:
u["Ti"]
Out[12]:
[(0.0, 350), (15.0, 350)]

Simulating system for three cases of inputs, and extracting results

In [13]:
sr_org.setInputs(Tc=300)
sr_org.simulate()
tm_org0, T_org0, Tc_org0, cA_org0 = sr_org.getSolutions("time","T","Tc","cA")
sr_org.setInputs(Tc=300+5)
sr_org.simulate()
tm_org0p, T_org0p, Tc_org0p, cA_org0p = sr_org.getSolutions("time","T","Tc","cA")
sr_org.setInputs(Tc=300-10)
sr_org.simulate()
tm_org0m, T_org0m, Tc_org0m, cA_org0m = sr_org.getSolutions("time","T","Tc","cA")

Plotting results

In [14]:
plt.plot(tm_org0,T_org0-273.15,LW=LW1,color=Cg1,label=r"$T_{\mathrm{c,org}}^{\ast}$")
plt.plot(tm_org0p,T_org0p-273.15,LW=LW1,color=Cr1,label=r"$T_{\mathrm{c,org}}^{\ast}+5$")
plt.plot(tm_org0m,T_org0m-273.15,LW=LW1,color=Cb1,label=r"$T_{\mathrm{c,org}}^{\ast}-10$")
plt.legend(ncol=3)
plt.title("Reactor temperature")
plt.xlabel(r"time $t$ [min]")
plt.ylabel(r"$T$ [${}^\circ \mathrm{C}$]")
plt.grid()
plt.xlim(0,15)
figfile = "tempSeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)
In [15]:
plt.plot(tm_org0,cA_org0,LW=LW1,color=Cg1,label=r"$T_{\mathrm{c,org}}^{\ast}$")
plt.plot(tm_org0p,cA_org0p,LW=LW1,color=Cr1,label=r"$T_{\mathrm{c,org}}^{\ast}+5$")
plt.plot(tm_org0m,cA_org0m,LW=LW1,color=Cb1,label=r"$T_{\mathrm{c,org}}^{\ast}-10$")
plt.legend(ncol=3,loc=7)
plt.title("Concentration of $\mathrm{A}$")
plt.xlabel(r"time $t$ [min]")
plt.ylabel(r"$c_{\mathrm{A}}$ [mol/L]")
plt.grid()
plt.xlim(0,15)
figfile = "concSeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)
In [16]:
plt.plot(tm_org0,Tc_org0-273.15,LW=LW1, color=Cg1, label=r"$T_\mathrm{c}^{\ast}$")
plt.plot(tm_org0p,Tc_org0p-273.15,LW=LW1,color=Cr1, label=r"$T_\mathrm{c}^{\ast}+5$")
plt.plot(tm_org0m,Tc_org0m-273.15,LW=LW1,color=Cb1, label=r"$T_\mathrm{c}^{\ast}-10$")
plt.title("Cooling liquid temperature")
plt.xlabel(r"time $t$ [min]")
plt.ylabel(r"$T_\mathrm{c}$ [${}^\circ \mathrm{C}$]")
plt.grid()
plt.axis(ymin=15,ymax=35)
plt.xlim(0,15)
plt.legend(ncol=3,loc=7)
figfile = "coolSeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)

Analysis: why oscillations?

Steady state model: steady concentration gives

dcAdt=00=V˙iV(cA,icA)kcAcA(T)=V˙i/Vk0exp(E/RT)+V˙i/VcA,i

while steady temperature gives

CpdTdt=00=V˙iVCp(TiT)+(ΔrH~)rV+Q˙0=V˙iVCp(TiT)+UA(TcT)=Q˙rm(T)+(ΔrH~)k0exp(E/RT)cA(T)V=Q˙gen(T)

We now plot "removed heat rate" Q˙rm(T) and "generated heat rate" Q˙gen(T): at equilibrium, they must be equal.

In [17]:
EdR = 8750
k0 = np.exp(EdR/350)
VdidV = 1
cAi = 1
def cA(T):
    return VdidV/(k0*np.exp(-EdR/T) + VdidV)*cAi
#
rho = 1e3
V = 100
cph = 0.239
Cp = rho*V*cph
DrHt = -5e4
UA = 5e4
Ti = 350
Tc = 300
def Qdrm(T,Tc):
    return -(VdidV*Cp*(Ti-T) + UA*(Tc-T))

def Qdgen(T):
    return (-DrHt)*k0*np.exp(-EdR/T)*cA(T)*V
In [18]:
TT = np.linspace(50,100)
plt.plot(TT, Qdgen(TT+273.15)/1e6,LW=LW1, color=Cr1, label=r"$\dot{Q}_\mathrm{gen}$")
plt.plot(TT, Qdrm(TT+273.15,Tc)/1e6,LW=LW1, color=Cb1, label=r"$\dot{Q}_\mathrm{rm}$")
plt.plot((350-273.15)*np.array([1,1]), [0,5],LW=LW1,color=Cg1,LS=LS2)
plt.plot([50,100], Qdrm(350,Tc)/1e6*np.array([1,1]),LW=LW1,color=Cg1,LS=LS2)
plt.plot((50.8)*np.array([1,1]), [0,5],LW=LW2,color=Cg2,LS=LS2)
plt.plot([50,100], Qdrm(50.8+273.15,Tc)/1e6*np.array([1,1]),LW=LW2,color=Cg2,LS=LS2)
plt.plot((97.2)*np.array([1,1]), [0,5],LW=LW2,color=Cg2,LS=LS2)
plt.plot([50,100], Qdrm(97.2+273.15,Tc)/1e6*np.array([1,1]),LW=LW2,color=Cg2,LS=LS2)
plt.title("Generated vs. removed heat rates")
plt.xlabel(r"temperature $T$ [C]")
plt.ylabel(r"$\dot{Q}$ [MJ/min]")
plt.grid()
plt.xlim(50,100)
plt.ylim(np.min(Qdrm(TT+273.15,Tc)/1e6), np.max(Qdrm(TT+273.15,Tc)/1e6))
plt.legend()
figfile = "equilSeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)

Observe that the system has 3 equilibrium points. For the Seborg reactor, we start the system at T=350273.15 [C], which is at the intersection of the thick green, dotted lines. If the reactor is operated so that the temperature becomes higher than 350273.15 [C], a higher heat rate is generated than removed, and the temperature increases. If, on the other hand, the temperature is operated so that the temperature becomes lower than 350273.15 [C], then the removal heat rate is higher than the generated heat rate, and the reactor is cooled. In summary: the equilibrium point given by the green cross is unstable. Similar arguments leads to the conclusion that the other two equilibrium points (at ca. T=50.8 [C] and ca. T=97.2 [C]) are stable. This analysis does not explain oscillations, but it seems reasonable that oscillations may occur:

  • If the temperature increases, generated heat rate (Q˙gen) will first increase. However, concentration of species A drops and eventually the reaction rate will draw to a halt - and generated heat rate will drop. Thus, the reactor temperature will drop.
  • As the reactor cools down, species A will increase in concentration due to addition of fresh A in the influent. Then, eventually, the reaction rate will "ignite" again, and the generated heat rate will again increase, leading to increasing temperature.
  • This will take place cyclically.
In [19]:
TT = np.linspace(35,110)
plt.plot(TT, Qdgen(TT+273.15)/1e6,LW=LW1, color=Cr1, label=r"$\dot{Q}_\mathrm{gen}$")
plt.plot(TT, Qdrm(TT+273.15,Tc)/1e6,LW=LW2, color="blue", LS = LS2, label=r"$\dot{Q}_\mathrm{rm}$, $T_\mathrm{c}$")
plt.plot(TT, Qdrm(TT+273.15,Tc+5)/1e6,LW=LW1, color=Cb2, label=r"$\dot{Q}_\mathrm{rm}$, $T_\mathrm{c}+5$")
plt.plot(TT, Qdrm(TT+273.15,Tc-10)/1e6,LW=LW1, color=Cb1, label=r"$\dot{Q}_\mathrm{rm}$, $T_\mathrm{c}-10$")
#
plt.plot((350-273.15)*np.array([1,1]), [-1,6],LW=LW1,color=Cg1,LS=LS2)
plt.plot([30,110], Qdrm(350,Tc)/1e6*np.array([1,1]),LW=LW1,color=Cg1,LS=LS2)
plt.plot((50.8)*np.array([1,1]), [-1,6],LW=LW2,color=Cg2,LS=LS2)
plt.plot([30,110], Qdrm(50.8+273.15,Tc)/1e6*np.array([1,1]),LW=LW2,color=Cg2,LS=LS2)
plt.plot((97.2)*np.array([1,1]), [-1,6],LW=LW2,color=Cg2,LS=LS2)
plt.plot([30,110], Qdrm(97.2+273.15,Tc)/1e6*np.array([1,1]),LW=LW2,color=Cg2,LS=LS2)
#
plt.title("Generated vs. removed heat rates")
plt.xlabel(r"temperature $T$ [C]")
plt.ylabel(r"$\dot{Q}$ [MJ/min]")
plt.grid()
plt.xlim(35,110)
plt.ylim(np.min(Qdrm(TT+273.15,Tc+5)/1e6), np.max(Qdrm(TT+273.15,Tc-10)/1e6))
plt.legend()
figfile = "equil_1_SeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)

Observe that for cooling temperature Tc10, the initial temperature leads to reduced temperature until the single equilibrium point at ca. Tc=40 [C] is reached. For the warmer cooling temperature (= less cooling) of Tc+5, the initial temperature leads to increased temperature. While it may seem like there is a single, stable equilibrium at ca. Tc=103 [C] (and possibly a point which is difficult to interpret at ca. Tc=[60,65] [C], it is important to realize that this analysis is purely based on steady state analysis, and that the dynamics may (and: in fact does!) influence the behavior.

Time varying inputs

In [20]:
uTc = [(0,300),(2,300),(2,300-10),(5,300-10),(5,300),(7,300),(7,300+5),(30,300+5)]
uTi = [(0,350),(15,350),(15,340),(20,340),(20,360),(30,360)]
In [21]:
sr_org.setInputs(Tc=uTc)
sr_org.setInputs(Ti=uTi)
In [22]:
sr_org.setSimulationOptions(stopTime=30)
sr_org.simulate()
In [23]:
tm, T, Tc, Ti = sr_org.getSolutions("time","T","Tc","Ti")
plt.plot(tm,Tc-273.15,LW=LW1,color=Cg1,label=r"$T_{\mathrm{c}}$")
plt.plot(tm,Ti-273.15,LW=LW1,color=Cb1,label=r"$T_{\mathrm{i}}$")
plt.plot(tm,T-273.15,LW=LW1,color=Cr1,label=r"$T$")
plt.title("Temperatures, original model")
plt.xlabel(r"time $t$ [min]")
plt.ylabel(r"$T$ [${}^\circ \mathrm{C}$]")
plt.grid()
plt.xlim(0,30)
plt.legend()
figfile = "varInputSeborgCSTRorg.pdf"
plt.savefig(figpath+figfile)