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/"
We consider the original Seborg et al. model, written as:
Here:
In the original Seborg et al. model, , while is constant, , and is constant.
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;
sr_org = ModelicaSystem("SeborgCSTR.mo","SeborgCSTR.ModSeborgCSTRorg")
pd.DataFrame(sr_org.getQuantities())
sr_org.getSimulationOptions()
sr_org.setSimulationOptions(stopTime=15,stepSize=0.05)
u = sr_org.getInputs()
u
u.keys()
u.values()
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()
u["Ti"]
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")
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)
Steady state model: steady concentration gives
while steady temperature gives
We now plot "removed heat rate" and "generated heat rate" : 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 [C], which is at the intersection of the thick green, dotted lines. If the reactor is operated so that the temperature becomes higher than [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 [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. [C] and ca. [C]) are stable. This analysis does not explain oscillations, but it seems reasonable that oscillations may occur:
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 , the initial temperature leads to reduced temperature until the single equilibrium point at ca. [C] is reached. For the warmer cooling temperature (= less cooling) of , the initial temperature leads to increased temperature. While it may seem like there is a single, stable equilibrium at ca. [C] (and possibly a point which is difficult to interpret at ca. [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.
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)