Dynamic Simulation Tutorial with DWSIM and Python, Part 2: Building the Dynamic Model

From DWSIM - Open Source Chemical Process Simulator
Jump to navigation Jump to search
Dialog-warning.png As of DWSIM v6.0, which include native Dynamic Simulation capabilities, this tutorial is obsolete.
Dialog-warning.png This tutorial requires advanced or above average Python programming skills.
Dialog-information.png You'll need at least DWSIM v5.7 (Cross-Platform UI) or newer on Windows, Linux or macOS to follow/reproduce the tasks within this tutorial.


Introduction

In order to make our model dynamic, we need to upgrade the valve model so it can set the cold water flow rate based on a given opening. Also, we need to add some 'delay' to heat exchanger model so the hot water doesn't instantaneously cools down upon a flow rate change. This is what happens in the real world, since both fluids have a residence time inside the exchanger. Any change in temperature and/or flow rate at the inlets take some time to "appear" at the outlets.

Creating the Dynamic Model

Go to the Script Manager and add a new script. Name it "InitVars". Add the following lines to the top of it: 

import System
from System import *

We need to add a few new variables to our flowsheet so we can keep track of the current time step during the dynamic simulation, record the time where a perturbation was added to the system, how much we will perturb the hot water flow rate and at which time we will do it: 

# Flowsheet
 
Flowsheet.SupressMessages = False
 
Flowsheet.ExtraProperties.TimeStep = 1.0
Flowsheet.ExtraProperties.CurrentTime = 0.0
Flowsheet.ExtraProperties.LastPerturbationTime = 0.0
Flowsheet.ExtraProperties.MaxTime = 400.0
Flowsheet.ExtraProperties.HotWaterPerturbationTime = 50.0
 
Flowsheet.ExtraProperties.PerturbationValue = 13.0

We will run our model from 0 to 400 s, at 1s-time steps. The hot water flow rate will be perturbed at 50 s, and its value will be set to 13 kg/s and remain at this value until the end of the simulation.

Valve

Our dynamic valve model will convert an opening value (from 0 to 100) to the actual flow rate through it. At the end of the calculation, we will set the dummy_in mass flow rate to the calculated value.

To convert the opening to flow rate, we will need a few new variables (add the following to the 'InitVars' script): 


# Control Valve
 
valve = Flowsheet.GetFlowsheetSimulationObject("FV-01")
 
valve.ExtraProperties.MaxCv = 3102.78
valve.ExtraProperties.Cv = 3102.78
valve.ExtraProperties.OpeningPct = 50.0
 
valve.ExtraProperties.CalculatedFlow = 0.0

For the actual dynamic valve model, create a new script and name it 'UpdateValve'. Associate this script to the 'Object Calculation Finished' event of the 'FV-01' object. This is the complete model: 

import System
from System import *
 
valve = Flowsheet.GetFlowsheetSimulationObject("FV-01")
inlet = Flowsheet.GetFlowsheetSimulationObject("dummy_in")
 
props = valve.ExtraProperties
 
x = props.OpeningPct / 100.0
 
props.Cv = props.MaxCv * 2.0 * Math.Exp(x * Math.Log(50.0)) / 100.0
 
# Cv = 11.6 Q (SG / dp)^0.5
# dp = SG/(Cv/11.6Q)^2
# where
# q = water flow (m3/hr)
# SG = specific gravity (1 for water)
# dp = pressure drop (kPa)
 
DP = valve.DeltaP / 1000
 
SG = inlet.Phases[0].Properties.density / 1000
 
Q = props.Cv / 11.6 / (SG/DP) ** 0.5 / 60 / 60 # m3/s
 
props.CalculatedFlow = Q # m3/s
 
inlet.Phases[0].Properties.massflow = Q * inlet.Phases[0].Properties.density

The above model will convert a given opening to the actual flow coefficient by using the equal percentage conversion equation.

Heat Exchanger

Our dynamic heat exchanger model will take into account the fluid residence time on both sides (cold and hot) and the calculated steady-state temperature value, applying a quadratic delay to the steady-state calculated value which will only be achieved after the fluid has exited the exchanger.

For our dynamic model to work, we need to keep track of the current time and the last 'perturbation' time. We need to add the following variables to the 'InitVars' script: 

# Heat Exchanger
 
hx = Flowsheet.GetFlowsheetSimulationObject("HX-01")
 
hx.ExtraProperties.VolumeForHotFluid = 0.0
hx.ExtraProperties.HotFluidResidenceTime = 0.0
 
hx.ExtraProperties.VolumeForColdFluid = 0.0
hx.ExtraProperties.ColdFluidResidenceTime = 0.0
 
hx.ExtraProperties.ColdFluidPrePerturbationOutletTemperature = 0.0
hx.ExtraProperties.HotFluidSPrePerturbationOutletTemperature = 0.0

For the actual dynamic heat exchanger model, create a new script and name it 'UpdateExchanger'. Associate this script to the 'Object Calculation Finished' event of the 'HX-01' object. This is the complete model: 

import System
from System import *
 
hx = Flowsheet.GetFlowsheetSimulationObject("HX-01")
 
t = Flowsheet.ExtraProperties.CurrentTime
t0 = Flowsheet.ExtraProperties.LastPerturbationTime
 
deltat = t - t0
 
hxhf = Flowsheet.GetFlowsheetSimulationObject("hot_water")
hxcf = Flowsheet.GetFlowsheetSimulationObject("flow_2")
 
hx.ExtraProperties.VolumeForHotFluid = hx.STProperties.Tube_NumberPerShell * hx.STProperties.Tube_Length * Math.PI * (hx.STProperties.Tube_Di) ** 2 / 4 # m3
hx.ExtraProperties.HotFluidResidenceTime = hx.ExtraProperties.VolumeForHotFluid / hxhf.Phases[0].Properties.volumetric_flow # s
 
hx.ExtraProperties.VolumeForColdFluid = hx.STProperties.Tube_Length * Math.PI * (hx.STProperties.Shell_Di) ** 2 / 4 - hx.ExtraProperties.VolumeForHotFluid # m3
hx.ExtraProperties.ColdFluidResidenceTime = hx.ExtraProperties.VolumeForColdFluid / hxcf.Phases[0].Properties.volumetric_flow # s
 
delay1 = hx.ExtraProperties.ColdFluidResidenceTime
delay2 = hx.ExtraProperties.HotFluidResidenceTime
 
if (deltat <= delay1 and t0 > 0.0):
    hx.ColdSideOutletTemperature = hx.ExtraProperties.ColdFluidPrePerturbationOutletTemperature + (deltat/delay1)**2 * (hx.ColdSideOutletTemperature - hx.ExtraProperties.ColdFluidPrePerturbationOutletTemperature)
 
if (deltat <= delay2 and t0 > 0.0):
    hx.HotSideOutletTemperature = hx.ExtraProperties.HotFluidPrePerturbationOutletTemperature + (deltat/delay2)**2 * (hx.HotSideOutletTemperature - hx.ExtraProperties.HotFluidPrePerturbationOutletTemperature)

Tanks

For the source and sink tanks, we will use simple dynamic models to keep track of the current water level based on the inlet and outlet flows.

Add to 'InitVars': 

# Tanks
 
source_tank = Flowsheet.GetFlowsheetSimulationObject("T-01")
 
source_tank.ExtraProperties.Diameter = 15.0
source_tank.ExtraProperties.CurrentLiquidHeight = 12.0
 
sink_tank = Flowsheet.GetFlowsheetSimulationObject("T-02")
 
sink_tank.ExtraProperties.Diameter = 15.0
sink_tank.ExtraProperties.CurrentLiquidHeight = 0.0

Helper Variables

Add the following to 'InitVars': 

# Cooling Water
 
cooling_water = Flowsheet.GetFlowsheetSimulationObject("dummy_in")
cooling_water.Phases[0].Properties.massflow = 23.6155
 
# Hot Water
 
hot_water = Flowsheet.GetFlowsheetSimulationObject("hot_water")
hot_water.Phases[0].Properties.massflow = 7.46

Helper Functions

Create a new script and name it 'Functions'. Add the following content to it: 

def SetValveOpening(x):
    # x in %   
    valve = Flowsheet.GetFlowsheetSimulationObject("FV-01")
    valve.ExtraProperties.OpeningPct = x
 
def GetValveOpening():
    # in % 
    valve = Flowsheet.GetFlowsheetSimulationObject("FV-01")
    return valve.ExtraProperties.OpeningPct
     
def SetHotWaterMassFlow(x):
    # x in kg/s
    hot_water = Flowsheet.GetFlowsheetSimulationObject("hot_water")
    hot_water.Phases[0].Properties.massflow = x
     
def GetHotWaterMassFlow():
    hot_water = Flowsheet.GetFlowsheetSimulationObject("hot_water")
    return hot_water.Phases[0].Properties.massflow
 
def GetCoolingWaterMassFlow():
    cold_water = Flowsheet.GetFlowsheetSimulationObject("dummy_in")
    return cold_water.Phases[0].Properties.massflow
     
def GetHotWaterOutletTemperature():
    hx = Flowsheet.GetFlowsheetSimulationObject("HX-01")
    return hx.HotSideOutletTemperature
 
def GetCoolingWaterOutletTemperature():
    hx = Flowsheet.GetFlowsheetSimulationObject("HX-01")
    return hx.ColdSideOutletTemperature
 
def StorePrePerturbationVariables():
    hx = Flowsheet.GetFlowsheetSimulationObject("HX-01")
    hx.ExtraProperties.ColdFluidPrePerturbationOutletTemperature = hx.ColdSideOutletTemperature
    hx.ExtraProperties.HotFluidPrePerturbationOutletTemperature = hx.HotSideOutletTemperature
     
def CalcSourceTankLevel():
    tank = Flowsheet.GetFlowsheetSimulationObject("T-01")
    valve = Flowsheet.GetFlowsheetSimulationObject("FV-01")
    tstep = Flowsheet.ExtraProperties.TimeStep
    voldelta = valve.ExtraProperties.CalculatedFlow * tstep
    hdelta = voldelta / Math.PI / tank.ExtraProperties.Diameter ** 2
    tank.ExtraProperties.CurrentLiquidHeight -= hdelta  
 
def GetSourceTankLevel():
    tank = Flowsheet.GetFlowsheetSimulationObject("T-01")
    return tank.ExtraProperties.CurrentLiquidHeight
     
def CalcSinkTankLevel():
    tank = Flowsheet.GetFlowsheetSimulationObject("T-02")
    valve = Flowsheet.GetFlowsheetSimulationObject("FV-01")
    tstep = Flowsheet.ExtraProperties.TimeStep
    voldelta = valve.ExtraProperties.CalculatedFlow * tstep
    hdelta = voldelta / Math.PI / tank.ExtraProperties.Diameter ** 2
    tank.ExtraProperties.CurrentLiquidHeight += hdelta  
 
def GetSinkTankLevel():
    tank = Flowsheet.GetFlowsheetSimulationObject("T-02")
    return tank.ExtraProperties.CurrentLiquidHeight
     
def GetSinkTankLevel():
    tank = Flowsheet.GetFlowsheetSimulationObject("T-02")
    return tank.ExtraProperties.CurrentLiquidHeight

Running the Dynamic Model

Create a new script and name it 'RunDynamicProcess_OpenLoop'. Add the following content to it: 

import clr
import System
from System import *
from System.Threading import *
 
clr.AddReference('System.Core')
clr.AddReference('DWSIM.GlobalSettings')
clr.ImportExtensions(System.Linq)
 
from DWSIM.GlobalSettings import *
 
source = Flowsheet.Scripts.Values.Where(lambda x: x.Title == 'Functions').FirstOrDefault().ScriptText.replace('\r', '')
exec(source)
 
initvars = Flowsheet.Scripts.Values.Where(lambda x: x.Title == 'InitVars').FirstOrDefault().ScriptText.replace('\r', '')
exec(initvars)
 
maxtime = Flowsheet.ExtraProperties.MaxTime
length = maxtime / Flowsheet.ExtraProperties.TimeStep
 
time = [0.0 + x*(maxtime - 0.0)/length for x in range(int(length))]

With the above commands, we get what we typed in the 'Functions' and 'InitVars' scripts using a LINQ/lambda command (Where(lambda x:...)) and append it to our current script, so they can be executed sequentially. This will initialize our variables and add our helpers functions to the memory. We also set our time variable, which will range from 0 to 400 at unit steps.

Run the Script

Our actual full dynamic process model is: 

source_level = {}
sink_level = {}
hot_water_flow = {}
hot_water_temp = {}
cooling_water_flow = {}
cooling_water_temp = {}
valve_opening = {}
 
perturbed_hotwater = False
 
for t in time:
    Flowsheet.SupressMessages = False
    Flowsheet.WriteMessage("Time Step: " + str(int(t+1)) + "/" + str(int(maxtime)))
    Flowsheet.SupressMessages = True
    Flowsheet.ExtraProperties.CurrentTime = t
    if (t >= Flowsheet.ExtraProperties.HotWaterPerturbationTime and not perturbed_hotwater): 
        StorePrePerturbationVariables()
        SetHotWaterMassFlow(Flowsheet.ExtraProperties.PerturbationValue)
        Flowsheet.ExtraProperties.LastPerturbationTime = t
        perturbed_hotwater = True
    Flowsheet.SolveFlowsheet2()
    CalcSourceTankLevel()
    CalcSinkTankLevel()
    source_level[t] = GetSourceTankLevel()
    sink_level[t] = GetSinkTankLevel()
    hot_water_flow[t] = GetHotWaterMassFlow()
    hot_water_temp[t] = GetHotWaterOutletTemperature()
    cooling_water_flow[t] = GetCoolingWaterMassFlow() 
    cooling_water_temp[t] = GetCoolingWaterOutletTemperature()
    valve_opening[t] = GetValveOpening()
 
Flowsheet.SupressMessages = False

Add the above code to 'RunDynamicProcess_OpenLoop' and click on "Update and Run (Async)".

This will run our simulation at each time step, keeping track of the calculated values and storing it on python lists. The Valve and Heat Exchanger dynamic models will be ran automatically by the flowsheet solver, since you've linked them and told the solver to run them after the corresponding objects are calculated.

Viewing Results in Excel

To view the results in a Excel spreadsheet, add a new script and name it 'SaveToExcel'. Add the following content: 

import clr
import sys
 
clr.AddReferenceByName('Microsoft.Office.Interop.Excel, Version=11.0.0.0, Culture=neutral, PublicKeyToken=71e9bce111e9429c')
 
from Microsoft.Office.Interop import Excel
 
ex = Excel.ApplicationClass()   
 
ex.Visible = False
ex.DisplayAlerts = True  
 
workbook = ex.Workbooks.Add()
workbook.Worksheets.Add()
ws1 = workbook.Worksheets[1]
 
ws1.UsedRange.Cells[1, 1].Value2 = "time"
ws1.UsedRange.Cells[1, 2].Value2 = "source_level"
ws1.UsedRange.Cells[1, 3].Value2 = "sink_level"
ws1.UsedRange.Cells[1, 4].Value2 = "hot_water_flow"
ws1.UsedRange.Cells[1, 5].Value2 = "hot_water_temp"
ws1.UsedRange.Cells[1, 6].Value2 = "cooling_water_flow"
ws1.UsedRange.Cells[1, 7].Value2 = "cooling_water_temp"
ws1.UsedRange.Cells[1, 8].Value2 = "valve_opening"
 
k = 2
     
for t in time:
    ws1.UsedRange.Cells[k, 1].Value2 = t
    ws1.UsedRange.Cells[k, 2].Value2 = source_level[t]
    ws1.UsedRange.Cells[k, 3].Value2 = sink_level[t]
    ws1.UsedRange.Cells[k, 4].Value2 = hot_water_flow[t]
    ws1.UsedRange.Cells[k, 5].Value2 = hot_water_temp[t]
    ws1.UsedRange.Cells[k, 6].Value2 = cooling_water_flow[t]
    ws1.UsedRange.Cells[k, 7].Value2 = cooling_water_temp[t]
    ws1.UsedRange.Cells[k, 8].Value2 = valve_opening[t]
    k += 1
 
ex.Visible = True

Add the following to 'RunDynamicProcess_OpenLoop' and run it again: 

outputresults = Flowsheet.Scripts.Values.Where(lambda x: x.Title == 'SaveToExcel').FirstOrDefault().ScriptText.replace('\r', '')
exec(outputresults)

You'll get this:

Dynamic excel.jpg

Viewing Results in Charts

To view the results in built-in DWSIM charts, add a new script and name it 'GenerateCharts'. Add the following content: 

import clr
clr.AddReference("Eto.OxyPlot")
clr.AddReference("OxyPlot")
clr.AddReference("DWSIM.ExtensionMethods.Eto")
clr.AddReference("DWSIM.ExtensionMethods")
import DWSIM.ExtensionMethods
clr.ImportExtensions(DWSIM.ExtensionMethods)
import OxyPlot
 
from Eto.OxyPlot import Plot
from DWSIM.UI.Shared import *
from System import Array
 
def CreatePlots():
 
    chart1 = Plot()
    chart1.Model = Common.CreatePlotModel(Array[float](time), Array[float](valve_opening.values()), "Valve Opening", "", "time (s)", "Opening (%)") 
    chart1.Model.AddYAxis("Flows (kg/s)", "flows", OxyPlot.Axes.AxisPosition.Right, 1) 
    chart1.Model.AddLineSeriesWithKeys(Array[float](time), Array[float](hot_water_flow.values()), "Hot Water Flow", "flows", "x")
    chart1.Model.AddLineSeriesWithKeys(Array[float](time), Array[float](cooling_water_flow.values()), "Cooling Water Flow", "flows", "x")
    chart1.Model.InvalidatePlot(True)
     
    form1 = Common.GetDefaultEditorForm("Valve Opening", 600, 600, chart1, False)
    form1.Show()
     
    chart2 = Plot()
    chart2.Model = Common.CreatePlotModel(Array[float](time), Array[float](hot_water_temp.values()), Array[float](cooling_water_temp.values()),"Outlet Temperatures", "", "time (s)", "Temperature (K)", "Hot Water", "Cooling Water") 
     
    form2 = Common.GetDefaultEditorForm("Fluid Temperatures", 600, 600, chart2, False)
    form2.Show()
     
Application.Invoke(lambda: CreatePlots())

Add the following to 'RunDynamicProcess_OpenLoop' and run it again: 

outputresults = Flowsheet.Scripts.Values.Where(lambda x: x.Title == 'GenerateCharts').FirstOrDefault().ScriptText.replace('\r', '')
exec(outputresults)

You'll get this:

Dynamic charts.jpg

Download File

Download the simulation file with what has been done so far: dynamic_part2.dwxmz

Return to Dynamic Simulation Tutorial with DWSIM and Python, Part 1: Concepts and Steady-State Model

Proceed to Dynamic Simulation Tutorial with DWSIM and Python, Part 3: Adding a PID Controller

Retrieved from "https://dwsim.org/wiki/index.php?title=Dynamic_Simulation_Tutorial_with_DWSIM_and_Python,_Part_2:_Building_the_Dynamic_Model&oldid=421"