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¶

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:

\frac{d c_{A}}{d t} = \frac{\dot{V_{i}}}{V} (c_{A, i} - c_{A}) - a \cdot r C_{p} \frac{d T}{d t} = \frac{{\dot{V}}_{i}}{V} C_{p, i} (T_{i} - T) + (- Δ_{r} \tilde{H}) r V + \dot{Q}

$\frac{dc_{\mathrm{A}}}{dt} = \frac{\dot{V_{\mathrm{i}}}}{V}\left(c_{\mathrm{A,i}}-c_{\mathrm{A}}\right)-a\cdot r\\ C_{p}\frac{dT}{dt} = \frac{\dot{V}_{\mathrm{i}}}{V}C_{p,\mathrm{i}}\left(T_{\mathrm{i}}-T\right)+\left(-\Delta_{\mathrm{r}}\tilde{H}\right)rV+\dot{Q}$

Here:

r = k c_{A}^{a} k = k_{0} \exp (- \frac{E}{R T}) \dot{Q} = U A (T_{c} - T)

$r = kc_{\mathrm{A}}^{a} \\ k = k_{0}\exp\left(-\frac{E}{RT}\right) \\ \dot{Q} = {\cal U}A\left(T_{\mathrm{c}}-T\right)$

In the original Seborg et al. model, $a = 1$ , while $\Delta_{\mathrm{r}}$ is constant, $\dot{V}_\mathrm{e} = \dot{V}_\mathrm{i}$ , and $C_{p} = C_{p,\mathrm{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¶

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¶

pd.DataFrame(sr_org.getQuantities())

Setting simulation length¶

sr_org.getSimulationOptions()

{'solver': 'dassl',
 'startTime': 0.0,
 'stepSize': 0.002,
 'stopTime': 1.0,
 'tolerance': 1e-06}

sr_org.setSimulationOptions(stopTime=15,stepSize=0.05)

Checking input¶

u = sr_org.getInputs()

u

{'Tc': None, 'Ti': None, 'Vdi': None, 'cAi': None}

u.keys()

['Vdi', 'cAi', 'Tc', 'Ti']

u.values()

[None, None, None, None]

Setting nominal inputs¶

sr_org.setInputs(Tc=300)
sr_org.setInputs(Ti=350)
sr_org.setInputs(Vdi=100)
sr_org.setInputs(cAi=1)

u = sr_org.getInputs()
u.keys()

['Vdi', 'cAi', 'Tc', 'Ti']

u["Ti"]

[(0.0, 350), (15.0, 350)]

Simulating system for three cases of inputs, and extracting results¶

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¶

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)

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)

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

\frac{d c_{A}}{d t} = 0 ⇓ 0 = \frac{{\dot{V}}_{i}}{V} (c_{A, i} - c_{A}) - k c_{A} ⇓ c_{A} (T) = \frac{{\dot{V}}_{i} / V}{k_{0} \exp (- \frac{E / R}{T}) + {\dot{V}}_{i} / V} c_{A, i}

$\frac{dc_{\mathrm{A}}}{dt} = 0\\ \Downarrow\\ 0 = \frac{\dot{V}_{\mathrm{i}}}{V}\left(c_{\mathrm{A,i}}-c_{\mathrm{A}}\right)-kc_{\mathrm{A}}\\ \Downarrow\\ c_{\mathrm{A}}(T) = \frac{\dot{V}_{\mathrm{i}}/V}{k_{0}\exp\left(-\frac{E/R}{T}\right)+\dot{V}_{\mathrm{i}}/V}c_{\mathrm{A,i}}$

while steady temperature gives

C_{p} \frac{d T}{d t} = 0 ⇓ 0 = \frac{{\dot{V}}_{i}}{V} C_{p} (T_{i} - T) + (- Δ_{r} \tilde{H}) r V + \dot{Q} ⇓ 0 = \underset{= - {\dot{Q}}_{r m} (T)}{\underset{⏟}{\frac{{\dot{V}}_{i}}{V} C_{p} (T_{i} - T) + U A (T_{c} - T)}} + \underset{= {\dot{Q}}_{g e n} (T)}{\underset{⏟}{(- Δ_{r} \tilde{H}) k_{0} \exp (- \frac{E / R}{T}) c_{A} (T) \cdot V}}

$C_{p}\frac{dT}{dt} = 0\\ \Downarrow\\ 0 = \frac{\dot{V}_{\mathrm{i}}}{V}C_p\left(T_{\mathrm{i}}-T\right)+\left(-\Delta_{\mathrm{r}}\tilde{H}\right)rV+\dot{Q}\\ \Downarrow\\ 0 = \underset{=-\dot{Q}_{\mathrm{rm}}\left(T\right)}{\underbrace{\frac{\dot{V}_{\mathrm{i}}}{V}C_p\left(T_{\mathrm{i}}-T\right)+{\cal U}A\left(T_{\mathrm{c}}-T\right)}}+\underset{=\dot{Q}_{\mathrm{gen}}\left(T\right)}{\underbrace{\left(-\Delta_{\mathrm{r}}\tilde{H}\right)k_{0}\exp\left(-\frac{E/R}{T}\right)c_{\mathrm{A}}(T)\cdot V}}$

We now plot "removed heat rate" $\dot{Q}_{\mathrm{rm}}(T)$ and "generated heat rate" $\dot{Q}_{\mathrm{gen}}(T)$ : at equilibrium, they must be equal.

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

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 = 350-273.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 $350-273.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 $350-273.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 ( $\dot{Q}_\mathrm{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.

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 $T_\mathrm{c} - 10$ , the initial temperature leads to reduced temperature until the single equilibrium point at ca. $T_\mathrm{c} = 40$ [C] is reached. For the warmer cooling temperature (= less cooling) of $T_\mathrm{c} + 5$ , the initial temperature leads to increased temperature. While it may seem like there is a single, stable equilibrium at ca. $T_\mathrm{c} = 103$ [C] (and possibly a point which is difficult to interpret at ca. $T_\mathrm{c} = [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¶

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)]

sr_org.setInputs(Tc=uTc)
sr_org.setInputs(Ti=uTi)

sr_org.setSimulationOptions(stopTime=30)
sr_org.simulate()

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)

	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