##############################################################################
# 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', '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 = 'total'),
**'none'** - exclude energy balances,
**'total'** - single energy balance for all phases.}"""))
CONFIG.declare("momentum_balance_type", ConfigValue(
default='use_parent_value',
domain=In(['use_parent_value', 'none', 'total']),
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 = 'total'),
**'none'** - exclude momentum balances,
**'total'** - single momentum balance for all phases.}"""))
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': 'total',
'momentum_balance_type': 'total',
'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 == "total":
self._make_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 == "total":
self._make_momentum_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 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 == "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_energy_balance(self):
"""
This method is called when energy_balance_type = '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_momentum_balance(self):
"""
This method is called when momentum_balance_type = 'total', and
constructs the momentum 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)"""))
[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()
# 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 == "total":
self._make_energy_balance()
# Build momentum balances if necessary
if self.config.momentum_balance_type == "total":
self._make_momentum_balance()
# Transform axial domain - not needed if inherited domain provided.
if self.config.inherited_length_domain is None:
self._transform_domain()
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']:
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 == "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")
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 == "total":
self.pressure = Var(self.time,
self.ldomain,
doc="Pressure [{}]".format(units['pressure']))
self.pressure_dx = DerivativeVar(self.pressure, wrt=self.ldomain)
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.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']:
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 [{}/{}]"
.format(units['holdup'], units['time']))
# 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 [{}/{}]"
.format(units['holdup'], units['time']))
# 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 [{}/{}]"
.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.ldomain,
self.phase_equilibrium_idx,
domain=Reals,
doc="Amount of generation in unit by phase "
"equilibria [{}/{}]"
.format(units['holdup'], units['time']))
# 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 [{}/{}]"
.format(units['holdup'], units['time']))
# 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 [{}/{}]"
.format(units['holdup'], units['time']))
@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 [{}/{}]"
.format(units['holdup'], units['time']))
@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']:
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 [{}/{}]"
.format(units['holdup'], units['time']))
self.elemental_diffusive_flux = Var(
self.time,
self.ldomain,
self.element_list,
domain=Reals,
doc="Axial flux of element due to diffusion [{}/{}]"
.format(units['amount'], units['time']))
# 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 [{}/{}]"
.format(units['amount'], units['time']))
# 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_energy_balance(self):
"""
This method is called when energy_balance_type = '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] = '-'
# 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 [{}/{}]"
.format(units['energy'], units['time']))
# Heat transfer term
if self.config.has_heat_transfer:
self.heat = Var(self.time,
self.ldomain,
domain=Reals,
doc="Heat transfered in unit [{}/{}]"
.format(units['energy'], units['time']))
# Work transfer term
if self.config.has_work_transfer:
self.work = Var(self.time,
self.ldomain,
domain=Reals,
doc="Work transfered in unit [{}/{}]"
.format(units['energy'], units['time']))
# 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_momentum_balance(self):
"""
This method is called when momentum_balance_type = 'total', and
constructs the momentum balance equations for the holdup block.
Args:
None
Returns:
None
"""
# Get units from property package
try:
punits = (self.config.property_package
.get_package_units()['pressure'])
except KeyError:
punits = '-'
# 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(punits))
# 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 == "total":
self._make_energy_balance()
# Build momentum balances if necessary
if self.config.momentum_balance_type == "total":
# 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 == "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_energy_balance(self):
"""
This method is called when energy_balance_type = '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 = 'total'
CONFIG_Base.get("momentum_balance_type")._default = 'total'
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', 'total'])
CONFIG_Base.get("momentum_balance_type")._domain = In(
['none', 'total'])
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])