© 2021 Jasper van Baten, AmsterCHEM

Examples:

Table of contents

Examples

The basics

Let's start with the basics. Each Unit Operation requires a Configure function and a Calculate function. Each of these functions takes a unit argument. During Configure, one can add parameters, ports and reports, using add_parameter, add_port and add_report.

During Calculate, the ports can be obtained from unit.ports[<name>], that is, the objects connected to the ports are obtained like that. If a ports is not connected, it will not appear in the unit.ports dictionary. Ports can only be left unconnected, if optional=True was specified as argument to add_port. Values can be obtained from FEED ports, while PRODUCT ports must be fully specified.

Similarly, parameters are available during Calculate from unit.parameters[<name>]. The value of an INPUT parameter is known at the start of Calculate, whereas the value of an OUTPUT parameter must be set during Calculate.

Reports are available from unit.reports[<name>]; reports can subsequently be written using e.g. print with the file argument set to the report.

Example 1: a simple heater.

def Configure(unit): #add a feed and product unit.add_port("Feed",unit.PortType.MATERIAL,unit.PortDirection.FEED) unit.add_port("Product",unit.PortType.MATERIAL,unit.PortDirection.PRODUCT) #add a heat duty parameter WATT=unit.Dimension.KILOGRAM*unit.Dimension.METER**2/unit.Dimension.SECOND**3 unit.add_parameter("Heat Duty",unit.ParameterType.REAL,unit.ParameterMode.INPUT,default=0,dimension=WATT) def Calculate(unit): #copy Feed to Product, apply heat duty fd=unit.ports["Feed"] duty=unit.parameters["Heat Duty"].value if fd.flow_rate==0 and duty!=0: raise RuntimeError('Zero feed flow rate; cannot apply duty') h=fd.get_property('enthalpy') #molar enthalpy if (duty!=0): h+=duty/fd.flow_rate unit.ports["Product"].set(x=fd.x,p=fd.p,flow_rate=fd.flow_rate,enthalpy=h)

During configure, two MATERIAL ports are added; a FEED called "Feed", and a PRODUCT called "Product".

Calculating thermo-physical properties

The example uses fd.flow_rate, fd.x and fd.t to access total flow rate, mol/s, overall composition, mol/mol, and temperature, K. Material streams, or rather phase equilibria in general, are defined by composition, temperature, pressure, phase presence, and for each phase, composition, temperature, pressure and phase fraction. These properties can be used directly.

Other properties can be calculated using get_property; for a FEED MATERIAL and PhaseEquilibrium, get_property delivers equilibrium properties; for a Phase, get_property delivers single phase properties, for a PhasePair, get_property delivers two-phase properties. The lists of equilibrium properties, single phase properties and two phase properties can be obtained from a MATERIAL stream using phase_equilibrium_properties, phase_properties and phase_pair_properties. Alternatively, a Phase returns single the list of available single phase properties as properties, and a phase pair returns the list of two phase properties as properties.

Pure compound constant properties can be obtained from a Compound; the list of compounds is exposed by a MATERIAL stream's compounds member. Constant properties can be obtained from a Compound using get_constant, and temperature dependent properties can be obtained using get_temperature_dependent_property; the lists of these properties are available from constant_properties and temperature_dependent_properties for a Compound, or, equivalently, constant_compound_properties and temperature_dependent_compound_properties from a MATERIAL stream. One can directly obtain name, molecular_weight, formula and cas_registry_number from a compound.

Example 2: thermo-physical property calculations.

#Python CAPE-OPEN Unit Operation script def Configure(unit): #example: add a feed unit.add_port("Feed",unit.PortType.MATERIAL,unit.PortDirection.FEED) def Calculate(unit): #example: calculate physical properties fd=unit.ports["Feed"] #some compound properties compound=fd.compounds[0]; #first compound print("For compound",compound.name) print("Molecular weight",compound.molecular_weight,"gr/mol") print("Chemical formula",compound.formula) print("Critical temperature",compound.get_constant('criticalTemperature')) for t in [300,350,400]: print("Vapor pressure at",t,"K = ",compound.get_temperature_dependent_property('vaporPressure',t),"Pa") #one can get multiple properties as a tuple by specifying a list of property names [Psat,dPsatdT]=compound.get_temperature_dependent_property(['vaporPressure','vaporPressure.Dtemperature'],300) #equilibrium properties for the feed stream print("Feed temperature",fd.t,"K") print("Feed pressure",fd.p,"Pa") print("Feed composition",fd.x,"mol/mol") print("Feed vapor fraction",fd.vapor_fraction,"mol/mol") #and rate print("Feed rate",fd.flow_rate,"mol/s") #other properties require calculation print("Feed enthalpy",fd.get_property('enthalpy'),"J/mol") #multiple properties can be calculated in one call, returned as tuple [h,s,v]=fd.get_property(['enthalpy','entropy','volume']) #the same calculations can be performed on any phase equilibrium eq=fd.create_phase_equilibrium(fd.x,'pressure',2e5,'enthalpy',h) #adiabatic re-flash of feed at 2 bar eq=fd.create_phase_equilibrium(x=fd.x,p=2e5,enthalpy=h) #same thing, alternative syntax print("Equilibrium temperature",eq.t,"K") print("Equilibrium pressure",eq.p,"Pa") print("Equilibrium composition",eq.x,"mol/mol") print("Equilibrium vapor fraction",eq.vapor_fraction,"mol/mol") #other properties require calculation print("Equilibrium entropy",eq.get_property('entropy'),"J/mol/K") #multiple properties can be calculated in one call, returned as tuple [h,s,v]=eq.get_property(['enthalpy','entropy','volume']) #one can iterate over phases of feed or calculated equilibrium and # obtain single phase properties: for p in fd.phases: print(p.name,'temperature',p.t,'K') print(p.name,'pressure',p.p,'Pa') print(p.name,'composition',p.x,'mol/mol') print(p.name,'phase fraction',p.phase_fraction,'mol/mol') #other single phase properties must be calculated print(p.name,'enthalpy',p.get_property('enthalpy'),'J/mol') #again mutiple properties can be calculated in one go [visc,density]=p.get_property(['viscosity','density']) #for two phase properties, two phases can be combined in a # PhasePair object eq=fd.create_phase_equilibrium(x=fd.x,p=1e5,vaporFraction=0.5) #flash results vapor and liquid vapor_liquid=fd.create_phase_pair(eq.phases[0],eq.phases[1]) print('surface tension',vapor_liquid.get_property('surfaceTension'),'N/m') #non-equilibrium phases can be created from a PhaseType or directly vapor1=fd.phase_types[0].create_phase(t=300,p=2e5,x=fd.x) vapor2=fd.create_vapor_phase(t=300,p=2e5,x=fd.x) print('vapor 1 density',vapor1.get_property('density'),'mol/m3') print('vapor 2 density',vapor2.get_property('density'),'mol/m3')

Units of measure

In example 1, an INPUT, REAL parameter "Heat duty" is added that has the dimension WATT.

WATT is constructed from the base dimensions as kg*m2/s3. All real numbers in CAPE-OPEN are exchanged in pure SI units. For real parameters and information ports you get to decide on the dimensionality of the number; for properties on MATERIAL ports, the dimension is predefined. In many cases CAPE-OPEN allows a choice between specific and molar units; the Python Unit Operation always selects molar units. So flow rates are in mol/s, compositions are in mol/mol, density is in mol/m3, enthalpy in J/mol, etc.

For temperature differences and pressure differences, conversions from guage pressure and Celcius or Fahrenheit should not take into account the offset. For such quantities, pass difference=True in addition to the dimension argument.

Example 3 shows some ways to setup units of measure. A Calculate method is not provided in example 3.

Example 3: Dimension and units of measure.

def Configure(unit): #some unit of measure examples JOULE=unit.Dimension.KILOGRAM*unit.Dimension.METER**2/unit.Dimension.SECOND**2 WATT=JOULE/unit.Dimension.SECOND HERTZ=1/unit.Dimension.SECOND NEWTON=unit.Dimension.KILOGRAM*unit.Dimension.METER/unit.Dimension.SECOND**2 PASCAL=NEWTON/unit.Dimension.METER**2 #a heat duty parameter unit.add_parameter("Heat Duty",unit.ParameterType.REAL,unit.ParameterMode.INPUT,default=0,dimension=WATT) #a pressure drop parameter unit.add_parameter("Pressure drop",unit.ParameterType.REAL,unit.ParameterMode.INPUT,default=0,dimension=PASCAL,difference=True) #dimensions are also passed to information streams unit.add_port("Frequency",unit.PortType.INFORMATION,unit.PortDirection.PRODUCT,dimension=HERTZ)

Validation

CAPE-OPEN unit operations are validated prior to calculation, and the PME is required to validate a Unit Operation after any change in configuration (not including changes to values on feed or product streams, provided the configuration of the streams does not change).

Unit operations do not need to be validated in case Calculate is called multiple times and the configuration is not changing, for example if the Unit Operation is in a recycle or during a parametric study.

Python Unit Operations may optionally implement a Validate function. A Validate function can raise an exception in case the configuration of the Unit Operation is not valid, or in case the streams connected to the Unit Operation have a configuration that is not suitable for the Unit Operation. In addition, the Unit Operation can store data during validation, that pertains to the configuration of the streams. After all, the PME must revalidate the Unit Operation in case this configuration changes. Such data typically includes information on phase or compound mapping.

During Validate, data can be stored in the validation_data dictionary, and then the data will be available during Calculate. The Unit Operation must not store any reference to external objects, such as streams, compounds, etc, directly or indirectly. Such objects cannot be referenced between calculations or validations. If the Unit Operation stores external references, behaviour is unpredictable.

The following example shows a water separator. This unit copies the feed stream to the product stream after taking out any water. The water is placed on a separate product stream called "Water". The example features a Validate function that locates water, first by CAS registry number, then by chemical formula, and if all fails, by name. Exception handling is used in case compounds are present that have no CAS registry number or chemical formula. If water is not found, validation fails. Otherwise the index of water in the compound slate is stored in validation_data and used during Calculate. For vector math, numpy is used. Note that no special care is taken in this example for zero flow rate of the Product stream, which will lead to a division by zero.

Example 4: Water separator.

import numpy as np #Example water separator def Configure(unit): unit.add_port("Feed",unit.PortType.MATERIAL,unit.PortDirection.FEED) unit.add_port("Product",unit.PortType.MATERIAL,unit.PortDirection.PRODUCT) unit.add_port("Water",unit.PortType.MATERIAL,unit.PortDirection.PRODUCT) def Validate(unit): water_index=-1 compounds=unit.ports["Feed"].compounds #try by cas number for i in range(0,len(compounds)): try: #in case of compound without CAS number if compounds[i].cas_registry_number=='7732-18-5': water_index=i break except: pass #try by chemical formula if water_index==-1: for i in range(0,len(compounds)): try: #in case of compound without formula if compounds[i].formula=='H2O': water_index=i break except: pass #try by name if water_index==-1: for i in range(0,len(compounds)): if compounds[i].name.lower()=='water': water_index=i break if water_index==-1: raise RuntimeError('water not found') #store water index for Calculate unit.validation_data['water_index']=water_index def Calculate(unit): fd=unit.ports["Feed"] #compound flow rates comp_rates=np.array(fd.x)*fd.flow_rate #take water out water_index=unit.validation_data['water_index'] water_rate=comp_rates[water_index] comp_rates[water_index]=0 #set product, at feed t and p product_rate=comp_rates.sum() unit.ports["Product"].set(flow_rate=product_rate,x=comp_rates/product_rate,p=fd.p,t=fd.t) #set water stream at feed t and p x_water=np.zeros(len(comp_rates)) x_water[water_index]=1 unit.ports["Water"].set(flow_rate=water_rate,x=x_water,p=fd.p,t=fd.t)

Your custom implementation of Validate is called after all built-in validation has already passed. Built in validation checks that all ports that are not marked as optional=True are connected, and that all input parameters not marked with optional=True are specified, and that the values are valid (e.g. not violating imposed minimum and maximum values). Then it will check that the compounds found on each of the connected material streams are consistent. For some PMEs this is always the case. Other PMEs allow you to define streams with different thermodynamic subsystems, and therefore different compounds, and will allow you to connect such streams to a unit operation. By default, the Python Unit Operation will not validate in such a case.

Sometimes it can be useful to allow for this. Consider a heat exchanger with a hot and a cold stream, or more generally, a stream 1 and stream 2. There is no particular reason that stream 1 and stream 2 should have the same compound slates. One side of your heat exchanger may be a process stream, while the other side may be steam for example. Consider the following configuration:

Example 5: Heat exchanger configuration.

def Configure(unit): unit.add_port("Feed 1",unit.PortType.MATERIAL,unit.PortDirection.FEED) unit.add_port("Product 1",unit.PortType.MATERIAL,unit.PortDirection.PRODUCT) unit.add_port("Feed 2",unit.PortType.MATERIAL,unit.PortDirection.FEED) unit.add_port("Product 2",unit.PortType.MATERIAL,unit.PortDirection.PRODUCT)

Now, connecting Feed 1 consistently with Product 1, and Feed 2 consistently with Product 2, but Feed 1 and Feed 2 inconsistently, will fail to validate with the following message:

material objects connected to ports 'Feed 2' and 'Feed 1' do not have same compounds

This behaviour can be changed by setting unit.allow_inconsistent_compounds=True in Configure. At this point you should take care of validation of consistent compound slates where required yourself.

Example 6: Counter-current heat exchanger with maximum heat exchange.

#Heat Exchanger, Counter-current, Max heat exchange def Configure(unit): #material ports unit.add_port("Feed 1",unit.PortType.MATERIAL,unit.PortDirection.FEED) unit.add_port("Product 1",unit.PortType.MATERIAL,unit.PortDirection.PRODUCT) unit.add_port("Feed 2",unit.PortType.MATERIAL,unit.PortDirection.FEED) unit.add_port("Product 2",unit.PortType.MATERIAL,unit.PortDirection.PRODUCT) #allow inconsistent compound slates unit.allow_inconsistent_compounds=True def Validate(unit): compoundSlate1=[comp.name for comp in unit.ports["Feed 1"].compounds] compoundSlate2=[comp.name for comp in unit.ports["Product 1"].compounds] if compoundSlate1!=compoundSlate2: raise RuntimeError('Feed 1 and Product 1 do not have the same compounds') compoundSlate1=[comp.name for comp in unit.ports["Feed 2"].compounds] compoundSlate2=[comp.name for comp in unit.ports["Product 2"].compounds] if compoundSlate1!=compoundSlate2: raise RuntimeError('Feed 2 and Product 2 do not have the same compounds') def Calculate(unit): fd1=unit.ports["Feed 1"] fd2=unit.ports["Feed 2"] pd1=unit.ports["Product 1"] pd2=unit.ports["Product 2"] if fd1.flow_rate==0: #zero flow case, no heat transfer possible; bring stream 1 to feed 2 temperature pd1.set(x=fd1.x,p=fd1.p,flow_rate=fd1.flow_rate,t=fd2.t) # do a PH flash in case stream 2 is a pure compound pd2.set(x=fd2.x,p=fd2.p,flow_rate=fd2.flow_rate,h=fd2.get_property('enthalpy')) elif fd2.flow_rate==0: #zero flow case, no heat transfer possible; bring stream 2 to feed 1 temperature pd2.set(x=fd2.x,p=fd2.p,flow_rate=fd2.flow_rate,t=fd1.t) # do a PH flash in case stream 1 is a pure compound pd1.set(x=fd1.x,p=fd1.p,flow_rate=fd1.flow_rate,h=fd1.get_property('enthalpy')) else: #work required to bring stream 2 to feed temperature 1 q1=(pd2.create_phase_equilibrium(x=fd2.x,t=fd1.t,p=fd2.p).get_property('enthalpy')-fd2.get_property('enthalpy'))*fd2.flow_rate #work required to bring stream 1 to feed temperature 2 q2=(pd1.create_phase_equilibrium(x=fd1.x,t=fd2.t,p=fd1.p).get_property('enthalpy')-fd1.get_property('enthalpy'))*fd1.flow_rate #take the smaller absolute heat duty, which is the max. possible heat duty # take the sign such that we apply it possitively to stream 1 and negatively to stream 2 q=-q1 if abs(q1)<abs(q2) else q2 #apply to both streams pd1.set(x=fd1.x,p=fd1.p,flow_rate=fd1.flow_rate,enthalpy=fd1.get_property('enthalpy')+q/fd1.flow_rate) pd2.set(x=fd2.x,p=fd2.p,flow_rate=fd2.flow_rate,enthalpy=fd2.get_property('enthalpy')-q/fd2.flow_rate)

Parameters

The above heat exchanger example is starting to look like an actual unit operation. But it is not very flexible. Let's add some INPUT parameters to allow further configuration of the unit operation.

Start with adding two pressure drop parameters, one for stream 1 and one for stream 2. These are created as INPUT, REAL parameters, with dimension=PASCAL, with Pa = N/m2 = (kg*m/s)/m2. Also, as these are pressure drops, in case of a unit of measure of e.g. barg, the offset in the unit of measure should be ignored. This is indicated with difference=True.

Next, an option is added that allows to switch between counter- and co-current operation. For counter-current operation, the maximum heat transfer is as in Example 6. For co-current heat transer, the maximum heat transfer is that for which the temperature of Product 1 equals the temperature of Product 2. This requires a root finder, and the toms748 routine from scipy is used. The limits are given by 0, and the maximum heat transfer for the counter-current case.

Finally, instead of applying the maximum possible heat transfer, an Efficiency parameter allows specifying a fraction of the maximum possible heat tranfer to be applied. For dimensionless parameters, by CAPE-OPEN convention, difference=True indicates a fraction.

Parameters cannot be used only to allow the user to configure the Unit Operation; they can also be used to report solution values. These types of parameters are OUTPUT parameters, and the Unit Operation reports out the maximum and actual heat transfer, both of which have the dimension of WATT, with W = kg*m2/s3.

For INPUT parameters, we can use the value; for OUTPUT parameters, we must set the value.

Example 7: Added options to the heat exchanger of Example 6.

from scipy import optimize #Heat Exchanger def Configure(unit): #material ports unit.add_port("Feed 1",unit.PortType.MATERIAL,unit.PortDirection.FEED) unit.add_port("Product 1",unit.PortType.MATERIAL,unit.PortDirection.PRODUCT) unit.add_port("Feed 2",unit.PortType.MATERIAL,unit.PortDirection.FEED) unit.add_port("Product 2",unit.PortType.MATERIAL,unit.PortDirection.PRODUCT) #allow inconsistent compound slates unit.allow_inconsistent_compounds=True #input parameters PASCAL=(unit.Dimension.KILOGRAM*unit.Dimension.METER/unit.Dimension.SECOND**2)/unit.Dimension.METER**2 unit.add_parameter("Pressure drop 1",unit.ParameterType.REAL,unit.ParameterMode.INPUT,default=0,min=0,dimension=PASCAL,difference=True) unit.add_parameter("Pressure drop 2",unit.ParameterType.REAL,unit.ParameterMode.INPUT,default=0,min=0,dimension=PASCAL,difference=True) unit.add_parameter("Mode",unit.ParameterType.STRING,unit.ParameterMode.INPUT,default="Counter-current",options=["Counter-current","Co-current"]) unit.add_parameter("Effectiveness",unit.ParameterType.REAL,unit.ParameterMode.INPUT,default=1,min=0,max=1,difference=True) #difference is True for dimensionless implies fraction #output parameters WATT=unit.Dimension.KILOGRAM*unit.Dimension.METER**2/unit.Dimension.SECOND**3 unit.add_parameter("Max. heat transfer",unit.ParameterType.REAL,unit.ParameterMode.OUTPUT,dimension=WATT) unit.add_parameter("Actual heat transfer",unit.ParameterType.REAL,unit.ParameterMode.OUTPUT,dimension=WATT) def Validate(unit): compoundSlate1=[comp.name for comp in unit.ports["Feed 1"].compounds] compoundSlate2=[comp.name for comp in unit.ports["Product 1"].compounds] if compoundSlate1!=compoundSlate2: raise RuntimeError('Feed 1 and Product 1 do not have the same compounds') compoundSlate1=[comp.name for comp in unit.ports["Feed 2"].compounds] compoundSlate2=[comp.name for comp in unit.ports["Product 2"].compounds] if compoundSlate1!=compoundSlate2: raise RuntimeError('Feed 2 and Product 2 do not have the same compounds') def Calculate(unit): fd1=unit.ports["Feed 1"] fd2=unit.ports["Feed 2"] pd1=unit.ports["Product 1"] pd2=unit.ports["Product 2"] p1=fd1.p-unit.parameters["Pressure drop 1"].value p2=fd2.p-unit.parameters["Pressure drop 2"].value if fd1.flow_rate==0: #zero flow case, no heat transfer possible; bring stream 1 to feed 2 temperature pd1.set(x=fd1.x,p=p1,flow_rate=fd1.flow_rate,t=fd2.t) # do a PH flash in case stream 2 is a pure compound pd2.set(x=fd2.x,p=p2,flow_rate=fd2.flow_rate,h=fd2.get_property('enthalpy')) #report zero heat transfer unit.parameters["Max. heat transfer"].value=0 unit.parameters["Actual heat transfer"].value=0 elif fd2.flow_rate==0: #zero flow case, no heat transfer possible; bring stream 2 to feed 1 temperature pd2.set(x=fd2.x,p=p2,flow_rate=fd2.flow_rate,t=fd1.t) # do a PH flash in case stream 1 is a pure compound pd1.set(x=fd1.x,p=p1,flow_rate=fd1.flow_rate,h=fd1.get_property('enthalpy')) #report zero heat transfer unit.parameters["Max. heat transfer"].value=0 unit.parameters["Actual heat transfer"].value=0 else: #store feed enthalpies: h_fd1=fd1.get_property('enthalpy') h_fd2=fd2.get_property('enthalpy') #work required to bring stream 2 to feed temperature 1 q1=(pd2.create_phase_equilibrium(x=fd2.x,t=fd1.t,p=p2).get_property('enthalpy')-fd2.get_property('enthalpy'))*fd2.flow_rate #work required to bring stream 1 to feed temperature 2 q2=(pd1.create_phase_equilibrium(x=fd1.x,t=fd2.t,p=p1).get_property('enthalpy')-fd1.get_property('enthalpy'))*fd1.flow_rate #take the smaller absolute heat duty, which is the max. possible heat duty q=-q1 if abs(q1)<abs(q2) else q2 if unit.parameters["Mode"]=="Counter-current": #Counter-current: the max amount of heat transfer taking place is q pass else: #Co-current: the amount of heat transfer taking place is between 0 and q, # solve for the actual max heat transfer so that exit temperatures are equal if q<0: qmin=q qmax=0 else: qmin=0 qmax=q #solution is bracketed, use a 1-dimensional root finder def f(q): #return temperature difference, given q return pd2.create_phase_equilibrium(x=fd2.x,enthalpy=h_fd2-q/fd2.flow_rate,p=p2).t-pd1.create_phase_equilibrium(x=fd1.x,enthalpy=h_fd1+q/fd1.flow_rate,p=p1).t q=optimize.toms748(f,qmin,qmax) #we found max q, report and calculate actual q unit.parameters["Max. heat transfer"].value=q q=q*unit.parameters["Effectiveness"].value unit.parameters["Actual heat transfer"].value=q #apply to both streams pd1.set(x=fd1.x,p=p1,flow_rate=fd1.flow_rate,enthalpy=h_fd1+q/fd1.flow_rate) pd2.set(x=fd2.x,p=p2,flow_rate=fd2.flow_rate,enthalpy=h_fd2-q/fd2.flow_rate)

In the above example, a STRING parameter was used with a predefined set of options. The control to select such an option is typically rendered as a drop down list by the PME. If no option list is specified, a STRING parameter can take any value, and is typically configured using a text box.

In the Python Unit Operation edit window, one can edit an INPUT parameter value by clicking on it. Any Python expression can be entered here:

Python Unit Operation editor

In addition to REAL and STRING parameters, the Python Unit Operation supports INTEGER, BOOLEAN and REAL_ARRAY parameters.

Information streams

As an alternative to parameters, one can use information streams. Most PMEs will allow Unit Operation parameters to be used in a computational manner in any case, but sometimes it makes sense to explicitly state that a specific piece of information must flow from one Unit Operation to another, and should therefore appear on an information stream. An information stream typically contains a single scalar real parameter, and that is exactly what is supported by the Python Unit Operation.

This example demonstrates a Unit Operation that measures the temperature difference between the actual stream temperature and the dew point temperature, for the purpose of exposing this information via an information stream to a controller.

Example 8: Expose temperature deviation from dew point via information port.

#Expose temperature deviation from dew point temperature def Configure(unit): unit.add_port("Feed",unit.PortType.MATERIAL,unit.PortDirection.FEED) unit.add_port("Product",unit.PortType.MATERIAL,unit.PortDirection.PRODUCT,optional=True) unit.add_port("DewPointTemperatureDeviation",unit.PortType.INFORMATION,unit.PortDirection.PRODUCT,dimension=unit.Dimension.KELVIN,difference=True) def Calculate(unit): #copy feed to product fd=unit.ports["Feed"] if "Product" in unit.ports: #else not connected, it is optional unit.ports["Product"].copy_from(fd) #dew point dew_point_t=fd.create_phase_equilibrium(x=fd.x,p=fd.p,vaporFraction=1).t unit.ports["DewPointTemperatureDeviation"].value=fd.t-dew_point_t

Keep in mind that not all PMEs support information streams.

Energy streams

Our simple heater in Example 1 used an INPUT parameter to specify the heat duty. Alternatively, an ENERGY port can be used.

Generic ENERGY energy streams transfer work, J/s.

An ENERGY FEED stream is an input specification. This does not mean that energy must be consumed; work with a negative sign means the input specification is on how much energy is produced. For a ENERGY PRODUCT stream, a positive sign means energy produced, a negative sign means energy consumed.

A thermal=True ENERGY port transfers heat. It does so within temperature_range which is available from a FEED and must be specified for a PRODUCT.

An axle=True ENERGY port transfers work via a rotating axle. It does so with axis_frequency which is available from a FEED and must be specified for a PRODUCT.

Here is a revision of Example 1, using an ENERGY port:

Example 9: Heater, with thermal energy feed port.

import warnings import math def Configure(unit): #add a feed and product unit.add_port("Feed",unit.PortType.MATERIAL,unit.PortDirection.FEED) unit.add_port("Product",unit.PortType.MATERIAL,unit.PortDirection.PRODUCT) #add a thermal energy inlet port WATT=unit.Dimension.KILOGRAM*unit.Dimension.METER**2/unit.Dimension.SECOND**3 unit.add_port("Heat Duty",unit.PortType.ENERGY,unit.PortDirection.FEED,thermal=True) def Calculate(unit): #copy Feed to Product, apply heat duty fd=unit.ports["Feed"] e=unit.ports["Heat Duty"] duty=e.work if fd.flow_rate==0 and duty!=0: raise RuntimeError('Zero feed flow rate; cannot apply duty') h=fd.get_property('enthalpy') #molar enthalpy if (duty!=0): h+=duty/fd.flow_rate pr=unit.ports["Product"] pr.set(x=fd.x,p=fd.p,flow_rate=fd.flow_rate,enthalpy=h) #check temperature range (t_min_e,t_max_e)=e.temperature_range if not math.isnan(t_min_e): #check against operating range t_out=pr.create_phase_equilibrium(x=fd.x,p=fd.p,enthalpy=h).t if duty>0: #heating if t_out>t_max_e: warnings.warn('heating up product stream to '+str(t_out)+' K is not feasible with heat stream with max. temperature of '+str(t_max_e)+' K') else: #cooling if t_out<t_min_e: warnings.warn('cooling down product stream to '+str(t_out)+' K is not feasible with heat stream with min. temperature of '+str(t_min_e)+' K')

Note that in addition to the heating or cooling of the product stream, the unit checks whether the temperature range on the ENERGY FEED stream is suitable for performing the heating or cooling. If the temperature range is not feasible, the unit produces a warning. Note that warnings issues via warnings.warn are automatically forwarded to the simulation environment:

Calculate warning displayed by COFE

An alternative arrangement is that the Product stream is heated to a specified temperature, and the required heat is reported as ENERGY PRODUCT stream.

Example 10: Heater, with thermal energy feed port.

def Configure(unit): #add a feed and product unit.add_port("Feed",unit.PortType.MATERIAL,unit.PortDirection.FEED) unit.add_port("Product",unit.PortType.MATERIAL,unit.PortDirection.PRODUCT) #outlet temperature parameter unit.add_parameter("Temperature",unit.ParameterType.REAL,unit.ParameterMode.INPUT,dimension=unit.Dimension.KELVIN) #add a thermal energy outlet port WATT=unit.Dimension.KILOGRAM*unit.Dimension.METER**2/unit.Dimension.SECOND**3 unit.add_port("Heat Duty",unit.PortType.ENERGY,unit.PortDirection.PRODUCT,thermal=True,optional=True) def Calculate(unit): #copy Feed to Product, apply temperature fd=unit.ports["Feed"] pr=unit.ports["Product"] pr.set(x=fd.x,flow_rate=fd.flow_rate,p=fd.p,t=unit.parameters['Temperature'].value) #check optional heat stream if "Heat Duty" in unit.ports: e=unit.ports["Heat Duty"] e.work=-(pr.create_phase_equilibrium(x=fd.x,p=fd.p,t=unit.parameters['Temperature'].value).get_property('enthalpy')-fd.get_property('enthalpy'))*fd.flow_rate if e.work<0: #note sign on product work #heating e.temperature_range=(fd.t,unit.parameters['Temperature'].value) else: #cooling e.temperature_range=(unit.parameters['Temperature'].value,fd.t)

Text reports

Reports are added using add_report during Configure.

The most convenient way to write to a report is to use built-in print function, specifying the report in the file argument:

Example 11: Writing a report.

def Configure(unit): #reports unit.add_report("Test Report") def Calculate(unit): #write to report rep=unit.reports['Test Report'] print("""Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum.""",file=rep)

Putting it all together

The last example is a bit more complicated, and uses material streams, parameters, private storage, custom validation, reporting and the abilty of displaying graphics. Here is the full listing of a multi-stream heat exchanger.

Example 12: Multi-stream shell and tube heat exchanger.

import math from scipy import interpolate from scipy import optimize import numpy as np import matplotlib import matplotlib.pyplot as plt def Configure(unit): #number of streams - not exposed as parameter number_of_tubes=3 unit.private_data['number_of_tubes']=number_of_tubes #ports unit.add_port("Shell Feed",unit.PortType.MATERIAL,unit.PortDirection.FEED) unit.add_port("Shell Product",unit.PortType.MATERIAL,unit.PortDirection.PRODUCT) for stream_index in range(0,number_of_tubes): tube_name="Tube "+str(stream_index+1) unit.add_port(tube_name+" Feed",unit.PortType.MATERIAL,unit.PortDirection.FEED) unit.add_port(tube_name+" Product",unit.PortType.MATERIAL,unit.PortDirection.PRODUCT) #allow inconsistent compound slates unit.allow_inconsistent_compounds=True #operational parameters WATT=unit.Dimension.KILOGRAM*unit.Dimension.METER**2/unit.Dimension.SECOND**3 for stream_index in range(0,number_of_tubes): tube_name="Tube "+str(stream_index+1) unit.add_parameter(tube_name+" Direction",unit.ParameterType.STRING,unit.ParameterMode.INPUT,default="Counter-current",options=['Co-current','Counter-current']) for stream_index in range(0,number_of_tubes): tube_name="Tube "+str(stream_index+1) unit.add_parameter(tube_name+" UA",unit.ParameterType.REAL,unit.ParameterMode.INPUT,default=5000,min=0,dimension=WATT/unit.Dimension.KELVIN) #numerical parameters unit.add_parameter("Interpolation temperature interval",unit.ParameterType.REAL,unit.ParameterMode.INPUT,default=5,min=0.1,max=20,dimension=unit.Dimension.KELVIN,difference=True,description="Max. interval between two points in enthalpty table") unit.add_parameter("Number of grid cells",unit.ParameterType.INTEGER,unit.ParameterMode.INPUT,default=20,min=1,max=200,description='Number of grid cells along the length of the heat exchanger') unit.add_parameter("Tolerance",unit.ParameterType.REAL,unit.ParameterMode.INPUT,default=1e-8,min=1e-15,max=1e-2,description='Relative convergence tolerance') #reports unit.add_report('Temperature profiles') unit.add_image('Temperature graph','CreateTemperatureGraph') def Validate(unit): compoundSlate1=[comp.name for comp in unit.ports["Shell Feed"].compounds] compoundSlate2=[comp.name for comp in unit.ports["Shell Product"].compounds] if compoundSlate1!=compoundSlate2: raise RuntimeError('Shell Feed and Shell Product do not have the same compounds') number_of_tubes=unit.private_data['number_of_tubes'] for stream_index in range(0,number_of_tubes): tube_name="Tube "+str(stream_index+1) compoundSlate1=[comp.name for comp in unit.ports[tube_name+" Feed"].compounds] compoundSlate2=[comp.name for comp in unit.ports[tube_name+" Product"].compounds] if compoundSlate1!=compoundSlate2: raise RuntimeError(tube_name+" Feed and "+tube_name+" Product do not have the same compounds") if 'enthalpy_tables' in unit.volatile_data: del unit.volatile_data['enthalpy_tables'] #thermodynamic data may have changed class enthalpy_interpolator: #A simple spline-based enthalpy interpolation table at constant pressure and composition. def __init__(self,feed,t_min,t_max,delta_temperature): #store these to check that table is up-to-date self.x=feed.x self.p=feed.p self.delta_temperature=delta_temperature #our content exists of an array of [temperature table temperature ... temperature] self.t_start=t_min-10 #allow some change in temperature before we need to rebuild self.t_end=t_max+10 #build list of phase boundaries so that we do not interpolate over a phase boundary phase_boundaries=[] try: b=feed.create_phase_equilibrium(x=feed.x,p=feed.p,vapor_fraction=0) if (b.t>self.t_start) and (b.t<self.t_end): phase_boundaries.append(b) except: pass try: b=feed.create_phase_equilibrium(x=feed.x,p=feed.p,vapor_fraction=1) if (b.t>t_start) and (b.t<t_end): phase_boundaries.append(b) except: pass if len(phase_boundaries)==2 and (phase_boundaries[0].t>phase_boundaries[1].t): phase_boundaries=[phase_boundaries[1],phase_boundaries[0]] #build tables to interpolate between limits and phase boundaries, with specified interval width self.content=[self.t_start] for interval_index in range(0,len(phase_boundaries)+1): t0=self.t_start if interval_index==0 else phase_boundaries[interval_index-1].t t1=self.t_end if interval_index==len(phase_boundaries) else phase_boundaries[interval_index].t n_interval=math.ceil((t1-t0)/delta_temperature) delta_t=(t1-t0)/n_interval t_list=[t0+i*delta_t for i in range(0,n_interval+1)] i_first=0 i_last=len(t_list) h_list=[0 for i in range(i_first,i_last)] if interval_index>0: h_list[0]=phase_boundaries[interval_index-1].get_property('enthalpy') i_first=1 if interval_index<len(phase_boundaries): i_last=i_last-1 h_list[i_last]=phase_boundaries[interval_index].get_property('enthalpy') for i in range(i_first,i_last): h_list[i]=feed.create_phase_equilibrium(x=feed.x,p=feed.p,t=t_list[i]).get_property('enthalpy') order=min(3,len(t_list)-1) self.content.append(interpolate.splrep(t_list,h_list,k=order)) self.content.append(t1) def needs_update(self,feed,t_min,t_max,delta_temperature): #rebuild the table in case any of the inputs changed or t range is not large enough return feed.x!=self.x or self.p!=feed.p or t_min<self.t_start or t_max>self.t_end or delta_temperature!=self.delta_temperature def enthalpy(self,t): #find proper interpolation region index=1 while (index<len(self.content)-2): if (self.content[index+1]>t): break index=index+2 #return interpolated value return interpolate.splev(t,self.content[index]) def d_enthalpy_d_t(self,t): #find proper interpolation region index=1 while index<len(self.content)-2: if self.content[index+1]>t: break index=index+2 #return derivative of interpolated value return interpolate.splev(t,self.content[index],der=1) def WriteProfileReport(profile,stream_info,f): #write a report of the temperature profiles vs dimensionless length width=13+16*len(stream_info) print(''.ljust(width,'-'),file=f) print('Distance'.rjust(13),*[s['name'].rjust(15) for s in stream_info],file=f) print(''.ljust(width,'-'),file=f) print(''.ljust(13),*[('product' if s['counter_current'] else 'feed').rjust(15) for s in stream_info],file=f) inx=1.0/(len(profile[0])-1) r=range(0,len(stream_info)) for i in range(0,len(profile[0])): print('{:13.4f}'.format(inx*i),*['{:15.3f}'.format(profile[j][i]) for j in r],'K',file=f) print(''.ljust(13),*[('feed' if s['counter_current'] else 'product').rjust(15) for s in stream_info],file=f) print(''.ljust(width,'-'),file=f) def Calculate(unit): #get stream pairs streams=[(unit.ports["Shell Feed"],unit.ports["Shell Product"])] number_of_tubes=unit.private_data['number_of_tubes'] number_of_streams=number_of_tubes+1 #including shell for stream_index in range(0,number_of_tubes): tube_name="Tube "+str(stream_index+1) streams.append((unit.ports[tube_name+" Feed"],unit.ports[tube_name+" Product"])) #operate at constant pressure p=[str[0].p for str in streams] #Pa #get operating temperature range t_feed=[str[0].t for str in streams] #K t_min=min(t_feed) t_max=max(t_feed) #set up interpolations delta_temperature=unit.parameters["Interpolation temperature interval"].value enthalpy_tables=unit.volatile_data['enthalpy_tables'] if 'enthalpy_tables' in unit.volatile_data else [] for i in range(0,number_of_streams): fd=streams[i][0] nm='Shell' if i==0 else 'Tube '+str(i) if (len(enthalpy_tables)>i) and not enthalpy_tables[i].needs_update(fd,t_min,t_max,delta_temperature): #does not require update print('Re-using enthalpy tables for '+nm) else: #make a new table print('Updating enthalpy tables for '+nm) t=enthalpy_interpolator(fd,t_min,t_max,delta_temperature) if (len(enthalpy_tables)>i): enthalpy_tables[i]=t else: enthalpy_tables.append(t) enthalpy_tables=enthalpy_tables[0:number_of_streams] #in case we had more streams before unit.volatile_data['enthalpy_tables']=enthalpy_tables #save the tables - they are re-used unless no longer valid #system of equations n_grid=unit.parameters["Number of grid cells"].value n_point=n_grid+1 #include feed point n_eq=number_of_streams*n_grid #pre-compile some information stream_info=[None for i in range(0,number_of_streams)] stream_info[0]={ 'name':"Shell", #name 'counter_current':False, #True for counter-current w.r.t. Shell 'equation_offset':0, #offset of first equation in equation system 't_indices':np.array([n_eq]+[i for i in range(0,n_grid)]), #indices of temperatures in [x]+[t_feed] 'direction':1 #1 for co-current, -1 for counter-current. } offset=n_grid for i in range(0,number_of_tubes): tube_name="Tube "+str(i+1) counter_current=unit.parameters[tube_name+" Direction"].value=='Counter-current' stream_info[i+1]={ 'name':tube_name, #name 'counter_current':counter_current, #True for counter-current w.r.t. Shell 'equation_offset':offset, #offset of first equation in equation system 't_indices': np.array([i+offset for i in range(0,n_grid)]+[n_eq+1+i]) if counter_current else np.array([n_eq+1+i]+[i+offset for i in range(0,n_grid)]), #indices of temperatures in [x]+[t_feed] 'half_cell_ua':0.5*unit.parameters[tube_name+" UA"].value/n_grid, 'ua':unit.parameters[tube_name+" UA"].value, #UA for tube (not available for Shell) 'direction':-1 if counter_current else 1 #1 for co-current, -1 for counter-current. } offset+=n_grid rate=[str[0].flow_rate for str in streams] #flow rates, mol/s #function of which to find root # q in each cell is given by ua times average temperature difference # in the limit of small intervals this matches LMTD, but the average # temperature does not have stability problems that LMTD has and # allows for cross-overs def fun(x): #append feed t to solved for t, so that we can index it all_t=np.append(x,t_feed) #t_indices are indices into this array f=np.ndarray(shape=(n_eq),dtype=float) shell_info=stream_info[0] t_shell=all_t[shell_info['t_indices']] #shell temperatures h=np.array([enthalpy_tables[0].enthalpy(t) for t in t_shell]) #shell enthalpy f_shell=(h[0:n_grid]-h[1:n_point])*rate[0] #shell direction is always 1, subtract q to all tubes f_tubes=np.array([]) #append tube f for each tube for i_tube in range(0,number_of_tubes): index=i_tube+1 tube_info=stream_info[index] t_tube=all_t[tube_info['t_indices']] #tube temperatures h=np.array([enthalpy_tables[index].enthalpy(t) for t in t_tube]) #tube enthalpy q_node=tube_info['half_cell_ua']*(t_shell-t_tube) #the contribution of each grid node to each of its neighbour cells q_tube=q_node[0:n_grid]+q_node[1:n_point] #q in each cell, W f_tube=(h[0:n_grid]-h[1:n_point])*rate[index]*tube_info['direction']+q_tube f_shell-=q_tube f_tubes=np.append(f_tubes,f_tube) return np.append(f_shell,f_tubes) #this is the jacobian # the structure of the jacobian is rather sparse - a sparse jacobian may speed # things up, but makes for a more complex example. A dense jacobian is used. def jac(x): j=np.zeros(shape=(n_eq,n_eq),dtype=float) #append feed t to solved for t, so that we can index it all_t=np.append(x,t_feed) #t_indices are indices into this array shell_info=stream_info[0] t_shell_indices=shell_info['t_indices'] #indices into all_t t_shell=all_t[t_shell_indices] #shell temperatures d_h_d_t=np.array([enthalpy_tables[0].d_enthalpy_d_t(t) for t in t_shell]) #shell d enthalpy /d t #f_shell=(h[1:n_point]-h[0:n_grid])*rate[0] #shell direction is always 1, subtract q to all tubes r=rate[0] for i in range(0,n_grid): t_index=t_shell_indices[i] if t_index<n_eq: j[i,t_index]=r*d_h_d_t[i] #derivative of shell balance w.r.t. shell temperature t_index=t_shell_indices[i+1] if t_index<n_eq: j[i,t_index]=-r*d_h_d_t[i+1] #derivative of shell balance w.r.t. shell temperature for i_tube in range(0,number_of_tubes): index=i_tube+1 tube_info=stream_info[index] t_tube_indices=tube_info['t_indices'] #indices into all_t t_tube=all_t[t_tube_indices] #tube temperatures d_h_d_t=np.array([enthalpy_tables[index].d_enthalpy_d_t(t) for t in t_tube]) #tube d enthalpy /d t #q_node=tube_info['half_cell_ua']*(t_shell-t_tube) #the contribution of each grid node to each of its neighbour cells #q_tube=q_node[0:n_grid]+q_node[1:n_point] #q in each cell, W #f_tube=(h[1:n_point]-h[0:n_grid])*rate[index]*tube_info['direction']+q_tube #f_shell-=q_tube eq_index=tube_info['equation_offset'] r=tube_info['direction']*rate[index] halfua=tube_info['half_cell_ua'] for i in range(0,n_grid): t_index=t_tube_indices[i] if t_index<n_eq: j[eq_index,t_index]=r*d_h_d_t[i]-halfua #derivative of tube balance w.r.t. tube temperature j[i,t_index]+=halfua #derivative of shell balance w.r.t. tube temperature t_index=t_tube_indices[i+1] if t_index<n_eq: j[eq_index,t_index]=-r*d_h_d_t[i+1]-halfua #derivative of tube balance w.r.t. tube temperature j[i,t_index]+=halfua #derivative of shell balance w.r.t. tube temperature t_index=t_shell_indices[i] if t_index<n_eq: j[eq_index,t_index]=halfua #derivative of tube balance w.r.t. shell temperature j[i,t_index]-=halfua #derivative of shell balance w.r.t. shell temperature t_index=t_shell_indices[i+1] if t_index<n_eq: j[eq_index,t_index]=halfua #derivative of tube balance w.r.t. shell temperature j[i,t_index]-=halfua #derivative of shell balance w.r.t. shell temperature eq_index+=1 return j #user specified tolerance on norm(f) # note that additional (systematic) error arises from limited accuracy of the temperature # interpolation which can be lowed by lowering the "Interpolation temperature interval" parameter tol=unit.parameters["Tolerance"].value #go in two rouns - use existing profiles as initial guess, or generate a new gues for use_initial_guess in [True,False]: if use_initial_guess: #use saved profile as initial guess if not 't_profiles' in unit.private_data: continue #initial guess not available profiles=unit.private_data['t_profiles'] if len(profiles)!=number_of_streams or len(profiles[0])!=n_point: continue #previous profiles do not 'fit' current problem size print("Using previous profiles as initial guess") x0=np.ndarray(shape=[n_eq+number_of_streams],dtype=float) #the 'extended' array for i in range(0,number_of_streams): x0[stream_info[i]['t_indices']]=profiles[i] x0=np.delete(x0,range(n_eq,n_eq+number_of_streams),0) #remove the trailing feed temperatures else: #generate new initial guess counter_current_streams=[] #we have a guess routine for mixed co-/counter-current, and one for co-current only co_current_streams=[] for i in range(0,number_of_streams): if stream_info[i]['counter_current']: counter_current_streams.append(i) else: co_current_streams.append(i) if counter_current_streams: #rough guess function co-/counter-current, # solve for Q=UA*<average temperature difference> for each stream # and close the heat balance print("calculating initial guess for counter-current flow") h0=[enthalpy_tables[i].enthalpy(t_feed[i]) for i in range(0,number_of_streams)] def f_init(t_out): f=np.ndarray(shape=(number_of_streams),dtype=float) jac=np.zeros(shape=(number_of_streams,number_of_streams),dtype=float) qtot=0 for i in range(0,number_of_tubes): halfua=stream_info[i+1]['ua']*0.5 q=halfua*(t_feed[0]+t_out[0]-t_feed[i+1]-t_out[i+1]) qtot+=q f[i+1]=q+rate[i+1]*(h0[i+1]-enthalpy_tables[i+1].enthalpy(t_out[i+1])) jac[i+1][i+1]=-halfua-rate[i+1]*enthalpy_tables[i+1].d_enthalpy_d_t(t_out[i+1]) jac[i+1][0]=halfua jac[0][i+1]=-halfua jac[0][0]+=halfua f[0]=qtot+rate[0]*(enthalpy_tables[0].enthalpy(t_out[0])-h0[0]) jac[0,0]+=rate[0]*enthalpy_tables[0].d_enthalpy_d_t(t_out[0]) return (f,jac) res=optimize.root(f_init,t_feed,jac=True,tol=1e-2) if not res.success: print("initial guess failed:",res.message.replace("\n"," ").replace(" "," ")) t0_prod=res.x #no need for accurate solution else: #rough guess function for co-current only # find common exit temperature for # which the heat balance closes. This # is the max heat transfer limit for all # co-current print("calculating initial guess for co-current flow") h_in=0 for i in range(0,number_of_streams): h_in+=rate[i]*enthalpy_tables[i].enthalpy(t_feed[i]) def f_heat_balance(t): balance=-h_in for i in range(0,number_of_streams): balance+=rate[i]*enthalpy_tables[i].enthalpy(t) return balance equal_t=optimize.toms748(f_heat_balance,t_min,t_max,xtol=1e-1) #no need for accurate solution t0_prod=[equal_t for i in range(0,number_of_streams)] #assume linear t profile, initialize t from guessed t0_prod x0=np.array([]) for i in range(0,number_of_streams): if stream_info[i]['counter_current']: dt=(t_feed[i]-t0_prod[i])/n_grid t0=t0_prod[i] else: dt=(t0_prod[i]-t_feed[i])/n_grid t0=t_feed[i]+dt x0=np.append(x0,[t0+dt*i for i in range(0,n_grid)]) #solve print("Solving for",n_eq,"degrees of freedom...") res=optimize.root(fun,x0,jac=jac,tol=tol,options={'xtol':tol}) if (res.success): break #no need for a second round x=res.x #print pofile report all_t=np.append(x,t_feed) profile=[all_t[s['t_indices']] for s in stream_info] WriteProfileReport(profile,stream_info,unit.reports['Temperature profiles']) unit.private_data['t_profiles']=profile #save solution for plotting and as initial guess #set the product streams for i in range(0,number_of_streams): s=streams[i] t_prod=profile[i][0 if stream_info[i]['counter_current'] else n_grid] h_prod=enthalpy_tables[i].enthalpy(t_prod) s[1].set(flow_rate=rate[i],x=s[0].x,p=p[i],enthalpy=h_prod) #check for solver error if not res.success: #print(res) raise RuntimeError(res.message.replace("\n"," ").replace(" "," ")) print("Converged using",res.nfev,"function evaluations,",res.njev,"Jacobian evaluations.") #prepare matplotlib matplotlib.use('Agg') #avoid gtk dependency issues - we don't need gtk here color=['r','g','b','c','m','y','b','orange'] #series cycle though these colors def CreateTemperatureGraph(unit,filename,width,height): #create a png file with temperature graph at requested width and height if not 't_profiles' in unit.private_data: raise RuntimeError('Temperature profiles are not available') profiles=unit.private_data['t_profiles'] #new matplotlib figure p=plt.figure(figsize=(width*1e-2,height*1e-2),dpi=100) #at 100 dpi this matches specified size in pixels try: ngrid=len(profiles[0]) #number of grid cells x=[i/(ngrid-1.0) for i in range(0,ngrid)] #dimensionless distance for index in range(0,len(profiles)): #temperature plot series plt.plot(x,profiles[index],linestyle='solid', marker='o',color=color[index%len(color)], label='Shell' if index==0 else 'Tube '+str(index)) plt.xlabel('Relative length') plt.ylabel('Temperature / [K]') plt.legend() plt.tight_layout() plt.savefig(filename) finally: plt.close(p) #clean up

The example depends on numpy and scipy for the numerics.

In the Configure function, a configuration parameter is defined: number_of_tubes. Such a parameter cannot be exposed as input parameter, as this changes the collection of ports and parameters; CAPE-OPEN allows such collections to change only while editing the unit operation.

Subsequently, this value is stored in private_date so that it is available during Validate and Calculate. The value is also used during configuration; two ports, feed and product, and two parameters, UA and flow direction, are added for each tube.

The unit is allowed to connect to streams with inconsistent compound slates, allow_inconsistent_compounds=True, and the validation of compounds on the connected streams is analogous to Example 5 and Example 6.

To speed up the calculations, enthalpy interpolation tables are created. For the simplicity of the example, zero pressure drop is assumed, so that enthalpy is a function only of temperature, at constant composition and pressure. Enthalpy is interpolated using spline interpolation from scipy.interpolate for enthalpy calculated at phase equilibrium, at user specified maximum temperature intervals of Interpolation temperature interval. Enthalpy at phase equilibrium vs. temperature is only 0th order continuous at phase boundaries, so as to prevent improper interpolation over the phase boundaries, the dew and bubble points are found using a pressure-vapor fraction flash. If within the specified temperature region, separate interpolation tables are produced on either side of the phase boundary.

Enthalpy tables can be re-used between successive calculations, provided the overall composition did not change, the pressure did not change, and the new temperange range is within bounds of the current table. The needs_update method checks this. In addition, the tables should not be used if the thermodynamic configuration of any of the materials has changed. If this happens, the PME must call Validate.

The enthalpy tables are stored in volatile_data. As opposed to private_data, volatile_data is not persisted between sessions. There would be no point to save the data in persistent storage, as after restoring, Validate will be called as thermdynamic configuration may have changed. At Validate, the tables are removed from volatile_data using del.

The calculation simulates heat transer between the shell and each of the tubes, under the assumption that the heat transfer coefficients are constant. The total U*A for each tube is specified as input parameter. In addition, each tube can be co- or counter-current to the shell, which is also specified as input parameter. Discretizing the heat balances over dimensionless length leads to the following equations for cell j:

initial guess equation

initial guess equation

initial guess equation

initial guess equation

Here, Q and each equation applies to a grid cell, where T applies to each grid node; the number of nodes is the number of cells plus 1. However feed temperatures are not calculated, as they are constant, and are excluded from the variables, so that the number of variables and equations match. All temperatures are organized in the direction of shell feed to shell product, so shell and co-current tubes have forward indexing with the feed temperature at the start, and counter-current tubes have reverse indexing, with feed temperature at the end. Temperature indexing is pre-calculated and stored as t_indices in temporary storage stream_info. The temperatures can then quickly be extracted in the order required for the equations, using t_indices on an array of [x]+[t_feed] where x is all temperatures solved for, and t_feed contains the constant feed temperatures.

The functions are implemented in fun; the Jacobian is rather sparse, but to not make the example overly complex, a dense Jacobian is implemented in jac. The equations are solved using scipy.optimize.root with a tolerance that the user provides as input parameter. Note that this tolerance, applies only to solving the equations. Some additional systemic error rises from the use of interpolation tables. This error can be reduced by choosing a smaller value for Interpolation temperature interval.

Any profiles stored from an initial guess can of course be used as an initial guess. The solution therefore proceeds in two attempts, due to the for loop over use_initial_guess. In the first attempt, if profiles are available, and if the shape of the profiles (number of tubes, number of grid cells) match the current problem, the initial guess is backed out from the stored profiles.

If this fails, or if appropriate profiles are not available, an initial guess is calculated. In the special case that all tubes are co-current with the shell, the initial guess exit temperatures are solved as one common temperature for which the heat balance closes. This temperature is bracketed between the min. and max. value of all inlet temperatures; these equations are implemented in f_heat_balance and solved using scipy.optimize.toms748.

If at least one tube if counter-current, a rough initial guess is calculated by solving

initial guess equation

initial guess equation

initial guess equation

initial guess equation

The equations are implemented in f_init along with the full Jacobian, and solved using scipy.optimize.root.

After estimating the outlet temperatures, the initial guess is constructed from linear interpolation between inlet and outlet temperatures.

The WriteProfileReport demonstrates how to use the print statement with the file argument set to a report; WriteProfileReport is called with unit.reports['Temperature profiles'] as argument. The reports are not only exposed to the PME, but also shown in the Reports tab:

Profile report in the editor

Finally, matplotlib is used to create a graph showing the temperature profiles as a function of dimensionless distance. Images can be added just like reports, using add_image in Configure. The script must supply a function that creates a PNG file for the image at requested width and height, see CreateTemperatureGraph. Note that, in contrast to creation of reports, images are not created as part of Calculate, but on request, with specified size. The unit should save relevant data to create the image. The image is displayed on the Report tab of the edit window:

Profile report in the editor

Building a library

So far, the module that defines the Unit Operation, its configuration and its calculation, has been typed directly in the editor. This script will be saved as part of the Unit Operation.

Sometimes you want the same Unit Operation definition to be used in multiple places. One way to do this is to copy and paste the script. Which works fine. But it is not very maintenance friendly. If you discover a bug in the script, or want to modify the calculation, all Unit Operations that are defined by the same script must be edited separately to correct for the same issues.

As an alternative, the Unit Operation can be defined from an external module. To demonstrate how this works, the following steps are followed.

Example 13: importing an external unit operation definition.

import os import sys def GetDocumentFolder(): #returns the My Documents folder for the current user import ctypes.wintypes pathbuffer=ctypes.create_unicode_buffer(ctypes.wintypes.MAX_PATH) ctypes.windll.shell32.SHGetFolderPathW(None,5,None,0,pathbuffer) return pathbuffer.value def GetUnitOperationModelFolder(): #assumes the My Documents folder has a sub-folder Python Unit Operation Models return os.path.join(GetDocumentFolder(),"Python Unit Operation Models") def Configure(unit): sys.path.append(GetUnitOperationModelFolder()) #add Python Unit Operation Models to the path import multistreamhtx #load the module of interest unit.define_from(multistreamhtx) #define the unit operation from this module

Of course if the folder that contains the Unit Operation models is already in the Python path, the script becomes a lot simpler:

Example 14: importing an external unit operation definition in the path.

import multistreamhtx def Configure(unit): unit.define_from(multistreamhtx)