Adding support for inherent rxns to mixer

idaes/models/unit_models/mixer.py

@@ -18,6 +18,7 @@
 from pyomo.environ import (
     Block,
     check_optimal_termination,
+    Constraint,
     Param,
     PositiveReals,
     Reals,
@@ -635,6 +636,50 @@ def add_material_mixing_equations(self, inlet_blocks, mixed_block, mb_type):
             # Let this pass for now with no units
             flow_units = None
 
+        if mixed_block.include_inherent_reactions:
+            if mb_type == MaterialBalanceType.total:
+                raise ConfigurationError(
+                    "Cannot do total flow mixing with inherent reaction; "
+                    "problem is under-constrained. Please use a different "
+                    "mixing type."
+                )
+
+            # Add extents of reaction and stoichiometric constraints
+            self.inherent_reaction_extent = Var(
+                self.flowsheet().time,
+                mixed_block.params.inherent_reaction_idx,
+                domain=Reals,
+                initialize=0.0,
+                doc=f"Extent of inherent reactions in outlet",
+                units=flow_units,
+            )
+
+            self.inherent_reaction_generation = Var(
+                self.flowsheet().time,
+                pc_set,
+                domain=Reals,
+                initialize=0.0,
+                doc=f"Generation due to inherent reactions in outlet",
+                units=flow_units,
+            )
+
+            @self.Constraint(
+                self.flowsheet().time,
+                pc_set,
+            )
+            def inherent_reaction_constraint(b, t, p, j):
+                if (p, j) in pc_set:
+                    return b.inherent_reaction_generation[t, p, j] == (
+                        sum(
+                            mixed_block[t].params.inherent_reaction_stoichiometry[
+                                r, p, j
+                            ]
+                            * self.inherent_reaction_extent[t, r]
+                            for r in mixed_block[t].params.inherent_reaction_idx
+                        )
+                    )
+                return Constraint.Skip
+
         if mb_type == MaterialBalanceType.componentPhase:
             # Create equilibrium generation term and constraints if required
             if self.config.has_phase_equilibrium is True:
@@ -678,6 +723,9 @@ def material_mixing_equations(b, t, p, j):
                         and pp.phase_equilibrium_list[r][1][1] == p
                     )
 
+                if mixed_block.include_inherent_reactions:
+                    rhs += b.inherent_reaction_generation[t, p, j]
+
                 return 0 == rhs
 
         elif mb_type == MaterialBalanceType.componentTotal:
@@ -688,7 +736,7 @@ def material_mixing_equations(b, t, p, j):
                 doc="Material mixing equations",
             )
             def material_mixing_equations(b, t, j):
-                return 0 == sum(
+                rhs = sum(
                     sum(
                         inlet_blocks[i][t].get_material_flow_terms(p, j)
                         for i in range(len(inlet_blocks))
@@ -698,11 +746,20 @@ def material_mixing_equations(b, t, j):
                     if (p, j) in pc_set
                 )
 
+                if mixed_block.include_inherent_reactions:
+                    rhs += sum(
+                        b.inherent_reaction_generation[t, p, j]
+                        for p in mixed_block.phase_list
+                        if (p, j) in pc_set
+                    )
+
+                return 0 == rhs
+
         elif mb_type == MaterialBalanceType.total:
             # Write phase-component balances
             @self.Constraint(self.flowsheet().time, doc="Material mixing equations")
             def material_mixing_equations(b, t):
-                return 0 == sum(
+                rhs = sum(
                     sum(
                         sum(
                             inlet_blocks[i][t].get_material_flow_terms(p, j)
@@ -715,6 +772,8 @@ def material_mixing_equations(b, t):
                     for p in mixed_block.phase_list
                 )
 
+                return 0 == rhs
+
         elif mb_type == MaterialBalanceType.elementTotal:
             raise ConfigurationError(
                 "{} Mixers do not support elemental "