Source code for idaes.models.pressure_changer

##############################################################################
# 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("has_phase_equilibrium")._default = True 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', 'adiabatic']), 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() elif self.config.thermodynamic_assumption == "adiabatic": self._make_adiabatic()
[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_adiabatic(self): """ Add constraints for adiabatic assumption. Args: None Returns: None """ # Isothermal constraint @self.Constraint(self.time, doc="Equate inlet and outlet enthalpy") def adiabatic(b, t): return b.holdup.properties_in[t].enth_mol == \ b.holdup.properties_out[t].enth_mol 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 if isinstance(blk.properties_isentropic[blk.time[1]].temperature, Var): 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)) for t in blk.time: blk.properties_isentropic[t].temperature.unfix() blk.isentropic.activate() elif outlvl > 0: logger.info('{} Initialisation Step 3 Skipped.'.format(blk.name)) # --------------------------------------------------------------------- # Solve unit 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))