From 68e2c4f93149c9466c11d78cff99ea1ee8fcd8ab Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Thu, 3 Aug 2023 15:48:59 -0400 Subject: [PATCH] Add support for heterogeneous reactions in MSContactor (#1229) * Fixing bugs in initialization * Working on heterogeneous rxns * Working on heterogeneous reactions * Working on tests * More testing * Updating docs * Fixing typos * Fixing f-strings * Fixing typo in docs * Fixing more typos in docs --- .../generic/unit_models/mscontactor.rst | 214 ++++++++++++++---- idaes/models/unit_models/mscontactor.py | 155 ++++++++++++- .../unit_models/tests/test_mscontactor.py | 212 ++++++++++++++++- 3 files changed, 525 insertions(+), 56 deletions(-) diff --git a/docs/reference_guides/model_libraries/generic/unit_models/mscontactor.rst b/docs/reference_guides/model_libraries/generic/unit_models/mscontactor.rst index 0b05b4f1f8..b50a52c08e 100644 --- a/docs/reference_guides/model_libraries/generic/unit_models/mscontactor.rst +++ b/docs/reference_guides/model_libraries/generic/unit_models/mscontactor.rst @@ -39,16 +39,19 @@ Degrees of Freedom As a general purpose model, the degrees of freedom of the multi-stream contactor models depend on the options chosen by the user. The potential degrees of freedom are: * states for feed blocks for each stream, -* material transfer terms (finite elements :math:`\times` interacting streams :math:`\times` common components), -* energy transfer terms if included (finite elements :math:`\times` interacting streams) -* pressure change terms if included (finite elements :math:`\times` streams with pressure change) +* material transfer terms (time points :math:`\times` finite elements :math:`\times` interacting streams :math:`\times` common components), +* energy transfer terms if included (time points :math:`\times` finite elements :math:`\times` interacting streams) +* pressure change terms if included (time points :math:`\times` finite elements :math:`\times` streams with pressure change) +* reaction extent terms for rate based and heterogeneous reactions if included (time points :math:`\times` finite elements :math:`\times` number of reactions) Model Structure --------------- Due to the custom nature of multi-stream contactors, this model does not make use of control volumes. Instead, a set of StateBlocks (named using the name given in the ``streams`` configuration dictionary) are created for each stream indexed by time and the set of finite elements, with an additional StateBlock for the feed state (named using the stream name appended with ``_inlet_state``) indexed only by time (unless ``has_feed`` is set to ``False`` for that stream). For streams with side streams (feed or draw), an additional set of indexed StateBlocks (named using the stream name appended with ``_side_stream_state``) is created for the side states which are indexed by time and the set of side states for that stream. -If reactions are required for a given stream, a set of indexed ReactionBlocks (named using the stream name appended with ``_reactions``) are created indexed by time and the set of finite elements. +If reactions are required for a given stream, a set of indexed ReactionBlocks (named using the stream name appended with ``_reactions``) are created indexed by time and the set of finite elements. The ``MSContactor`` model also supports the concept of heterogeneous reactions which involve species across multiple streams; in these cases users can declare a single "heterogeneous reaction block" for the system in which case a Heterogeneous ReactionBlock will be created for each finite element (see later for more details on defining heterogeneous reaction blocks). + +**Note:** due to the custom nature of multi-stream contactors, the ``MSContactor`` model **does not** define a constraint for the extent of reaction as the basis for this is not well defined (e.g. volume can be defined either by the volume of each stream/phase in a finite element or the total volume in each element. Users are required to add the necessary constraints appropriate to their system. All other variables and constraints are written at the unit model level. @@ -57,41 +60,43 @@ Variables The multi-stream contactor creates the following variables. Here ``t`` indicates the time domain and ``x`` indicates finite element. -============================ =========================================== ============================================================================================================ ========================================== -Variable Name Description Notes -============================ =========================================== ============================================================================================================ ========================================== -:math:`M_{t,x,s1,s2,j}` material_transfer_term Material transfer term for component ``j`` between stream ``s1`` and ``s2`` at ``x`` and ``t`` -:math:`E_{t,x,s1,s2}` energy_transfer_term Energy transfer term between stream ``s1`` and ``s2`` at ``x`` and ``t`` Only if energy balances included -:math:`Q_{t,x,s}` stream + "_heat" External heat transfer into stream ``s`` at ``x`` and ``t`` Only if ``has_heat_transfer`` for stream -:math:`\Delta P_{t,x,s}` stream + "_deltaP" Pressure change in stream ``s`` at ``x`` and ``t`` Only if ``has_pressure_change`` for stream -:math:`G_{rate,t,x,s,p,j}` stream + "_rate_reaction_generation" Generation of component ``j`` in phase ``p`` due to rate reactions in stream ``s`` at ``x`` ``t`` Only if rate reactions present for stream -:math:`G_{equil,t,x,s,p,j}` stream + "_equilibrium_reaction_generation" Generation of component ``j`` in phase ``p`` due to equilibrium reactions in stream ``s`` at ``x`` and ``t`` Only if equilibrium reactions present for stream -:math:`G_{inher,t,x,s,p,j}` stream + "_inherent_reaction_generation" Generation of component ``j`` in phase ``p`` due to inherent reactions in stream ``s`` at ``x`` and ``t`` Only if inherent reactions present for stream -:math:`X_{rate,t,x,s,r}` stream + "_rate_reaction_extent" Extent of rate reaction ``r`` in stream ``s`` at ``x`` and ``t`` Only if rate reactions present for stream -:math:`X_{equil,t,x,s,r}` stream + "_equilibrium_reaction_extent" Extent of equilibrium reaction ``r`` in stream ``s`` at ``x`` and ``t`` Only if equilibrium reactions present for stream -:math:`X_{inher,t,x,s,r}` stream + "_inherent_reaction_extent" Extent of inherent reaction ``r`` in stream ``s`` at ``x`` and ``t`` Only if inherent reactions present for stream -============================ =========================================== ============================================================================================================ ========================================== +============================= ============================================== ============================================================================================================= ========================================= +Variable Name Description Notes +============================= ============================================== ============================================================================================================= ========================================= +:math:`M_{t,x,s1,s2,j}` material_transfer_term Material transfer term for component ``j`` between stream ``s1`` and ``s2`` at ``x`` and ``t`` +:math:`E_{t,x,s1,s2}` energy_transfer_term Energy transfer term between stream ``s1`` and ``s2`` at ``x`` and ``t`` Only if energy balances included +:math:`Q_{t,x,s}` stream + "_heat" External heat transfer into stream ``s`` at ``x`` and ``t`` Only if ``has_heat_transfer`` for stream +:math:`\Delta P_{t,x,s}` stream + "_deltaP" Pressure change in stream ``s`` at ``x`` and ``t`` Only if ``has_pressure_change`` for stream +:math:`G_{rate,t,x,s,p,j}` stream + "_rate_reaction_generation" Generation of component ``j`` in phase ``p`` due to rate reactions in stream ``s`` at ``x`` ``t`` Only if rate reactions present for stream +:math:`G_{equil,t,x,s,p,j}` stream + "_equilibrium_reaction_generation" Generation of component ``j`` in phase ``p`` due to equilibrium reactions in stream ``s`` at ``x`` and ``t`` Only if equilibrium reactions present for stream +:math:`G_{inher,t,x,s,p,j}` stream + "_inherent_reaction_generation" Generation of component ``j`` in phase ``p`` due to inherent reactions in stream ``s`` at ``x`` and ``t`` Only if inherent reactions present for stream +:math:`G_{hetero,t,x,s,p,j}` stream + "_heterogeneous_reaction_generation" Generation of component ``j`` in phase ``p`` due to heterogeneous reactions in stream ``s`` at ``x`` ``t`` Only if heterogeneous reactions present +:math:`X_{rate,t,x,s,r}` stream + "_rate_reaction_extent" Extent of rate reaction ``r`` in stream ``s`` at ``x`` and ``t`` Only if rate reactions present for stream +:math:`X_{equil,t,x,s,r}` stream + "_equilibrium_reaction_extent" Extent of equilibrium reaction ``r`` in stream ``s`` at ``x`` and ``t`` Only if equilibrium reactions present for stream +:math:`X_{inher,t,x,s,r}` stream + "_inherent_reaction_extent" Extent of inherent reaction ``r`` in stream ``s`` at ``x`` and ``t`` Only if inherent reactions present for stream +:math:`X_{hetero,t,x,r}` heterogeneous_reaction_extent Extent of heterogeneous reaction ``r`` at ``x`` and ``t`` Only if heterogeneous reactions present +============================= ============================================== ============================================================================================================= ========================================= Constraints ----------- In all cases, the multi-stage contactor model writes a set of material balances for each stream in the model. For component ``j`` in stream ``s`` the following constraint, named ``stream + "_material_balance"``, is written for all finite elements ``x``: -.. math:: 0 = \sum_p{F_{t,x-,s,p,j}} - \sum_p{F_{t,x,s,p,j}} + \left[ \sum_p{F_{side,t,x,s,p,j}} \right] + \sum_o{M_{t,x,s,o,j}} + \left[ \sum_p{G_{rate,t,x,s,p,j}} + \sum_p{G_{equil,t,x,s,p,j}} + \sum_p{G_{inher,t,x,s,p,j}} \right] +.. math:: 0 = \sum_p{F_{t,x-,s,p,j}} - \sum_p{F_{t,x,s,p,j}} + \left[ \sum_p{F_{side,t,x,s,p,j}} \right] + \sum_o{M_{t,x,s,o,j}} + \left[ \sum_p{G_{rate,t,x,s,p,j}} + \sum_p{G_{equil,t,x,s,p,j}} + \sum_p{G_{inher,t,x,s,p,j}} + \sum_p{G_{hetero,t,x,s,p,j}} \right] -where ``F`` is the material flow term, ``x-`` represents the the previous finite element (``x-1`` in the case of co-current flow and ``x+1`` in the case of counter-current flow), ``F_side`` is the material flow term for a side stream (if present) and ``o`` represents all other streams in the model (for cases where ``s`` is the second index (i.e., M_{t,x,o,s,j}) the term is multiplied by -1). The reaction generation terms are only included if the appropriate reaction type is supported by the reaction or property package for the stream. +where ``F`` is the material flow term, ``x-`` represents the previous finite element (``x-1`` in the case of co-current flow and ``x+1`` in the case of counter-current flow), ``F_side`` is the material flow term for a side stream (if present) and ``o`` represents all other streams in the model (for cases where ``s`` is the second index (i.e., M_{t,x,o,s,j}) the term is multiplied by -1). The reaction generation terms are only included if the appropriate reaction type is supported by the reaction or property package for the stream. For systems including rate reactions, the following constraint, names ``stream + "_rate_reaction_constraint"``, is written to relate the generation of component ``j`` in phase ``p`` to the extent of each rate reaction as shown below where :math:`\alpha_{r,p,j}` is the stoichiometric coefficient for component ``j`` in phase ``p`` for reaction ``r``. .. math:: G_{rate,t,x,s,p,j} = \sum_r{\alpha_{rate_r,p,j} \times X_{rate,t,x,s,r}} -Equivalent constraints are written for equilibrium and inherent reactions as necessary. +Equivalent constraints are written for equilibrium, inherent and heterogeneous reactions as necessary. For streams including energy balances (``has_energy_balance = True``) the following constraint (named ``stream + "_energy_balance"``) is written at each finite element: .. math:: 0 = \sum_p{H_{t,x-,s,p}} - \sum_p{H_{t,x,s,p}} + \biggl[ \sum_p{H_{side,t,x,s,p}} \biggr] + \sum_o{E_{t,x,s,o}} + \biggl[ Q_{t,x,s} \biggr] + \biggl[ \sum_{rate}{\Delta H_{rxn,r} \times X_{rate,t,x,s,r}} + \sum_{equil}{\Delta H_{rxn,r} \times X_{equil,t,x,s,r}} + \sum_{inher}{\Delta H_{rxn,r} \times X_{inher,t,x,s,r}} \biggr] -where ``H`` represent enthalpy flow terms and :math:`\Delta H_{rxn}` represents heat of reaction. The heat of reaction terms are only included if a reaction package is provided for the stream AND the configuration option ``has_heat_of_reaction = True`` is set for the stream. +where ``H`` represent enthalpy flow terms and :math:`\Delta H_{rxn}` represents heat of reaction. The heat of reaction terms are only included if a reaction package is provided for the stream AND the configuration option ``has_heat_of_reaction = True`` is set for the stream. **Note** heterogeneous reactions **do not** support heat of reaction terms as it is uncertain which stream/phase the heat should be added too. For streams including pressure balances (``has_pressure_balance = True``) the following constraint (named ``stream + "_pressure_balance"``) is written at each finite element: @@ -116,21 +121,150 @@ MSContactor Class Stream Configuration Options ---------------------------- -========================== ========================= ======== ============================================================================================== -Argument Type Default Description -========================== ========================= ======== ============================================================================================== -property_package PropertyParameter Block None Property package associated with stream -property_package_args dict None Configuration arguments for State Blocks -reaction_package Reaction Parameter Block None Reaction package associated with stream -reaction_package_args dict None Configuration arguments for Reaction Blocks -flow_direction FlowDirection Enum forward Direction of flow for stream -has_feed bool True Whether stream has a feed Port and inlet state, or if all flow is provided via mass transfer. -has_rate_reactions bool False Whether rate-based reactions occur in stream. -has_equilibrium_reactions bool False Whether equilibrium-based reactions occur in stream. -has_energy_balance bool True Whether to include energy balance for stream. -has_heat_transfer bool False Whether to include external heat transfer terms in energy balance for stream. -has_heat_of_reaction bool False Whether heat of reaction terms should be included in energy balance for stream. -has_pressure_balance bool True Whether to include pressure balance for stream. -has_pressure_change bool False Whether to include :math:`\Delta P` terms in pressure balance for stream. -side_streams list None Finite elements at which a side stream should be included. -========================== ========================= ======== ============================================================================================== +============================= ============================= ======== ============================================================================================== +Argument Type Default Description +============================= ============================= ======== ============================================================================================== +property_package PropertyParameter Block None Property package associated with stream +property_package_args dict None Configuration arguments for State Blocks +reaction_package Reaction Parameter Block None Reaction package associated with stream +reaction_package_args dict None Configuration arguments for Reaction Blocks +flow_direction FlowDirection Enum forward Direction of flow for stream +has_feed bool True Whether stream has a feed Port and inlet state, or if all flow is provided via mass transfer. +has_rate_reactions bool False Whether rate-based reactions occur in stream. +has_equilibrium_reactions bool False Whether equilibrium-based reactions occur in stream. +has_energy_balance bool True Whether to include energy balance for stream. +has_heat_transfer bool False Whether to include external heat transfer terms in energy balance for stream. +has_heat_of_reaction bool False Whether heat of reaction terms should be included in energy balance for stream. +has_pressure_balance bool True Whether to include pressure balance for stream. +has_pressure_change bool False Whether to include :math:`\Delta P` terms in pressure balance for stream. +side_streams list None Finite elements at which a side stream should be included. +heterogeneous_reactions heterogeneous reaction block None Heterogeneous reaction package to use for system +heterogeneous_reactions_args dict None Configuration arguments for heterogeneous reaction blocks +============================= ============================= ======== ============================================================================================== + +Heterogeneous Reaction Blocks +----------------------------- + +Heterogeneous reaction blocks are a new feature in IDAES, and are currently still in beta development. Due to this, there is no base class for heterogeneous reaction blocks yet as the API is not yet finalized. + +Currently, the requirements for are similar to those for ReactionBlocks and are demonstrated in the example outline below: + +.. code-block:: + + from pyomo.environ import Constraint, Set, units, Var + from pyomo.common.config import ConfigValue + + from idaes.core import ( + declare_process_block_class, + ProcessBlockData, + ProcessBlock, + ) + from idaes.core.base import property_meta + from idaes.core.util.misc import add_object_reference + + # ----------------------------------------------------------------------------- + # Heterogeneous Reaction Parameter Block + @declare_process_block_class("MyHeterogeneousReactionParameters") + class MyHeterogeneousReactionParametersData(ProcessBlockData, property_meta.HasPropertyClassMetadata): + def build(self): + super().build() + + self._reaction_block_class = MyHeterogeneousReactionsBlock + + self.reaction_idx = Set( + initialize=[ + "reaction1", + "reaction2", + ... + ] + ) + + self.reaction_stoichiometry = { + ("reaction1", "phase1", "component1"): stoichiometric_coefficient, + ("reaction1", "phase1", "component2"): stoichiometric_coefficient, + ... + } + + # Add any global parameters here, such as activation energies and Arrhenius constants + + # Define base units for reactions + @classmethod + def define_metadata(cls, obj): + obj.add_default_units( + { + "time": units.hour, + "length": units.m, + "mass": units.kg, + "amount": units.mol, + "temperature": units.K, + } + ) + + # The next few lines are boilerplate - you should be able to just copy these + @property + def reaction_block_class(self): + return self._reaction_block_class + + def build_reaction_block(self, *args, **kwargs): + """ + Methods to construct a ReactionBlock associated with this + ReactionParameterBlock. This will automatically set the parameters + construction argument for the ReactionBlock. + + Returns: + ReactionBlock + + """ + default = kwargs.pop("default", {}) + initialize = kwargs.pop("initialize", {}) + + if initialize == {}: + default["parameters"] = self + else: + for i in initialize.keys(): + initialize[i]["parameters"] = self + + return self.reaction_block_class( # pylint: disable=not-callable + *args, **kwargs, **default, initialize=initialize + ) + + + # Define the heterogenous ReactionBlock + @declare_process_block_class( + "MyHeterogeneousReactionsBlock", block_class=ProcessBlock + ) + class MyHeterogeneousReactionsData(ProcessBlockData): + # Create Class ConfigBlock - this needs to be here + CONFIG = ProcessBlockData.CONFIG() + CONFIG.declare( + "parameters", + ConfigValue( + description="""A reference to an instance of the Heterogeneous Reaction Parameter + Block associated with this property package.""", + ), + ) + # Add any additional configuration arguments you want + + def build(self): + super().build() + + # This creates an easy link back to the parameters + add_object_reference(self, "_params", self.config.parameters) + + # Need to define a reaction rate variable - make sure units are correct + self.reaction_rate = Var( + self.params.reaction_idx, + initialize=0, + units=#units + ) + + # Define rule for calculating reaction rate + def rule_reaction_rate_eq(b, r): + return b.reaction_rate[r] == #rate expression + + self.reaction_rate_eq = Constraint(self.params.reaction_idx, rule=rule_reaction_rate_eq) + + # This is boilerplate + @property + def params(self): + return self._params diff --git a/idaes/models/unit_models/mscontactor.py b/idaes/models/unit_models/mscontactor.py index ee11414945..611928b912 100644 --- a/idaes/models/unit_models/mscontactor.py +++ b/idaes/models/unit_models/mscontactor.py @@ -43,7 +43,6 @@ __author__ = "Andrew Lee" -# TODO: Initializer object # TODO: Could look at using Pyomo DAE for the length domain, but this would make # it harder to do side feeds. @@ -114,13 +113,28 @@ def initialization_routine( # First, build list of names for known constraints const_names = [] for s in model.config.streams.keys(): - const_names.append(s + "_rate_reaction_constraint") - const_names.append(s + "_equilibrium_reaction_constraint") - const_names.append(s + "_inherent_reaction_constraint") - const_names.append(s + "_material_balance") - const_names.append(s + "_energy_balance") - const_names.append(s + "_pressure_balance") - const_names.append(s + "_side_stream_pressure_balance") + const_names.append(model.name + "." + s + "_rate_reaction_constraint") + const_names.append( + model.name + "." + s + "_equilibrium_reaction_constraint" + ) + const_names.append(model.name + "." + s + "_inherent_reaction_constraint") + const_names.append( + model.name + "." + s + "_heterogeneous_reaction_constraint" + ) + const_names.append(model.name + "." + s + "_material_balance") + const_names.append(model.name + "." + s + "_energy_balance") + const_names.append(model.name + "." + s + "_pressure_balance") + const_names.append(model.name + "." + s + "_side_stream_pressure_balance") + + try: + # If has rate reactions ,fi extent to 0 for first pass + getattr(model, s + "_rate_reaction_extent").fix(0) + except AttributeError: + pass + + # Fix extents for heterogeneous reactions to 0 for first pass if present + if hasattr(model, "heterogeneous_reaction_extent"): + model.heterogeneous_reaction_extent.fix(0) # Iterate through all constraints attached to model - do not search sub-blocks for c in model.component_objects(Constraint, descend_into=False): @@ -322,6 +336,26 @@ class MSContactorData(UnitModelBlockData): doc="List of interacting stream pairs as 2-tuples ('stream1', 'stream2').", ), ) + # TODO: Consider a base call for heterogeneous reactions and set domain + CONFIG.declare( + "heterogeneous_reactions", + ConfigValue( + default=None, + # domain=list, + description="Heterogeneous reaction package to use in contactor.", + doc="Heterogeneous reaction package to use in contactor. Heterogeneous " + "reaction packages are expected to have a certain structure and methods; " + "please refer to the documentation for more details.", + ), + ) + CONFIG.declare( + "heterogeneous_reactions_args", + ConfigDict( + implicit=True, + description="Arguments to use for constructing reaction packages", + doc="ConfigBlock with arguments to be passed to heterogeneous reaction block(s)", + ), + ) def build(self): """ @@ -338,6 +372,10 @@ def build(self): self._verify_inputs() flow_basis, uom = self._build_state_blocks() + + if self.config.heterogeneous_reactions is not None: + self._build_heterogeneous_reaction_blocks() + self._build_material_balance_constraints(flow_basis, uom) self._build_energy_balance_constraints(uom) self._build_pressure_balance_constraints(uom) @@ -386,10 +424,14 @@ def _verify_inputs(self): ): # Common component, assume interaction self.stream_component_interactions.add((stream1, stream2, j)) - if len(self.stream_component_interactions) == 0: + if ( + len(self.stream_component_interactions) == 0 + and self.config.heterogeneous_reactions is None + ): raise ConfigurationError( - "No common components found in property packages. MSContactor model assumes " - "mass transfer occurs between components with the same name in different streams." + "No common components found in property packages and no heterogeneous reactions " + "specified. The MSContactor model assumes that mass transfer occurs between " + "components with the same name in different streams or due to heterogeneous reactions." ) # Check that reaction block was provided if reactions requested @@ -487,6 +529,30 @@ def _build_state_blocks(self): return flow_basis, uom + def _build_heterogeneous_reaction_blocks(self): + rpack = self.config.heterogeneous_reactions + rpack_args = self.config.heterogeneous_reactions_args + + try: + self.heterogeneous_reactions = rpack.build_reaction_block( + self.flowsheet().time, + self.elements, + doc="Heterogeneous reaction block for contactor.", + **rpack_args, + ) + except AttributeError: + raise ConfigurationError( + "Heterogeneous reaction package has not implemented a " + "build_reaction_block method. Please ensure that your " + "reaction block conforms to the required standards." + ) + + if not hasattr(self.config.heterogeneous_reactions, "reaction_idx"): + raise PropertyNotSupportedError( + "Heterogeneous reaction package does not contain a list of " + "reactions (reaction_idx)." + ) + def _build_material_balance_constraints(self, flow_basis, uom): # Get units for transfer terms if flow_basis is MaterialFlowBasis.molar: @@ -509,6 +575,19 @@ def _build_material_balance_constraints(self, flow_basis, uom): doc="Inter-stream mass transfer term", ) + if hasattr(self, "heterogeneous_reactions"): + # Add extents of reaction and stoichiometric constraints + # We will assume the user will define how extent will be calculated + self.heterogeneous_reaction_extent = Var( + self.flowsheet().time, + self.elements, + self.config.heterogeneous_reactions.reaction_idx, + domain=Reals, + initialize=0.0, + doc="Extent of heterogeneous reactions", + units=mb_units, + ) + # Build balance equations for stream, sconfig in self.config.streams.items(): state_block = getattr(self, stream) @@ -520,7 +599,7 @@ def _build_material_balance_constraints(self, flow_basis, uom): if hasattr(self, stream + "_reactions"): reaction_block = getattr(self, stream + "_reactions") - # Add equilibrium reaction terms (if required) + # Add homogeneous rate reaction terms (if required) if sconfig.has_rate_reactions: if not hasattr(sconfig.reaction_package, "rate_reaction_idx"): raise PropertyNotSupportedError( @@ -701,6 +780,51 @@ def inherent_reaction_rule(b, t, s, p, j): inherent_reaction_constraint, ) + # Add heterogeneous reaction terms (if required) + if hasattr(self, "heterogeneous_reactions"): + heterogeneous_reactions_generation = Var( + self.flowsheet().time, + self.elements, + pc_set, + domain=Reals, + initialize=0.0, + doc="Generation due to heterogeneous reactions", + units=mb_units, + ) + self.add_component( + stream + "_heterogeneous_reactions_generation", + heterogeneous_reactions_generation, + ) + + def heterogeneous_reaction_rule(b, t, s, p, j): + if (p, j) in pc_set: + return heterogeneous_reactions_generation[t, s, p, j] == ( + sum( + self.heterogeneous_reactions[ + t, s + ].params.reaction_stoichiometry[r, p, j] + * b.heterogeneous_reaction_extent[t, s, r] + for r in self.config.heterogeneous_reactions.reaction_idx + if (r, p, j) + in self.heterogeneous_reactions[ + t, s + ].params.reaction_stoichiometry + ) + ) + return Constraint.Skip + + heterogeneous_reaction_constraint = Constraint( + self.flowsheet().time, + self.elements, + pc_set, + doc="Heterogeneous reaction stoichiometry constraint", + rule=heterogeneous_reaction_rule, + ) + self.add_component( + stream + "_heterogeneous_reaction_constraint", + heterogeneous_reaction_constraint, + ) + # Material balance for stream def material_balance_rule(b, t, s, j): in_state, out_state, side_state = _get_state_blocks(b, t, s, stream) @@ -760,6 +884,13 @@ def material_balance_rule(b, t, s, j): inherent_reaction_generation[t, s, p, j] for p in phase_list ) + # Add heterogeneous reactions (if required) + if self.config.heterogeneous_reactions is not None: + rhs += sum( + heterogeneous_reactions_generation[t, s, p, j] + for p in phase_list + ) + return 0 == rhs mbal = Constraint( diff --git a/idaes/models/unit_models/tests/test_mscontactor.py b/idaes/models/unit_models/tests/test_mscontactor.py index b28b087616..66fcebca86 100644 --- a/idaes/models/unit_models/tests/test_mscontactor.py +++ b/idaes/models/unit_models/tests/test_mscontactor.py @@ -16,9 +16,11 @@ """ import pytest +from types import MethodType from pyomo.environ import ( assert_optimal_termination, + Block, ConcreteModel, Constraint, log, @@ -53,8 +55,9 @@ MSContactorInitializer, ) from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.util.misc import add_object_reference from idaes.core.solvers import get_solver -from idaes.core.util.exceptions import ConfigurationError +from idaes.core.util.exceptions import ConfigurationError, PropertyNotSupportedError from idaes.core.util.initialization import ( propagate_state, fix_state_vars, @@ -389,12 +392,16 @@ def test_verify_inputs_no_common_components(self): with pytest.raises( ConfigurationError, - match="No common components found in property packages. MSContactor " - "model assumes mass transfer occurs between components with the " - "same name in different streams.", + match="No common components found in property packages and no heterogeneous reactions " + "specified. The MSContactor model assumes that mass transfer occurs between " + "components with the same name in different streams or due to heterogeneous reactions.", ): m.fs.unit._verify_inputs() + # Should pass if heterogeneous reaction argument provided + m.fs.unit.config.heterogeneous_reactions = True + m.fs.unit._verify_inputs() + @pytest.mark.unit def test_verify_inputs_reactions_with_no_package(self): m = ConcreteModel() @@ -1834,6 +1841,203 @@ def test_equilibrium_reactions(self, model): ) ) + @pytest.mark.unit + def test_heterogeneous_reactions_no_build_method(self, model): + model.fs.unit.config.heterogeneous_reactions = True + + model.fs.unit._verify_inputs() + _, _ = model.fs.unit._build_state_blocks() + + with pytest.raises( + ConfigurationError, + match="Heterogeneous reaction package has not implemented a " + "build_reaction_block method. Please ensure that your " + "reaction block conforms to the required standards.", + ): + model.fs.unit._build_heterogeneous_reaction_blocks() + + @pytest.mark.unit + def test_heterogeneous_reactions_no_rxn_index(self, model): + model.hetero_dummy = Block() + + def build_reaction_block(*args, **kwargs): + pass + + model.hetero_dummy.build_reaction_block = MethodType( + build_reaction_block, model.hetero_dummy + ) + + model.fs.unit.config.heterogeneous_reactions = model.hetero_dummy + + model.fs.unit._verify_inputs() + _, _ = model.fs.unit._build_state_blocks() + + with pytest.raises( + PropertyNotSupportedError, + match="Heterogeneous reaction package does not contain a list of " + "reactions \(reaction_idx\).", + ): + model.fs.unit._build_heterogeneous_reaction_blocks() + + @pytest.mark.unit + def test_heterogeneous_reactions(self, model): + model.fs.hetero_dummy = Block() + model.fs.hetero_dummy.reaction_idx = Set(initialize=["R1", "R2", "R3", "R4"]) + model.fs.hetero_dummy.reaction_stoichiometry = { + ("R1", "p1", "c1"): 1, + ("R2", "p1", "c2"): 1, + ("R3", "p2", "c1"): 1, + ("R4", "p2", "c2"): 1, + } + + model.fs.unit.config.heterogeneous_reactions = model.fs.hetero_dummy + + model.fs.unit._verify_inputs() + flow_basis, uom = model.fs.unit._build_state_blocks() + + model.fs.unit.heterogeneous_reactions = Block( + model.fs.time, + model.fs.unit.elements, + ) + for e in model.fs.unit.elements: + add_object_reference( + model.fs.unit.heterogeneous_reactions[0, e], + "params", + model.fs.hetero_dummy, + ) + + model.fs.unit._build_material_balance_constraints(flow_basis, uom) + + assert isinstance(model.fs.unit.heterogeneous_reaction_extent, Var) + for k in model.fs.unit.heterogeneous_reaction_extent.keys(): + assert k in [ + (0, 1, "R1"), + (0, 1, "R2"), + (0, 1, "R3"), + (0, 1, "R4"), + (0, 2, "R1"), + (0, 2, "R2"), + (0, 2, "R3"), + (0, 2, "R4"), + ] + + for s in ["stream1", "stream2"]: + gen = getattr(model.fs.unit, s + "_heterogeneous_reactions_generation") + assert isinstance(gen, Var) + for k in gen: + assert k in [ + (0, 1, "p1", "c1"), + (0, 1, "p1", "c2"), + (0, 1, "p2", "c1"), + (0, 1, "p2", "c2"), + (0, 2, "p1", "c1"), + (0, 2, "p1", "c2"), + (0, 2, "p2", "c1"), + (0, 2, "p2", "c2"), + ] + + con = getattr(model.fs.unit, s + "_heterogeneous_reaction_constraint") + assert isinstance(con, Constraint) + for k, c in con.items(): + assert k in [ + (0, 1, "p1", "c1"), + (0, 1, "p1", "c2"), + (0, 1, "p2", "c1"), + (0, 1, "p2", "c2"), + (0, 2, "p1", "c1"), + (0, 2, "p1", "c2"), + (0, 2, "p2", "c1"), + (0, 2, "p2", "c2"), + ] + + if k[2] == "p1" and k[3] == "c1": + r = "R1" + elif k[2] == "p1" and k[3] == "c2": + r = "R2" + elif k[2] == "p2" and k[3] == "c1": + r = "R3" + else: + r = "R4" + + expr = str( + gen[k] - model.fs.unit.heterogeneous_reaction_extent[0, k[1], r] + ) + assert str(c.body) == expr + + for j in [ + "c1", + "c2", + ]: # has +ve mass transfer, forward flow, heterogeneous reactions + assert str(model.fs.unit.stream1_material_balance[0, 1, j]._expr) == str( + 0 + == sum( + model.fs.unit.stream1_inlet_state[0].get_material_flow_terms(p, j) + for p in ["p1", "p2"] + ) + - sum( + model.fs.unit.stream1[0, 1].get_material_flow_terms(p, j) + for p in ["p1", "p2"] + ) + + model.fs.unit.material_transfer_term[0, 1, "stream1", "stream2", j] + + sum( + model.fs.unit.stream1_heterogeneous_reactions_generation[0, 1, p, j] + for p in ["p1", "p2"] + ) + ) + assert str(model.fs.unit.stream1_material_balance[0, 2, j]._expr) == str( + 0 + == sum( + model.fs.unit.stream1[0, 1].get_material_flow_terms(p, j) + for p in ["p1", "p2"] + ) + - sum( + model.fs.unit.stream1[0, 2].get_material_flow_terms(p, j) + for p in ["p1", "p2"] + ) + + model.fs.unit.material_transfer_term[0, 2, "stream1", "stream2", j] + + sum( + model.fs.unit.stream1_heterogeneous_reactions_generation[0, 2, p, j] + for p in ["p1", "p2"] + ) + ) + + for j in [ + "c1", + "c2", + ]: # has -ve mass transfer, forward flow, heterogeneous reactions + assert str(model.fs.unit.stream2_material_balance[0, 2, j]._expr) == str( + 0 + == sum( + model.fs.unit.stream2_inlet_state[0].get_material_flow_terms(p, j) + for p in ["p1", "p2"] + ) + - sum( + model.fs.unit.stream2[0, 2].get_material_flow_terms(p, j) + for p in ["p1", "p2"] + ) + - model.fs.unit.material_transfer_term[0, 2, "stream1", "stream2", j] + + sum( + model.fs.unit.stream2_heterogeneous_reactions_generation[0, 2, p, j] + for p in ["p1", "p2"] + ) + ) + assert str(model.fs.unit.stream2_material_balance[0, 1, j]._expr) == str( + 0 + == sum( + model.fs.unit.stream2[0, 2].get_material_flow_terms(p, j) + for p in ["p1", "p2"] + ) + - sum( + model.fs.unit.stream2[0, 1].get_material_flow_terms(p, j) + for p in ["p1", "p2"] + ) + - model.fs.unit.material_transfer_term[0, 1, "stream1", "stream2", j] + + sum( + model.fs.unit.stream2_heterogeneous_reactions_generation[0, 1, p, j] + for p in ["p1", "p2"] + ) + ) + @pytest.mark.unit def test_rate_reactions(self, model): model.fs.unit.config.streams["stream2"].reaction_package = model.fs.reactions