Merge branch 'main' into pysmo_custom_sampling

docs/reference_guides/model_libraries/generic/unit_models/mscontactor.rst

@@ -75,14 +75,21 @@ Variable                      Name                                           Des
 :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
+:math:`V_x`                   volume                                         Total volume of element ``x``                                                                                 Only if ``has_holdup``
+:math:`f_{t,x,s}`             volume_frac_stream                             Volume fraction of stream ``s`` in element ``x`` at time ``t``                                                Only if ``has_holdup``
+:math:`\phi_{t,x,s,p}`        stream + "_phase_fraction"                     Volume fraction of phase ``p`` in stream ``s`` in element ``x`` at time ``t``                                 Only if ``has_holdup``
+:math:`N_{t,x,s,p,j}`         stream + "_material_holdup"                    Holdup of component ``j`` in phase ``p`` for stream ``s`` at ``x`` and ``t``                                  Only if ``has_holdup``
+:math:`dN/dt_{t,x,s,p,j}`     stream + "_material_accumulation"              Accumulation of component ``j`` in phase ``p`` for stream ``s`` at ``x`` and ``t``                            Only if ``dynamic``
+:math:`U_{t,x,s,p}`           stream + "_energy_holdup"                      Holdup of energy in phase ``p`` for stream ``s`` at ``x`` and ``t``                                           Only if ``has_holdup``
+:math:`dU/dt_{t,x,s,p}`       stream + "_energy_accumulation"                Accumulation of energy in phase ``p`` for stream ``s`` at ``x`` and ``t``                                     Only if ``dynamic``
 ============================= ============================================== ============================================================================================================= =========================================
 
 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}} + \sum_p{G_{hetero,t,x,s,p,j}} \right]
+.. math:: dN/dt_{t,x,s,p,j} = \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 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.
 
@@ -94,18 +101,38 @@ Equivalent constraints are written for equilibrium, inherent and heterogeneous r
 
 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]
+.. math:: \sum_p{dU/dt_{t,x,s,p}} = \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. **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:
 
 .. math:: 0 = P_{t,x-,s} - P_{t,x,s} + \biggl[ \Delta P_{t,x,s} \biggr]
 
-where ``P`` represents pressure. For streams with side streams, the following pressure equality constraint (names ``stream + "_side_stream_pressure_balance"``) is also written:
+where ``P`` represents pressure. For streams with side streams, the following pressure equality constraint (named ``stream + "_side_stream_pressure_balance"``) is also written:
 
 .. math:: P_{t,x,s} = P_{side,t,x,s}
 
+If ``has_holdup`` is true, the following additional constraints are included to calculate holdup terms. First, ``sum_volume_frac`` constrains the sum of all volume fractions to be 1.
+
+.. math:: 1 = \sum_s{f_{t,x,s}}
+
+Additionally, constraints are written for the sum of phase fractions in each stream (named ``stream + "_sum_phase_fractions"``):
+
+.. math:: 1 = \sum_p{\phi_{t,x,s,p}}
+
+The material holdup is defined by the following constraint (named ``stream + "_material_holdup_constraint"``):
+
+.. math:: N_{t,x,s,p,j} = V \times f_{t,x,s} \phi_{t,x,s,p} \times \C_{t,x,s,p,j}
+
+where :math:`C_{t,x,s,p,j}` is the concentration of component ``j`` in phase ``p`` for stream ``s`` at ``x`` and ``t``.
+
+The energy holdup is defined by the following constraint (named ``stream + "_energy_holdup_constraint"``):
+
+.. math:: U_{t,x,s,p} = V*f_{t,x,s} \times \phi_{t,x,s,p} \times \u_{t,x,s,p}
+
+where :math:`u_{t,x,s,p}` is the internal energy density of phase ``p`` for stream ``s`` at ``x`` and ``t``.
+
 Initialization
 --------------
 

docs/reference_guides/model_libraries/generic/unit_models/solid_liquid/thickener0d.rst

@@ -1,6 +1,9 @@
 Thickener (0D)
 ==============
 
+.. warning::
+    The Thickener model is currently in beta status and will likely change in the next release as a more predictive version is developed.
+
 The ``Thickener0D`` unit model is an extension of the :ref:`SLSeparator <reference_guides/model_libraries/generic/unit_models/solid_liquid/sl_separator:Generic Solid-Liquid Separator>` model which adds constraints to estimate the area and height of a vessel required to achieve the desired separation of solid and liquid based on experimental measurements of the settling velocity. This model is based on correlations described in:
 
 [1] Coulson & Richardson's Chemical Engineering, Volume 2 Particle Technology & Separation Processes (4th Ed.), Butterworth-Heinemann (2001)

idaes/core/dmf/datasets.py

@@ -39,7 +39,7 @@
 NAME = "idaes-pse"
 
 # Keep this the same as constant DMF_DATA_ROOT in setup.py
-DMF_DATA_ROOT = "data"
+DMF_DATA_ROOT = "dmf_data"
 
 _log = logging.getLogger(__name__)
 

idaes/models/unit_models/mscontactor.py

@@ -18,9 +18,19 @@
 from pandas import DataFrame
 
 # Import Pyomo libraries
-from pyomo.environ import Block, Constraint, RangeSet, Reals, Set, units, Var
+from pyomo.environ import (
+    Block,
+    Constraint,
+    Expression,
+    RangeSet,
+    Reals,
+    Set,
+    units,
+    Var,
+)
 from pyomo.common.config import ConfigDict, ConfigValue, Bool, In
 from pyomo.contrib.incidence_analysis import solve_strongly_connected_components
+from pyomo.dae import DerivativeVar
 
 # Import IDAES cores
 from idaes.core import (
@@ -381,6 +391,8 @@ def build(self):
         if self.config.heterogeneous_reactions is not None:
             self._build_heterogeneous_reaction_blocks()
 
+        self._add_geometry(uom)
+
         self._build_material_balance_constraints(flow_basis, uom)
         self._build_energy_balance_constraints(uom)
         self._build_pressure_balance_constraints(uom)
@@ -394,17 +406,16 @@ def _verify_inputs(self):
                 f"{list(self.config.streams.keys())}"
             )
 
-        if self.config.dynamic:
-            raise NotImplementedError(
-                "MSContactor model does not support dynamics yet."
-            )
-
         # Build indexing sets
         self.elements = RangeSet(
             1,
             self.config.number_of_finite_elements,
             doc="Set of finite elements in cascade (1 to number of elements)",
         )
+        self.streams = Set(
+            initialize=[k for k in self.config.streams.keys()],
+            doc="Set of streams in unit",
+        )
 
         interacting_streams = self.config.interacting_streams
         # If user did not provide interacting streams list, assume all steams interact
@@ -558,15 +569,49 @@ def _build_heterogeneous_reaction_blocks(self):
                 "reactions (reaction_idx)."
             )
 
+    def _add_geometry(self, uom):
+        if self.config.has_holdup:
+            # Add volume for each element
+            # TODO: Assuming constant volume for now
+            self.volume = Var(
+                self.elements,
+                initialize=1,
+                units=uom.VOLUME,
+                doc="Volume of element",
+            )
+            self.volume_frac_stream = Var(
+                self.flowsheet().time,
+                self.elements,
+                self.streams,
+                initialize=1 / len(self.streams),
+                units=units.dimensionless,
+                doc="Volume fraction of each stream in element",
+            )
+
+            @self.Constraint(
+                self.flowsheet().time,
+                self.elements,
+                doc="Sum of volume fractions constraint",
+            )
+            def sum_volume_frac(b, t, e):
+                return 1 == sum(b.volume_frac_stream[t, e, s] for s in b.streams)
+
+            for stream in self.config.streams.keys():
+                phase_list = getattr(self, stream).phase_list
+                _add_phase_fractions(self, stream, phase_list)
+
     def _build_material_balance_constraints(self, flow_basis, uom):
         # Get units for transfer terms
         if flow_basis is MaterialFlowBasis.molar:
             mb_units = uom.FLOW_MOLE
+            hu_units = uom.AMOUNT
         elif flow_basis is MaterialFlowBasis.mass:
             mb_units = uom.FLOW_MASS
+            hu_units = uom.MASS
         else:
             # Flow type other, so cannot determine units
             mb_units = None
+            hu_units = None
 
         # Material transfer terms are indexed by stream pairs and components.
         # Convention is that a positive material flow term indicates flow into
@@ -603,6 +648,49 @@ def _build_material_balance_constraints(self, flow_basis, uom):
             if hasattr(self, stream + "_reactions"):
                 reaction_block = getattr(self, stream + "_reactions")
 
+            # Material holdup and accumulation
+            if self.config.has_holdup:
+                material_holdup = Var(
+                    self.flowsheet().time,
+                    self.elements,
+                    pc_set,
+                    domain=Reals,
+                    initialize=1.0,
+                    doc="Material holdup of stream in element",
+                    units=hu_units,
+                )
+                self.add_component(
+                    stream + "_material_holdup",
+                    material_holdup,
+                )
+
+                holdup_eq = Constraint(
+                    self.flowsheet().time,
+                    self.elements,
+                    pc_set,
+                    doc=f"Holdup constraint for stream {stream}",
+                    rule=partial(
+                        _holdup_rule,
+                        stream=stream,
+                    ),
+                )
+                self.add_component(
+                    stream + "_material_holdup_constraint",
+                    holdup_eq,
+                )
+
+            if self.config.dynamic:
+                material_accumulation = DerivativeVar(
+                    material_holdup,
+                    wrt=self.flowsheet().time,
+                    doc="Material accumulation for in element",
+                    units=mb_units,
+                )
+                self.add_component(
+                    stream + "_material_accumulation",
+                    material_accumulation,
+                )
+
             # Add homogeneous rate reaction terms (if required)
             if sconfig.has_rate_reactions:
                 if not hasattr(sconfig.reaction_package, "rate_reaction_idx"):
@@ -812,6 +900,52 @@ def _build_energy_balance_constraints(self, uom):
 
         for stream, pconfig in self.config.streams.items():
             if pconfig.has_energy_balance:
+                state_block = getattr(self, stream)
+                phase_list = state_block.phase_list
+
+                # Material holdup and accumulation
+                if self.config.has_holdup:
+                    energy_holdup = Var(
+                        self.flowsheet().time,
+                        self.elements,
+                        phase_list,
+                        domain=Reals,
+                        initialize=1.0,
+                        doc="Energy holdup of stream in element",
+                        units=uom.ENERGY,
+                    )
+                    self.add_component(
+                        stream + "_energy_holdup",
+                        energy_holdup,
+                    )
+
+                    energy_holdup_eq = Constraint(
+                        self.flowsheet().time,
+                        self.elements,
+                        phase_list,
+                        doc=f"Energy holdup constraint for stream {stream}",
+                        rule=partial(
+                            _energy_holdup_rule,
+                            stream=stream,
+                        ),
+                    )
+                    self.add_component(
+                        stream + "_energy_holdup_constraint",
+                        energy_holdup_eq,
+                    )
+
+                if self.config.dynamic:
+                    energy_accumulation = DerivativeVar(
+                        energy_holdup,
+                        wrt=self.flowsheet().time,
+                        doc="Energy accumulation for in element",
+                        units=uom.POWER,
+                    )
+                    self.add_component(
+                        stream + "_energy_accumulation",
+                        energy_accumulation,
+                    )
+
                 if pconfig.has_heat_transfer:
                     heat = Var(
                         self.flowsheet().time,
@@ -1136,7 +1270,15 @@ def _material_balance_rule(blk, t, s, j, stream, mb_units):
             for p in phase_list
         )
 
-    return 0 == rhs
+    if not blk.config.dynamic:
+        lhs = 0 * mb_units
+    else:
+        acc = getattr(blk, stream + "_material_accumulation")
+        lhs = sum(acc[t, s, p, j] for p in phase_list)
+        if mb_units is not None:
+            lhs = units.convert(lhs, mb_units)
+
+    return lhs == rhs
 
 
 def _get_energy_transfer_term(blk, uom):
@@ -1225,7 +1367,13 @@ def _energy_balance_rule(blk, t, s, stream, uom):
                 for e in pconfig.reaction_package.equilibrium_reaction_idx
             )
 
-    return 0 == rhs
+    if not blk.config.dynamic:
+        lhs = 0 * uom.POWER
+    else:
+        acc = getattr(blk, stream + "_energy_accumulation")
+        lhs = units.convert(sum(acc[t, s, p] for p in phase_list), uom.POWER)
+
+    return lhs == rhs
 
 
 def _pressure_balance_rule(blk, t, s, stream, uom):
@@ -1246,11 +1394,78 @@ def _pressure_balance_rule(blk, t, s, stream, uom):
     if pconfig.has_pressure_change:
         rhs += getattr(blk, stream + "_deltaP")[t, s]
 
-    return 0 == rhs
+    return 0 * uom.PRESSURE == rhs
 
 
 def _side_stream_pressure_rule(b, t, s, stream):
     stage_state = getattr(b, stream)[t, s]
     side_state = getattr(b, stream + "_side_stream_state")[t, s]
 
     return stage_state.pressure == side_state.pressure
+
+
+def _holdup_rule(b, t, e, p, j, stream):
+    holdup = getattr(b, stream + "_material_holdup")
+    stage_state = getattr(b, stream)[t, e]
+    phase_frac = getattr(b, stream + "_phase_fraction")
+
+    return holdup[t, e, p, j] == (
+        b.volume[e]
+        * b.volume_frac_stream[t, e, stream]
+        * phase_frac[t, e, p]
+        * stage_state.get_material_density_terms(p, j)
+    )
+
+
+def _energy_holdup_rule(b, t, e, p, stream):
+    holdup = getattr(b, stream + "_energy_holdup")
+    stage_state = getattr(b, stream)[t, e]
+    phase_frac = getattr(b, stream + "_phase_fraction")
+
+    return holdup[t, e, p] == (
+        b.volume[e]
+        * b.volume_frac_stream[t, e, stream]
+        * phase_frac[t, e, p]
+        * stage_state.get_energy_density_terms(p)
+    )
+
+
+def _add_phase_fractions(b, stream, phase_list):
+    if len(phase_list) > 1:
+        phase_fraction = Var(
+            b.flowsheet().time,
+            b.elements,
+            phase_list,
+            initialize=1 / len(phase_list),
+            doc=f"Volume fraction of holdup by phase in stream {stream}",
+        )
+        b.add_component(stream + "_phase_fraction", phase_fraction)
+
+        sum_of_phase_fractions = Constraint(
+            b.flowsheet().time,
+            b.elements,
+            rule=partial(
+                _sum_phase_frac_rule, phase_frac=phase_fraction, phase_list=phase_list
+            ),
+            doc=f"Sum of phase fractions constraint for stream {stream}",
+        )
+        b.add_component(stream + "_sum_phase_fractions", sum_of_phase_fractions)
+
+    else:
+
+        def phase_frac_rule(b, t, x, p):
+            return 1
+
+        phase_fraction = Expression(
+            b.flowsheet().time,
+            b.elements,
+            phase_list,
+            rule=phase_frac_rule,
+            doc=f"Volume fraction of holdup by phase in stream {stream}",
+        )
+
+        b.add_component(stream + "_phase_fraction", phase_fraction)
+
+
+def _sum_phase_frac_rule(b, t, x, phase_frac, phase_list):
+    return 1 == sum(phase_frac[t, x, p] for p in phase_list)