Source code for idaes.core.holdup

##############################################################################
# Institute for the Design of Advanced Energy Systems Process Systems
# Engineering Framework (IDAES PSE Framework) Copyright (c) 2018, by the
# software owners: The Regents of the University of California, through
# Lawrence Berkeley National Laboratory,  National Technology & Engineering
# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia
# University Research Corporation, et al. All rights reserved.
# 
# Please see the files COPYRIGHT.txt and LICENSE.txt for full copyright and
# license information, respectively. Both files are also available online
# at the URL "https://github.com/IDAES/idaes".
##############################################################################
"""
Core IDAES holdup models.
"""
from __future__ import division

# Import Python libraries
import logging

# Import Pyomo libraries
from pyomo.environ import Constraint, Param, Reals, TransformationFactory, Var
from pyomo.dae import ContinuousSet, DerivativeVar
from pyomo.common.config import ConfigValue, In

# Import IDAES cores
from idaes.core import ProcessBlockData, declare_process_block_class
from idaes.core.util.config import is_parameter_block, list_of_floats
from idaes.core.util.misc import add_object_ref

__author__ = "Andrew Lee, Jaffer Ghouse"

# TODO
# Add support for other forms of pressure (momentum) and energy balances
# Support non-constant volume for Holdup1D

# Set up logger
logger = logging.getLogger('idaes.unit_model.holdup')


# Set up general ConfigBlock for Holdups and Unit Models based on Holdups
CONFIG = ProcessBlockData.CONFIG()
CONFIG.declare("dynamic", ConfigValue(
        default='use_parent_value',
        domain=In(['use_parent_value', True, False]),
        description="Dynamic model flag",
        doc="""Indicates whether this model will be dynamic,
**default** - 'use_parent_value'.  **Valid values:** {
**'use_parent_value'** - get flag from parent,
**True** - set as a dynamic model,
**False** - set as a steady-state model}"""))
CONFIG.declare("include_holdup", ConfigValue(
        default='use_parent_value',
        domain=In(['use_parent_value', True, False]),
        description="Holdup construction flag",
        doc="""Indicates whether holdup terms should be constructed or not.
Must be True if dynamic = True, **default** - 'use_parent_value'.
**Valid values:** {
**'use_parent_value`** - get flag from parent (default = True),
**True** - construct holdup terms,
**False** - do not construct holdup terms}"""))
CONFIG.declare("material_balance_type", ConfigValue(
        default='use_parent_value',
        domain=In(['use_parent_value', 'none', 'component_phase',
                   'component_total', 'element_total']),
        description="Material balance construction flag",
        doc="""Indicates what type of mass balance should be constructed,
**default** - 'use_parent_value'. **Valid values:** {
**'use_parent_value'** - get flag from parent (default = 'component'),
**'none'** - exclude material balances,
**'component_phase' - use phase component balances,
**'component_total' - use total component balances,
**'element_total' - use total element balances.}"""))
CONFIG.declare("energy_balance_type", ConfigValue(
        default='use_parent_value',
        domain=In(['use_parent_value', 'none', 'enthalpy_total']),
        description="Energy balance construction flag",
        doc="""Indicates what type of energy balance should be constructed,
**default** - 'use_parent_value'. **Valid values:** {
**'use_parent_value'** - get flag from parent (default = 'enthalpy_total'),
**'none'** - exclude energy balances,
**'enthalpy_total'** - single ethalpy balance for material.}"""))
CONFIG.declare("momentum_balance_type", ConfigValue(
        default='use_parent_value',
        domain=In(['use_parent_value', 'none', 'pressure']),
        description="Momentum balance construction flag",
        doc="""Indicates what type of momentum balance should be constructed,
**default** - 'use_parent_value'. **Valid values:** {
**'use_parent_value'** - get flag from parent (default = 'pressure'),
**'none'** - exclude momentum balances,
**'pressure'** - single pressure balance for material.}"""))
CONFIG.declare("has_rate_reactions", ConfigValue(
        default='use_parent_value',
        domain=In(['use_parent_value', True, False]),
        description="Rate reaction construction flag",
        doc="""Indicates whether terms for rate controlled reactions should be
constructed, **default** - 'use_parent_value'. **Valid values:** {
**'use_parent_value'** - get flag from parent (default = False),
**True** - include kinetic reaction terms,
**False** - exclude kinetic reaction terms.}"""))
CONFIG.declare("has_equilibrium_reactions", ConfigValue(
        default='use_parent_value',
        domain=In(['use_parent_value', True, False]),
        description="Equilibrium reaction construction flag",
        doc="""Indicates whether terms for equilibrium controlled reactions
should be constructed, **default** - 'use_parent_value'. **Valid values:** {
**'use_parent_value'** - get flag from parent (default = False),
**True** - include equilibrium reaction terms,
**False** - exclude equilibrium reaction terms.}"""))
CONFIG.declare("has_phase_equilibrium", ConfigValue(
        default='use_parent_value',
        domain=In(['use_parent_value', True, False]),
        description="Phase equilibrium construction flag",
        doc="""Indicates whether terms for phase equilibrium should be
            constructed (default = 'use_parent_value')
                - 'use_parent_value' - get flag from parent (default = False)
                - True - include phase equilibrium terms
                - False - exclude phase equilibrium terms"""))
CONFIG.declare("has_mass_transfer", ConfigValue(
        default='use_parent_value',
        domain=In(['use_parent_value', True, False]),
        description="Mass transfer term construction flag",
        doc="""Indicates whether terms for mass transfer should be constructed,
**default** - 'use_parent_value'. **Valid values:** {
**'use_parent_value'** - get flag from parent (default = False),
**True** - include mass transfer terms,
**False** - exclude mass transfer terms.}"""))
CONFIG.declare("has_heat_transfer", ConfigValue(
        default='use_parent_value',
        domain=In(['use_parent_value', True, False]),
        description="Heat transfer term construction flag",
        doc="""Indicates whether terms for heat transfer should be constructed,
**default** - 'use_parent_value'. **Valid values:** {
**'use_parent_value'** - get flag from parent (default = False),
**True** - include heat transfer terms,
**False** - exclude heat transfer terms.}"""))
CONFIG.declare("has_work_transfer", ConfigValue(
        default='use_parent_value',
        domain=In(['use_parent_value', True, False]),
        description="Work transfer term construction flag",
        doc="""Indicates whether terms for work transfer should be constructed,
**default** - 'use_parent_value'. **Valid values** {
**'use_parent_value'** - get flag from parent (default = False),
**True** - include work transfer terms,
**False** - exclude work transfer terms.}"""))
CONFIG.declare("has_pressure_change", ConfigValue(
        default='use_parent_value',
        domain=In(['use_parent_value', True, False]),
        description="Pressure change term construction flag",
        doc="""Indicates whether terms for pressure change should be
constructed, **default** - 'use_parent_value'. **Valid values:** {
**'use_parent_value'** - get flag from parent (default = False),
**True** - include pressure change terms,
**False** - exclude pressure change terms.}"""))


[docs]@declare_process_block_class("Holdup", doc=""" Holdup is a specialized Pyomo block for IDAES holdup blocks, and contains instances of HoldupData. In most cases, users will want to use one of the derived holdup classes (Holdup0D, Holdup1D or HoldupStatic for their models, however this class can be used to create the framework of a holdup block to which the user can then add custom constraints.""") class HoldupData(ProcessBlockData): """ The HoldupData Class forms the base class for all IDAES holdup models. The purpose of this class is to automate the tasks common to all holdup blockss and ensure that the necessary attributes of a holdup block are present. The most signfiicant role of the Holdup class is to set up the build arguments for the holdup block, automatically link to the time domain of the parent block, and to get the information about the property package. """ CONFIG = CONFIG() CONFIG.declare("property_package", ConfigValue( default='use_parent_value', domain=is_parameter_block, description="Property package to use for holdup", doc="""Property parameter object used to define property calculations, **default** - 'use_parent_value'. **Valid values:** { **'use_parent_value'** - get package from parent (default = None), **a ParameterBlock object**.}""")) CONFIG.declare("property_package_args", ConfigValue( default='use_parent_value', description="Arguments to use for constructing property packages", doc="""A dict of arguments to be passed to a property block and used when constructing these, **default** - 'use_parent_value'. **Valid values:** { **'use_parent_value'** - get package from parent (default = None), **dict** - see property package for documentation.}"""))
[docs] def build(self): """ General build method for Holdup blocks. This method calls a number of sub-methods which automate the construction of expected attributes of all Holdup blocks. Inheriting models should call `super().build`. Args: None Returns: None """ # Check construction flags self._inherit_default_build_arguments() # Call UnitBlockData.build to setup dynamics self._setup_dynamics() # Get property pacakge details self.get_property_package()
def _inherit_default_build_arguments(self): """ Method to collect build arguments from parent blocks as required. If an argument is set as 'use_parent_value', this method attempts to set the argument to that of the parent model, otherwise a default value is used. Args: None Returns: None """ parent = self.parent_block() build_arguments = {'property_package': None, 'property_package_args': {}, 'include_holdup': True, 'material_balance_type': 'component_phase', 'energy_balance_type': 'enthalpy_total', 'momentum_balance_type': 'pressure', 'has_rate_reactions': False, 'has_equilibrium_reactions': False, 'has_phase_equilibrium': False, 'has_mass_transfer': False, 'has_heat_transfer': False, 'has_work_transfer': False, 'has_pressure_change': False } for arg in build_arguments: # If argument is 'use_parent_value', try to get from parent if getattr(self.config, arg) is 'use_parent_value': try: setattr(self.config, arg, getattr(parent.config, arg)) except AttributeError: # If parent does not have flag, resort to defaults self.config[arg] = build_arguments[arg] def _setup_dynamics(self): """ This method automates the setting of the dynamic flag and time domain for holdup blocks. If dynamic flag is 'use_parent_value', method attempts to get the value of the dynamic flag from the parent model, otherwise the local vlaue is used. The time domain is aloways collected from the parent model. Finally, if dynamic = True, the include_holdup flag is checked to ensure it is also True. Args: None Returns: None """ # Check the dynamic flag, and retrieve if necessary if self.config.dynamic == 'use_parent_value': # Get dynamic flag from parent try: self.config.dynamic = self.parent_block().config.dynamic except AttributeError: # If parent does not have dynamic flag, raise Exception raise AttributeError('{} has a parent model ' 'with no dynamic attribute.' .format(self.name)) # Try to get reference to time object from parent try: object.__setattr__(self, "time", self.parent_block().time) except AttributeError: raise AttributeError('{} has a parent model ' 'with no time domain'.format(self.name)) # Check include_holdup, if present if self.config.dynamic: if not self.config.include_holdup: # Dynamic model must have include_holdup = True logger.warning('{} Dynamic holdup blocks must have ' 'include_holdup = True. ' 'Overwritting argument.' .format(self.name)) self.config.include_holdup = True
[docs] def get_property_package(self): """ This method gathers the necessary information about the property package to be used in the holdup block. If a property package has not been provided by the user, the method searches up the model tree until it finds an object with the 'default_property_package' attribute and uses this package for the holdup block. The method also gathers any default construction arguments specified for the property package and combines these with any arguments specified by the user for the holdup block (user specified arguments take priority over defaults). Args: None Returns: None """ # Get property_package block if not provided in arguments parent = self.parent_block() if self.config.property_package in (None, 'use_parent_value'): # Try to get property_package from parent if hasattr(parent.config, "property_package"): if parent.config.property_package is None: parent.config.property_package = \ self._get_default_prop_pack() self.config.property_package = parent.config.property_package else: self.config.property_package = self._get_default_prop_pack() # Get module of property package self.property_module = self.config.property_package.property_module # Check for any flowsheet level build arguments try: # If flowsheet arguments exist, merge with local arguments # Local arguments take precedence arg_dict = self.config.property_package.config.default_arguments arg_dict.update(self.config.property_package_args) self.config.property_package_args = arg_dict except AttributeError: # Otherwise, just use local arguments pass
def _get_default_prop_pack(self): """ This method is used to find a default property package defined at the flowsheet level if a package is not provided as an argument when instantiating the holdup block. Args: None Returns: None """ parent = self.parent_block() while True: if hasattr(parent.config, "default_property_package"): break else: if parent.parent_block() is None: raise AttributeError( '{} no property package provided and ' 'no default defined. Found end of ' 'parent tree.'.format(self.name)) elif parent.parent_block() == parent: raise ValueError( '{} no property package provided and ' 'no default defined. Found recursive ' 'loop in parent tree.'.format(self.name)) parent = parent.parent_block() logger.info('{} Using default property package' .format(self.name)) if parent.config.default_property_package is None: raise ValueError('{} no default property package has been ' 'specified at flowsheet level (' 'default_property_package = None)' .format(self.name)) return parent.config.default_property_package def _get_indexing_sets(self): """ This method collects all necessary indexing sets from property parameter block and makes references to these for use within the holdup block. Collected indexing sets are: phase_list, component_list, rate_reaction_idx, equilibrium_reaction_idx, and element_list. Args: None Returns: None """ # Get outlet property block object try: prop_block = self.properties_out[0] except AttributeError: try: prop_block = self.properties[0, self.ldomain.last()] except AttributeError: prop_block = self.properties[0] # Get phase and component list(s) try: add_object_ref(self, "phase_list", prop_block.phase_list) add_object_ref(self, "component_list", prop_block.component_list) except AttributeError: raise AttributeError('{} property_package provided does not appear' ' to be a valid property package (does not ' 'contain a component_list and/or phase_list)' .format(self.name)) # Get reaction indices and stoichiometry if self.config.has_rate_reactions: try: add_object_ref(self, "rate_reaction_idx", prop_block.rate_reaction_idx) add_object_ref(self, "rate_reaction_stoichiometry", prop_block.rate_reaction_stoichiometry) except AttributeError: self.config.has_rate_reactions = False logger.info('{} property package does not support kinetic ' 'reactions. has_rate_reactions set to False.' .format(self.name)) if self.config.has_equilibrium_reactions: try: add_object_ref(self, "equilibrium_reaction_idx", prop_block.equilibrium_reaction_idx) add_object_ref(self, "equilibrium_reaction_stoichiometry", prop_block.equilibrium_reaction_stoichiometry) except AttributeError: self.config.has_equilibrium_reactions = False logger.info('{} property package does not support equilibrium ' 'reactions. has_equilibrium_reactions set to ' 'False.'.format(self.name)) if self.config.has_phase_equilibrium: try: add_object_ref(self, "phase_equilibrium_idx", prop_block.phase_equilibrium_idx) add_object_ref(self, "phase_equilibrium_list", prop_block.phase_equilibrium_list) except AttributeError: self.config.has_phase_equilibrium = False logger.info('{} property package does not support phase ' 'equilibrium has_phase_equilibrium set to ' 'False.'.format(self.name)) # If element balances, check properties for list of elements if self.config.material_balance_type == 'element_total': try: add_object_ref(self, "element_list", prop_block.element_list) except AttributeError: raise AttributeError("{} Selected property package does not " "support elemental mass balances" .format(self.name))
[docs]@declare_process_block_class("Holdup0D", doc=""" Holdup0D is a specialized Pyomo block for IDAES non-discretized holdup blocks, and contains instances of Holdup0dData. Holdup0D should be used for any holdup with a defined volume and distinct inlets and outlets which does not require spatial discretization. This encompases most basic unit models used in process modeling.""") class Holdup0dData(HoldupData): """ 0-Dimensional (Non-Discretised) Holdup Class This class forms the core of all non-discretized IDAES models. It builds property blocks and adds mass, energy and momentum balances. The form of the terms used in these constraints is specified in the chosen property package. """
[docs] def build(self): """ Build method for Holdup0D blocks. This method calls submethods to setup the necessary property blocks, distributed variables, material, energy and momentum balances based on the arguments provided by the user. Args: None Returns: None """ # Call build method from base class super(Holdup0dData, self).build() # Construct Property Blocks self._make_prop_blocks() # Get indexing sets from property parameters block self._get_indexing_sets() # Build Common Variables self._make_vars() # Make phase fraction variables if necessary if self.config.include_holdup: self._make_phase_frac() # Build material balances if necessary if self.config.material_balance_type in \ ["component_phase", "component_total"]: self._make_component_balances() elif self.config.material_balance_type == "element_total": self._make_element_total_balances() # Build energy balances if necessary if self.config.energy_balance_type == "enthalpy_total": self._make_total_energy_balance() elif self.config.energy_balance_type is not "none": raise ValueError("{} Unrecognised energy_balance_type. This is " "probably due to a failure of the Holdup block " "to inherit a value from its parent. Current " "value of energy_balance_type = {}" .format(self.name, self.config.energy_balance_type)) # Build momentum balances if necessary if self.config.momentum_balance_type == "pressure": self._make_pressure_balance() elif self.config.momentum_balance_type is not "none": raise ValueError("{} Unrecognised momentum_balance_type. This is " "probably due to a failure of the Holdup block " "to inherit a value from its parent. Current " "value of momentum_balance_type = {}" .format(self.name, self.config.momentum_balance_type))
def _make_prop_blocks(self): """ This method constructs the inlet and outlet property blocks for the holdup block. Args: None Returns: None """ self.properties_in = self.property_module.PropertyBlock( self.time, doc="Material properties at inlet", has_sum_fractions=False, calculate_equilibrium_reactions= self.config.has_equilibrium_reactions, calculate_phase_equilibrium=self.config.has_phase_equilibrium, parameters=self.config.property_package, **self.config.property_package_args) self.properties_out = self.property_module.PropertyBlock( self.time, doc="Material properties at outlet", has_sum_fractions=True, calculate_equilibrium_reactions= self.config.has_equilibrium_reactions, calculate_phase_equilibrium=self.config.has_phase_equilibrium, parameters=self.config.property_package, **self.config.property_package_args) def _make_vars(self): """ This method constructs all the temporal derivative variables in the holdup block, and their associated distributed variables. Args: None Returns: None """ # Get units from property package units = {} for u in ['length', 'holdup', 'amount', 'time', 'energy']: try: units[u] = self.config.property_package.get_package_units()[u] except KeyError: units[u] = '-' if self.config.include_holdup: self.volume = Var(self.time, initialize=1.0, doc='Holdup Volume [{}^3]' .format(units['length'])) # Material holdup and accumulation if (self.config.material_balance_type in ['component_phase', 'component_total']): if self.config.include_holdup: self.material_holdup = Var(self.time, self.phase_list, self.component_list, domain=Reals, doc="Material holdup in unit [{}]" .format(units['holdup'])) if self.config.dynamic is True: self.material_accumulation = DerivativeVar( self.material_holdup, wrt=self.time, doc="Material accumulation in unit [{}/{}]" .format(units['holdup'], units['time'])) elif self.config.material_balance_type == 'element_total': if self.config.include_holdup: self.element_holdup = Var( self.time, self.element_list, domain=Reals, doc="Elemental holdup in unit [{}]" .format(units['amount'])) if self.config.dynamic is True: self.element_accumulation = DerivativeVar( self.element_holdup, wrt=self.time, doc="Elemental accumulation in unit [{}/{}]" .format(units['amount'], units['time'])) # Energy holdup and accumulation if self.config.energy_balance_type == "enthalpy_total": if self.config.include_holdup: self.energy_holdup = Var( self.time, self.phase_list, domain=Reals, doc="Energy holdup in unit [{}]" .format(units['energy'])) if self.config.dynamic is True: self.energy_accumulation = DerivativeVar( self.energy_holdup, wrt=self.time, doc="Energy holdup in unit [{}/{}]" .format(units['energy'], units['time'])) def _make_phase_frac(self): """ This method constructs the phase_fraction variables for the holdup block, and the associated constraint o nthe sum of phase_fraction == 1. For systems with only one phase, phase_fraction is created as a Pyomo Expression with a value of 1. Args: None Returns: None """ if len(self.phase_list) > 1: self.phase_fraction = Var( self.time, self.phase_list, initialize=1/len(self.phase_list), doc='Volume fraction of holdup by phase') @self.Constraint(self.time, doc='Sum of phase fractions == 1') def sum_of_phase_fractions(self, t): return 1 == sum(self.phase_fraction[t, p] for p in self.phase_list) else: @self.Expression(self.time, self.phase_list, doc='Volume fraction of holdup by phase') def phase_fraction(self, t, p): return 1 def _make_component_balances(self): """ This method is called when material_balance_type calls for component balances, and constructs the component balance equations for the holdup block.. Args: None Returns: None """ # Get units from property package units = {} for u in ['holdup', 'time']: try: units[u] = self.config.property_package.get_package_units()[u] except KeyError: units[u] = '-' # Get phase component list(s) if hasattr(self.config.property_package, "phase_component_list"): phase_component_list = ( self.config.property_package.phase_component_list) else: # Otherwise assume all components in all phases phase_component_list = {} for p in self.phase_list: phase_component_list[p] = self.component_list # Create material balance terms as required # Kinetic reaction generation if self.config.has_rate_reactions: self.rate_reaction_generation = Var( self.time, self.phase_list, self.component_list, domain=Reals, doc="Amount of component generated in " "unit by kinetic reactions [{}/{}]" .format(units['holdup'], units['time'])) # Equilibrium reaction generation if self.config.has_equilibrium_reactions: self.equilibrium_reaction_generation = Var( self.time, self.phase_list, self.component_list, domain=Reals, doc="Amount of component generated in unit " "by equilibrium reactions [{}/{}]" .format(units['holdup'], units['time'])) # Phase equilibrium generation if (self.config.has_phase_equilibrium and self.config.material_balance_type != 'component_total'): self.phase_equilibrium_generation = Var( self.time, self.phase_equilibrium_idx, domain=Reals, doc="Amount of generation in unit by phase " "equilibria [{}/{}]" .format(units['holdup'], units['time'])) # Material transfer term if self.config.has_mass_transfer: self.mass_transfer_term = Var( self.time, self.phase_list, self.component_list, domain=Reals, doc="Component material transfer into unit [{}/{}]" .format(units['holdup'], units['time'])) # Create rules to substitute material balance terms # Accumulation term def accumulation_term(b, t, p, j): return b.material_accumulation[t, p, j] if b.config.dynamic else 0 def kinetic_term(b, t, p, j): return (b.rate_reaction_generation[t, p, j] if b.config.has_rate_reactions else 0) def equilibrium_term(b, t, p, j): if self.config.has_equilibrium_reactions: return b.equilibrium_reaction_generation[t, p, j] else: return 0 def phase_equilibrium_term(b, t, p, j): if (self.config.has_phase_equilibrium and self.config.material_balance_type != 'component_total'): sd = {} for r in b.phase_equilibrium_idx: if self.phase_equilibrium_list[r][0] == j: if self.phase_equilibrium_list[r][1][0] == p: sd[r] = 1 elif self.phase_equilibrium_list[r][1][1] == p: sd[r] = -1 else: sd[r] = 0 else: sd[r] = 0 return sum(b.phase_equilibrium_generation[t, r]*sd[r] for r in b.phase_equilibrium_idx) else: return 0 def transfer_term(b, t, p, j): if self.config.has_mass_transfer: return b.mass_transfer_term[t, p, j] else: return 0 # Component balances if self.config.material_balance_type == 'component_phase': @self.Constraint(self.time, self.phase_list, self.component_list, doc="Material balances") def material_balance(b, t, p, j): if j in phase_component_list[p]: return accumulation_term(b, t, p, j) == ( b.properties_in[t].material_balance_term[p, j] - b.properties_out[t].material_balance_term[p, j] + kinetic_term(b, t, p, j) + equilibrium_term(b, t, p, j) + phase_equilibrium_term(b, t, p, j) + transfer_term(b, t, p, j)) else: return Constraint.Skip elif self.config.material_balance_type == 'component_total': @self.Constraint(self.time, self.component_list, doc="Material balances") def material_balance(b, t, j): cplist = [] for p in self.phase_list: if j in phase_component_list[p]: cplist.append(p) return ( sum(accumulation_term(b, t, p, j) for p in cplist) == sum(b.properties_in[t].material_balance_term[p, j] for p in cplist) - sum(b.properties_out[t].material_balance_term[p, j] for p in cplist) + sum(kinetic_term(b, t, p, j) for p in cplist) + sum(equilibrium_term(b, t, p, j) for p in cplist) + sum(transfer_term(b, t, p, j) for p in cplist)) else: # This should never trigger. raise ValueError('Unrecognised value for argument ' 'material_balance_type') # TODO: Need to set material_holdup = 0 for non-present component-phase # pairs. Not ideal, but needed to close DoF. Is there a better way? # Material Holdup if self.config.include_holdup: @self.Constraint(self.time, self.phase_list, self.component_list, doc="Material holdup calculations") def material_holdup_calculation(b, t, p, j): if j in phase_component_list[p]: return b.material_holdup[t, p, j] == ( b.volume[t]*self.phase_fraction[t, p] * b.properties_out[t].material_density_term[p, j]) else: return b.material_holdup[t, p, j] == 0 if self.config.has_rate_reactions: # Add extents of reaction and stoichiometric constraints self.rate_reaction_extent = Var( self.time, self.rate_reaction_idx, domain=Reals, doc="Extent of kinetic reactions[{}/{}]" .format(units['holdup'], units['time'])) @self.Constraint(self.time, self.phase_list, self.component_list, doc="Kinetic reaction stoichiometry constraint") def rate_reaction_stoichiometry_constraint(b, t, p, j): if j in phase_component_list[p]: return b.rate_reaction_generation[t, p, j] == ( sum(b.rate_reaction_stoichiometry[r, p, j] * b.rate_reaction_extent[t, r] for r in b.rate_reaction_idx)) else: return Constraint.Skip if self.config.has_equilibrium_reactions: # Add extents of reaction and stoichiometric constraints self.equilibrium_reaction_extent = Var( self.time, self.equilibrium_reaction_idx, domain=Reals, doc="Extent of equilibrium reactions[{}/{}]" .format(units['holdup'], units['time'])) @self.Constraint(self.time, self.phase_list, self.component_list, doc="Equilibrium reaction stoichiometry") def equilibrium_reaction_stoichiometry_constraint(b, t, p, j): if j in phase_component_list[p]: return b.equilibrium_reaction_generation[t, p, j] == ( sum(b.equilibrium_reaction_stoichiometry[r, p, j] * b.equilibrium_reaction_extent[t, r] for r in b.equilibrium_reaction_idx)) else: return Constraint.Skip def _make_element_total_balances(self): """ This method is called when material_balance_type = 'element', and constructs the element balance equations for the holdup block (indexed by element). Args: None Returns: None """ # Get units from property package units = {} for u in ['amount', 'time']: try: units[u] = self.config.property_package.get_package_units()[u] except KeyError: units[u] = '-' # Add Material Balance terms @self.Expression(self.time, self.phase_list, self.element_list, doc="Inlet elemental flow terms [{}/{}]" .format(units['amount'], units['time'])) def elemental_flow_in(b, t, p, e): return sum(b.properties_in[t].material_balance_term[p, j] * b.properties_out[t].element_comp[j][e] for j in b.component_list) @self.Expression(self.time, self.phase_list, self.element_list, doc="Outlet elemental flow terms [{}/{}]" .format(units['amount'], units['time'])) def elemental_flow_out(b, t, p, e): return sum(b.properties_out[t].material_balance_term[p, j] * b.properties_out[t].element_comp[j][e] for j in b.component_list) # Create material balance terms as needed if self.config.has_mass_transfer: self.elemental_mass_transfer = Var( self.time, self.element_list, domain=Reals, doc="Element material transfer into unit [{}/{}]" .format(units['amount'], units['time'])) # Create rules to substitute material balance terms # Accumulation term def accumulation_term(b, t, e): return b.element_accumulation[t, e] if b.config.dynamic else 0 # Mass transfer term def transfer_term(b, t, e): return (b.elemental_mass_transfer[t, e] if b.config.has_mass_transfer else 0) # Element balances @self.Constraint(self.time, self.element_list, doc="Elemental material balances") def element_balance(b, t, e): return accumulation_term(b, t, e) == ( sum(b.elemental_flow_in[t, p, e] for p in b.phase_list) - sum(b.elemental_flow_out[t, p, e] for p in b.phase_list) + transfer_term(b, t, e)) # Elemental Holdup if self.config.include_holdup: @self.Constraint(self.time, self.element_list, doc="Elemental holdup calcuation") def elemental_holdup_calculation(b, t, e): return b.element_holdup[t, e] == ( b.volume[t] * sum(b.phase_fraction[t, p] * b.properties_out[t].material_density_term[p, j] * b.properties_out[t].element_comp[j][e] for p in b.phase_list for j in b.component_list)) def _make_total_energy_balance(self): """ This method is called when energy_balance_type = 'enthalpy_total', and constructs the energy balance equations for the holdup block. Args: None Returns: None """ # Get units from property package units = {} for u in ['energy', 'time']: try: units[u] = self.config.property_package.get_package_units()[u] except KeyError: units[u] = '-' # Create scaling factor self.scaling_factor_energy = Param( default=1e-6, mutable=True, doc='Energy balance scaling parameter') # Create energy balance terms as needed # Heat transfer term if self.config.has_heat_transfer: self.heat = Var(self.time, domain=Reals, initialize=0.0, doc="Heat transfered in unit [{}/{}]" .format(units['energy'], units['time'])) # Work transfer if self.config.has_work_transfer: self.work = Var(self.time, domain=Reals, initialize=0.0, doc="Work transfered in unit [{}/{}]" .format(units['energy'], units['time'])) # Create rules to substitute energy balance terms # Accumulation term def accumulation_term(b, t, p): return b.energy_accumulation[t, p] if b.config.dynamic else 0 def heat_term(b, t): return b.heat[t] if b.config.has_heat_transfer else 0 def work_term(b, t): return b.work[t] if b.config.has_work_transfer else 0 # Energy balance equation @self.Constraint(self.time, doc="Energy balances") def energy_balance(b, t): return (sum(accumulation_term(b, t, p) for p in b.phase_list) * b.scaling_factor_energy) == ( sum(b.properties_in[t].energy_balance_term[p] for p in b.phase_list) * b.scaling_factor_energy - sum(self.properties_out[t].energy_balance_term[p] for p in b.phase_list) * b.scaling_factor_energy + heat_term(b, t)*b.scaling_factor_energy + work_term(b, t)*b.scaling_factor_energy) # Energy Holdup if self.config.include_holdup: @self.Constraint(self.time, self.phase_list, doc="Energy holdup constraint") def energy_holdup_calculation(b, t, p): return b.energy_holdup[t, p] == ( b.volume[t]*self.phase_fraction[t, p] * b.properties_out[t].energy_density_term[p]) def _make_pressure_balance(self): """ This method is called when momentum_balance_type = 'pressure', and constructs the pressure balance equations for the holdup block. Args: None Returns: None """ # Get units from property package try: p_units = \ self.config.property_package.get_package_units()['pressure'] except KeyError: p_units = '-' # Add Momentum Balance Variables as necessary if self.config.has_pressure_change: self.deltaP = Var(self.time, domain=Reals, doc="Pressure difference across unit [{}]" .format(p_units)) # Create rules to substitute energy balance terms # Pressure change term def deltaP_term(b, t): return b.deltaP[t] if b.config.has_pressure_change else 0 # Create scaling factor self.scaling_factor_pressure = Param( default=1e-4, mutable=True, doc='Momentum balance scaling parameter') # Momentum balance equation @self.Constraint(self.time, doc='Momentum balance') def momentum_balance(b, t): return 0 == (b.properties_in[t].pressure * b.scaling_factor_pressure - b.properties_out[t].pressure * b.scaling_factor_pressure + deltaP_term(b, t)*b.scaling_factor_pressure)
[docs] def model_check(blk): """ This method exectues the model_check methods on the associated property blocks (if they exist). This method is generally called by a unit model as part of the unit's model_check method. Args: None Returns: None """ # Try property block model check for t in blk.time: try: blk.properties_in[t].model_check() except AttributeError: logger.warning('{} Holdup inlet property block has no model ' 'check. To correct this, add a model_check ' 'method to the associated PropertyBlock class.' .format(blk.name)) try: blk.properties_out[t].model_check() except AttributeError: logger.warning('{} Holdup outlet property block has no ' 'model check. To correct this, add a ' 'model_check method to the associated ' 'PropertyBlock class.'.format(blk.name))
[docs] def initialize(blk, state_args=None, outlvl=0, optarg=None, solver='ipopt', hold_state=True): ''' Initialisation routine for holdup (default solver ipopt) Keyword Arguments: state_args : a dict of arguments to be passed to the property package(s) to provide an initial state for initialization (see documentation of the specific property package) (default = {}). outlvl : sets output level of initialisation routine. **Valid values:** **0** - no output (default), **1** - return solver state for each step in routine, **2** - include solver output infomation (tee=True) optarg : solver options dictionary object (default=None) solver : str indicating whcih solver to use during initialization (default = 'ipopt') hold_state : flag indicating whether the initialization routine should unfix any state variables fixed during initialization, **default** - True. **Valid values:** **True** - states varaibles are not unfixed, and a dict of returned containing flags for which states were fixed during initialization, **False** - state variables are unfixed after initialization by calling the relase_state method. Returns: If hold_states is True, returns a dict containing flags for which states were fixed during initialization. ''' # Check if holdup has associated inlet_mixer, and initialize if present # Check for associated inlet_mixer block try: # Call inlet_mixer initialization routine mix_flags = blk.inlet_mixer.initialize(outlvl=outlvl-1, optarg=optarg, solver=solver, hold_state=hold_state, state_args=state_args) mixer_flag = True hold_state_in = False except AttributeError: mixer_flag = False hold_state_in = hold_state # Get inlet state if not provided if state_args is None: state_args = {} state_dict = blk.properties_in[0].declare_port_members() for k in state_dict.keys(): if state_dict[k].is_indexed(): state_args[k] = {} for m in state_dict[k].keys(): state_args[k][m] = state_dict[k][m].value else: state_args[k] = state_dict[k].value # Initialize property blocks in_flags = blk.properties_in.initialize(outlvl=outlvl-1, optarg=optarg, solver=solver, hold_state=hold_state_in, **state_args) blk.properties_out.initialize(outlvl=outlvl-1, optarg=optarg, solver=solver, hold_state=False, **state_args) # Check for assoicated outlet block try: # Call outlet_splitter initialization routine blk.outlet_splitter.initialize(outlvl=outlvl-1, optarg=optarg, solver=solver, hold_state=False, state_args=state_args) except AttributeError: pass if outlvl > 0: logger.info('{} Initialisation Complete'.format(blk.name)) if mixer_flag: return mix_flags else: return in_flags
[docs] def release_state(blk, flags, outlvl=0): ''' Method to relase state variables fixed during initialisation. Keyword Arguments: flags : dict containing information of which state variables were fixed during initialization, and should now be unfixed. This dict is returned by initialize if hold_state = True. outlvl : sets output level of of logging Returns: None ''' try: blk.inlet_mixer.release_state(flags, outlvl-1) except AttributeError: blk.properties_in.release_state(flags, outlvl=outlvl-1)
[docs]@declare_process_block_class("Holdup1D") class Holdup1dData(HoldupData): """ 1-Dimensional Holdup Class This class is designed to be the core of all 1D discretized IDAES models. It builds property blocks, inlet/outlet ports and adds mass, energy and momentum balances. The form of the terms used in these constraints is specified in the chosen property package. Assumes constant reactor dimensions """ CONFIG = HoldupData.CONFIG() CONFIG.declare("inherited_length_domain", ConfigValue( default=None, description="Length domain inherited from parent", doc="""A Pyomo ContinuousSet to use as the length domain in the holdup. """ )) CONFIG.declare("length_domain_set", ConfigValue( default=[0.0, 1.0], domain=list_of_floats, description="Set to use to initialize length domain", doc="""List of floats between 0 and 1 to used to initialize length domain, if inherited_length_domain not set (default = [0.0, 1.0] """ )) CONFIG.declare("flow_direction", ConfigValue( default='forward', domain=In(['forward', 'backward']), description="Flow direction flag", doc="""Flag indicating the direction of flow within the length domain (default = `forward`). - 'forward' - flow from 0 to 1 - 'backward' - flow from 1 to 0""")) CONFIG.declare("discretization_method", ConfigValue( default='OCLR', domain=In(['OCLR', 'OCLL', 'BFD', 'FFD']), description="Discretization method to apply to length domain", doc="""Method to be used by DAE transformation when discretizing length domain (default = `OCLR`). - 'OCLR' - orthogonal collocation (Radau roots) - 'OCLL' - orthogonal collocation (Legendre roots) - 'BFD' - backwards finite difference (1st order) - 'FFD' - forwards finite difference (1st order)""")) CONFIG.declare("finite_elements", ConfigValue( default=20, domain=int, description="Number of finite elements length domain", doc="""Number of finite elements to use when discretizing length domain (default=20)""")) CONFIG.declare("collocation_points", ConfigValue( default=3, domain=int, description="Number of collocation points per finite element", doc="""Number of collocation points to use per finite element when discretizing length domain (default=3)""")) CONFIG.declare("has_mass_diffusion", ConfigValue( default=False, domain=In([True, False]), description="Mass diffusion flag", doc="""Flag indicating whether mass diffusion/dispersion should be included in material balance equations (default=False)""")) CONFIG.declare("has_energy_diffusion", ConfigValue( default=False, domain=In([True, False]), description="Energy diffusion flag", doc="""Flag indicating whether energy diffusion/dispersion should be included in energy balance equations (default=False)""")) CONFIG.declare("velocity_type", ConfigValue( default='use_parent_value', domain=In(['use_parent_value', 'none', 'mixture', 'phase']), description="Velocity calculation flag", doc="""Flag indicating whether velocity should be include in the model **default** - 'use_parent_value'. **Valid values:** { **'none'** - do not calculate velocity, **'mixture'** - calculate a single velocity for all phases, **'phase'** - calculate separate velocities for each phase}"""))
[docs] def build(self): """ Build method for Holdup1D blocks. This method calls submethods to setup the necessary property blocks, distributed variables, material, energy and momentum balances based on the arguments provided by the user. Args: None Returns: None """ # Call build method from base class super(Holdup1dData, self).build() # Build the axial domain self._make_domain() # Construct Property Blocks self._make_prop_blocks() # Get indexing sets from property parameters block self._get_indexing_sets() # Build Common Variables self._make_vars() if self.config.velocity_type != 'none': self._make_velocity() # Make phase fraction variables if necessary if self.config.include_holdup: self._make_phase_frac() # Build material balances if necessary if self.config.material_balance_type in \ ["component_phase", "component_total"]: self._make_component_balances() elif self.config.material_balance_type == "element_total": self._make_element_total_balances() # Build energy balances if necessary if self.config.energy_balance_type == "enthalpy_total": self._make_total_energy_balance() # Build momentum balances if necessary if self.config.momentum_balance_type == "pressure": self._make_pressure_balance() # Transform axial domain - not needed if inherited domain provided. if self.config.inherited_length_domain is None: self._transform_domain()
def _inherit_default_build_arguments(self): """ Overloading base method to collect build arguments from parent blocks as required to include additional arguments. If an argument is set as 'use_parent_value', this method attempts to set the argument to that of the parent model, otherwise a default value is used. Args: None Returns: None """ parent = self.parent_block() build_arguments = {'property_package': None, 'property_package_args': {}, 'include_holdup': True, 'material_balance_type': 'component_phase', 'energy_balance_type': 'enthalpy_total', 'momentum_balance_type': 'pressure', 'has_rate_reactions': False, 'has_equilibrium_reactions': False, 'has_phase_equilibrium': False, 'has_mass_transfer': False, 'has_heat_transfer': False, 'has_work_transfer': False, 'has_pressure_change': False, 'velocity_type': 'none', } for arg in build_arguments: # If argument is 'use_parent_value', try to get from parent if getattr(self.config, arg) is 'use_parent_value': try: setattr(self.config, arg, getattr(parent.config, arg)) except AttributeError: # If parent does not have flag, resort to defaults self.config[arg] = build_arguments[arg] def _make_domain(self): ''' Build the axial domain ''' # Get units from property package try: lunits = self.config.property_package.get_package_units()['length'] except KeyError: lunits = '-' # Check type of length_domain argument if self.config.inherited_length_domain is not None: # If ContinuousSet, make ldomain a reference if isinstance(self.config.inherited_length_domain, ContinuousSet): add_object_ref(self, "ldomain", self.config.inherited_length_domain) else: # TODO : This should be validated through config block raise TypeError("{} inherited_length_domain must be a Pyomo " "ContinuousSet".format(self.name)) else: self.ldomain = ContinuousSet( bounds=(0.0, 1.0), initialize=self.config.length_domain_set, doc="Length domain [{}]".format(lunits)) def _make_prop_blocks(self): """ This method constructs the property blocks for the holdup block. Args: None Returns: None """ # Need to set calculate_equilibrium flag at inlet to False # Create new rule for constructing PropertyBlock which changes flag def property_rule(b, t, x): if ((self.config.flow_direction is "forward" and x == self.ldomain.first()) or (self.config.flow_direction is "backward" and x == self.ldomain.last())): # If first element, set flag to False self.properties[t, x].config.has_sum_fractions = False else: # Otherwise, set flag based on construction arguments self.properties[t, x].config.has_sum_fractions = True # Call PropertyBlock build method self.properties[t, x].build() self.properties = self.property_module.PropertyBlock( self.time, self.ldomain, doc="Material properties in holdup", rule=property_rule, calculate_equilibrium_reactions= self.config.has_equilibrium_reactions, calculate_phase_equilibrium=self.config.has_phase_equilibrium, parameters=self.config.property_package, **self.config.property_package_args) def _make_vars(self): ''' Build discretized variables ''' # Get units from property package units = {} for u in ['length', 'holdup', 'amount', 'time', 'energy', 'pressure', 'temperature']: try: units[u] = self.config.property_package.get_package_units()[u] except KeyError: units[u] = '-' # Reactor Design self.area = Var(domain=Reals, initialize=1.0, doc="Reactor area [{}^2]".format(units['length'])) self.length = Var(domain=Reals, initialize=1.0, doc="Reactor length [{}]".format(units['length'])) self.volume = Var(domain=Reals, initialize=1.0, doc="Reactor volume [{}^3]".format(units['length'])) def rule_geometry(b): return b.volume == b.length*b.area self.eq_geometry = Constraint(rule=rule_geometry) # Material flow, holdup and accumulation if self.config.material_balance_type in ['component_phase', 'component_total']: self.material_flow = Var( self.time, self.ldomain, self.phase_list, self.component_list, doc="Axial material flowrate [{}/{}]" .format(units['holdup'], units['time'])) self.material_flow_dx = DerivativeVar(self.material_flow, wrt=self.ldomain) if self.config.include_holdup: self.material_holdup = Var( self.time, self.ldomain, self.phase_list, self.component_list, domain=Reals, doc="Material holdup in unit [{}]" .format(units['holdup'])) if self.config.dynamic is True: self.material_accumulation = DerivativeVar( self.material_holdup, wrt=self.time, doc="Material accumulation in unit [{}/{}]" .format(units['holdup'], units['time'])) elif self.config.material_balance_type == 'element_total': self.elemental_flow = Var( self.time, self.ldomain, self.phase_list, self.element_list, doc="Axial elemental flowrate [{}/{}]" .format(units['amount'], units['time'])) self.elemental_flow_dx = DerivativeVar(self.elemental_flow, wrt=self.ldomain) if self.config.include_holdup: self.element_holdup = Var( self.time, self.ldomain, self.element_list, domain=Reals, doc="Elemental holdup in unit [{}]" .format(units['amount'])) if self.config.dynamic is True: self.element_accumulation = DerivativeVar( self.element_holdup, wrt=self.time, doc="Elemental accumulation in unit [{}/{}]" .format(units['amount'], units['time'])) # Mass diffusion variables if self.config.has_mass_diffusion: # Check that properties support diffusion if (hasattr(self.properties[0, 1], "material_concentration_term") and hasattr(self.properties[0, 1], "diffusivity")): self.material_concentration_term = Var( self.time, self.ldomain, self.phase_list, self.component_list, domain=Reals, doc="Concentration variable for diffusion") self.material_concentration_dx2 = DerivativeVar( self.material_concentration_term, wrt=(self.ldomain, self.ldomain)) else: self.config.has_mass_diffusion = False logger.info('{} property package does not support mass ' 'diffusion. has_mass_diffusion set to False.' .format(self.name)) # Energy holdup and accumulation if self.config.energy_balance_type == "enthalpy_total": self.energy_flow = Var( self.time, self.ldomain, self.phase_list, doc="Axial energy flowrate [{}/{}]" .format(units['energy'], units['time'])) self.energy_flow_dx = DerivativeVar(self.energy_flow, wrt=self.ldomain) if self.config.include_holdup: self.energy_holdup = Var( self.time, self.ldomain, self.phase_list, domain=Reals, doc="Energy holdup in unit [{}]" .format(units['energy'])) if self.config.dynamic is True: self.energy_accumulation = DerivativeVar( self.energy_holdup, wrt=self.time, doc="Energy accumulation in unit [{}/{}]" .format(units['energy'], units['time'])) # Energy diffusion variables if self.config.has_energy_diffusion is True: # Check that properties support diffusion if (hasattr(self.properties[0, 1], "temperature") and hasattr(self.properties[0, 1], "conductivity")): # Add diffusion transfer variables self.temperature = Var(self.time, self.ldomain, domain=Reals, doc="Temperature [{}]" .format(units['temperature'])) self.temperature_dx2 = DerivativeVar( self.temperature, wrt=(self.ldomain, self.ldomain)) else: self.config.has_energy_diffusion = False logger.info('{} property package does not support energy ' 'diffusion. has_energy_diffusion set to False.' .format(self.name)) # Pressure if self.config.momentum_balance_type == "pressure": self.pressure = Var(self.time, self.ldomain, doc="Pressure [{}]".format(units['pressure'])) self.pressure_dx = DerivativeVar(self.pressure, wrt=self.ldomain) def _make_velocity(self): """ This method constructs the velocity variables and constraints as required. Args: None Returns: None """ # Get units from property package units = {} for u in ['length', 'time', ]: try: units[u] = self.config.property_package.get_package_units()[u] except KeyError: units[u] = '-' if self.config.velocity_type == 'mixture': self.velocity = Var(self.time, self.ldomain, initialize=1, doc="Mixture superficial velocity within " "reactor [{}/{}]" .format(units['length'], units['time'])) try: @self.Constraint(self.time, self.ldomain, doc="Velocity calculation") def velocity_calculation(b, t, x): return (b.velocity[t, x]*b.area == b.properties[t, x].flow_vol) except AttributeError: raise AttributeError( '{} could not construct velocity constraint.' ' This is likely due to the property package' ' not supporting flow_vol.'.format(self.name)) elif self.config.velocity_type == 'phase': self.velocity = Var(self.time, self.phase_list, self.ldomain, initialize=1, doc="Superficial phase velocity within reactor" " [{}/{}]".format(units['length'], units['time'])) try: @self.Constraint(self.time, self.phase_list, self.ldomain, doc="Velocity calculation") def velocity_calculation(b, t, p, x): return (b.velocity[t, p, x]*b.area == b.properties[t, x].flow_vol_phase[p]) except AttributeError: raise AttributeError( '{} could not construct velocity constraint.' ' This is likely due to the property package' ' not supporting flow_vol_phase.'.format(self.name)) else: raise ValueError('{} unrecognised value for argument ' 'velocity_type {}'.format( self.name, self.config.velocity_type)) def _make_phase_frac(self): """ This method constructs the phase_fraction variables for the holdup block, and the associated constraint on the sum of phase_fraction == 1. For systems with only one phase, phase_fraction is created as a Pyomo Expression with a value of 1. Args: None Returns: None """ if len(self.phase_list) > 1: self.phase_fraction = Var(self.time, self.ldomain, self.phase_list, initialize=1/len(self.phase_list), doc='Volume fraction of holdup by phase') @self.Constraint(self.time, self.ldomain, doc="Sum of phase fractions == 1") def sum_of_phase_fractions(b, t, x): return 1 == sum(b.phase_fraction[t, x, p] for p in b.phase_list) else: @self.Expression(self.time, self.ldomain, self.phase_list, doc='Volume fraction of holdup by phase') def phase_fraction(self, t, x, p): return 1 def _flow_direction_term(self): """ This method sets the sign for convective terms according to flow direction. Args: None Returns: None """ if self.config.flow_direction is "forward": return -1 else: return 1 def _make_component_balances(self): """ This method is called when material_balance_type calls for component balances, and constructs the component balance equations for the holdup block.. Args: None Returns: None """ # Get units from property package units = {} for u in ['holdup', 'time', 'length']: try: units[u] = self.config.property_package.get_package_units()[u] except KeyError: units[u] = '-' # Get phase component list(s) if hasattr(self.config.property_package, "phase_component_list"): phase_component_list = ( self.config.property_package.phase_component_list) else: # Otherwise assume all components in all phases phase_component_list = {} for p in self.phase_list: phase_component_list[p] = self.component_list # Link flowrates to property blocks @self.Constraint(self.time, self.ldomain, self.phase_list, self.component_list, doc="Link material flow terms to property packages") def link_material_flows(b, t, x, p, j): return b.material_flow[t, x, p, j] == ( b.properties[t, x].material_balance_term[p, j]) # Add Material Balance Variables as necessary # Mass diffusion term if self.config.has_mass_diffusion: self.material_diffusive_flux = Var( self.time, self.ldomain, self.phase_list, self.component_list, domain=Reals, doc="Axial flux of component due to diffusion per unit " "length [{}/{}.{}]".format(units['holdup'], units['time'], units['length'])) # Kinetic reaction generation term if self.config.has_rate_reactions: self.rate_reaction_generation = Var( self.time, self.ldomain, self.phase_list, self.component_list, domain=Reals, doc="Amount of component generated in " "unit by kinetic reactions per unit length " "[{}/{}.{}]".format(units['holdup'], units['time'], units['length'])) # Equilibrium reaction term if self.config.has_equilibrium_reactions: self.equilibrium_reaction_generation = Var( self.time, self.ldomain, self.phase_list, self.component_list, domain=Reals, doc="Amount of component generated in unit " "by equilibrium reactions per unit length " "[{}/{}.{}]".format(units['holdup'], units['time'], units['length'])) # Phase equilibrium generation if (self.config.has_phase_equilibrium and self.config.material_balance_type != 'component_total'): self.phase_equilibrium_generation = Var( self.time, self.ldomain, self.phase_equilibrium_idx, domain=Reals, doc="Amount of generation in unit by phase " "equilibria per unit length " "[{}/{}.{}]".format(units['holdup'], units['time'], units['length'])) # Mass transfer term if self.config.has_mass_transfer: self.mass_transfer_term = Var( self.time, self.ldomain, self.phase_list, self.component_list, domain=Reals, doc="Component material transfer into unit per unit " "length [{}/{}.{}]".format(units['holdup'], units['time'], units['length'])) # Create rules to substitute material balance terms # Accumulation term def accumulation_term(b, t, x, p, j): return (b.material_accumulation[t, x, p, j] if b.config.dynamic else 0) # Mass diffusion def diffusion_term(b, t, x, p, j): return (b.material_diffusive_flux[t, x, p, j] if b.config.has_mass_diffusion else 0) # Kineitc reaction term def kinetic_term(b, t, x, p, j): return (b.rate_reaction_generation[t, x, p, j] if b.config.has_rate_reactions else 0) # Equilibrium reaction term def equilibrium_term(b, t, x, p, j): return (b.equilibrium_reaction_generation[t, x, p, j] if b.config.has_equilibrium_reactions else 0) def phase_equilibrium_term(b, t, x, p, j): if (self.config.has_phase_equilibrium and self.config.material_balance_type != 'component_total'): sd = {} for r in b.phase_equilibrium_idx: if self.phase_equilibrium_list[r][0] == j: if self.phase_equilibrium_list[r][1][0] == p: sd[r] = 1 elif self.phase_equilibrium_list[r][1][1] == p: sd[r] = -1 else: sd[r] = 0 else: sd[r] = 0 return sum(b.phase_equilibrium_generation[t, x, r]*sd[r] for r in b.phase_equilibrium_idx) else: return 0 # Mass transfer term def transfer_term(b, t, x, p, j): return (b.mass_transfer_term[t, x, p, j] if b.config.has_mass_transfer else 0) # Component balances if self.config.material_balance_type == 'component_phase': @self.Constraint(self.time, self.ldomain, self.phase_list, self.component_list, doc="Material balances") def material_balance(b, t, x, p, j): if ((b.config.flow_direction is "forward" and x == b.ldomain.first()) or (b.config.flow_direction is "backward" and x == b.ldomain.last())): return Constraint.Skip else: if j in phase_component_list[p]: return b.length*accumulation_term(b, t, x, p, j) == ( b._flow_direction_term() * b.material_flow_dx[t, x, p, j] + b.length*kinetic_term(b, t, x, p, j) + b.length*equilibrium_term(b, t, x, p, j) + b.length*phase_equilibrium_term(b, t, x, p, j) + b.length*transfer_term(b, t, x, p, j) + b.area*diffusion_term(b, t, x, p, j)/b.length) else: return Constraint.Skip elif self.config.material_balance_type == 'component_total': @self.Constraint(self.time, self.ldomain, self.component_list, doc="Material balances") def material_balance(b, t, x, j): cplist = [] for p in b.phase_list: if j in phase_component_list[p]: cplist.append(p) if ((b.config.flow_direction is "forward" and x == b.ldomain.first()) or (b.config.flow_direction is "backward" and x == b.ldomain.last())): return Constraint.Skip else: return ( b.length*sum(accumulation_term(b, t, x, p, j) for p in cplist) == b._flow_direction_term() * sum(b.material_flow_dx[t, x, p, j] for p in cplist) + b.length*sum(kinetic_term(b, t, x, p, j) for p in cplist) + b.length*sum(equilibrium_term(b, t, x, p, j) for p in cplist) + b.length*sum(transfer_term(b, t, x, p, j) for p in cplist) + b.area*sum(diffusion_term(b, t, x, p, j) for p in cplist)/b.length) else: # This should never trigger. raise ValueError('Unrecognised value for argument ' 'material_balance_type') # TODO: Need to set material_holdup = 0 for non-present component-phase # pairs. Not ideal, but needed to close DoF. Is there a better way? # Material Holdup if self.config.include_holdup: @self.Constraint(self.time, self.ldomain, self.phase_list, self.component_list, doc="Material holdup calculation") def material_holdup_calculation(b, t, x, p, j): if j in phase_component_list[p]: return b.material_holdup[t, x, p, j] == ( b.area*b.phase_fraction[t, x, p] * b.properties[t, x].material_density_term[p, j]) else: return b.material_holdup[t, x, p, j] == 0 if self.config.has_rate_reactions: # Add extents of reaction and stoichiometric constraints self.rate_reaction_extent = Var( self.time, self.ldomain, self.rate_reaction_idx, domain=Reals, doc="Extent of kinetic reactions per unit length " "[{}/{}.{}]".format(units['holdup'], units['time'], units['length'])) @self.Constraint(self.time, self.ldomain, self.phase_list, self.component_list, doc="Kineitc reaction stoichiometry constraint") def rate_reaction_stoichiometry_constraint(b, t, x, p, j): if ((b.config.flow_direction is "forward" and x == b.ldomain.first()) or (b.config.flow_direction is "backward" and x == b.ldomain.last())): return Constraint.Skip else: if j in phase_component_list[p]: return b.rate_reaction_generation[t, x, p, j] == ( sum(b.rate_reaction_stoichiometry[r, p, j] * b.rate_reaction_extent[t, x, r] for r in b.rate_reaction_idx)) else: return Constraint.Skip if self.config.has_equilibrium_reactions: # Add extents of reaction and stoichiometric constraints self.equilibrium_reaction_extent = Var( self.time, self.ldomain, self.equilibrium_reaction_idx, domain=Reals, doc="Extent of equilibrium reactions per unit " "length [{}/{}.{}]".format(units['holdup'], units['time'], units['length'])) @self.Constraint(self.time, self.ldomain, self.phase_list, self.component_list, doc="Equilibrium reaction stoichiometry " "constraint") def equilibrium_reaction_stoichiometry_constraint(b, t, x, p, j): if ((b.config.flow_direction is "forward" and x == b.ldomain.first()) or (b.config.flow_direction is "backward" and x == b.ldomain.last())): return Constraint.Skip else: if j in phase_component_list[p]: return b.equilibrium_reaction_generation[t, x, p, j] \ == sum(b.equilibrium_reaction_stoichiometry[r, p, j] * b.equilibrium_reaction_extent[t, x, r] for r in b.equilibrium_reaction_idx) else: return Constraint.Skip if self.config.has_mass_diffusion: # Add diffusion transfer constraints @self.Constraint(self.time, self.ldomain, self.phase_list, self.component_list, doc="Link material concentration term to " "properties") def link_material_concentration(b, t, x, p, j): return b.material_concentration_term[t, x, p, j] == ( b.properties[t, x].material_concentration_term[p, j]) @self.Constraint(self.time, self.ldomain, self.phase_list, self.component_list, doc="Material diffusive flux calculation") def material_flux_calculation(b, t, x, p, j): if ((b.config.flow_direction is "forward" and x == b.ldomain.first()) or (b.config.flow_direction is "backward" and x == b.ldomain.last())): return Constraint.Skip else: return b.material_diffusive_flux[t, x, p, j] == ( -b.material_concentration_dx2[t, x, p, j] * b.properties[t, x].diffusivity[p, j]) def _make_element_total_balances(self): """ This method is called when material_balance_type = 'element', and constructs the element balance equations for the holdup block (indexed by element). Args: None Returns: None """ # Get units from property package units = {} for u in ['amount', 'holdup', 'time', 'length']: try: units[u] = self.config.property_package.get_package_units()[u] except KeyError: units[u] = '-' # Link flowrates to property blocks @self.Constraint(self.time, self.ldomain, self.phase_list, self.element_list, doc="Elemental flow calcuations") def elemental_flow_calculation(b, t, x, p, e): return b.elemental_flow[t, x, p, e] == ( sum(b.properties[t, x].material_balance_term[p, j] * b.properties[t, x].element_comp[j][e] for j in b.component_list)) # Add Material Balance Variables as necessary # Elemental diffusion term if self.config.has_mass_diffusion: self.material_diffusive_flux = Var( self.time, self.ldomain, self.phase_list, self.component_list, domain=Reals, doc="Axial flux of component due to diffusion per " "unit length [{}/{}.{}]".format(units['holdup'], units['time'], units['length'])) self.elemental_diffusive_flux = Var( self.time, self.ldomain, self.element_list, domain=Reals, doc="Axial flux of element due to diffusion per unit " "length [{}/{}.{}]".format(units['amount'], units['time'], units['length'])) # Mass transfer term if self.config.has_mass_transfer: self.elemental_mass_transfer = Var( self.time, self.ldomain, self.element_list, domain=Reals, doc="Element material transfer into unit per unit " "length [{}/{}.{}]".format(units['amount'], units['time'], units['length'])) # Create rules to substitute material balance terms # Accumulation term def accumulation_term(b, t, x, e): return b.element_accumulation[t, x, e] if b.config.dynamic else 0 # Elemental diffusive flux term def diffusion_term(b, t, x, e): return (b.elemental_diffusive_flux[t, x, e] if b.config.has_mass_diffusion else 0) # Elemental mass transfer term def transfer_term(b, t, x, e): return (b.elemental_mass_transfer[t, x, e] if b.config.has_mass_transfer else 0) # Element balances @self.Constraint(self.time, self.ldomain, self.element_list, doc="Element balances") def element_balance(b, t, x, e): if ((b.config.flow_direction is "forward" and x == b.ldomain.first()) or (b.config.flow_direction is "backward" and x == b.ldomain.last())): return Constraint.Skip else: return b.length*accumulation_term(b, t, x, e) == ( b._flow_direction_term() * sum(b.elemental_flow_dx[t, x, p, e] for p in b.phase_list) + b.length*transfer_term(b, t, x, e) + b.area*diffusion_term(b, t, x, e)/b.length) # Elemental Holdup if self.config.include_holdup: @self.Constraint(self.time, self.ldomain, self.element_list, doc="Elemental holdup calculation") def elemental_holdup_calculation(b, t, x, e): return b.element_holdup[t, x, e] == ( b.volume * sum(b.phase_fraction[t, x, p] * b.properties[t, x].material_density_term[p, j] * b.properties[t, x].element_comp[j][e] for p in b.phase_list for j in b.component_list)) if self.config.has_mass_diffusion: # Add diffusion transfer constraints @self.Constraint(self.time, self.ldomain, self.phase_list, self.component_list, doc="Link material concentration term to " "properties") def link_material_concentration(b, t, x, p, j): return b.material_concentration_term[t, x, p, j] == ( b.properties[t, x].material_concentration_term[p, j]) @self.Constraint(self.time, self.ldomain, self.phase_list, self.component_list, doc="Material diffusive flux calculation") def material_flux_calculation(b, t, x, p, j): if ((b.config.flow_direction is "forward" and x == b.ldomain.first()) or (b.config.flow_direction is "backward" and x == b.ldomain.last())): return Constraint.Skip else: return b.material_diffusive_flux[t, x, p, j] == ( -b.material_concentration_dx2[t, x, p, j] * b.properties[t, x].diffusivity[p, j]) # Calculate elemental flux terms @self.Constraint(self.time, self.ldomain, self.element_list, doc="Elemental diffusive flux calcuations") def elemental_diffusion_calculation(b, t, x, e): return b.elemental_diffusive_flux[t, x, e] == ( sum(b.material_diffusive_flux[t, x, p, j] * b.properties[t, x].element_comp[j][e] for j in b.component_list for p in b.phase_list)) def _make_total_energy_balance(self): """ This method is called when energy_balance_type = 'enthalpy_total', and constructs the energy balance equations for the holdup block. Args: None Returns: None """ # Get units from property package units = {} for u in ['energy', 'time', 'length']: try: units[u] = self.config.property_package.get_package_units()[u] except KeyError: units[u] = '-' # Link flow terms to property blocks @self.Constraint(self.time, self.ldomain, self.phase_list, doc="Link energy flow term to properties") def link_energy_flow(b, t, x, p,): return b.energy_flow[t, x, p] == ( b.properties[t, x].energy_balance_term[p]) # Create scaling factor self.scaling_factor_energy = Param( default=1e-6, mutable=True, doc='Energy balance scaling parameter') # Add energy balance variables as necessary # Energy diffusion term if self.config.has_energy_diffusion: self.energy_conduction_term = Var( self.time, self.ldomain, domain=Reals, doc="Axial conduction of heat per unit length " "[{}/{}.{}]".format(units['energy'], units['time'], units['length'])) # Heat transfer term if self.config.has_heat_transfer: self.heat = Var(self.time, self.ldomain, domain=Reals, doc="Heat transfered in unit per unit length " "[{}/{}.{}]".format(units['energy'], units['time'], units['length'])) # Work transfer term if self.config.has_work_transfer: self.work = Var(self.time, self.ldomain, domain=Reals, doc="Work transfered in unit per unit length " "[{}/{}.{}]".format(units['energy'], units['time'], units['length'])) # Create rules to substitute material balance terms # Accumulation term def accumulation_term(b, t, x, p): return b.energy_accumulation[t, x, p] if b.config.dynamic else 0 # Energy diffusion term def diffusion_term(b, t, x): return (b.energy_conduction_term[t, x] if b.config.has_energy_diffusion else 0) # Heat transfer term def heat_term(b, t, x): return b.heat[t, x] if b.config.has_heat_transfer else 0 # Work transfer term def work_term(b, t, x): return b.work[t, x] if b.config.has_work_transfer else 0 # Energy balance equation @self.Constraint(self.time, self.ldomain, doc="Energy balance") def energy_balance(b, t, x): if ((b.config.flow_direction is "forward" and x == b.ldomain.first()) or (b.config.flow_direction is "backward" and x == b.ldomain.last())): return Constraint.Skip else: return (b.length*sum(accumulation_term(b, t, x, p) for p in b.phase_list) * b.scaling_factor_energy) == ( b._flow_direction_term() * sum(b.energy_flow_dx[t, x, p] for p in b.phase_list) * b.scaling_factor_energy + b.length*heat_term(b, t, x) * b.scaling_factor_energy + b.length*work_term(b, t, x) * b.scaling_factor_energy - b.area*diffusion_term(b, t, x) * b.scaling_factor_energy/b.length) # Energy Holdup if self.config.include_holdup: @self.Constraint(self.time, self.ldomain, self.phase_list, doc="Energy holdup calcualtion") def energy_holdup_calculation(b, t, x, p): return b.energy_holdup[t, x, p] == ( b.area*b.phase_fraction[t, x, p] * b.properties[t, x].energy_density_term[p]) if self.config.has_energy_diffusion: # Link Temperature to property blocks @self.Constraint(self.time, self.ldomain, doc="Link temperature to properties") def link_temperature(b, t, x): return b.temperature[t, x] == b.properties[t, x].temperature @self.Constraint(self.time, self.ldomain, doc="Thermal diffusion calculation") def heat_conduction(b, t, x): if ((b.config.flow_direction is "forward" and x == b.ldomain.first()) or (b.config.flow_direction is "backward" and x == b.ldomain.last())): return Constraint.Skip else: return b.energy_conduction_term[t, x] == ( -b.temperature_dx2[t, x] * b.properties[t, x].conductivity) def _make_pressure_balance(self): """ This method is called when momentum_balance_type = 'pressure', and constructs the pressure balance equations for the holdup block. Args: None Returns: None """ # Get units from property package units = {} for u in ['pressure', 'length']: try: units[u] = self.config.property_package.get_package_units()[u] except KeyError: units[u] = '-' # Link pressure to property blocks @self.Constraint(self.time, self.ldomain, doc="Link pressure to properties") def link_pressure(b, t, x): return b.pressure[t, x] == b.properties[t, x].pressure # Add momentum balance variables as necessary if self.config.has_pressure_change: self.deltaP = Var( self.time, self.ldomain, domain=Reals, doc="Axial rate of change of pressure per unit length " "[{}/{}]".format(units['pressure'], units['length'])) # Create rules to substitute material balance terms # Pressure change term def deltaP_term(b, t, x): return b.deltaP[t, x] if b.config.has_pressure_change else 0 # Create scaling factor self.scaling_factor_pressure = Param( default=1e-4, mutable=True, doc='Momentum balance scaling parameter') # Momentum balance equation @self.Constraint(self.time, self.ldomain, doc="Pressure balance") def momentum_balance(b, t, x): if ((b.config.flow_direction is "forward" and x == b.ldomain.first()) or (b.config.flow_direction is "backward" and x == b.ldomain.last())): return Constraint.Skip else: return 0 == (b._flow_direction_term()*b.pressure_dx[t, x] * b.scaling_factor_pressure + b.length*deltaP_term(b, t, x) * b.scaling_factor_pressure) def _transform_domain(self): ''' This method applies the DAE transformation to the spatial domain of the Holdup block. Args: None Returns: None ''' # Apply DAE Transformation if self.config.discretization_method == 'BFD': self.discretizer = TransformationFactory('dae.finite_difference') self.discretizer.apply_to(self, nfe=self.config.finite_elements, wrt=self.ldomain, scheme='BACKWARD') elif self.config.discretization_method == 'FFD': self.discretizer = TransformationFactory('dae.finite_difference') self.discretizer.apply_to(self, nfe=self.config.finite_elements, wrt=self.ldomain, scheme='FORWARD') elif self.config.discretization_method == 'OCLL': self.discretizer = TransformationFactory('dae.collocation') self.discretizer.apply_to( self, wrt=self.ldomain, nfe=self.config.finite_elements, collocation_points=self.config.collocation_points, scheme='LAGRANGE-LEGENDRE') else: self.discretizer = TransformationFactory('dae.collocation') self.discretizer.apply_to( self, wrt=self.ldomain, nfe=self.config.finite_elements, collocation_points=self.config.collocation_points, scheme='LAGRANGE-RADAU')
[docs] def model_check(blk): """ This method exectues the model_check methods on the associated property blocks (if they exist). This method is generally called by a unit model as part of the unit's model_check method. Args: None Returns: None """ for t in blk.time: for x in blk.ldomain: # Try property block model check try: blk.properties[t, x].model_check() except AttributeError: logger.warning('{} Holdup property block has no model ' 'check. To correct this, add a model_check ' 'method to the associated PropertyBlock ' 'class.'.format(blk.name))
[docs] def initialize(blk, state_args=None, outlvl=0, hold_state=True, solver='ipopt', optarg=None): ''' Initialisation routine for holdup (default solver ipopt) Keyword Arguments: state_args : a dict of arguments to be passed to the property package(s) to provide an initial state for initialization (see documentation of the specific property package) (default = {}). outlvl : sets output level of initialisation routine * 0 = no output * 1 = return solver state for each step in routine * 2 = include solver output infomation (tee=True) optarg : solver options dictionary object (default=None) solver : str indicating whcih solver to use during initialization (default = 'ipopt') hold_state : flag indicating whether the initialization routine should unfix any state variables fixed during initialization (default=True). - True - states varaibles are not unfixed, and a dict of returned containing flags for which states were fixed during initialization. - False - state variables are unfixed after initialization by calling the relase_state method Returns: If hold_states is True, returns a dict containing flags for which states were fixed during initialization. ''' # Check if holdup has associated inlet_mixer, and initialize if present # Check for assoicated inlet_mixer block try: # Call inlet_mixer initialization routine mix_flags = blk.inlet_mixer.initialize(outlvl=outlvl-1, optarg=optarg, solver=solver, hold_state=hold_state, state_args=state_args) mixer_flag = True except AttributeError: mixer_flag = False # Get inlet state if not provided if state_args is None: state_args = {} if blk.config.flow_direction == "forward": inlet_key = blk.ldomain.first() else: inlet_key = blk.ldomain.last() state_dict = blk.properties[0, inlet_key].declare_port_members() for k in state_dict.keys(): if state_dict[k].is_indexed(): state_args[k] = {} for m in state_dict[k].keys(): state_args[k][m] = state_dict[k][m].value else: state_args[k] = state_dict[k].value # Initialize property blocks blk.properties.initialize(outlvl=outlvl-1, optarg=optarg, solver=solver, **state_args) if outlvl > 0: logger.info('{} Initialisation Complete'.format(blk.name)) if mixer_flag: return mix_flags
[docs] def release_state(blk, flags, outlvl=0): ''' Method to relase state variables fixed during initialisation. Keyword Arguments: flags : dict containing information of which state variables were fixed during initialization, and should now be unfixed. This dict is returned by initialize if hold_state = True. outlvl : sets output level of of logging Returns: None ''' try: blk.inlet_mixer.release_state(flags, outlvl-1) except AttributeError: pass
[docs]@declare_process_block_class("HoldupStatic") class HoldupStaticData(HoldupData): """ Static Holdup Class This class is designed to be used for unit operations zero volume or holdups with no through flow (such as dead zones). This type of holdup has only a single PropertyBlock index by time (Holdup0D has two). """
[docs] def build(self): """ Build method for HoldupStatic blocks. This method calls submethods to setup the necessary property blocks, distributed variables, material, energy and momentum balances based on the arguments provided by the user. Args: None Returns: None """ # Call build method from base class super(HoldupStaticData, self).build() # Construct Property Blocks self._make_prop_blocks() # Get indexing sets from property parameters block self._get_indexing_sets() # Build Common Variables self._make_vars() # Make phase fraction variables if necessary if self.config.include_holdup: self._make_phase_frac() # Build material balances if necessary if (self.config.has_rate_reactions or self.config.has_equilibrium_reactions or self.config.has_mass_transfer): if self.config.material_balance_type in \ ["component_phase", "component_total"]: self._make_component_balances() elif self.config.material_balance_type == "element_total": self._make_element_total_balances() # Build energy balances if necessary if (self.config.has_heat_transfer or self.config.has_work_transfer): if self.config.energy_balance_type == "enthalpy_total": self._make_total_energy_balance() # Build momentum balances if necessary if self.config.momentum_balance_type == "pressure": # No momentum balance in Holdup0D # However, flag needed for mixers/splitters pass
def _make_prop_blocks(self): """ This method constructs the property block for the holdup block. Args: None Returns: None """ self.properties = self.property_module.PropertyBlock( self.time, doc="Material properties in holdup", has_sum_fractions=False, calculate_equilibrium_reactions= self.config.has_equilibrium_reactions, calculate_phase_equilibrium=self.config.has_phase_equilibrium, parameters=self.config.property_package, **self.config.property_package_args) def _make_vars(self): """ This method constructs all the temporal derivative variables in the holdup block, and their associated dsitributed variables. Args: None Returns: None """ # Get units from property package units = {} for u in ['length', 'holdup', 'amount', 'time', 'energy']: try: units[u] = self.config.property_package.get_package_units()[u] except KeyError: units[u] = '-' if self.config.include_holdup: self.volume = Var(self.time, initialize=1.0, doc='Holdup Volume [{}^3]' .format(units['length'])) # Material holdup and accumulation if self.config.material_balance_type in ['component_phase', 'component_total']: if self.config.include_holdup: self.material_holdup = Var(self.time, self.phase_list, self.component_list, domain=Reals, doc="Material holdup in unit [{}]" .format(units['holdup'])) if self.config.dynamic is True: self.material_accumulation = DerivativeVar( self.material_holdup, wrt=self.time, doc="Material accumulation in unit[{}/{}]" .format(units['holdup'], units['time'])) elif self.config.material_balance_type == 'element_total': if self.config.include_holdup: self.element_holdup = Var( self.time, self.element_list, domain=Reals, doc="Elemental holdup in unit [{}]" .format(units['amount'])) if self.config.dynamic is True: self.element_accumulation = DerivativeVar( self.element_holdup, wrt=self.time, doc="Elemental accumulation in unit [{}/{}]" .format(units['amount'], units['time'])) # Energy holdup and accumulation if self.config.energy_balance_type == "enthalpy_total": if self.config.include_holdup: self.energy_holdup = Var( self.time, self.phase_list, domain=Reals, doc="Energy holdup in unit [{}]" .format(units['energy'])) if self.config.dynamic is True: self.energy_accumulation = DerivativeVar( self.energy_holdup, wrt=self.time, doc="Energy accumulation in unit [{}/{}]" .format(units['energy'], units['time'])) def _make_phase_frac(self): """ This method constructs the phase_fraction variables for the holdup block, and the associated constraint o nthe sum of phase_fraction == 1. For systems with only one phase, phase_fraction is created as a Pyomo Expression with a value of 1. Args: None Returns: None """ if len(self.phase_list) > 1: self.phase_fraction = Var( self.time, self.phase_list, initialize=1/len(self.phase_list), doc='Volume fraction of holdup by phase') @self.Constraint(self.time, doc='Sum of phase fractions == 1') def sum_of_phase_fractions(self, t): return 1 == sum(self.phase_fraction[t, p] for p in self.phase_list) else: @self.Expression(self.time, self.phase_list, doc='Volume fraction of holdup by phase') def phase_fraction(self, t, p): return 1 def _make_component_balances(self): """ This method is called when material_balance_type calls for component balances, and constructs the component balance equations for the holdup block. Args: None Returns: None """ # Get units from property package units = {} for u in ['holdup', 'time']: try: units[u] = self.config.property_package.get_package_units()[u] except KeyError: units[u] = '-' # Get phase component list(s) if hasattr(self.config.property_package, "phase_component_list"): phase_component_list = ( self.config.property_package.phase_component_list) else: # Otherwise assume all components in all phases phase_component_list = {} for p in self.phase_list: phase_component_list[p] = self.component_list # Create material balance terms as required # Kinetic reaction generation if self.config.has_rate_reactions: self.rate_reaction_generation = Var( self.time, self.phase_list, self.component_list, domain=Reals, doc="Amount of component generated in " "unit by kinetic reactions [{}/{}]" .format(units['holdup'], units['time'])) # Equilibrium reation generation if self.config.has_equilibrium_reactions: self.equilibrium_reaction_generation = Var( self.time, self.phase_list, self.component_list, domain=Reals, doc="Amount of component generated in unit " "by equilibrium reactions [{}/{}]" .format(units['holdup'], units['time'])) # Phase equilibrium generation if (self.config.has_phase_equilibrium and self.config.material_balance_type != 'component_total'): self.phase_equilibrium_generation = Var( self.time, self.phase_equilibrium_idx, domain=Reals, doc="Amount of generation in unit by phase " "equilibria [{}/{}]" .format(units['holdup'], units['time'])) # Material transfer term if self.config.has_mass_transfer: self.mass_transfer_term = Var( self.time, self.phase_list, self.component_list, domain=Reals, doc="Component material transfer into unit [{}/{}]" .format(units['holdup'], units['time'])) # Create rules to substitute material balance terms # Accumulation term def accumulation_term(b, t, p, j): return b.material_accumulation[t, p, j] if b.config.dynamic else 0 def kinetic_term(b, t, p, j): return (b.rate_reaction_generation[t, p, j] if b.config.has_rate_reactions else 0) def equilibrium_term(b, t, p, j): return (b.equilibrium_reaction_generation[t, p, j] if b.config.has_equilibrium_reactions else 0) def phase_equilibrium_term(b, t, p, j): if (self.config.has_phase_equilibrium and self.config.material_balance_type != 'component_total'): sd = {} for r in b.phase_equilibrium_idx: if self.phase_equilibrium_list[r][0] == j: if self.phase_equilibrium_list[r][1][0] == p: sd[r] = 1 elif self.phase_equilibrium_list[r][1][1] == p: sd[r] = -1 else: sd[r] = 0 else: sd[r] = 0 return sum(b.phase_equilibrium_generation[t, r]*sd[r] for r in b.phase_equilibrium_idx) else: return 0 def transfer_term(b, t, p, j): return (b.mass_transfer_term[t, p, j] if b.config.has_mass_transfer else 0) # Component balances if self.config.material_balance_type == 'component_phase': @self.Constraint(self.time, self.phase_list, self.component_list, doc="Material balances") def material_balance(b, t, p, j): if j in phase_component_list[p]: return accumulation_term(b, t, p, j) == ( kinetic_term(b, t, p, j) + equilibrium_term(b, t, p, j) + phase_equilibrium_term(b, t, p, j) + transfer_term(b, t, p, j)) else: return Constraint.Skip elif self.config.material_balance_type == 'component_total': @self.Constraint(self.time, self.component_list, doc="Material balances") def material_balance(b, t, j): cplist = [] for p in phase_component_list: if j in phase_component_list[p]: cplist.append(p) return ( sum(accumulation_term(b, t, p, j) for p in cplist) == sum(kinetic_term(b, t, p, j) for p in cplist) + sum(equilibrium_term(b, t, p, j) for p in cplist) + sum(transfer_term(b, t, p, j) for p in cplist)) else: # This should never trigger (should be caught earlier) raise ValueError('Unrecognised value for argument ' 'material_balance_type') # TODO: Need to set material_holdup = 0 for non-present component-phase # pairs. Not ideal, but needed to close DoF. Is there a better way? # Material Holdup if self.config.include_holdup: @self.Constraint(self.time, self.phase_list, self.component_list, doc="Material holdup calculations") def material_holdup_calculation(b, t, p, j): if j in phase_component_list[p]: return b.material_holdup[t, p, j] == ( b.volume[t]*self.phase_fraction[t, p] * b.properties[t].material_density_term[p, j]) else: return b.material_holdup[t, p, j] == 0 if self.config.has_rate_reactions: # Add extents of reaction and stoichiometric constraints self.rate_reaction_extent = Var( self.time, self.rate_reaction_idx, domain=Reals, doc="Extent of kinetic reactions [{}/{}]" .format(units['holdup'], units['time'])) @self.Constraint(self.time, self.phase_list, self.component_list, doc="Kinetic reaction stoichiometry constraint") def rate_reaction_stoichiometry_constraint(b, t, p, j): if j in phase_component_list[p]: return b.rate_reaction_generation[t, p, j] == ( sum(b.rate_reaction_stoichiometry[r, p, j] * b.rate_reaction_extent[t, r] for r in b.rate_reaction_idx)) else: return Constraint.Skip if self.config.has_equilibrium_reactions: # Add extents of reaction and stoichiometric constraints self.equilibrium_reaction_extent = Var( self.time, self.equilibrium_reaction_idx, domain=Reals, doc="Extent of equilibrium reactions [{}/{}]" .format(units['holdup'], units['time'])) @self.Constraint(self.time, self.phase_list, self.component_list, doc="Equilibrium reaction stoichiometry") def equilibrium_reaction_stoichiometry_constraint(b, t, p, j): if j in phase_component_list[p]: return b.equilibrium_reaction_generation[t, p, j] == ( sum(b.equilibrium_reaction_stoichiometry[r, p, j] * b.equilibrium_reaction_extent[t, r] for r in b.equilibrium_reaction_idx)) else: return Constraint.Skip def _make_element_total_balances(self): """ This method is called when material_balance_type = 'element', and constructs the element balance equations for the holdup block (indexed by element). Args: None Returns: None """ # Get units from property package units = {} for u in ['amount', 'time']: try: units[u] = self.config.property_package.get_package_units()[u] except KeyError: units[u] = '-' # Create material balance terms as needed if self.config.has_mass_transfer: self.elemental_mass_transfer = Var( self.time, self.element_list, domain=Reals, doc="Element material transfer into unit [{}/{}]" .format(units['amount'], units['time'])) # Create rules to substitute material balance terms # Accumulation term def accumulation_term(b, t, e): return b.element_accumulation[t, e] if b.config.dynamic else 0 # Mass transfer term def transfer_term(b, t, e): return (b.elemental_mass_transfer[t, e] if b.config.has_mass_transfer else 0) # Element balances @self.Constraint(self.time, self.element_list, doc="Elemental material balances") def element_balance(b, t, e): return accumulation_term(b, t, e) == transfer_term(b, t, e) # Elemental Holdup if self.config.include_holdup: @self.Constraint(self.time, self.element_list, doc="Elemental holdup calcuation") def elemental_holdup_calculation(b, t, e): return b.element_holdup[t, e] == ( b.volume[t] * sum(b.phase_fraction[t, p] * b.properties[t].material_density_term[p, j] * b.properties[t].element_comp[j][e] for p in b.phase_list for j in b.component_list)) def _make_total_energy_balance(self): """ This method is called when energy_balance_type = 'enthalpy_total', and constructs the energy balance equations for the holdup block. Args: None Returns: None """ # Get units from property package units = {} for u in ['energy', 'time']: try: units[u] = self.config.property_package.get_package_units()[u] except KeyError: units[u] = '-' # Create scaling factor self.scaling_factor_energy = Param( default=1e-6, mutable=True, doc='Energy balance scaling parameter') # Create energy balance terms as needed # Heat transfer term if self.config.has_heat_transfer: self.heat = Var(self.time, domain=Reals, doc="Heat transfered in unit [{}/{}]" .format(units['energy'], units['time'])) # Work transfer if self.config.has_work_transfer: self.work = Var(self.time, domain=Reals, doc="Work transfered in unit [{}/{}]" .format(units['energy'], units['time'])) # Create rules to substitute energy balance terms # Accumulation term def accumulation_term(b, t, p): return b.energy_accumulation[t, p] if b.config.dynamic else 0 def heat_term(b, t): return b.heat[t] if b.config.has_heat_transfer else 0 def work_term(b, t): return b.work[t] if b.config.has_work_transfer else 0 # Energy balance equation @self.Constraint(self.time, doc="Energy balances") def energy_balance(b, t): return (sum(accumulation_term(b, t, p) for p in b.phase_list) * b.scaling_factor_energy) == ( heat_term(b, t)*b.scaling_factor_energy + work_term(b, t)*b.scaling_factor_energy) # Energy Holdup if self.config.include_holdup: @self.Constraint(self.time, self.phase_list, doc="Energy holdup constraint") def energy_holdup_calculation(b, t, p): return b.energy_holdup[t, p] == ( b.volume[t]*self.phase_fraction[t, p] * b.properties[t].energy_density_term[p])
[docs] def model_check(blk): """ This method exectues the model_check methods on the associated property blocks (if they exist). This method is generally called by a unit model as part of the unit's model_check method. Args: None Returns: None """ # Try property block model check for t in blk.time: try: blk.properties[t].model_check() except AttributeError: logger.warning('{} Holdup property block has no model check. ' 'To correct this, add a model_check ' 'method to the associated PropertyBlock class.' .format(blk.name))
[docs] def initialize(blk, state_args=None, outlvl=0, optarg=None, solver='ipopt', hold_state=True): ''' Initialisation routine for holdup (default solver ipopt) Keyword Arguments: state_args : a dict of arguments to be passed to the property package(s) to provide an initial state for initialization (see documentation of the specific property package) (default = {}). outlvl : sets output level of initialisation routine * 0 = no output (default) * 1 = return solver state for each step in routine * 2 = include solver output infomation (tee=True) optarg : solver options dictionary object (default=None) solver : str indicating whcih solver to use during initialization (default = 'ipopt') hold_state : flag indicating whether the initialization routine should unfix any state variables fixed during initialization (default=True). - True - states varaibles are not unfixed, and a dict of returned containing flags for which states were fixed during initialization. - False - state variables are unfixed after initialization by calling the relase_state method Returns: If hold_states is True, returns a dict containing flags for which states were fixed during initialization. ''' # Check if holdup has associated inlet_mixer, and initialize if present # Check for assoicated inlet_mixer block try: # Call inlet_mixer initialization routine mix_flags = blk.inlet_mixer.initialize(outlvl=outlvl-1, optarg=optarg, solver=solver, hold_state=hold_state, state_args=state_args) mixer_flag = True hold_state_holdup = False except AttributeError: mixer_flag = False hold_state_holdup = hold_state # Check state arguments state_args = {} if state_args is None else state_args # Initialize property blocks flags = blk.properties.initialize(outlvl=outlvl-1, optarg=optarg, solver=solver, hold_state=hold_state_holdup, **state_args) # Check for assoicated outlet block try: # Call outlet_splitter initialization routine blk.outlet_splitter.initialize(outlvl=outlvl-1, optarg=optarg, solver=solver, hold_state=False, state_args=state_args) except AttributeError: pass if outlvl > 0: logger.info('{} Initialisation Complete'.format(blk.name)) if mixer_flag: return mix_flags else: return flags
[docs] def release_state(blk, flags, outlvl=0): ''' Method to relase state variables fixed during initialisation. Keyword Arguments: flags : dict containing information of which state variables were fixed during initialization, and should now be unfixed. This dict is returned by initialize if hold_state = True. outlvl : sets output level of of logging Returns: None ''' try: blk.inlet_mixer.release_state(flags, outlvl-1) except AttributeError: blk.properties.release_state(flags, outlvl=outlvl-1)
""" This section lays out a standard template for UnitModel Config Blocks. It is prepopulated with most of the standard construction arguments for Holdup Blocks (with the exception of property package arguments), along with preset default and domain settings. """ # Set up a standard ConfigBlock for models using Holdups CONFIG_Base = CONFIG() # Set default values CONFIG_Base.get('include_holdup')._default = False CONFIG_Base.get("material_balance_type")._default = 'component_phase' CONFIG_Base.get("energy_balance_type")._default = 'enthalpy_total' CONFIG_Base.get("momentum_balance_type")._default = 'pressure' CONFIG_Base.get("has_rate_reactions")._default = False CONFIG_Base.get("has_equilibrium_reactions")._default = False CONFIG_Base.get("has_phase_equilibrium")._default = False CONFIG_Base.get("has_mass_transfer")._default = False CONFIG_Base.get("has_heat_transfer")._default = False CONFIG_Base.get("has_work_transfer")._default = False CONFIG_Base.get("has_pressure_change")._default = False # Set domains CONFIG_Base.get('include_holdup')._domain = In([True, False]) CONFIG_Base.get("material_balance_type")._domain = In( ['none', 'component_phase', 'component_total', 'element_total']) CONFIG_Base.get("energy_balance_type")._domain = In( ['none', 'enthalpy_total']) CONFIG_Base.get("momentum_balance_type")._domain = In( ['none', 'pressure']) CONFIG_Base.get("has_rate_reactions")._domain = In([True, False]) CONFIG_Base.get("has_phase_equilibrium")._domain = In([True, False]) CONFIG_Base.get("has_equilibrium_reactions")._domain = In([True, False]) CONFIG_Base.get("has_mass_transfer")._domain = In([True, False]) CONFIG_Base.get("has_heat_transfer")._domain = In([True, False]) CONFIG_Base.get("has_work_transfer")._domain = In([True, False]) CONFIG_Base.get("has_pressure_change")._domain = In([True, False]) """ This section modifies CONFIG_Base with arguments specific to Holdup1D blocks. """ # Set up a standard ConfigBlock for models using Holdups CONFIG_Base_1D = CONFIG_Base() CONFIG_Base_1D.declare("velocity_type", ConfigValue( default='none', domain=In(['none', 'mixture', 'phase']), description="Velocity calculation flag", doc="""Flag indicating whether velocity should be include in the model **default** - 'use_parent_value'. **Valid values:** { **'none'** - do not calculate velocity, **'mixture'** - calculate a single velocity for all phases, **'phase'** - calculate separate velocities for each phase}"""))