##############################################################################
# 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".
##############################################################################
"""
Basic IDAES pressure changer models.
"""
from __future__ import division
# Import Python libraries
import logging
# Import Pyomo libraries
from pyomo.environ import SolverFactory, value, Var
from pyomo.opt import TerminationCondition
from pyomo.common.config import ConfigValue, In
# Import IDAES cores
from idaes.core import UnitBlockData, declare_process_block_class, \
Holdup0D, CONFIG_Base
from idaes.core.util.config import is_parameter_block, list_of_strings
from idaes.core.util.misc import add_object_ref
__author__ = "Andrew Lee"
# Set up logger
logger = logging.getLogger('idaes.unit_model')
# TODO : address equilibrium issues with setting isentropic property state
[docs]@declare_process_block_class("PressureChanger")
class PressureChangerData(UnitBlockData):
"""
Compressor/Expander Unit Class
"""
CONFIG = CONFIG_Base()
# Set default values of inherited attributes
CONFIG.get('include_holdup')._default = False
CONFIG.get("material_balance_type")._default = 'component_phase'
CONFIG.get("energy_balance_type")._default = 'total'
CONFIG.get("momentum_balance_type")._default = 'total'
CONFIG.get("has_rate_reactions")._default = False
CONFIG.get("has_equilibrium_reactions")._default = False
CONFIG.get("has_phase_equilibrium")._default = True
CONFIG.get("has_mass_transfer")._default = False
CONFIG.get("has_heat_transfer")._default = False
CONFIG.get("has_work_transfer")._default = True
CONFIG.get("has_pressure_change")._default = True
# Add unit model attributes
CONFIG.declare("property_package", ConfigValue(
default=None,
domain=is_parameter_block,
description="Property package to use for holdup",
doc="""Property parameter object used to define property calculations
(default = 'use_parent_value')
- 'use_parent_value' - get package from parent (default = None)
- a ParameterBlock object"""))
CONFIG.declare("property_package_args", ConfigValue(
default={},
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')
- 'use_parent_value' - get package from parent (default = None)
- a dict (see property package for documentation)"""))
CONFIG.declare("inlet_list", ConfigValue(
domain=list_of_strings,
description="List of inlet names",
doc="""A list containing names of inlets (default = None)
- None - default single inlet
- list - a list of names for inlets"""))
CONFIG.declare("num_inlets", ConfigValue(
domain=int,
description="Number of inlets to unit",
doc="""Argument indication number (int) of inlets to construct
(default = None). Not used if inlet_list arg is provided.
- None - use inlet_list arg instead
- int - Inlets will be named with sequential numbers from 1
to num_inlets"""))
CONFIG.declare("outlet_list", ConfigValue(
domain=list_of_strings,
description="List of outlet names",
doc="""A list containing names of outlets (default = None)
- None - default single outlet
- list - a list of names for outlets"""))
CONFIG.declare("num_outlets", ConfigValue(
domain=int,
description="Number of outlets to unit",
doc="""Argument indication number (int) of outlets to construct
(default = None). Not used if outlet_list arg is provided.
- None - use outlet_list arg instead
- int - Outlets will be named with sequential numbers from 1
to num_outlets"""))
CONFIG.declare("compressor", ConfigValue(
default=True,
domain=In([True, False]),
description="Compressor flag",
doc="""Indicates whether this unit should be considered a
compressor (True (default), pressure increase) or an expander
(False, pressure decrease)."""))
CONFIG.declare("thermodynamic_assumption", ConfigValue(
default='isothermal',
domain=In(['isothermal', 'isentropic', 'pump']),
description="Thermodynamic assumption to use",
doc="""Flag to set the thermodynamic assumption to use for the unit.
- 'isothermal' (default)
- 'isentropic'"""))
[docs] def build(self):
"""
Begin building model (pre-DAE transformation)
Args:
None
Returns:
None
"""
# Call UnitModel.build to setup dynamics
super(PressureChangerData, self).build()
# Build Holdup Block
self.holdup = Holdup0D()
# Set Unit Geometry and holdup Volume
self._set_geometry()
# Construct performance equations
self._make_performance()
# Construct equations for thermodynamic assumption
if self.config.thermodynamic_assumption == "isothermal":
self._make_isothermal()
elif self.config.thermodynamic_assumption == "isentropic":
self._make_isentropic()
elif self.config.thermodynamic_assumption == "pump":
self._make_pump()
[docs] def post_transform_build(self):
"""
Continue model construction after DAE transformation
Args:
None
Returns:
None
"""
# Construct Inlets
self.build_inlets(inlets=self.config.inlet_list,
num_inlets=self.config.num_inlets)
# Build Outlets
self.build_outlets(outlets=self.config.outlet_list,
num_outlets=self.config.num_outlets)
def _set_geometry(self):
"""
Define the geometry of the unit as necessary, and link to holdup volume
Args:
None
Returns:
None
"""
# For this case, just create a reference to holdup volume
if self.config.include_holdup is True:
add_object_ref(self, "volume", self.holdup.volume)
def _make_performance(self):
"""
Define constraints which describe the behaviour of the unit model.
Args:
None
Returns:
None
"""
# Set references to balance terms at unit level
add_object_ref(self, "work_mechanical", self.holdup.work)
add_object_ref(self, "deltaP", self.holdup.deltaP)
# Performance Variables
self.ratioP = Var(self.time, initialize=1.0,
doc="Pressure Ratio")
# Pressure Ratio
@self.Constraint(self.time, doc="Pressure ratio constraint")
def ratioP_calculation(b, t):
sf = b.holdup.scaling_factor_pressure
return sf*b.ratioP[t]*b.holdup.properties_in[t].pressure == \
sf*b.holdup.properties_out[t].pressure
def _make_pump(self):
"""
Add constraints for the incompresisble fluid assumption
Args:
None
Returns:
None
"""
self.work_fluid = Var(
self.time,
initialize=1.0,
doc="Work required to increase the pressure of the liquid")
self.efficiency_pump = Var(
self.time,
initialize=1.0,
doc="Pump efficiency")
@self.Constraint(self.time, doc="Pump fluid work constraint")
def fluid_work_calculation(b, t):
return b.work_fluid[t] == (
(b.holdup.properties_out[t].pressure -
b.holdup.properties_in[t].pressure) *
b.holdup.properties_out[t].flow_vol)
# Actual work
@self.Constraint(self.time, doc="Actual mechanical work calculation")
def actual_work(b, t):
sf = b.holdup.scaling_factor_energy
if b.config.compressor:
return sf*b.work_fluid[t] == sf*(
b.work_mechanical[t]*b.efficiency_pump[t])
else:
return sf*b.work_mechanical[t] == sf*(
b.work_fluid[t]*b.efficiency_pump[t])
def _make_isothermal(self):
"""
Add constraints for isothermal assumption.
Args:
None
Returns:
None
"""
# Isothermal constraint
@self.Constraint(self.time, doc="Equate inlet and outlet temperature")
def isothermal(b, t):
return b.holdup.properties_in[t].temperature == \
b.holdup.properties_out[t].temperature
def _make_isentropic(self):
"""
Add constraints for isentropic assumption.
Args:
None
Returns:
None
"""
# Get indexing sets from holdup block
add_object_ref(self, "phase_list", self.holdup.phase_list)
add_object_ref(self, "component_list", self.holdup.component_list)
# Add isentropic variables
self.efficiency_isentropic = Var(self.time,
initialize=0.8,
doc="Efficiency with respect to an "
"isentropic process [-]")
self.work_isentropic = Var(self.time,
initialize=0.0,
doc="Work input to unit if isentropic "
"process [-]")
# Build Isentropic Property block
self.properties_isentropic = \
self.holdup.property_module.PropertyBlock(
self.time,
parameters=self.holdup.config.property_package,
has_sum_fractions=True,
calculate_equilibrium_reactions=
self.config.has_equilibrium_reactions,
calculate_phase_equilibrium=self.config.has_phase_equilibrium,
**self.config.property_package_args)
# Connect isentropic property states
@self.Constraint(self.time, doc="Pressure for isentropic calculations")
def isentropic_pressure(b, t):
sf = b.holdup.scaling_factor_pressure
return sf*b.properties_isentropic[t].pressure == \
sf*b.ratioP[t]*b.holdup.properties_in[t].pressure
# TODO: This assumes isentropic composition is the same as outlet
@self.Constraint(self.time,
self.component_list,
doc="Material flows for isentropic properties")
def isentropic_material(b, t, j):
return b.properties_isentropic[t].flow_mol_comp[j] == \
b.holdup.properties_out[t].flow_mol_comp[j]
@self.Constraint(self.time, doc="Isentropic assumption")
def isentropic(b, t):
return b.properties_isentropic[t].entr_mol == \
b.holdup.properties_in[t].entr_mol
# Isentropic work
@self.Constraint(self.time, doc="Calculate work of isentropic process")
def isentropic_energy_balance(b, t):
sf = b.holdup.scaling_factor_energy
return sf*b.work_isentropic[t] == sf*(
sum(b.properties_isentropic[t].energy_balance_term[p]
for p in b.phase_list) -
sum(b.holdup.properties_in[t].energy_balance_term[p]
for p in b.phase_list))
# Actual work
@self.Constraint(self.time, doc="Actual mechanical work calculation")
def actual_work(b, t):
sf = b.holdup.scaling_factor_energy
if b.config.compressor:
return sf*b.work_isentropic[t] == sf*(
b.work_mechanical[t]*b.efficiency_isentropic[t])
else:
return sf*b.work_mechanical[t] == sf*(
b.work_isentropic[t]*b.efficiency_isentropic[t])
[docs] def model_check(blk):
"""
Check that pressure change matches with compressor argument (i.e. if
compressor = True, pressure should increase or work should be positive)
Args:
None
Returns:
None
"""
if blk.config.compressor:
# Compressor
# Check that pressure does not decrease
if any(blk.deltaP[t].fixed and
(value(blk.deltaP[t]) < 0.0) for t in blk.time):
logger.warning('{} Compressor set with negative deltaP.'
.format(blk.name))
if any(blk.ratioP[t].fixed and
(value(blk.ratioP[t]) < 1.0) for t in blk.time):
logger.warning('{} Compressor set with ratioP less than 1.'
.format(blk.name))
if any(blk.holdup.properties_out[t].pressure.fixed and
(value(blk.holdup.properties_in[t].pressure) >
value(blk.holdup.properties_out[t].pressure))
for t in blk.time):
logger.warning('{} Compressor set with pressure decrease.'
.format(blk.name))
# Check that work is not negative
if any(blk.work_mechanical[t].fixed and
(value(blk.work_mechanical[t]) < 0.0) for t in blk.time):
logger.warning('{} Compressor maybe set with negative work.'
.format(blk.name))
else:
# Expander
# Check that pressure does not increase
if any(blk.deltaP[t].fixed and
(value(blk.deltaP[t]) > 0.0) for t in blk.time):
logger.warning('{} Expander/turbine set with positive deltaP.'
.format(blk.name))
if any(blk.ratioP[t].fixed and
(value(blk.ratioP[t]) > 1.0) for t in blk.time):
logger.warning('{} Expander/turbine set with ratioP greater '
'than 1.'.format(blk.name))
if any(blk.holdup.properties_out[t].pressure.fixed and
(value(blk.holdup.properties_in[t].pressure) <
value(blk.holdup.properties_out[t].pressure))
for t in blk.time):
logger.warning('{} Expander/turbine maybe set with pressure ',
'increase.'.format(blk.name))
# Check that work is not positive
if any(blk.work_mechanical[t].fixed and
(value(blk.work_mechanical[t]) > 0.0) for t in blk.time):
logger.warning('{} Expander/turbine set with positive work.'
.format(blk.name))
# Run holdup block model checks
blk.holdup.model_check()
# Run model checks on isentropic property block
try:
for t in blk.time:
blk.properties_isentropic[t].model_check()
except AttributeError:
pass
[docs] def initialize(blk, state_args={}, routine=None, outlvl=0,
solver='ipopt', optarg={'tol': 1e-6}):
'''
General wrapper for pressure changer initialisation routines
Keyword Arguments:
routine : str stating which initialization routine to execute
* None - use routine matching thermodynamic_assumption
* 'isentropic' - use isentropic initialization routine
* 'isothermal' - use isothermal initialization routine
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 = return solver state for each step in subroutines
* 3 = include solver output infomation (tee=True)
optarg : solver options dictionary object (default={'tol': 1e-6})
solver : str indicating whcih solver to use during
initialization (default = 'ipopt')
Returns:
None
'''
if routine is None:
# Use routine for specific type of unit
routine = blk.config.thermodynamic_assumption
# Call initialisation routine
if routine is "isentropic":
blk.init_isentropic(state_args=state_args,
outlvl=outlvl,
solver=solver,
optarg=optarg)
else:
# Call the general initialization routine in UnitBlockData
super(PressureChangerData, blk).initialize(state_args=state_args,
outlvl=outlvl,
solver=solver,
optarg=optarg)
[docs] def init_isentropic(blk, state_args, outlvl, solver, optarg):
'''
Initialisation routine for unit (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 = return solver state for each step in subroutines
* 3 = include solver output infomation (tee=True)
optarg : solver options dictionary object (default={'tol': 1e-6})
solver : str indicating whcih solver to use during
initialization (default = 'ipopt')
Returns:
None
'''
# Set solver options
if outlvl > 3:
stee = True
else:
stee = False
opt = SolverFactory(solver)
opt.options = optarg
# ---------------------------------------------------------------------
# Initialize Isentropic block
blk.properties_isentropic.initialize(outlvl=outlvl-1,
optarg=optarg,
solver=solver,
**state_args)
if outlvl > 0:
logger.info('{} Initialisation Step 1 Complete.'.format(blk.name))
# ---------------------------------------------------------------------
# Initialize holdup block
flags = blk.holdup.initialize(outlvl=outlvl-1,
optarg=optarg,
solver=solver,
state_args=state_args)
if outlvl > 0:
logger.info('{} Initialisation Step 2 Complete.'.format(blk.name))
# ---------------------------------------------------------------------
# Solve for isothermal conditions
for t in blk.time:
blk.properties_isentropic[t].temperature.fix()
blk.isentropic.deactivate()
results = opt.solve(blk, tee=stee)
if outlvl > 0:
if results.solver.termination_condition == \
TerminationCondition.optimal:
logger.info('{} Initialisation Step 3 Complete.'
.format(blk.name))
else:
logger.warning('{} Initialisation Step 3 Failed.'
.format(blk.name))
# ---------------------------------------------------------------------
# Solve unit
for t in blk.time:
blk.properties_isentropic[t].temperature.unfix()
blk.isentropic.activate()
results = opt.solve(blk, tee=stee)
if outlvl > 0:
if results.solver.termination_condition == \
TerminationCondition.optimal:
logger.info('{} Initialisation Step 4 Complete.'
.format(blk.name))
else:
logger.warning('{} Initialisation Step 4 Failed.'
.format(blk.name))
# ---------------------------------------------------------------------
# Release Inlet state
blk.holdup.release_state(flags, outlvl-1)
if outlvl > 0:
logger.info('{} Initialisation Complete.'.format(blk.name))