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/"
The original model is based on several simplifying assumptions:
We maintain that the mixture density is constant, with constant volumetric flow rate through the reactor. However,
Complex model 1 will be given descriptive name related to , which indicates constant density (). The following is a DAE implementation in Modelica.
mo
model ModSeborgCSTRrho
// Model of Seborg CSTR when assuming dominant solvent and Constant flow rate Vd
// 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 m = rho*V "Reactor content mass, g";
parameter Real mS = m "Reactor solvent mass, g";
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 UA = 5e4 "Heat transfer parameter, J/(min.K)";
//
parameter Real HtA_o = 5e4 "Molar enthalpy of A at std state, J/mol";
parameter Real HtB_o = 0 "Molar enthalpy of B at std state, J/mol";
parameter Real HhS_o = 0 "Specific enthalpy of solvent at std state, J/g";
parameter Real cptA = 5 "Molar heat capacity of A, J/(mol.K)";
parameter Real cptB = cptA "Molar heat capacity of B, J/(mol.K)";
parameter Real cphS = 0.239 "Specific heat capacity of solvent, J/(g.K)";
parameter Real cph = cphS "Specific heat capacity of mixture, J/(g.K)";
parameter Real T_o = 293.15 "Temperature in std state, K";
parameter Real p_o = 1.01e5 "Pressure in std state, Pa";
// Initial state parameters
parameter Real cA0 = 0.5 "Initial concentration of A, mol/L";
parameter Real nA0 = cA0*V "Initial number of moles of A, mol";
parameter Real nB0 = 0 "Initial number of moles of B, mol";
parameter Real T0 = 350 "Initial temperature, K";
parameter Real HtA0 = HtA_o + cptA*(T0-T_o) "Initial pure molar enthalpy of A, J/mol";
parameter Real HtB0 = HtB_o + cptB*(T0-T_o) "Initial pure molar enthalpy of B, J/mol";
parameter Real HhS0 = HhS_o + cphS*(T0-T_o) "Initial pure specific enthalpy of solvent, J/g";
parameter Real HA0 = nA0*HtA0 "Initial pure enthalpy of A, J";
parameter Real HB0 = nB0*HtB0 "Initial pure enthalpy of B, J";
parameter Real HS0 = mS*HhS0 "Initial pure enthalpy of solvent, J";
parameter Real H0 = HA0 + HB0 + HS0 "Initial total enthalpy of ideal solution, J";
parameter Real U0 = H0 - p_o*V "Initial internal energy, J";
// Declaring variables
// -- states
Real nA(start = nA0, fixed = true) "Initializing amount of A in reactor, mol";
Real nB(start = nB0, fixed = true) "Initializing amount of B in reactor, mol";
Real U(start = U0, fixed = true) "Initializing internal energy in reactor, J";
// -- auxiliary variables
Real ndAi "Influent molar flow rate of A, mol/s";
Real ndAe "Effluent molar flow rate of A, mol/s";
Real ndAg "Molar rate of generation of A, mol/s";
Real ndBi "Influent molar flow rate of B, mol/s";
Real ndBe "Effluent molar flow rate of B, mol/s";
Real ndBg "Molar rate of generation of B, mol/s";
Real mdSi "Influent mass flow rate of solvent, g/min";
Real mdSe "Effluent mass flow rate of solvent, g/min";
Real Hdi "Influent enthalpy flow rate, J/min";
Real Hde "Effluent enthalpy flow rate, J/min";
Real cA "Molar concentration of A, mol/L";
Real cB "Molar concentration of B, mol/L";
Real r "Rate of reaction, mol/(L.s)";
Real T "Reactor temperature, K";
Real k "Reaction 'constant', ...";
Real H "Reactor enthalpy, J";
Real HA "Enthalpy of pure A, J";
Real HB "Enthalpy of pure B, J";
Real HS "Enthalpy of pure solvent, J";
Real HtA "Molar enthalpy of pure A, J/mol";
Real HtB "Molar enthalpy of pure B, J/mol";
Real HhS "Specific enthalpy of pure solvent, J/g";
Real HdAi "Influent enthalpy flow of pure A,J/min";
Real HdBi "Influent enthalpy flow of pure B, J/min";
Real HdSi "Influent enthalpy flow of pure solvent, J/min";
Real HtAi "Influent molar enthalpy of pure A, J/mol";
Real HtBi "Influent molar enthalpy of pure B, J/mol";
Real HhSi "Influent specific enthalpy of solvent, J/g";
Real HdAe "Effluent enthalpy flow of pure A, J/min";
Real HdBe "Effluent enthalpy flow of pure B, J/min";
Real HdSe "Effluent enthalpy flow of pure solvent, J/min";
Real Qd "Heat flow rate, J/min";
// -- input variables
input Real Vdi "Volumetric flow rate through reactor, L/s";
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";
output Real y_cB "Molar concentration of B, mol/L";
// Equations constituting the model
equation
// Differential equations
der(nA) = ndAi - ndAe + ndAg;
der(nB) = ndBi - ndBe + ndBg;
der(U) = Hdi - Hde + Qd;
// Algebraic equations
nA = cA*V;
nB = cB*V;
ndAi = cAi*Vdi;
ndBi = 0;
mdSi = rho*Vdi;
ndAe = cA*Vdi;
ndBe = cB*Vdi;
mdSe = mdSi;
ndAg = -a*r*V;
ndBg = r*V;
r = k*cA^a;
k = k0*exp(-EdR/T);
U = H-p_o*V;
H = HA + HB + HS;
HA = nA*HtA;
HB = nB*HtB;
HS = mS*HhS;
HtA = HtA_o + cptA*(T-T_o);
HtB = HtB_o + cptB*(T-T_o);
HhS = HhS_o + cphS*(T-T_o);
Hdi = HdAi + HdBi + HdSi;
HdAi = ndAi*HtAi;
HdBi = ndBi*HtBi;
HdSi = mdSi*HhSi;
HtAi = HtA_o + cptA*(Ti-T_o);
HtBi = 0;
HhSi = HhS_o + cphS*(Ti-T_o);
Hde = HdAe + HdBe + HdSe;
HdAe = ndAe*HtA;
HdBe = ndBe*HtB;
HdSe = mdSe*HhS;
Qd = UA*(Tc-T);
// Outputs
y_T = T;
y_cB = cB;
end ModSeborgCSTRrho;
In a more general model, we only assume ideal solution in the reaction mixture. The model is given a descriptive name related to .
For the ideal solution model:
It is relatively straightforward to develop a DAE model for this case. It is also possible to develop an ODE model: this time, we need three states (we need to include information about species B). The ODE formulation is simpler than the DAE model, but it is relatively complicated to develop the ODE model. The DAE model is as follows:
mo
model ModSeborgCSTRis
// Model of Seborg CSTR when assuming Ideal Solution, only
// author: Bernt Lie
// University of Southeast Norway
// November 7, 2017
//
// Parameters
parameter Real rhoS_o = 1e3 "Density of pure S, g/L";
parameter Real rhoA_o = 1.5e3 "Density of pure A, g/L";
parameter Real rhoB_o = 2.5e3 "Density of pure B, g/L";
parameter Real MA = 50 "Molar mass of A, g/mol";
parameter Real MB = MA*a "Molar mass of B, g/mol";
parameter Real V = 100 "Reactor volume, L";
parameter Real a = 1 "Stoichiometric constant, -";
parameter Real k0 = exp(EdR/350) "Pre-exponential factor, 1/min";
parameter Real EdR = 8750 "Activation temperature, K";
parameter Real p_o = 1.01e5 "Pressure in std state, Pa";
//
parameter Real HhS_o = 0 "Specific enthalpy of solvent at std state, J/g";
parameter Real HtA_o = 5e4 "Molar enthalpy of A at std state, J/mol";
parameter Real HtB_o = 0 "Molar enthalpy of B at std state, J/mol";
//
parameter Real cphS = 0.239 "Specific heat capacity of solvent, J/(g.K)";
parameter Real cptA = 5 "Molar heat capacity of A, J/(mol.K)";
parameter Real cptB = cptA "Molar heat capacity of B, J/(mol.K)";
//
parameter Real T_o = 293.15 "Temperature in std state, K";
//
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 nA0 = cA0*V "Initial number of moles of A, mol";
parameter Real nB0 = 0 "Initial number of moles of B, mol";
parameter Real mA0 = nA0*MA "Initial mass of A, g";
parameter Real mB0 = nB0*MB "Initial mass of B, g";
parameter Real VA0 = mA0/rhoA_o "Initial volume of A, L";
parameter Real VB0 = mB0/rhoB_o "Initial volume of B, L";
parameter Real VS0 = V - VA0 - VB0 "Initial volume of S, L";
parameter Real mS0 = VS0*rhoS_o "Initial mass of S, g";
parameter Real T0 = 350 "Initial temperature, K";
parameter Real HhS0 = HhS_o + cphS*(T0-T_o) "Initial pure specific enthalpy of solvent, J/g";
parameter Real HtA0 = HtA_o + cptA*(T0-T_o) "Initial pure molar enthalpy of A, J/mol";
parameter Real HtB0 = HtB_o + cptB*(T0-T_o) "Initial pure molar enthalpy of B, J/mol";
parameter Real HS0 = mS0*HhS0 "Initial pure enthalpy of solvent, J";
parameter Real HA0 = nA0*HtA0 "Initial pure enthalpy of A, J";
parameter Real HB0 = nB0*HtB0 "Initial pure enthalpy of B, J";
parameter Real H0 = HS0 + HA0 + HB0 "Initial total enthalpy of ideal solution, J";
parameter Real U0 = H0 - p_o*V "Initial internal energy, J";
// Declaring variables
// -- states
Real mS(start = mS0) "Initial amount of S in reactor, kg";
Real nA(start = nA0, fixed = true) "Initializing amount of A in reactor, mol";
Real nB(start = nB0, fixed = true) "Initializing amount of B in reactor, mol";
Real U(start = U0, fixed = true) "Initializing internal energy in reactor, J";
// -- auxiliary variables
// Real mS "Mass of S, g";
Real V_x;
Real mdSi "Influent mass flow rate of S, g/s";
Real mdSe "Effluent mass flow rate of S, g/s";
Real ndAi "Influent molar flow rate of A, mol/s";
Real ndAe "Effluent molar flow rate of A, mol/s";
Real ndAg "Molar rate of generation of A, mol/s";
Real ndBi "Influent molar flow rate of B, mol/s";
Real ndBe "Effluent molar flow rate of B, mol/s";
Real ndBg "Molar rate of generation of B, mol/s";
Real Hdi "Influent enthalpy flow rate, J/min";
Real Hde "Effluent enthalpy flow rate, J/min";
Real Qd "Heat flow rate, J/min";
//
Real VS "Volume of pure S, L";
Real VA "Volume of pure A, L";
Real VB "Volume of pure B, L";
Real mA "Mass of A, g";
Real mB "Mass of B, g";
//
Real Vde "Effluent volumetric flow rate, L/min";
//
Real rhoS "Specific concentration of S, g/L";
Real cA "Molar concentration of A, mol/L";
Real cB "Molar concentration of B, mol/L";
Real r "Rate of reaction, mol/(L.s)";
Real k "Reaction 'constant', ...";
Real T "Absolute temperature, K";
//
Real H "Reactor enthalpy, J";
Real HS "Enthalpy of pure solvent, J";
Real HA "Enthalpy of pure A, J";
Real HB "Enthalpy of pure B, J";
Real HhS "Specific enthalpy of pure solvent, J/g";
Real HtA "Molar enthalpy of pure A, J/mol";
Real HtB "Molar enthalpy of pure B, J/mol";
Real HdSi "Influent enthalpy flow of pure solvent, J/min";
Real HdAi "Influent enthalpy flow of pure A,J/min";
Real HdBi "Influent enthalpy flow of pure B, J/min";
Real HhSi "Influent specific enthalpy of solvent, J/g";
Real HtAi "Influent molar enthalpy of pure A, J/mol";
Real HtBi "Influent molar enthalpy of pure B, J/mol";
Real HdSe "Effluent enthalpy flow of pure solvent, J/min";
Real HdAe "Effluent enthalpy flow of pure A, J/min";
Real HdBe "Effluent enthalpy flow of pure B, J/min";
//
Real mdSi_x;
// -- input variables
input Real Vdi "Volumetric flow rate through reactor, L/s";
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";
output Real y_cB "Molar concentration of B, mol/L";
// Equations constituting the model
equation
// Differential equations
mdSi_x = mdSe + der(mS);
der(mS) = mdSi - mdSe;
der(nA) = ndAi - ndAe + ndAg;
der(nB) = ndBi - ndBe + ndBg;
der(U) = Hdi - Hde + Qd;
// Algebraic equations
V = VS + VA + VB;
VS = mS/rhoS_o;
VA = mA/rhoA_o;
VB = mB/rhoB_o;
mS = rhoS*V;
V_x = mS/rhoS_o + nA*MA/rhoA_o + nB*MB/rhoB_o;
mA = nA*MA;
mB = nB*MB;
nA = cA*V;
nB = cB*V;
//
mdSi = rhoS_o*(1-cAi*MA/rhoA_o)*Vdi;
ndAi = cAi*Vdi;
ndBi = 0;
mdSe = rhoS*Vde;
ndAe = cA*Vde;
ndBe = cB*Vde;
ndAg = -a*r*V;
ndBg = r*V;
r = k*cA^a;
k = k0*exp(-EdR/T);
//
U = H-p_o*V;
H = HS + HA + HB;
HS = mS*HhS;
HA = nA*HtA;
HB = nB*HtB;
//
HhS = HhS_o + cphS*(T-T_o);
HtA = HtA_o + cptA*(T-T_o);
HtB = HtB_o + cptB*(T-T_o);
//
Hdi = HdSi + HdAi + HdBi;
HdSi = mdSi*HhSi;
HdAi = ndAi*HtAi;
HdBi = ndBi*HtBi;
//
HhSi = HhS_o + cphS*(Ti-T_o);
HtAi = HtA_o + cptA*(Ti-T_o);
HtBi = 0;
//
Hde = HdAe + HdBe + HdSe;
HdSe = mdSe*HhS;
HdAe = ndAe*HtA;
HdBe = ndBe*HtB;
//
Qd = UA*(Tc-T);
// Outputs
y_T = T;
y_cB = cB;
end ModSeborgCSTRis;
sr_is = ModelicaSystem("SeborgCSTR.mo","SeborgCSTR.ModSeborgCSTRis")
sr_rho = ModelicaSystem("SeborgCSTR.mo","SeborgCSTR.ModSeborgCSTRrho")
sr_org = ModelicaSystem("SeborgCSTR.mo","SeborgCSTR.ModSeborgCSTRorg")
sr_is.setSimulationOptions(stopTime=15,stepSize=0.05)
sr_rho.setSimulationOptions(stopTime=15,stepSize=0.05)
sr_org.setSimulationOptions(stopTime=15,stepSize=0.05)
#
sr_is.setInputs(Tc=300)
sr_is.setInputs(Ti=350)
sr_is.setInputs(Vdi=100)
sr_is.setInputs(cAi=1)
sr_rho.setInputs(Tc=300)
sr_rho.setInputs(Ti=350)
sr_rho.setInputs(Vdi=100)
sr_rho.setInputs(cAi=1)
sr_org.setInputs(Tc=300)
sr_org.setInputs(Ti=350)
sr_org.setInputs(Vdi=100)
sr_org.setInputs(cAi=1)
sr_is.setInputs(Tc=300)
sr_is.simulate()
tm_is0, T_is0, Tc_is0, cA_is0 = sr_is.getSolutions("time","T","Tc","cA")
sr_is.setInputs(Tc=300+5)
sr_is.simulate()
tm_is0p, T_is0p, Tc_is0p, cA_is0p = sr_is.getSolutions("time","T","Tc","cA")
sr_is.setInputs(Tc=300-10)
sr_is.simulate()
tm_is0m, T_is0m, Tc_is0m, cA_is0m = sr_is.getSolutions("time","T","Tc","cA")
sr_rho.setInputs(Tc=300)
sr_rho.simulate()
tm_rho0, T_rho0, Tc_rho0, cA_rho0 = sr_rho.getSolutions("time","T","Tc","cA")
sr_rho.setInputs(Tc=300+5)
sr_rho.simulate()
tm_rho0p, T_rho0p, Tc_rho0p, cA_rho0p = sr_rho.getSolutions("time","T","Tc","cA")
sr_rho.setInputs(Tc=300-10)
sr_rho.simulate()
tm_rho0m, T_rho0m, Tc_rho0m, cA_rho0m = sr_rho.getSolutions("time","T","Tc","cA")
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.plot(tm_rho0,T_rho0-273.15,LW=LW1,LS=LS2,color=Cg2,label=r"$T_{\mathrm{c,\rho}}^{\ast}$")
plt.plot(tm_rho0p,T_rho0p-273.15,LW=LW1,LS=LS2,color=Cr2,label=r"$T_{\mathrm{c,\rho}}^{\ast}+5$")
plt.plot(tm_rho0m,T_rho0m-273.15,LW=LW1,LS=LS2,color=Cb2,label=r"$T_{\mathrm{c,\rho}}^{\ast}-10$")
plt.plot(tm_is0,T_is0-273.15,LW=LW1,LS=LS3,color=Cg2,label=r"$T_{\mathrm{c,is}}^{\ast}$")
plt.plot(tm_is0p,T_is0p-273.15,LW=LW1,LS=LS3,color=Cr2,label=r"$T_{\mathrm{c,is}}^{\ast}+5$")
plt.plot(tm_is0m,T_is0m-273.15,LW=LW1,LS=LS3,color=Cb2,label=r"$T_{\mathrm{c,is}}^{\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 = "tempSeborgCSTR.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.plot(tm_rho0,cA_rho0,LW=LW1,LS=LS2,color=Cg2,label=r"$T_{\mathrm{c,\rho}}^{\ast}$")
plt.plot(tm_rho0p,cA_rho0p,LW=LW1,LS=LS2,color=Cr2,label=r"$T_{\mathrm{c,\rho}}^{\ast}+5$")
plt.plot(tm_rho0m,cA_rho0m,LW=LW1,LS=LS2,color=Cb2,label=r"$T_{\mathrm{c,\rho}}^{\ast}-10$")
plt.plot(tm_is0,cA_is0,LW=LW1,LS=LS3,color=Cg2,label=r"$T_{\mathrm{c,is}}^{\ast}$")
plt.plot(tm_is0p,cA_is0p,LW=LW1,LS=LS3,color=Cr2,label=r"$T_{\mathrm{c,is}}^{\ast}+5$")
plt.plot(tm_is0m,cA_is0m,LW=LW1,LS=LS3,color=Cb2,label=r"$T_{\mathrm{c,is}}^{\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 = "concSeborgCSTR.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$ [s]")
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 = "coolSeborgCSTR.pdf"
plt.savefig(figpath+figfile)