Source code for idaes.models_extra.gas_solid_contactors.unit_models.fixed_bed_1D

#################################################################################
# The Institute for the Design of Advanced Energy Systems Integrated Platform
# Framework (IDAES IP) was produced under the DOE Institute for the
# Design of Advanced Energy Systems (IDAES).
#
# Copyright (c) 2018-2024 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.md and LICENSE.md
# for full copyright and license information.
#################################################################################
"""
IDAES 1D Fixed Bed model.

"""

# Import Pyomo libraries
from pyomo.environ import (
    Var,
    Param,
    Reals,
    value,
    TransformationFactory,
    Constraint,
    check_optimal_termination,
    units as pyunits,
)
from pyomo.common.config import ConfigBlock, ConfigValue, In, Bool
from pyomo.util.calc_var_value import calculate_variable_from_constraint
from pyomo.dae import ContinuousSet, DerivativeVar
from pyomo.util.subsystems import TemporarySubsystemManager
from pyomo.common.collections import ComponentSet
from pyomo.contrib.incidence_analysis import (
    IncidenceGraphInterface,
    solve_strongly_connected_components,
)

# Import IDAES cores
from idaes.core import (
    ControlVolume1DBlock,
    UnitModelBlockData,
    declare_process_block_class,
    MaterialBalanceType,
    EnergyBalanceType,
    MomentumBalanceType,
    FlowDirection,
    DistributedVars,
)
from idaes.core.util.config import (
    is_physical_parameter_block,
    is_reaction_parameter_block,
)
from idaes.core.util.exceptions import (
    ConfigurationError,
    BurntToast,
    InitializationError,
)
from idaes.core.util.initialization import fix_state_vars, revert_state_vars
from idaes.core.util.model_statistics import degrees_of_freedom
from idaes.core.util.dyn_utils import get_index_set_except
from idaes.core.util.tables import create_stream_table_dataframe
from idaes.core.util.constants import Constants as constants
from idaes.core.util.math import smooth_abs
import idaes.logger as idaeslog
from idaes.core.util import scaling as iscale
from idaes.core.solvers import get_solver

__author__ = "Chinedu Okoli"

# Set up logger
_log = idaeslog.getLogger(__name__)

# Assumptions:
# Only gas and solid phases, which are explicitly named as such
# Perfect mixing in gas phase
# Axially varying gas and solid phase
# Reaction/adsorption on the solid is the rate limiting step


[docs] @declare_process_block_class("FixedBed1D") class FixedBed1DData(UnitModelBlockData): """ 1D Fixed Bed Reactor Model Class """ # Create template for unit level config arguments CONFIG = UnitModelBlockData.CONFIG() # Unit level config arguments CONFIG.declare( "finite_elements", ConfigValue( default=10, domain=int, description="Number of finite elements length domain", doc="""Number of finite elements to use when discretizing length domain (default=20)""", ), ) CONFIG.declare( "length_domain_set", ConfigValue( default=[0.0, 1.0], domain=list, description="Number of finite elements length domain", doc="""length_domain_set - (optional) list of point to use to initialize a new ContinuousSet if length_domain is not provided (default = [0.0, 1.0])""", ), ) CONFIG.declare( "transformation_method", ConfigValue( default="dae.finite_difference", description="Method to use for DAE transformation", doc="""Method to use to transform domain. Must be a method recognized by the Pyomo TransformationFactory, **default** - "dae.finite_difference". **Valid values:** { **"dae.finite_difference"** - Use a finite difference transformation method, **"dae.collocation"** - use a collocation transformation method}""", ), ) CONFIG.declare( "transformation_scheme", ConfigValue( default=None, domain=In([None, "BACKWARD", "FORWARD", "LAGRANGE-RADAU"]), description="Scheme to use for DAE transformation", doc="""Scheme to use when transforming domain. See Pyomo documentation for supported schemes, **default** - None. **Valid values:** { **None** - defaults to "BACKWARD" for finite difference transformation method, and to "LAGRANGE-RADAU" for collocation transformation method, **"BACKWARD"** - Use a finite difference transformation method, **"FORWARD""** - use a finite difference transformation method, **"LAGRANGE-RADAU""** - use a collocation transformation method}""", ), ) 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( "flow_type", ConfigValue( default="forward_flow", domain=In(["forward_flow", "reverse_flow"]), description="Flow configuration of Fixed Bed", doc="""Flow configuration of Fixed Bed **default** - "forward_flow". **Valid values:** { **"forward_flow"** - gas flows from 0 to 1, **"reverse_flow"** - gas flows from 1 to 0.}""", ), ) CONFIG.declare( "material_balance_type", ConfigValue( default=MaterialBalanceType.componentTotal, domain=In(MaterialBalanceType), description="Material balance construction flag", doc="""Indicates what type of mass balance should be constructed, **default** - MaterialBalanceType.componentTotal. **Valid values:** { **MaterialBalanceType.none** - exclude material balances, **MaterialBalanceType.componentPhase** - use phase component balances, **MaterialBalanceType.componentTotal** - use total component balances, **MaterialBalanceType.elementTotal** - use total element balances, **MaterialBalanceType.total** - use total material balance.}""", ), ) CONFIG.declare( "energy_balance_type", ConfigValue( default=EnergyBalanceType.enthalpyTotal, domain=In(EnergyBalanceType), description="Energy balance construction flag", doc="""Indicates what type of energy balance should be constructed, **default** - EnergyBalanceType.enthalpyTotal. **Valid values:** { **EnergyBalanceType.none** - exclude energy balances, **EnergyBalanceType.enthalpyTotal** - single enthalpy balance for material, **EnergyBalanceType.enthalpyPhase** - enthalpy balances for each phase, **EnergyBalanceType.energyTotal** - single energy balance for material, **EnergyBalanceType.energyPhase** - energy balances for each phase.}""", ), ) CONFIG.declare( "momentum_balance_type", ConfigValue( default=MomentumBalanceType.pressureTotal, domain=In(MomentumBalanceType), description="Momentum balance construction flag", doc="""Indicates what type of momentum balance should be constructed, **default** - MomentumBalanceType.pressureTotal. **Valid values:** { **MomentumBalanceType.none** - exclude momentum balances, **MomentumBalanceType.pressureTotal** - single pressure balance for material, **MomentumBalanceType.pressurePhase** - pressure balances for each phase, **MomentumBalanceType.momentumTotal** - single momentum balance for material, **MomentumBalanceType.momentumPhase** - momentum balances for each phase.}""", ), ) CONFIG.declare( "has_pressure_change", ConfigValue( default=True, domain=Bool, description="Pressure change term construction flag", doc="""Indicates whether terms for pressure change should be constructed, **default** - False. **Valid values:** { **True** - include pressure change terms, **False** - exclude pressure change terms.}""", ), ) CONFIG.declare( "pressure_drop_type", ConfigValue( default="ergun_correlation", domain=In(["ergun_correlation", "simple_correlation"]), description="Construction flag for type of pressure drop", doc="""Indicates what type of pressure drop correlation should be used, **default** - "ergun_correlation". **Valid values:** { **"ergun_correlation"** - Use the Ergun equation, **"simple_correlation"** - Use a simplified pressure drop correlation.}""", ), ) # Create template for phase specific config arguments _PhaseTemplate = UnitModelBlockData.CONFIG() _PhaseTemplate.declare( "property_package", ConfigValue( default=None, domain=is_physical_parameter_block, description="Property package to use for control volume", 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""", ), ) _PhaseTemplate.declare( "property_package_args", ConfigValue( default={}, domain=dict, description="Arguments for constructing gas property package", doc="""A dict of arguments to be passed to the PropertyBlockData 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)""", ), ) _PhaseTemplate.declare( "reaction_package", ConfigValue( default=None, domain=is_reaction_parameter_block, description="Reaction package to use for control volume", doc="""Reaction parameter object used to define reaction calculations, **default** - None. **Valid values:** { **None** - no reaction package, **ReactionParameterBlock** - a ReactionParameterBlock object.}""", ), ) _PhaseTemplate.declare( "reaction_package_args", ConfigBlock( implicit=True, implicit_domain=ConfigBlock, description="Arguments to use for constructing reaction packages", doc="""A ConfigBlock with arguments to be passed to a reaction block(s) and used when constructing these, **default** - None. **Valid values:** { see reaction package for documentation.}""", ), ) _PhaseTemplate.declare( "has_equilibrium_reactions", ConfigValue( default=False, domain=Bool, description="Equilibrium reaction construction flag", doc="""Indicates whether terms for equilibrium controlled reactions should be constructed, **default** - True. **Valid values:** { **True** - include equilibrium reaction terms, **False** - exclude equilibrium reaction terms.}""", ), ) # Create individual config blocks for gas and solid sides CONFIG.declare("gas_phase_config", _PhaseTemplate(doc="gas phase config arguments")) CONFIG.declare( "solid_phase_config", _PhaseTemplate(doc="solid phase config arguments") ) # =========================================================================
[docs] def build(self): """ Begin building model (pre-DAE transformation). Args: None Returns: None """ # Call UnitModel.build to build default attributes super().build() # local aliases used to shorten object names solid_phase = self.config.solid_phase_config gas_phase = self.config.gas_phase_config # Get units meta data from property packages units_meta_solid = solid_phase.property_package.get_metadata().get_derived_units # Set flow direction for the gas control volume # Gas flows from 0 to 1 if self.config.flow_type == "forward_flow": set_direction_gas = FlowDirection.forward # Gas flows from 1 to 0 if self.config.flow_type == "reverse_flow": set_direction_gas = FlowDirection.backward # Consistency check for flow direction, transformation method and # transformation scheme if ( self.config.flow_type == "forward_flow" and self.config.transformation_method == "dae.finite_difference" and self.config.transformation_scheme is None ): self.config.transformation_scheme = "BACKWARD" elif ( self.config.flow_type == "reverse_flow" and self.config.transformation_method == "dae.finite_difference" and self.config.transformation_scheme is None ): self.config.transformation_scheme = "FORWARD" elif ( self.config.flow_type == "forward_flow" and self.config.transformation_method == "dae.collocation" and self.config.transformation_scheme is None ): self.config.transformation_scheme = "LAGRANGE-RADAU" elif ( self.config.flow_type == "reverse_flow" and self.config.transformation_method == "dae.collocation" ): raise ConfigurationError( "{} invalid value for " "transformation_method argument." "Must be " "dae.finite_difference " "if " "flow_type is" " " "reverse_flow" ".".format(self.name) ) elif ( self.config.flow_type == "forward_flow" and self.config.transformation_scheme == "FORWARD" ): raise ConfigurationError( "{} invalid value for " "transformation_scheme argument. " "Must be " "BACKWARD " "if flow_type is" " " "forward_flow" ".".format(self.name) ) elif ( self.config.flow_type == "reverse_flow" and self.config.transformation_scheme == "BACKWARD" ): raise ConfigurationError( "{} invalid value for " "transformation_scheme argument." "Must be " "FORWARD " "if " "flow_type is" " " "reverse_flow" ".".format(self.name) ) elif ( self.config.transformation_method == "dae.finite_difference" and self.config.transformation_scheme != "BACKWARD" and self.config.transformation_scheme != "FORWARD" ): raise ConfigurationError( "{} invalid value for " "transformation_scheme argument. " "Must be " "BACKWARD" " or " "FORWARD" " " "if transformation_method is" " " "dae.finite_difference" ".".format(self.name) ) elif ( self.config.transformation_method == "dae.collocation" and self.config.transformation_scheme != "LAGRANGE-RADAU" ): raise ConfigurationError( "{} invalid value for " "transformation_scheme argument." "Must be " "LAGRANGE-RADAU" " if " "transformation_method is" " " "dae.collocation" ".".format(self.name) ) # Set arguments for gas sides if homoogeneous reaction block if gas_phase.reaction_package is not None: has_rate_reaction_gas_phase = True else: has_rate_reaction_gas_phase = False # Set arguments for gas and solid sides if heterogeneous reaction block if solid_phase.reaction_package is not None: has_mass_transfer_gas_phase = True else: has_mass_transfer_gas_phase = False # Set heat transfer terms if self.config.energy_balance_type != EnergyBalanceType.none: has_heat_transfer = True else: has_heat_transfer = False # Set heat of reaction terms if ( self.config.energy_balance_type != EnergyBalanceType.none and gas_phase.reaction_package is not None ): has_heat_of_reaction_gas_phase = True else: has_heat_of_reaction_gas_phase = False # Create a unit model length domain self.length_domain = ContinuousSet( bounds=(0.0, 1.0), initialize=self.config.length_domain_set, doc="Normalized length domain", ) self.bed_height = Var( domain=Reals, initialize=1, doc="Bed length", units=units_meta_solid("length"), ) # ========================================================================= # Build Control volume 1D for gas phase and populate gas control volume self.gas_phase = ControlVolume1DBlock( transformation_method=self.config.transformation_method, transformation_scheme=self.config.transformation_scheme, finite_elements=self.config.finite_elements, collocation_points=self.config.collocation_points, dynamic=True, # Fixed beds must be dynamic has_holdup=True, # holdup must be True for fixed beds area_definition=DistributedVars.variant, property_package=gas_phase.property_package, property_package_args=gas_phase.property_package_args, reaction_package=gas_phase.reaction_package, reaction_package_args=gas_phase.reaction_package_args, ) self.gas_phase.add_geometry( length_domain=self.length_domain, length_domain_set=self.config.length_domain_set, length_var=self.bed_height, flow_direction=set_direction_gas, ) self.gas_phase.add_state_blocks( information_flow=set_direction_gas, has_phase_equilibrium=False ) if gas_phase.reaction_package is not None: self.gas_phase.add_reaction_blocks( has_equilibrium=gas_phase.has_equilibrium_reactions ) self.gas_phase.add_material_balances( balance_type=self.config.material_balance_type, has_phase_equilibrium=False, has_mass_transfer=has_mass_transfer_gas_phase, has_rate_reactions=has_rate_reaction_gas_phase, ) self.gas_phase.add_energy_balances( balance_type=self.config.energy_balance_type, has_heat_transfer=has_heat_transfer, has_heat_of_reaction=has_heat_of_reaction_gas_phase, ) self.gas_phase.add_momentum_balances( balance_type=self.config.momentum_balance_type, has_pressure_change=self.config.has_pressure_change, ) # ========================================================================= # Build Control volume 1D for solid phase and populate solid control volume # Build Solid Phase StateBlock # As there is no solid flow, there is not need for a control volume as a # set of indexed state blocks will be sufficient. # Defined state is set to True so that the "sum(mass_frac)=1" eqn in # the solid state block is deactivated. This is done here as there is # currently no way to deactivate the constraint at the initial time # for batch systems (i.e. no inlet or outlet ports). # The "sum(mass_frac)=1 for all t neq 0" eqn is written in the # unit model instead self.solid_properties = solid_phase.property_package.state_block_class( self.flowsheet().time, self.length_domain, parameters=solid_phase.property_package, defined_state=True, ) if solid_phase.reaction_package is not None: # TODO - a generalization of the heterogeneous reaction block # The heterogeneous reaction block does not use the # add_reaction_blocks in control volumes as control volumes are # currently setup to handle only homogeneous reaction properties. # Thus appending the heterogeneous reaction block to the # solid state block is currently hard coded here. tmp_dict = dict(**solid_phase.reaction_package_args) tmp_dict["gas_state_block"] = self.gas_phase.properties tmp_dict["solid_state_block"] = self.solid_properties tmp_dict["parameters"] = solid_phase.reaction_package self.solid_reactions = solid_phase.reaction_package.reaction_block_class( self.flowsheet().time, self.length_domain, doc="Reaction properties in control volume", **tmp_dict, ) # ========================================================================= # Add Ports for gas side self.add_inlet_port(name="gas_inlet", block=self.gas_phase) self.add_outlet_port(name="gas_outlet", block=self.gas_phase) # ========================================================================= # Add performance equation method self._apply_transformation() self._make_performance()
# ========================================================================= def _apply_transformation(self): """ Method to apply DAE transformation to the Control Volume length domain. Transformation applied will be based on the Control Volume configuration arguments. """ if self.config.finite_elements is None: raise ConfigurationError( "{} was not provided a value for the finite_elements" " configuration argument. Please provide a valid value.".format( self.name ) ) if self.config.transformation_method == "dae.finite_difference": self.discretizer = TransformationFactory(self.config.transformation_method) self.discretizer.apply_to( self, wrt=self.length_domain, nfe=self.config.finite_elements, scheme=self.config.transformation_scheme, ) elif self.config.transformation_method == "dae.collocation": self.discretizer = TransformationFactory(self.config.transformation_method) self.discretizer.apply_to( self, wrt=self.length_domain, nfe=self.config.finite_elements, ncp=self.config.collocation_points, scheme=self.config.transformation_scheme, ) def _make_performance(self): """ Constraints for unit model. Args: None Returns: None """ # local aliases used to shorten object names gas_phase = self.config.gas_phase_config solid_phase = self.config.solid_phase_config # Get units meta data from property packages units_meta_gas = gas_phase.property_package.get_metadata().get_derived_units units_meta_solid = solid_phase.property_package.get_metadata().get_derived_units # Declare Mutable Parameters self.eps = Param( mutable=True, default=1e-8, doc="Smoothing Factor for Smooth IF Statements", units=pyunits.dimensionless, ) # Unit Model variables self.bed_diameter = Var( initialize=1, doc="Reactor diameter", units=units_meta_solid("length") ) self.bed_area = Var( domain=Reals, initialize=1, doc="Reactor cross-sectional area", units=units_meta_solid("area"), ) # Phase specific variables self.velocity_superficial_gas = Var( self.flowsheet().time, self.length_domain, domain=Reals, initialize=0.05, doc="Gas superficial velocity", units=units_meta_gas("velocity"), ) # Dimensionless numbers, mass and heat transfer coefficients self.Re_particle = Var( self.flowsheet().time, self.length_domain, domain=Reals, initialize=1.0, doc="Particle Reynolds number", units=pyunits.dimensionless, ) self.Pr_particle = Var( self.flowsheet().time, self.length_domain, domain=Reals, initialize=1.0, doc="Prandtl number of gas in bed", units=pyunits.dimensionless, ) self.Nu_particle = Var( self.flowsheet().time, self.length_domain, domain=Reals, initialize=1.0, doc="Particle Nusselt number", units=pyunits.dimensionless, ) self.gas_solid_htc = Var( self.flowsheet().time, self.length_domain, domain=Reals, initialize=1.0, doc="Gas-solid heat transfer coefficient", units=units_meta_gas("heat_transfer_coefficient"), ) # Solid phase variables self.solid_phase_area = Var( self.flowsheet().time, self.length_domain, domain=Reals, initialize=1, doc="Solid phase area", units=units_meta_solid("area"), ) self.solid_material_holdup = Var( self.flowsheet().time, self.length_domain, solid_phase.property_package.component_list, domain=Reals, initialize=1.0, doc="Solid material holdup per unit length", units=units_meta_solid("mass") / units_meta_solid("length"), ) self.solid_material_accumulation = DerivativeVar( self.solid_material_holdup, wrt=self.flowsheet().config.time, doc="Solid material accumulation", units=units_meta_solid("mass") / units_meta_solid("length") / units_meta_solid("time"), ) if self.config.energy_balance_type != EnergyBalanceType.none: # Solid phase - gas to solid heat transfer self.solid_phase_heat = Var( self.flowsheet().time, self.length_domain, domain=Reals, initialize=0.0, doc="Heat transferred per unit length", units=units_meta_solid("power") / units_meta_solid("length"), ) self.solid_energy_holdup = Var( self.flowsheet().time, self.length_domain, initialize=1, doc="Solid phase energy holdup", units=units_meta_solid("energy") / units_meta_solid("length"), ) self.solid_energy_accumulation = DerivativeVar( self.solid_energy_holdup, initialize=0, wrt=self.flowsheet().config.time, doc="Solid energy accumulation", units=units_meta_solid("power") / units_meta_solid("length"), ) # ========================================================================= # Add performance equations # --------------------------------------------------------------------- # TODO - Deactivate/delete sum_component_eqns of initial conditions for # all space in gas_phase as currently no way to tell CV1D to not build # sum_eqns for initial conditions if self.gas_phase.config.dynamic: t0 = self.flowsheet().time.first() for x in self.length_domain: blk = self.gas_phase.properties[t0, x] if hasattr(blk, "sum_component_eqn"): blk.sum_component_eqn.deactivate() # --------------------------------------------------------------------- # Geometry constraints # Bed area @self.Constraint(doc="Bed area") def bed_area_eqn(b): return b.bed_area == (constants.pi * (0.5 * b.bed_diameter) ** 2) # Area of gas side, and solid side @self.Constraint(self.flowsheet().time, self.length_domain, doc="Gas side area") def gas_phase_area_constraint(b, t, x): return ( b.gas_phase.area[t, x] == pyunits.convert(b.bed_area, to_units=units_meta_gas("area")) * b.solid_properties[t, x].params.voidage ) @self.Constraint( self.flowsheet().time, self.length_domain, doc="Solid side area" ) def solid_phase_area_constraint(b, t, x): return b.solid_phase_area[t, x] == b.bed_area * ( 1 - b.solid_properties[t, x].params.voidage ) # --------------------------------------------------------------------- # Hydrodynamic constraints # Gas superficial velocity @self.Constraint( self.flowsheet().time, self.length_domain, doc="Gas superficial velocity" ) def gas_super_vel(b, t, x): return ( b.velocity_superficial_gas[t, x] * pyunits.convert(b.bed_area, to_units=units_meta_gas("area")) * b.gas_phase.properties[t, x].dens_mol == b.gas_phase.properties[t, x].flow_mol ) # Gas side pressure drop calculation if ( self.config.has_pressure_change and self.config.pressure_drop_type == "simple_correlation" ): # Simplified pressure drop @self.Constraint( self.flowsheet().time, self.length_domain, doc="Gas side pressure drop calculation -" "simplified pressure drop", ) def gas_phase_config_pressure_drop(b, t, x): # 0.2/s is a unitted constant in the correlation return pyunits.convert( b.gas_phase.deltaP[t, x], to_units=units_meta_solid("pressure") / units_meta_solid("length"), ) == -(0.2 / pyunits.s) * ( b.velocity_superficial_gas[t, x] * ( pyunits.convert( b.solid_properties[t, x].dens_mass_particle, to_units=units_meta_solid("density_mass"), ) - b.gas_phase.properties[t, x].dens_mass ) ) elif ( self.config.has_pressure_change and self.config.pressure_drop_type == "ergun_correlation" ): # Ergun equation @self.Constraint( self.flowsheet().time, self.length_domain, doc="Gas side pressure drop calculation -" "Ergun equation", ) def gas_phase_config_pressure_drop(b, t, x): return -pyunits.convert( b.gas_phase.deltaP[t, x], to_units=units_meta_solid("pressure") / units_meta_solid("length"), ) == ( (150 * pyunits.dimensionless) * (1 - b.solid_properties[t, x].params.voidage) ** 2 * b.gas_phase.properties[t, x].visc_d * b.velocity_superficial_gas[t, x] / ( pyunits.convert( b.solid_properties[t, x].params.particle_dia, to_units=units_meta_solid("length"), ) ** 2 * b.solid_properties[t, x].params.voidage ** 3 ) ) + ( (1.75 * pyunits.dimensionless) * b.gas_phase.properties[t, x].dens_mass * (1 - b.solid_properties[t, x].params.voidage) * b.velocity_superficial_gas[t, x] ** 2 / ( pyunits.convert( b.solid_properties[t, x].params.particle_dia, to_units=units_meta_solid("length"), ) * b.solid_properties[t, x].params.voidage ** 3 ) ) elif self.config.has_pressure_change is False: pass else: raise BurntToast( "{} encountered unrecognized argument for " "the pressure drop correlation. Please contact the IDAES" " developers with this bug.".format(self.name) ) # --------------------------------------------------------------------- # Reaction constraints # Build homogeneous reaction constraints if gas_phase.reaction_package is not None: # Gas side rate reaction extent @self.Constraint( self.flowsheet().time, self.length_domain, gas_phase.reaction_package.rate_reaction_idx, doc="Gas side rate reaction extent", ) def gas_phase_config_rxn_ext(b, t, x, r): return b.gas_phase.rate_reaction_extent[t, x, r] == ( b.gas_phase.reactions[t, x].reaction_rate[r] * b.gas_phase.area[t, x] ) # Build hetereogeneous reaction constraints if solid_phase.reaction_package is not None: # Gas side heterogeneous rate reaction generation @self.Constraint( self.flowsheet().time, self.length_domain, gas_phase.property_package.phase_list, gas_phase.property_package.component_list, doc="Gas side heterogeneous" "rate reaction generation", ) def gas_comp_hetero_rxn(b, t, x, p, j): return b.gas_phase.mass_transfer_term[t, x, p, j] == ( sum( b.solid_reactions[t, x].rate_reaction_stoichiometry[r, p, j] * b.solid_reactions[t, x].reaction_rate[r] for r in (solid_phase.reaction_package.rate_reaction_idx) ) * b.solid_phase_area[t, x] ) # --------------------------------------------------------------------- # Solid phase component balance # material holdup constraint @self.Constraint( self.flowsheet().config.time, self.length_domain, solid_phase.property_package.component_list, doc="Solid phase material holdup constraints", ) def solid_material_holdup_calculation(b, t, x, j): return b.solid_material_holdup[t, x, j] == ( b.solid_phase_area[t, x] * b.solid_properties[t, x].get_material_density_terms("Sol", j) ) # Add component balances @self.Constraint( self.flowsheet().time, self.length_domain, solid_phase.property_package.component_list, doc="Material balances", ) def solid_material_balances(b, t, x, j): if solid_phase.reaction_package is not None: return b.solid_material_accumulation[t, x, j] == ( b.solid_phase_area[t, x] * b.solid_properties[t, x].params.mw_comp[j] * sum( b.solid_reactions[t, x].reaction_rate[r] * solid_phase.reaction_package.rate_reaction_stoichiometry[ r, "Sol", j ] for r in solid_phase.reaction_package.rate_reaction_idx ) ) else: return b.solid_material_accumulation[t, x, j] == 0 # Sum of mass fractions at all time equals 1 @self.Constraint( self.flowsheet().config.time, self.length_domain, doc="Sum of mass fractions at all time", ) def solid_sum_component_eqn(b, t, x): if t == b.flowsheet().config.time.first(): return Constraint.Skip else: return 1 == sum( b.solid_properties[t, x].mass_frac_comp[j] for j in b.solid_properties[t, x].params.component_list ) # --------------------------------------------------------------------- if self.config.energy_balance_type != EnergyBalanceType.none: # Solid phase - gas to solid heat transfer @self.Constraint( self.flowsheet().time, self.length_domain, doc="Solid phase - gas to solid heat transfer", ) def solid_phase_heat_transfer(b, t, x): return pyunits.convert( b.solid_phase_heat[t, x], to_units=units_meta_gas("power") / units_meta_gas("length"), ) * pyunits.convert( b.solid_properties[t, x].params.particle_dia, to_units=units_meta_gas("length"), ) == 6 * b.gas_solid_htc[ t, x ] * ( b.gas_phase.properties[t, x].temperature - pyunits.convert( b.solid_properties[t, x].temperature, to_units=units_meta_gas("temperature"), ) ) * pyunits.convert( b.solid_phase_area[t, x], to_units=units_meta_gas("area") ) # Solid phase energy balance # Accumulation equal to heat transfer @self.Constraint( self.flowsheet().config.time, self.length_domain, doc="Solid phase energy holdup constraints", ) def solid_energy_holdup_calculation(b, t, x): return b.solid_energy_holdup[t, x] == ( b.solid_phase_area[t, x] * pyunits.convert( b.solid_properties[t, x].get_energy_density_terms("Sol"), to_units=units_meta_solid("energy") / units_meta_solid("volume"), ) ) @self.Constraint( self.flowsheet().config.time, self.length_domain, doc="Solid phase energy balances", ) def solid_enthalpy_balances(b, t, x): if solid_phase.reaction_package is not None: return b.solid_energy_accumulation[t, x] / units_meta_solid( "time" ) == b.solid_phase_heat[t, x] + -sum( b.solid_reactions[t, x].reaction_rate[r] * b.solid_phase_area[t, x] * pyunits.convert( b.solid_reactions[t, x].dh_rxn[r], to_units=units_meta_solid("energy_mole"), ) for r in solid_phase.reaction_package.rate_reaction_idx ) else: return b.solid_energy_accumulation[t, x] == b.solid_phase_heat[t, x] # Dimensionless numbers, mass and heat transfer coefficients # Particle Reynolds number @self.Constraint( self.flowsheet().time, self.length_domain, doc="Particle Reynolds number", ) def reynolds_number_particle(b, t, x): return ( b.Re_particle[t, x] * b.gas_phase.properties[t, x].visc_d == b.velocity_superficial_gas[t, x] * pyunits.convert( b.solid_properties[t, x].params.particle_dia, to_units=units_meta_gas("length"), ) * b.gas_phase.properties[t, x].dens_mass ) # Prandtl number @self.Constraint( self.flowsheet().time, self.length_domain, doc="Prandtl number of gas in bed", ) def prandtl_number(b, t, x): return ( b.Pr_particle[t, x] * b.gas_phase.properties[t, x].therm_cond == b.gas_phase.properties[t, x].cp_mass * b.gas_phase.properties[t, x].visc_d ) # Particle Nusselt number @self.Constraint( self.flowsheet().time, self.length_domain, doc="Particle Nusselt number" ) def nusselt_number_particle(b, t, x): return ( b.Nu_particle[t, x] ** 3 == ( 2.0 + 1.1 * (smooth_abs(b.Re_particle[t, x], b.eps) ** 0.6) ** 3 ) * b.Pr_particle[t, x] ) # Gas-solid heat transfer coefficient @self.Constraint( self.flowsheet().time, self.length_domain, doc="Gas-solid heat transfer coefficient", ) def gas_solid_htc_eqn(b, t, x): return ( b.gas_solid_htc[t, x] * pyunits.convert( b.solid_properties[t, x].params.particle_dia, to_units=units_meta_gas("length"), ) == b.Nu_particle[t, x] * b.gas_phase.properties[t, x].therm_cond ) # Gas phase - gas to solid heat transfer @self.Constraint( self.flowsheet().time, self.length_domain, doc="Gas phase - gas to solid heat transfer", ) def gas_phase_heat_transfer(b, t, x): return b.gas_phase.heat[t, x] * pyunits.convert( b.solid_properties[t, x].params.particle_dia, to_units=units_meta_gas("length"), ) == -6 * b.gas_solid_htc[t, x] * ( b.gas_phase.properties[t, x].temperature - pyunits.convert( b.solid_properties[t, x].temperature, to_units=units_meta_gas("temperature"), ) ) * pyunits.convert( b.solid_phase_area[t, x], to_units=units_meta_gas("area") ) elif self.config.energy_balance_type == EnergyBalanceType.none: # if self.config.energy_balance_type == EnergyBalanceType.none: # If energy balance is none fix gas and solid temperatures to # initial solid temperature @self.Constraint( self.flowsheet().time, self.length_domain, doc="Isothermal gas phase constraint", ) def isothermal_gas_phase(b, t, x): # Skip constraint at the initial and boundary points t0 = b.flowsheet().config.time.first() if b.config.flow_type == "forward_flow": x_inlet = b.length_domain.first() else: x_inlet = b.length_domain.last() if x == x_inlet or t == b.flowsheet().config.time.first(): return Constraint.Skip else: return ( b.gas_phase.properties[t, x].temperature == b.solid_properties[t0, x].temperature ) @self.Constraint( self.flowsheet().config.time, self.length_domain, doc="Isothermal solid phase constraint", ) def isothermal_solid_phase(b, t, x): # Skip constraint at the initial point t0 = b.flowsheet().config.time.first() if t == b.flowsheet().config.time.first(): return Constraint.Skip else: return ( b.solid_properties[t, x].temperature == b.solid_properties[t0, x].temperature ) # Set default scaling values for some CV1D vars (to avoid warnings) if hasattr(self.gas_phase, "heat"): for v in self.gas_phase.heat.values(): if iscale.get_scaling_factor(v) is None: iscale.set_scaling_factor(v, 1e-6) if hasattr(self.gas_phase, "area"): for v in self.gas_phase.area.values(): if iscale.get_scaling_factor(v) is None: iscale.set_scaling_factor(v, 1) # ========================================================================= # Model initialization routine
[docs] def initialize_build( blk, gas_phase_state_args=None, solid_phase_state_args=None, outlvl=idaeslog.NOTSET, solver=None, optarg=None, ): """ Initialization routine for 1DFixedBed unit. Keyword Arguments: gas_phase_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 = None). solid_phase_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 = None). outlvl : sets output level of initialization routine optarg : solver options dictionary object (default=None, use default solver options) solver : str indicating which solver to use during initialization (default = None, use default solver) Returns: None """ # Set up logger for initialization and solve init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="unit") solve_log = idaeslog.getSolveLogger(blk.name, outlvl, tag="unit") # Create solver opt = get_solver(solver, optarg) # --------------------------------------------------------------------- # local aliases used to shorten object names gas_phase = blk.config.gas_phase_config solid_phase = blk.config.solid_phase_config # Keep all unit model geometry constraints, derivative_var constraints, # and property block constraints active. Additionally, in control # volumes - keep conservation linking constraints and # holdup calculation (for dynamic flowsheets) constraints active geometry_constraints_terms = [ "bed_area_eqn", "solid_phase_area_constraint", "gas_phase_area_constraint", ] endswith_terms = ( "_disc_eq", "linking_constraint", "linking_constraints", "_holdup_calculation", ) startswith_terms = ("properties", "solid_properties") for c in blk.component_objects(Constraint, descend_into=True): if ( not c.parent_block().local_name.startswith(startswith_terms) and not c.local_name.endswith(endswith_terms) and c.local_name not in geometry_constraints_terms ): c.deactivate() # --------------------------------------------------------------------- # Initialize thermophysical property constraints init_log.info("Initialize Thermophysical Properties") # Initialize gas_phase block gas_phase_flags = blk.gas_phase.properties.initialize( state_args=gas_phase_state_args, outlvl=outlvl, optarg=optarg, solver=solver, hold_state=True, ) # Initialize solid_phase properties block solid_phase_flags = blk.solid_properties.initialize( state_args=solid_phase_state_args, outlvl=outlvl, optarg=optarg, solver=solver, hold_state=True, ) init_log.info_high("Initialization Step 1 Complete.") # --------------------------------------------------------------------- # Initialize hydrodynamics (gas velocity) for t in blk.flowsheet().time: for x in blk.length_domain: calculate_variable_from_constraint( blk.velocity_superficial_gas[t, x], blk.gas_super_vel[t, x] ) blk.gas_super_vel.activate() init_log.info("Initialize Hydrodynamics") with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: results = opt.solve(blk, tee=slc.tee) if check_optimal_termination(results): init_log.info_high( "Initialization Step 2 {}.".format(idaeslog.condition(results)) ) else: _log.warning("{} Initialization Step 2 Failed.".format(blk.name)) # --------------------------------------------------------------------- # Initialize mass balance - no reaction and no pressure drop # Unfix material balance state variables (including particle porosity) # but keep other states fixed blk.gas_phase.properties.release_state(gas_phase_flags) blk.solid_properties.release_state(solid_phase_flags) for t in blk.flowsheet().time: for x in blk.length_domain: blk.gas_phase.properties[t, x].pressure.fix() blk.gas_phase.properties[t, x].temperature.fix() blk.solid_properties[t, x].temperature.fix() if blk.config.has_holdup is True: # Fix initial conditions of flowrate t0 = blk.flowsheet().time.first() if blk.config.flow_type == "forward_flow": x_inlet = blk.length_domain.first() else: x_inlet = blk.length_domain.last() for x in blk.length_domain: if x != x_inlet: blk.gas_phase.properties[t0, x].flow_mol.fix() blk.gas_phase.properties[t0, x].sum_component_eqn.deactivate() blk.gas_phase.material_balances.activate() if gas_phase.reaction_package is not None: for t in blk.flowsheet().time: gas_rxn_gen = blk.gas_phase.rate_reaction_generation for x in blk.length_domain: for p in gas_phase.property_package.phase_list: for j in gas_phase.property_package.component_list: (gas_rxn_gen[t, x, p, j].fix(0.0)) blk.solid_material_balances.activate() blk.solid_sum_component_eqn.activate() if solid_phase.reaction_package is not None: for t in blk.flowsheet().time: for x in blk.length_domain: for r in solid_phase.reaction_package.rate_reaction_idx: blk.solid_reactions[t, x].reaction_rate[r].fix(0.0) # Gas side heterogeneous rate reaction generation for x in blk.length_domain: for p in gas_phase.property_package.phase_list: for j in gas_phase.property_package.component_list: (blk.gas_phase.mass_transfer_term[t, x, p, j].fix(0.0)) init_log.info("Initialize Mass Balances") init_log.info_high( "initialize mass balances - no reactions and no pressure drop" ) with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: results = opt.solve(blk, tee=slc.tee, symbolic_solver_labels=True) if check_optimal_termination(results): init_log.info_high( "Initialization Step 3a {}.".format(idaeslog.condition(results)) ) else: _log.warning("{} Initialization Step 3a Failed.".format(blk.name)) # Initialize mass balance - with reaction and no pressure drop if gas_phase.reaction_package is not None: # local aliases used to shorten object names gas_rxn_gen = blk.gas_phase.rate_reaction_generation gas_phase_stoichiometry_eqn = ( blk.gas_phase.rate_reaction_stoichiometry_constraint ) # Initialize reaction property package blk.gas_phase.reactions.activate() for t in blk.flowsheet().time: for x in blk.length_domain: obj = blk.gas_phase.reactions[t, x] for c in obj.component_objects(Constraint, descend_into=False): c.activate() blk.gas_phase.reactions.initialize( outlvl=outlvl, optarg=optarg, solver=solver ) for t in blk.flowsheet().time: for x in blk.length_domain: for r in gas_phase.reaction_package.rate_reaction_idx: calculate_variable_from_constraint( blk.gas_phase.rate_reaction_extent[t, x, r], blk.gas_phase_config_rxn_ext[t, x, r], ) for p in gas_phase.property_package.phase_list: for j in gas_phase.property_package.component_list: (gas_rxn_gen[t, x, p, j].unfix()) if not ( ( blk.gas_phase.config.transformation_scheme != "FORWARD" and x == blk.length_domain.first() ) or ( blk.gas_phase.config.transformation_scheme == "FORWARD" and x == blk.length_domain.last() ) ): calculate_variable_from_constraint( gas_rxn_gen[t, x, p, j], gas_phase_stoichiometry_eqn[t, x, p, j], ) gas_phase_stoichiometry_eqn.activate() blk.gas_phase_config_rxn_ext.activate() if solid_phase.reaction_package is not None: gas_mass_transfer_term = blk.gas_phase.mass_transfer_term # Initialize reaction property package blk.solid_reactions.activate() for t in blk.flowsheet().time: for x in blk.length_domain: blk.solid_reactions[t, x].reaction_rate.unfix() obj = blk.solid_reactions[t, x] for c in obj.component_objects(Constraint, descend_into=False): c.activate() blk.solid_reactions.initialize(outlvl=outlvl, optarg=optarg, solver=solver) for t in blk.flowsheet().time: for x in blk.length_domain: for p in gas_phase.property_package.phase_list: for j in gas_phase.property_package.component_list: (gas_mass_transfer_term[t, x, p, j].unfix()) calculate_variable_from_constraint( gas_mass_transfer_term[t, x, p, j], blk.gas_comp_hetero_rxn[t, x, p, j], ) blk.gas_comp_hetero_rxn.activate() if ( gas_phase.reaction_package is not None or solid_phase.reaction_package is not None ): init_log.info_high( "initialize mass balances - with reactions and no pressure drop" ) with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: results = opt.solve(blk, tee=slc.tee, symbolic_solver_labels=True) if check_optimal_termination(results): init_log.info_high( "Initialization Step 3b {}.".format(idaeslog.condition(results)) ) else: _log.warning("{} Initialization Step 3b Failed.".format(blk.name)) # Initialize mass balance - with pressure drop for t in blk.flowsheet().time: for x in blk.length_domain: # Unfix all pressure variables except at the inlet if blk.gas_phase.properties[t, x].config.defined_state is False: blk.gas_phase.properties[t, x].pressure.unfix() blk.gas_phase.pressure_balance.activate() if blk.config.has_pressure_change: blk.gas_phase_config_pressure_drop.activate() for t in blk.flowsheet().time: for x in blk.length_domain: calculate_variable_from_constraint( blk.gas_phase.deltaP[t, x], blk.gas_phase_config_pressure_drop[t, x], ) init_log.info_high( "initialize mass balances - with reactions and pressure drop" ) with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: results = opt.solve(blk, tee=slc.tee, symbolic_solver_labels=True) if check_optimal_termination(results): init_log.info_high( "Initialization Step 3c {}.".format(idaeslog.condition(results)) ) else: _log.warning("{} Initialization Step 3c Failed.".format(blk.name)) # --------------------------------------------------------------------- # Initialize energy balance calc_var_kwds = {"eps": 5e-6} if blk.config.energy_balance_type != EnergyBalanceType.none: # Unfix temperatures if blk.config.flow_type == "forward_flow": x_inlet = blk.length_domain.first() else: x_inlet = blk.length_domain.last() for t in blk.flowsheet().time: for x in blk.length_domain: if t != blk.flowsheet().config.time.first() and x != x_inlet: # Unfix gas temperature variables except at the inlet blk.gas_phase.properties[t, x].temperature.unfix() for t in blk.flowsheet().time: if t != blk.flowsheet().config.time.first(): for x in blk.length_domain: # Unfix solid temperature variables except at the inlet blk.solid_properties[t, x].temperature.unfix() # Initialize dimensionless numbers, # mass and heat transfer coefficients for t in blk.flowsheet().time: for x in blk.length_domain: # Initialize gas temperature to solid temperature values calculate_variable_from_constraint( blk.Re_particle[t, x], blk.reynolds_number_particle[t, x], **calc_var_kwds, ) calculate_variable_from_constraint( blk.Pr_particle[t, x], blk.prandtl_number[t, x], **calc_var_kwds ) calculate_variable_from_constraint( blk.Nu_particle[t, x], blk.nusselt_number_particle[t, x], **calc_var_kwds, ) calculate_variable_from_constraint( blk.gas_solid_htc[t, x], blk.gas_solid_htc_eqn[t, x], **calc_var_kwds, ) calculate_variable_from_constraint( blk.gas_phase.heat[t, x], blk.gas_phase_heat_transfer[t, x], **calc_var_kwds, ) calculate_variable_from_constraint( blk.solid_phase_heat[t, x], blk.solid_phase_heat_transfer[t, x], **calc_var_kwds, ) if t != blk.flowsheet().time.first(): calculate_variable_from_constraint( blk.gas_phase.energy_accumulation[t, x, "Vap"], blk.gas_phase.energy_accumulation_disc_eq[t, x, "Vap"], **calc_var_kwds, ) calculate_variable_from_constraint( blk.solid_energy_accumulation[t, x], blk.solid_energy_accumulation_disc_eq[t, x], **calc_var_kwds, ) blk.reynolds_number_particle.activate() blk.prandtl_number.activate() blk.nusselt_number_particle.activate() blk.gas_solid_htc_eqn.activate() # Activate gas phase energy balance equations blk.gas_phase_heat_transfer.activate() blk.gas_phase.enthalpy_balances.activate() # Activate solid phase energy balance equations blk.solid_phase_heat_transfer.activate() blk.solid_enthalpy_balances.activate() init_log.info("Initialize Energy Balances") with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: results = opt.solve(blk, tee=slc.tee, symbolic_solver_labels=True) if check_optimal_termination(results): init_log.info_high( "Initialization Step 4 {}.".format(idaeslog.condition(results)) ) else: raise InitializationError( f"{blk.name} failed to initialize successfully. Please " f"check the output logs for more information." ) # Initialize energy balance if blk.config.energy_balance_type == EnergyBalanceType.none: if blk.config.flow_type == "forward_flow": x_inlet = blk.length_domain.first() else: x_inlet = blk.length_domain.last() t0 = blk.flowsheet().config.time.first() for t in blk.flowsheet().time: for x in blk.length_domain: if t != blk.flowsheet().config.time.first() and x != x_inlet: # Unfix gas temperature variables except at the inlet blk.gas_phase.properties[t, x].temperature.unfix() blk.gas_phase.properties[t, x].temperature = value( blk.solid_properties[t0, x].temperature ) for t in blk.flowsheet().time: t0 = blk.flowsheet().config.time.first() if t != blk.flowsheet().config.time.first(): for x in blk.length_domain: # Unfix solid temperature variables except initial blk.solid_properties[t, x].temperature = value( blk.solid_properties[t0, x].temperature ) blk.solid_properties[t, x].temperature.unfix() blk.isothermal_gas_phase.activate() blk.isothermal_solid_phase.activate() init_log.info("Initialize Energy Balances") with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: results = opt.solve(blk, tee=slc.tee) if check_optimal_termination(results): init_log.info_high( "Initialization Step 4 {}.".format(idaeslog.condition(results)) ) else: raise InitializationError( f"{blk.name} failed to initialize successfully. Please " f"check the output logs for more information." )
[docs] def block_triangularization_initialize( blk, gas_phase_state_args=None, solid_phase_state_args=None, outlvl=idaeslog.NOTSET, solver=None, calc_var_kwds=None, ): """ Block triangularization (BT) initialization routine for 1D FixedBed unit. The BT initialization routine initializes the unit model in the following steps: Step 1: Convert the system of equations from a partial differential set of equations (PDAE) to an algebraic set of equations (AE) by deactivating the discretization equations and sum_component_equations (if present), and fixing the state variables to a guess (typically inlet conditions). Step 2: Decompose the AE into strongly connected components via block triangularization and solve individual blocks Step 3: revert the fixed variables and deactivated constraints to their original states. The initialization procedure is wrapped in the TemporarySubsystemManager utility to ensure that the deactivation of constraints and fixing of variables are only temporary and the final model structure is unchanged after initialization. Keyword Arguments: gas_phase_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 = None). solid_phase_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 = None). outlvl : sets output level of initialization routine optarg : solver options dictionary object (default=None, use default solver options) solver : str indicating which solver to use during initialization (default = None, use default solver) calc_var_kwds: Dictionary of keyword arguments for calculate_variable_from_constraint Returns: None """ # Set up logger for initialization and solve init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="unit") # ========================================================================= # Check that the original model has zero degrees of freedom if degrees_of_freedom(blk) != 0: msg = ( "Original model has nonzero degrees of freedom. This was " "unexpected. " + "\n" + "Degrees of freedom: " + str(degrees_of_freedom(blk)) ) init_log.error(msg) raise ValueError("Nonzero degrees of freedom.") # ========================================================================= # Create list of constraints to deactivate to_deactivate = [] # Deactivate gas/solid phase sum equations and discretization equations # (mass/energy flow, mass/energy accumulation, pressure) endswith_terms = ("disc_eq", "sum_component_eqn") for c in blk.component_objects(Constraint, descend_into=True): if c.local_name.endswith(endswith_terms): c.attribute_errors_generate_exceptions = False to_deactivate.append(c) # Additional equations to deactivate if energy balance type is none if hasattr(blk, "isothermal_gas_phase"): isothermal_gas_phase_eqn = blk.isothermal_gas_phase[...] isothermal_gas_phase_eqn.attribute_errors_generate_exceptions = False to_deactivate.extend(isothermal_gas_phase_eqn) if hasattr(blk, "isothermal_solid_phase"): isothermal_solid_phase_eqn = blk.isothermal_solid_phase[...] isothermal_solid_phase_eqn.attribute_errors_generate_exceptions = False to_deactivate.extend(isothermal_solid_phase_eqn) # ========================================================================= # Variables to fix during initialization # List of state vars to fix (typically fixed to boundary/initial conditions) gas_flags = fix_state_vars(blk.gas_phase.properties, gas_phase_state_args) solid_flags = fix_state_vars(blk.solid_properties, solid_phase_state_args) # List of additional variables to fix to_fix = [] # Flow derivative variables are fixed to ensure square AE set_value = 0 time_set = blk.flowsheet().time for t in blk.flowsheet().time: for var in ComponentSet(blk.component_objects(Var)): if var.name.endswith("_flow_dx"): n = var.index_set().dimen if n == 1: var.set_value(set_value) to_fix.append(var[t]) elif n >= 2: info = get_index_set_except(var, time_set) non_time_set = info["set_except"] index_getter = info["index_getter"] for non_time_index in non_time_set: index = index_getter(non_time_index, t) var[index].set_value(set_value) to_fix.append(var[index]) # ========================================================================= with TemporarySubsystemManager(to_fix=to_fix, to_deactivate=to_deactivate): # Check if the system is structurally singular igraph = IncidenceGraphInterface(blk) N = len(igraph.constraints) M = len(igraph.variables) matching = igraph.maximum_matching() if not (N == len(matching) and M == len(matching)): msg = ( "Model is structurally singular as maximum matching " "does not include all costraints and variables " "\n" "Number of constraints: " + str(N) + "\n" + "Number of variables: " + str(M) + "\n" + "Maximum matching: " + str(len(matching)) ) init_log.warning(msg) raise ValueError("Structural singularity") # Decompose AE system into strongly connected components and solve if calc_var_kwds is None: calc_var_kwds = {"eps": 1e-8} try: solve_strongly_connected_components( blk, calc_var_kwds=calc_var_kwds, solver=solver ) except RuntimeError: msg = ( "Initial value for variable results in a derivative value " "that is very close to zero. Please provide a different initial " "guess and/or adjust calc_var_kwds tolerance setting." ) init_log.warning(msg) raise # Revert the state vars to their original state revert_state_vars(blk.gas_phase.properties, gas_flags) revert_state_vars(blk.solid_properties, solid_flags)
def calculate_scaling_factors(self): super().calculate_scaling_factors() # local aliases used to shorten object names gas_phase = self.config.gas_phase_config solid_phase = self.config.solid_phase_config # Get units meta data from property packages units_meta_gas = gas_phase.property_package.get_metadata().get_derived_units # scale some variables if hasattr(self, "bed_diameter"): sf = 1 / value(self.bed_diameter) iscale.set_scaling_factor(self.bed_diameter, 1e-1 * sf) if hasattr(self, "bed_area"): sf = 1 / value(constants.pi * (0.5 * self.bed_diameter) ** 2) iscale.set_scaling_factor(self.bed_area, 1e-1 * sf) if hasattr(self, "solid_phase_area"): sf = iscale.get_scaling_factor(self.bed_area) iscale.set_scaling_factor(self.solid_phase_area, sf) if hasattr(self.gas_phase, "area"): sf = iscale.get_scaling_factor(self.bed_area) iscale.set_scaling_factor(self.gas_phase.area, sf) if self.config.energy_balance_type != EnergyBalanceType.none: if hasattr(self, "Re_particle"): for (t, x), v in self.Re_particle.items(): if iscale.get_scaling_factor(v) is None: sf1 = iscale.get_scaling_factor( self.gas_phase.properties[t, x].dens_mass ) sf2 = 1 / value( pyunits.convert( self.solid_properties[t, x].params.particle_dia, to_units=units_meta_gas("length"), ) ) iscale.set_scaling_factor(v, 1e-1 * sf1 * sf2) if hasattr(self, "solid_material_accumulation"): for (t, x, j), v in self.solid_material_accumulation.items(): if iscale.get_scaling_factor(v) is None: sf = iscale.get_scaling_factor( self.solid_properties[t, x].get_material_density_terms( "Sol", j ), default=1, warning=True, ) iscale.set_scaling_factor(v, sf) if self.config.energy_balance_type != EnergyBalanceType.none: if hasattr(self, "Pr_particle"): for (t, x), v in self.Pr_particle.items(): if iscale.get_scaling_factor(v) is None: sf1 = iscale.get_scaling_factor( self.gas_phase.properties[t, x].cp_mass ) sf2 = iscale.get_scaling_factor( self.gas_phase.properties[t, x].visc_d ) iscale.set_scaling_factor(v, sf1 * sf2) if hasattr(self, "Nu_particle"): for (t, x), v in self.Nu_particle.items(): if iscale.get_scaling_factor(v) is None: sf = 1 / ( value( ( ( 2.0 + 1.1 * ( smooth_abs(self.Re_particle[t, x], self.eps) ** 0.6 ) ** 3 ) * self.Pr_particle[t, x] ) ** (1 / 3) ) ) iscale.set_scaling_factor(v, 1e-2 * sf) if hasattr(self, "gas_solid_htc"): for (t, x), v in self.gas_solid_htc.items(): if iscale.get_scaling_factor(v) is None: sf1 = iscale.get_scaling_factor(self.Nu_particle[t, x]) sf2 = iscale.get_scaling_factor( self.gas_phase.properties[t, x].therm_cond ) iscale.set_scaling_factor(v, sf1 * sf2) if hasattr(self, "solid_energy_accumulation"): for (t, x), v in self.solid_energy_accumulation.items(): if iscale.get_scaling_factor(v) is None: sf = iscale.get_scaling_factor( self.solid_properties[t, x].get_energy_density_terms("Sol"), default=1, warning=True, ) iscale.set_scaling_factor(v, sf) if hasattr(self, "solid_energy_holdup"): for (t, x), v in self.solid_energy_holdup.items(): sf1 = iscale.get_scaling_factor(self.solid_phase_area[t, x]) sf2 = iscale.get_scaling_factor( self.solid_properties[t, x].get_energy_density_terms("Sol") ) iscale.set_scaling_factor(v, sf1 * sf2) # scale heat variables in CV1D for (t, x), v in self.gas_phase.heat.items(): sf1 = iscale.get_scaling_factor( self.gas_phase.properties[t, x].enth_mol ) sf2 = iscale.get_scaling_factor( self.gas_phase.properties[t, x].flow_mol ) iscale.set_scaling_factor(v, sf1 * sf2) for (t, x), v in self.solid_phase_heat.items(): sf1 = iscale.get_scaling_factor(self.solid_properties[t, x].enth_mass) sf2 = iscale.get_scaling_factor( self.solid_properties[t, x].dens_mass_particle ) sf3 = iscale.get_scaling_factor(self.solid_phase_area[t, x]) iscale.set_scaling_factor(v, sf1 * sf2 * sf3) # Scale some constraints if hasattr(self, "bed_area_eqn"): for c in self.bed_area_eqn.values(): iscale.constraint_scaling_transform( c, iscale.get_scaling_factor(self.bed_area), overwrite=False ) if hasattr(self, "gas_phase_area_constraint"): for (t, x), c in self.gas_phase_area_constraint.items(): iscale.constraint_scaling_transform( c, iscale.get_scaling_factor(self.bed_area), overwrite=False ) if hasattr(self, "solid_phase_area_constraint"): for (t, x), c in self.solid_phase_area_constraint.items(): iscale.constraint_scaling_transform( c, iscale.get_scaling_factor(self.bed_area), overwrite=False ) if hasattr(self, "gas_super_vel"): for (t, x), c in self.gas_super_vel.items(): iscale.constraint_scaling_transform( c, iscale.get_scaling_factor(self.gas_phase.properties[t, x].flow_mol), overwrite=False, ) if hasattr(self, "gas_phase_config_pressure_drop"): for (t, x), c in self.gas_phase_config_pressure_drop.items(): iscale.constraint_scaling_transform( c, iscale.get_scaling_factor(self.gas_phase.deltaP[t, x]), overwrite=False, ) if hasattr(self, "gas_phase_config_rxn_ext"): for (t, x, r), c in self.gas_phase_config_rxn_ext.items(): sf1 = iscale.get_scaling_factor( self.gas_phase.reactions[t, x].reaction_rate[r] ) sf2 = iscale.get_scaling_factor(self.bed_area) iscale.constraint_scaling_transform(c, sf1 * sf2, overwrite=False) if hasattr(self, "gas_comp_hetero_rxn"): for (t, x, p, j), c in self.gas_comp_hetero_rxn.items(): iscale.constraint_scaling_transform( c, 1e-3 * iscale.get_scaling_factor( self.gas_phase.mass_transfer_term[t, x, p, j] ), overwrite=False, ) if self.config.energy_balance_type != EnergyBalanceType.none: if hasattr(self, "reynolds_number_particle"): for (t, x), c in self.reynolds_number_particle.items(): iscale.constraint_scaling_transform( c, iscale.get_scaling_factor(self.Re_particle[t, x]), overwrite=False, ) if hasattr(self, "prandtl_number"): for (t, x), c in self.prandtl_number.items(): iscale.constraint_scaling_transform( c, iscale.get_scaling_factor(self.Pr_particle[t, x]), overwrite=False, ) if hasattr(self, "nusselt_number_particle"): for (t, x), c in self.nusselt_number_particle.items(): iscale.constraint_scaling_transform( c, iscale.get_scaling_factor(self.Nu_particle[t, x]), overwrite=False, ) if hasattr(self, "gas_solid_htc_eqn"): for (t, x), c in self.gas_solid_htc_eqn.items(): iscale.constraint_scaling_transform( c, iscale.get_scaling_factor(self.gas_solid_htc[t, x]), overwrite=False, ) if hasattr(self, "gas_phase_heat_transfer"): for (t, x), c in self.gas_phase_heat_transfer.items(): iscale.constraint_scaling_transform( c, 1e-2 * iscale.get_scaling_factor(self.gas_phase.heat[t, x]), overwrite=False, ) if hasattr(self, "solid_phase_heat_transfer"): for (t, x), c in self.solid_phase_heat_transfer.items(): iscale.constraint_scaling_transform( c, 1e-2 * iscale.get_scaling_factor(self.solid_phase_heat[t, x]), overwrite=False, ) elif self.config.energy_balance_type == EnergyBalanceType.none: if hasattr(self, "isothermal_gas_phase"): for (t, x), c in self.isothermal_gas_phase.items(): iscale.constraint_scaling_transform( c, iscale.get_scaling_factor( self.gas_phase.properties[t, x].temperature ), overwrite=False, ) if hasattr(self, "isothermal_solid_phase"): for (t, x), c in self.isothermal_solid_phase.items(): iscale.constraint_scaling_transform( c, iscale.get_scaling_factor( self.solid_properties[t, x].temperature ), overwrite=False, ) # Scale some constraints if hasattr(self, "solid_material_holdup_calculation"): for (t, x, j), c in self.solid_material_holdup_calculation.items(): sf1 = iscale.get_scaling_factor(self.solid_phase_area[t, x]) sf2 = iscale.get_scaling_factor( self.solid_properties[t, x].dens_mass_particle ) iscale.constraint_scaling_transform(c, sf1 * sf2, overwrite=False) if hasattr(self, "solid_material_accumulation_disc_eq"): for (t, x, j), c in self.solid_material_accumulation_disc_eq.items(): iscale.constraint_scaling_transform( c, 1e-6 * iscale.get_scaling_factor( self.solid_material_accumulation[t, x, j] ), overwrite=False, ) if hasattr(self, "solid_material_balances"): component_list = solid_phase.property_package.component_list # Get a single representative value in component list for j in component_list: break if solid_phase.reaction_package is not None: rate_reaction_idx = solid_phase.reaction_package.rate_reaction_idx # Get a single representative value in rate_reaction index list for r in rate_reaction_idx: break for (t, x, j), c in self.solid_material_balances.items(): sf1 = iscale.get_scaling_factor(self.solid_phase_area[t, x]) sf2 = iscale.get_scaling_factor( self.solid_properties[t, x].params.mw_comp[j], default=1e1, ) sf3 = iscale.get_scaling_factor( self.solid_reactions[t, x].reaction_rate[r] ) iscale.constraint_scaling_transform( c, sf1 * sf2 * sf3, overwrite=False ) else: for (t, x, j), c in self.solid_material_balances.items(): sf1 = iscale.get_scaling_factor(self.solid_phase_area[t, x]) sf2 = iscale.get_scaling_factor( self.solid_properties[t, x].params.mw_comp[j], default=1e1, ) iscale.constraint_scaling_transform(c, sf1 * sf2, overwrite=False) if hasattr(self, "solid_sum_component_eqn"): component_list = solid_phase.property_package.component_list # Get a single representative value in component list for j in component_list: break for (t, x), c in self.solid_sum_component_eqn.items(): iscale.constraint_scaling_transform( c, iscale.get_scaling_factor( self.solid_properties[t, x].mass_frac_comp[j] ), overwrite=False, ) if self.config.energy_balance_type != EnergyBalanceType.none: if hasattr(self, "solid_energy_holdup_calculation"): for (t, x), c in self.solid_energy_holdup_calculation.items(): sf1 = iscale.get_scaling_factor(self.solid_phase_area[t, x]) sf2 = iscale.get_scaling_factor( self.solid_properties[t, x].get_energy_density_terms("Sol") ) iscale.constraint_scaling_transform(c, sf1 * sf2, overwrite=False) if hasattr(self, "solid_enthalpy_balances"): for (t, x), c in self.solid_enthalpy_balances.items(): sf1 = iscale.get_scaling_factor(self.solid_phase_area[t, x]) sf2 = iscale.get_scaling_factor( self.solid_properties[t, x].enth_mass ) iscale.constraint_scaling_transform( c, 1e-2 * sf1 * sf2, overwrite=False ) if hasattr(self.gas_phase, "enthalpy_balances"): phase_list = self.gas_phase.properties.phase_list for (t, x), c in self.gas_phase.enthalpy_balances.items(): # This is OK sf = iscale.min_scaling_factor( [ # pylint: disable-next=protected-access self.gas_phase._enthalpy_flow[t, x, p] for p in phase_list ] ) iscale.constraint_scaling_transform(c, 1e-3 * sf, overwrite=True) if hasattr(self, "solid_energy_accumulation_disc_eq"): for (t, x), c in self.solid_energy_accumulation_disc_eq.items(): iscale.constraint_scaling_transform( c, 1e-3 * iscale.get_scaling_factor( self.solid_energy_accumulation[t, x] ), overwrite=False, ) if hasattr(self.gas_phase, "energy_accumulation_disc_eq"): for (t, x, p), c in self.gas_phase.energy_accumulation_disc_eq.items(): iscale.constraint_scaling_transform( c, 1e-3 * iscale.get_scaling_factor( self.gas_phase.energy_accumulation[t, x, p] ), overwrite=True, ) def _get_stream_table_contents(self, time_point=0): return create_stream_table_dataframe( {"Gas Inlet": self.gas_inlet, "Gas Outlet": self.gas_outlet}, time_point=time_point, ) def _get_performance_contents(self, time_point=0): var_dict = {} var_dict["Bed Height"] = self.bed_height var_dict["Bed Area"] = self.bed_area var_dict["Gas Inlet Velocity"] = self.velocity_superficial_gas[time_point, 0] var_dict["Gas Outlet Velocity"] = self.velocity_superficial_gas[time_point, 1] return {"vars": var_dict}