#!/usr/bin/env python # -*- coding: utf-8 -*- """******************************************************************************** tutorial_che_1.py DAE Tools: pyDAE module, www.daetools.com Copyright (C) Dragan Nikolic *********************************************************************************** DAE Tools is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License version 3 as published by the Free Software Foundation. DAE Tools is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with the DAE Tools software; if not, see <http://www.gnu.org/licenses/>. ********************************************************************************""" __doc__ = """ Continuously Stirred Tank Reactor with energy balance and Van de Vusse reactions: .. code-block:: none A -> B -> C 2A -> D Reference: G.A. Ridlehoover, R.C. Seagrave. Optimization of Van de Vusse Reaction Kinetics Using Semibatch Reactor Operation, Ind. Eng. Chem. Fundamen. 1973;12(4):444-447. `doi:10.1021/i160048a700 <https://doi.org/10.1021/i160048a700>`_ The concentrations plot: .. image:: _static/tutorial_che_1-results.png :width: 500px The temperatures plot: .. image:: _static/tutorial_che_1-results2.png :width: 500px """ import sys from daetools.pyDAE import * from time import localtime, strftime # Standard variable types are defined in variable_types.py from pyUnits import m, kg, s, K, Pa, mol, J, W, kJ, hour, l K_t = daeVariableType("K_t", s**(-1), 0, 1E20, 0, 1e-5) K2_t = daeVariableType("K2_t", m**3/(mol*s), 0, 1E20, 0, 1e-5) class modTutorial(daeModel): def __init__(self, Name, Parent = None, Description = ""): daeModel.__init__(self, Name, Parent, Description) # Parameters self.k10 = daeParameter("k_10", s**(-1), self, "A->B pre-exponential factor") self.k20 = daeParameter("k_20", s**(-1), self, "B->C pre-exponential factor") self.k30 = daeParameter("k_30", m**3/(mol*s), self, "2A->D pre-exponential factor") self.E1 = daeParameter("E_1", K, self, "A->B activation energy") self.E2 = daeParameter("E_2", K, self, "B->C activation energy") self.E3 = daeParameter("E_3", K, self, "2A->D activation energy") self.dHr1 = daeParameter("dHr1", J/mol, self, "A->B heat of reaction") self.dHr2 = daeParameter("dHr2", J/mol, self, "B->C heat of reaction") self.dHr3 = daeParameter("dHr3", J/mol, self, "2A->D heat of reaction") self.rho = daeParameter("ρ", kg/(m**3), self, "Density") self.cp = daeParameter("c_p", J/(kg*K), self, "Heat capacity of reactants") self.kw = daeParameter("k_w", J/(K*s*(m**2)), self, "Heat transfer coefficient") self.AR = daeParameter("A_r", m**2, self, "Area of jacket cooling") self.mK = daeParameter("m_K", kg, self, "Mass of the cooling fluid") self.cpK = daeParameter("c_pk", J/(kg*K), self, "Heat capacity of the cooling fluid") # Degrees of freedom (for optimisation) self.VR = daeVariable("V_r", volume_t, self, "Reactor volume") self.Qk = daeVariable("Q_k", power_t, self, "Jacket cooling rate") self.Qf = daeVariable("Qf", volume_flowrate_t, self, "Feed flowrate") self.Caf = daeVariable("Caf", molar_concentration_t, self, "Feed concentration") self.Tf = daeVariable("Tf", temperature_t, self, "Feed temperature") # Variables self.Ca = daeVariable("Ca", molar_concentration_t, self, "Concentration of A") self.Cb = daeVariable("Cb", molar_concentration_t, self, "Concentration of B") self.Cc = daeVariable("Cc", molar_concentration_t, self, "Concentration of C") self.Cd = daeVariable("Cd", molar_concentration_t, self, "Concentration of D") self.T = daeVariable("T", temperature_t, self, "Temperature in the reactor") self.Tk = daeVariable("T_k", temperature_t, self, "Temperature of the cooling jacket") self.k1 = daeVariable("k_1", K_t, self, "Reaction A->B rate constant") self.k2 = daeVariable("k_2", K_t, self, "Reaction B->C rate constant") self.k3 = daeVariable("k_3", K2_t, self, "Reaction 2A->D rate constant") self.V0 = daeVariable("V0", volume_t, self, "Minimum reactor volume") self.Vm = daeVariable("Vm", volume_t, self, "Maximum reactor volume") self.V0_star = daeVariable("V0_star", no_t, self, "Dimensionless minimum volume") self.Ca_star = daeVariable("Ca_star", no_t, self, "Dimensionless concentration of A") self.Cb_star = daeVariable("Cb_star", no_t, self, "Dimensionless concentration of B") self.P1 = daeVariable("P1", no_t, self, "Dimensionless A->B reaction rate constant") self.P2 = daeVariable("P2", no_t, self, "Dimensionless B->C reaction rate constant") self.P3 = daeVariable("P3", no_t, self, "Dimensionless 2A->D reaction rate constant") self.theta = daeVariable("theta", no_t, self, "Dimensionless time") #self.tau_ref = daeVariable("tau_ref", time_t, self, "Mean residence time of the reference reactor") self.sigma_f = daeVariable("sigma_f", no_t, self, "The proportion of the semibatch cycle occupied by filling") self.sigma_b = daeVariable("sigma_b", no_t, self, "The proportion of the semibatch cycle occupied by batch") self.sigma_e = daeVariable("sigma_e", no_t, self, "The proportion of the semibatch cycle occupied by emptying") self.time_f = daeVariable("time_f", time_t, self, "The time duration of filling phase") self.time_b = daeVariable("time_b", time_t, self, "The time duration of batch phase") self.time_e = daeVariable("time_e", time_t, self, "The time duration of emptying phase") def DeclareEquations(self): # Create adouble objects to make equations more readable rho = self.rho() cp = self.cp() kw = self.kw() AR = self.AR() mK = self.mK() cpK = self.cpK() Qk = self.Qk() dHr1 = self.dHr1() dHr2 = self.dHr2() dHr3 = self.dHr3() k10 = self.k10() k20 = self.k20() k30 = self.k30() E1 = self.E1() E2 = self.E2() E3 = self.E3() Qf = self.Qf() VR = self.VR() Tf = self.Tf() Caf = self.Caf() # Variables k1 = self.k1() k2 = self.k2() k3 = self.k3() T = self.T() Tk = self.Tk() Ca = self.Ca() Cb = self.Cb() Cc = self.Cc() Cd = self.Cd() V0 = self.V0() Vm = self.Vm() V0_star = self.V0_star() Ca_star = self.Ca_star() Cb_star = self.Cb_star() P1 = self.P1() P2 = self.P2() P3 = self.P3() theta = self.theta() #tau_ref = self.tau_ref() sigma_f = self.sigma_f() sigma_b = self.sigma_b() sigma_e = self.sigma_e() time_f = self.time_f() time_b = self.time_b() time_e = self.time_e() # Derivatives dVr_dt = dt(self.VR()) dVrCa_dt = dt(self.VR() * self.Ca()) dVrCb_dt = dt(self.VR() * self.Cb()) dVrCc_dt = dt(self.VR() * self.Cc()) dVrCd_dt = dt(self.VR() * self.Cd()) dVrT_dt = dt(self.VR() * self.T()) dTk_dt = dt(self.Tk()) # Intermediates r1 = k1 * VR * Ca r2 = k2 * VR * Cb r3 = k3 * VR * (Ca**2) ra = -r1 - 2*r3 + Qf*(Caf-Ca) rb = r1 - r2 - Qf*Cb rc = r2 - Qf*Cc rd = r3 - Qf*Cd # Volume (constant in this case) eq = self.CreateEquation("k1", "") eq.Residual = dVr_dt #- Qf # Reaction rate constants eq = self.CreateEquation("k1", "") eq.Residual = k1 - k10 * Exp(-E1 / T) eq = self.CreateEquation("k2", "") eq.Residual = k2 - k20 * Exp(-E2 / T) eq = self.CreateEquation("k3", "") eq.Residual = k3 - k30 * Exp(-E3 / T) # Mass balance eq = self.CreateEquation("Ca", "") eq.Residual = dVrCa_dt - ra eq = self.CreateEquation("Cb", "") eq.Residual = dVrCb_dt - rb eq = self.CreateEquation("Cc", "") eq.Residual = dVrCc_dt - rc eq = self.CreateEquation("Cd", "") eq.Residual = dVrCd_dt - rd # Energy balance - reactor eq = self.CreateEquation("EnergyBalanceReactor", "") eq.Residual = rho * cp * dVrT_dt - ( Qf * rho * cp * (Tf - T) \ - r1 * dHr1 \ - r2 * dHr2 \ - r3 * dHr3 \ + kw * AR * (Tk - T) ) # Energy balance - cooling fluid eq = self.CreateEquation("EnergyBalanceCooling", "") eq.Residual = mK * cpK * dTk_dt - (Qk + kw * AR * (T - Tk)) # V0* eq = self.CreateEquation("V0_star", "") eq.Residual = V0_star - V0/Vm # Ca* eq = self.CreateEquation("Ca_star", "") eq.Residual = Ca_star * Caf - Ca # Cb* eq = self.CreateEquation("Cb_star", "") eq.Residual = Cb_star * Caf - Cb # P1 eq = self.CreateEquation("P1", "") eq.Residual = P1 * Qf - k1*Vm eq.CheckUnitsConsistency = False # P2 eq = self.CreateEquation("P2", "") eq.Residual = P2 * Qf - k2*Vm eq.CheckUnitsConsistency = False # P3 eq = self.CreateEquation("P3", "") eq.Residual = P3 * Qf - k3*Vm eq.CheckUnitsConsistency = False # theta eq = self.CreateEquation("theta", "") eq.Residual = theta - Time() * Qf / Vm # tau_ref #eq = self.CreateEquation("tau_ref", "") #eq.Residual = tau_ref * (sigma_f * Qf) - Vm # sigma_f eq = self.CreateEquation("sigma_f", "") eq.Residual = sigma_f * (time_f + time_b + time_e) - time_f # sigma_b eq = self.CreateEquation("sigma_b", "") eq.Residual = sigma_b * (time_f + time_b + time_e) - time_b # sigma_e eq = self.CreateEquation("sigma_e", "") eq.Residual = sigma_e * (time_f + time_b + time_e) - time_e class simTutorial(daeSimulation): def __init__(self): daeSimulation.__init__(self) self.m = modTutorial("tutorial_che_1") self.m.Description = __doc__ def SetUpParametersAndDomains(self): self.m.k10.SetValue(1.287e10 * 1/hour) self.m.k20.SetValue(1.287e10 * 1/hour) self.m.k30.SetValue(9.043e9 * l/(mol*hour)) self.m.E1.SetValue(9758.3 * K) self.m.E2.SetValue(9758.3 * K) self.m.E3.SetValue(8560 * K) self.m.dHr1.SetValue(4.2 * kJ/mol) self.m.dHr2.SetValue(-11 * kJ/mol) self.m.dHr3.SetValue(-41.85 * kJ/mol) self.m.rho.SetValue(0.9342 * kg/l) self.m.cp.SetValue(3.01 * kJ/(kg*K)) self.m.kw.SetValue(4032 * kJ/(K*hour*(m**2))) self.m.AR.SetValue(0.215 * m**2) self.m.mK.SetValue(5 * kg) self.m.cpK.SetValue(2 * kJ/(kg*K)) def SetUpVariables(self): self.m.time_f.AssignValue( 1.0 / 16.0) self.m.time_b.AssignValue(14.0 / 16.0) self.m.time_e.AssignValue( 1.0 / 16.0) self.m.V0.AssignValue(1 * l) self.m.Vm.AssignValue(100 * l) self.m.Qf.AssignValue(14.19 * l/hour) self.m.Qk.AssignValue(-1579.5 * kJ/hour) self.m.Caf.AssignValue(5.1 * mol/l) self.m.Tf.AssignValue((273.15 + 104.9) * K) self.m.VR.SetInitialCondition(10.0 * l) self.m.Ca.SetInitialCondition(2.2291 * mol/l) self.m.Cb.SetInitialCondition(1.0417 * mol/l) self.m.Cc.SetInitialCondition(0.91397 * mol/l) self.m.Cd.SetInitialCondition(0.91520 * mol/l) self.m.T.SetInitialCondition((273.15 + 79.591) * K) self.m.Tk.SetInitialCondition((273.15 + 77.69) * K) # Variable types can be: # cnAlgebraic = algebraic state variables # cnDifferential = differential state variables # cnAssigned = assigned variables (degrees of freedom) # Distributed variables can also be of a mixed type - i.e. some points can be differential, # while the others can be algebraic or assigned: # cnSomePointsAssigned = algebraic + some assigned # cnSomePointsDifferential = algebraic + some differential # cnMixedAlgebraicAssignedDifferential = algebraic + assigned + differential varTypes = {cnAlgebraic : [], cnDifferential : [], cnAssigned : []} for var in self.m.Variables: varTypes[var.Type].append(var.Name) for varType, varNames in varTypes.items(): if varType == cnAlgebraic: print('Algebraic: %s' % varNames) elif varType == cnDifferential: print('Differential: %s' % varNames) elif varType == cnAssigned: print('Assigned: %s' % varNames) else: print('Mixed type: %s' % varNames) def run(**kwargs): simulation = simTutorial() return daeActivity.simulate(simulation, reportingInterval = 600, # 10 min timeHorizon = 3*60*60, # 3h **kwargs) def runAsWebService(): # The loader function should have the argument 'initializeAndReturn' set to True. # Optionally, the argument 'datareporter' should be set to daeNoOpDataReporter() # to allow the web service clients to get the results. def loaderFunction(**kwargs): print('Arguments sent by the web service client: %s' % kwargs) return run(initializeAndReturn = True, datareporter = daeNoOpDataReporter(), **kwargs) # Add the simulation name to the 'availableSimulations' dictionary. # The value is a Python callable objects that returns an initialised simulation object. availableSimulations = {} availableSimulations['tutorial_che_1'] = loaderFunction daeSimulationWebService.runSimulationsAsWebService(availableSimulations) if __name__ == "__main__": guiRun = False if (len(sys.argv) > 1 and sys.argv[1] == 'console') else True run(guiRun = guiRun)