Reformulate solubility_products constraints (#1197)

* Modified log_solubility_product function to accomodate the comparison betweem s and Q. Modified bounds of mole_frac_phase_comp to be consistent with those of log_mole_frac_phase_comp * Remove comments * run black * Update state_difinitions vars lower bouds * Add mutable parameters for reformated solubility_product * Update solubility_product tests * test_plate_heat_exchanger.py Remove initialization opts after updating bounds in state_difinitions * run black * Reverse changes for lower bounds in state_definition * Add assertion for default s_norm and a_scale checking * Typo corrected * Add documentation for solubility_product * Update doc string * Equation corrected for sulubility product doc * Switch to set_value for default s_norm * Scaling applied in subsaturated cases (low concentrations) to use default solver options * Run black --------- Co-authored-by: Xinhong Liu <[email protected]> Co-authored-by: Bernard Knueven <[email protected]> Co-authored-by: Andrew Lee <[email protected]> Co-authored-by: Ludovico Bianchi <[email protected]>
IDAES · Jul 27, 2023 · 52d7d5c · 52d7d5c
1 parent ad5b8ad
commit 52d7d5c
Show file tree

Hide file tree

Showing 4 changed files with 192 additions and 44 deletions.
diff --git a/docs/explanations/components/property_package/general_reactions/equil_form.rst b/docs/explanations/components/property_package/general_reactions/equil_form.rst
@@ -7,7 +7,7 @@ Equilibrium Reaction Forms
 Power Law (power_law_equil)
 ---------------------------
 
-The method uses a power law form using the concentration form provided  to calculate the reaction rate.
+The method uses a power law form using the concentration form provided to calculate the reaction equilibrium.
 
 .. math:: k_{eq} = \prod_{(p, j)}{x_{(p,j)}^{O_{(p,j)}}}
 
@@ -18,14 +18,14 @@ The method uses a power law form using the concentration form provided  to calcu
 
    ":math:`O`", "reaction_order", "phase, component", "Reaction order"
 
-Providing a `reaction_order` dict is optional. If one is not provided, it will be assumed that this is an elementary reaction and that the reaction order is equal to the stoichiometric coefficient for all component in non-solid phases (the contribution of solid phases is assumed to be constant and included in the equilibrium constant, thus an order of zero is assumed).
+Providing a `reaction_order` dict is optional. If one is not provided, it will be assumed that this is an elementary reaction and that the reaction order is equal to the stoichiometric coefficient for all components in non-solid phases (the contribution of solid phases is assumed to be constant and included in the equilibrium constant, thus an order of zero is assumed).
 
 Log Power Law (log_power_law_equil)
 -----------------------------------
 
-The method uses a log form of a power law using the concentration form provided  to calculate the reaction rate.
+The method uses a log form of a power law using the concentration form provided to calculate the reaction equilibrium.
 
-.. math:: log(k_{eq}) = \sum_{(p, j)}{O_{(p,j)}*log(x_{(p,j))}
+.. math:: log(k_{eq}) = \sum_{(p, j)}{O_{(p,j)}*log(x_{(p,j)})}
 
 **Parameters**
 
@@ -34,4 +34,77 @@ The method uses a log form of a power law using the concentration form provided
 
    ":math:`O`", "reaction_order", "phase, component", "Reaction order"
 
+Providing a `reaction_order` dict is optional. If one is not provided, it will be assumed that this is an elementary reaction and that the reaction order is equal to the stoichiometric coefficient for all components in non-solid phases (the contribution of solid phases is assumed to be constant and included in the equilibrium constant, thus an order of zero is assumed).
+
+Solubility Product (solubility_product)
+---------------------------------------
+
+The method uses a complementarity formulation for solid precipitation using solubility product form.
+
+For precipitation at equilibrium, 
+
+.. math:: Q = k_{sp} - \prod_{(liq, j)}{x_{(liq,j)}^{O_{(liq,j)}}}
+
+* :math:`Q = 0` only holds if solids (S) are present in the system, :math:`S \geq 0`.
+* :math:`Q \geq 0` if the system is subsaturated and no solids (S) are present, :math:`S = 0`.
+
+:math:`S` represents the solid presentation and is normalized and scaled as:
+
+.. math:: S = C*\frac{s}{s+N}
+
+where :math:`s` is assumed to be the sum of the flowrates of any solids formed in the reaction. 
+
+Only one of :math:`Q` and :math:`S` can be greater than zero at any time, which can be written in the form of a complementarity constraint as:
+
+.. math:: Q - max(0, Q-S) == 0
+
+The :math:`\max()` function is provided as an IDAES utility which provides a smooth max expression with mutable smoothing parameter, :math:`\epsilon`.
+
+**Parameters**
+
+.. csv-table::
+   :header: "Symbol", "Parameter Name", "Indices", "Default", "Description"
+
+   ":math:`O`", "reaction_order", "liq, component", "\-", "Reaction order"    
+    ":math:`C`", "s_scale", "\-", "10", "Scaling factor for the solid term"    
+    ":math:`N`", "s_norm", "\-", "1e-4 (`k_eq_ref` if provided)", "Normalization parameter for the solid term" 
+    ":math:`\epsilon`", "eps", "\-", "1e-4", "Smoothing factor"
+
+Providing a `reaction_order` dict is optional. If one is not provided, it will be assumed that this is an elementary reaction and that the reaction order is equal to the stoichiometric coefficient for all components in non-solid phases (the contribution of solid phases is assumed to be constant and included in the equilibrium constant, thus an order of zero is assumed).
+
+Log Solubility Product (log_solubility_product)
+-----------------------------------------------
+
+The method uses a complementarity formulation for solid precipitation using a log form of solubility product.
+
+For precipitation at equilibrium, 
+
+.. math:: Q = log(k_{sp}) - \sum_{(liq, j)}{O_{(liq,j)}*log(x_{(liq,j)})}
+
+* :math:`Q = 0` only holds if solids (S) are present in the system, :math:`S \geq 0`.
+* :math:`Q \geq 0` if the system is subsaturated and no solids (S) are present, :math:`S = 0`.
+
+:math:`S` represents the solid presentation and is normalized and scaled as:
+
+.. math:: S = C\frac{s}{s+N}
+
+where :math:`s` is assumed to be the sum of the flowrates of any solids formed in the reaction. 
+
+Only one of :math:`Q` and :math:`S` can be greater than zero at any time, which can be written in the form of a complementarity constraint as:
+
+.. math:: Q - max(0, Q-S) == 0
+
+The :math:`\max()` function is provided as an IDAES utility which provides a smooth max expression with mutable smoothing parameter, :math:`\epsilon`.
+
+**Parameters**
+
+.. csv-table::
+   :header: "Symbol", "Parameter Name", "Indices", "Default", "Description"
+
+   ":math:`O`", "reaction_order", "liq, component", "\-", "Reaction order"    
+    ":math:`C`", "s_scale", "\-", "1", "Scaling factor for the solid term"    
+    ":math:`N`", "s_norm", "\-", "1e-4 (`k_eq_ref` if provided)", "Normalization parameter for the solid term" 
+    ":math:`\epsilon`", "eps", "\-", "1e-4", "Smoothing factor"
+
 Providing a `reaction_order` dict is optional. If one is not provided, it will be assumed that this is an elementary reaction and that the reaction order is equal to the stoichiometric coefficient for all component in non-solid phases (the contribution of solid phases is assumed to be constant and included in the equilibrium constant, thus an order of zero is assumed).
+
diff --git a/idaes/models/properties/modular_properties/reactions/equilibrium_forms.py b/idaes/models/properties/modular_properties/reactions/equilibrium_forms.py
@@ -19,7 +19,7 @@
 # TODO: Look into protected access issues
 # pylint: disable=protected-access
 
-from pyomo.environ import Param, units as pyunits
+from pyomo.environ import Param, units as pyunits, value
 
 from idaes.core.util.math import smooth_max
 from idaes.core.util.exceptions import ConfigurationError
@@ -96,7 +96,7 @@ def calculate_scaling_factors(b, sf_keq):
 # ----------------------------------------------------------------------------
 class solubility_product:
     """
-    Complementariity formulation for solid precipitation
+    Complementarity formulation for solid precipitation
     Thanks to Larry Biegler for the formulation
 
     Like any phase equilibrium, solid precipitation is complicated by
@@ -113,7 +113,7 @@ class solubility_product:
     Thus, only one of S and Q can be greater than zero at any time.
     This can be written in the form of a complementarity constraint as:
 
-        * S - MAX(0, S-Q) == 0
+        * Q - MAX(0, Q-S) == 0
 
     where S is assumed to be the sum of the flowrates any solids formed in the
     reaction. This allows for multiple solid products, and only applies the
@@ -131,6 +131,20 @@ def build_parameters(rblock, config):
             mutable=True, initialize=1e-4, doc="Smoothing parameter for smooth maximum"
         )
 
+        rblock.s_norm = Param(
+            mutable=True,
+            initialize=1e-4,
+            doc="Normalizing factor for solid precipitation term",
+        )
+        if hasattr(rblock, "k_eq_ref"):
+            rblock.s_norm.set_value(value(rblock.k_eq_ref))
+
+        rblock.s_scale = Param(
+            mutable=True,
+            initialize=1,
+            doc="Scaling factor for solid precipitation term w.r.t saturated status Q = Ksp - f(C)",
+        )
+
     @staticmethod
     def return_expression(b, rblock, r_idx, T):
         e = None
@@ -181,6 +195,8 @@ def return_expression(b, rblock, r_idx, T):
             if sunits is not None:
                 s = s / sunits
 
+        s = rblock.s_scale * s / (s + rblock.s_norm)
+
         Q = b.k_eq[r_idx] - e
 
         # Need to remove units again
@@ -189,7 +205,7 @@ def return_expression(b, rblock, r_idx, T):
         if Qunits is not None:
             Q = Q / Qunits
 
-        return s - smooth_max(0, s - Q, rblock.eps) == 0
+        return Q - smooth_max(0, Q - s, rblock.eps) == 0
 
     @staticmethod
     def calculate_scaling_factors(b, sf_keq):
@@ -198,7 +214,7 @@ def calculate_scaling_factors(b, sf_keq):
 
 class log_solubility_product:
     """
-    Complementariity formulation for solid precipitation
+    Complementarity formulation for solid precipitation
     Thanks to Larry Biegler for the formulation
 
     Like any phase equilibrium, solid precipitation is complicated by
@@ -216,7 +232,7 @@ class uses a log form to represent to equilibrium constraint.
     Thus, only one of S and Q can be greater than zero at any time.
     This can be written in the form of a complementarity constraint as:
 
-        * S - MAX(0, S-Q) == 0
+        * Q - MAX(0, Q-S) == 0
 
     where S is assumed to be the sum of the flowrates any solids formed in the
     reaction. This allows for multiple solid products, and only applies the
@@ -234,6 +250,20 @@ def build_parameters(rblock, config):
             mutable=True, initialize=1e-4, doc="Smoothing parameter for smooth maximum"
         )
 
+        rblock.s_norm = Param(
+            mutable=True,
+            initialize=1e-4,
+            doc="Normalizing factor for solid precipitation term",
+        )
+        if hasattr(rblock, "k_eq_ref"):
+            rblock.s_norm.set_value(value(rblock.k_eq_ref))
+
+        rblock.s_scale = Param(
+            mutable=True,
+            initialize=10,
+            doc="Scaling factor for solid precipitation term w.r.t saturated status Q = ln(Ksp) - ln(f(C))",
+        )
+
     @staticmethod
     def return_expression(b, rblock, r_idx, T):
         e = None
@@ -284,10 +314,12 @@ def return_expression(b, rblock, r_idx, T):
             if sunits is not None:
                 s = s / sunits
 
+        s = rblock.s_scale * s / (s + rblock.s_norm)
+
         Q = b.log_k_eq[r_idx] - e
         # Q should be unitless due to log form
 
-        return s - smooth_max(0, s - Q, rblock.eps) == 0
+        return Q - smooth_max(0, Q - s, rblock.eps) == 0
 
     @staticmethod
     def calculate_scaling_factors(b, sf_keq):

diff --git a/idaes/models/properties/modular_properties/reactions/tests/test_equilibrium_forms.py b/idaes/models/properties/modular_properties/reactions/tests/test_equilibrium_forms.py
@@ -414,6 +414,10 @@ def test_solubility_no_order():
     # Check parameter construction
     assert isinstance(m.rparams.reaction_r1.eps, Param)
     assert value(m.rparams.reaction_r1.eps) == 1e-4
+    assert isinstance(m.rparams.reaction_r1.s_norm, Param)
+    assert value(m.rparams.reaction_r1.s_norm) == 1e-4
+    assert isinstance(m.rparams.reaction_r1.s_norm, Param)
+    assert value(m.rparams.reaction_r1.s_scale) == 1
     assert isinstance(m.rparams.reaction_r1.reaction_order, Var)
     assert len(m.rparams.reaction_r1.reaction_order) == 6
 
@@ -434,16 +438,19 @@ def test_solubility_no_order():
         m.thermo[1].flow_mol_phase_comp["sol", "c1"]
         + m.thermo[1].flow_mol_phase_comp["sol", "c2"]
     ) / (pyunits.mol / pyunits.s)
-    Q = m.rxn[1].k_eq["r1"] - (
-        m.thermo[1].mole_frac_phase_comp["p1", "c1"]
-        ** m.rparams.reaction_r1.reaction_order["p1", "c1"]
-        * m.thermo[1].mole_frac_phase_comp["p1", "c2"]
-        ** m.rparams.reaction_r1.reaction_order["p1", "c2"]
-    )
+    s = m.rparams.reaction_r1.s_scale * s / (s + m.rparams.reaction_r1.s_norm)
 
-    assert str(rform) == str(
-        s - smooth_max(0, s - Q / pyunits.dimensionless, m.rparams.reaction_r1.eps) == 0
-    )
+    Q = (
+        m.rxn[1].k_eq["r1"]
+        - (
+            m.thermo[1].mole_frac_phase_comp["p1", "c1"]
+            ** m.rparams.reaction_r1.reaction_order["p1", "c1"]
+            * m.thermo[1].mole_frac_phase_comp["p1", "c2"]
+            ** m.rparams.reaction_r1.reaction_order["p1", "c2"]
+        )
+    ) / (pyunits.dimensionless)
+
+    assert str(rform) == str(Q - smooth_max(0, Q - s, m.rparams.reaction_r1.eps) == 0)
 
 
 @pytest.mark.unit
@@ -505,6 +512,10 @@ def test_solubility_product_with_order():
     # Check parameter construction
     assert isinstance(m.rparams.reaction_r1.eps, Param)
     assert value(m.rparams.reaction_r1.eps) == 1e-4
+    assert isinstance(m.rparams.reaction_r1.s_norm, Param)
+    assert value(m.rparams.reaction_r1.s_norm) == 1e-4
+    assert isinstance(m.rparams.reaction_r1.s_norm, Param)
+    assert value(m.rparams.reaction_r1.s_scale) == 1
     assert isinstance(m.rparams.reaction_r1.reaction_order, Var)
     assert len(m.rparams.reaction_r1.reaction_order) == 6
     assert m.rparams.reaction_r1.reaction_order["p1", "c1"].value == 1
@@ -523,24 +534,27 @@ def test_solubility_product_with_order():
         m.thermo[1].flow_mol_phase_comp["sol", "c1"]
         + m.thermo[1].flow_mol_phase_comp["sol", "c2"]
     ) / (pyunits.mol / pyunits.s)
-    Q = m.rxn[1].k_eq["r1"] - (
-        m.thermo[1].mole_frac_phase_comp["p1", "c1"]
-        ** m.rparams.reaction_r1.reaction_order["p1", "c1"]
-        * m.thermo[1].mole_frac_phase_comp["p1", "c2"]
-        ** m.rparams.reaction_r1.reaction_order["p1", "c2"]
-        * m.thermo[1].mole_frac_phase_comp["p2", "c1"]
-        ** m.rparams.reaction_r1.reaction_order["p2", "c1"]
-        * m.thermo[1].mole_frac_phase_comp["p2", "c2"]
-        ** m.rparams.reaction_r1.reaction_order["p2", "c2"]
-        * m.thermo[1].mole_frac_phase_comp["sol", "c1"]
-        ** m.rparams.reaction_r1.reaction_order["sol", "c1"]
-        * m.thermo[1].mole_frac_phase_comp["sol", "c2"]
-        ** m.rparams.reaction_r1.reaction_order["sol", "c2"]
-    )
+    s = m.rparams.reaction_r1.s_scale * s / (s + m.rparams.reaction_r1.s_norm)
 
-    assert str(rform) == str(
-        s - smooth_max(0, s - Q / pyunits.dimensionless, m.rparams.reaction_r1.eps) == 0
-    )
+    Q = (
+        m.rxn[1].k_eq["r1"]
+        - (
+            m.thermo[1].mole_frac_phase_comp["p1", "c1"]
+            ** m.rparams.reaction_r1.reaction_order["p1", "c1"]
+            * m.thermo[1].mole_frac_phase_comp["p1", "c2"]
+            ** m.rparams.reaction_r1.reaction_order["p1", "c2"]
+            * m.thermo[1].mole_frac_phase_comp["p2", "c1"]
+            ** m.rparams.reaction_r1.reaction_order["p2", "c1"]
+            * m.thermo[1].mole_frac_phase_comp["p2", "c2"]
+            ** m.rparams.reaction_r1.reaction_order["p2", "c2"]
+            * m.thermo[1].mole_frac_phase_comp["sol", "c1"]
+            ** m.rparams.reaction_r1.reaction_order["sol", "c1"]
+            * m.thermo[1].mole_frac_phase_comp["sol", "c2"]
+            ** m.rparams.reaction_r1.reaction_order["sol", "c2"]
+        )
+    ) / (pyunits.dimensionless)
+
+    assert str(rform) == str(Q - smooth_max(0, Q - s, m.rparams.reaction_r1.eps) == 0)
 
 
 @pytest.mark.unit
@@ -647,6 +661,10 @@ def test_log_solubility_no_order():
     # Check parameter construction
     assert isinstance(m.rparams.reaction_r1.eps, Param)
     assert value(m.rparams.reaction_r1.eps) == 1e-4
+    assert isinstance(m.rparams.reaction_r1.s_norm, Param)
+    assert value(m.rparams.reaction_r1.s_norm) == 1e-4
+    assert isinstance(m.rparams.reaction_r1.s_norm, Param)
+    assert value(m.rparams.reaction_r1.s_scale) == 10
     assert isinstance(m.rparams.reaction_r1.reaction_order, Var)
     assert len(m.rparams.reaction_r1.reaction_order) == 6
 
@@ -667,14 +685,16 @@ def test_log_solubility_no_order():
         m.thermo[1].flow_mol_phase_comp["sol", "c1"]
         + m.thermo[1].flow_mol_phase_comp["sol", "c2"]
     ) / (pyunits.mol / pyunits.s)
+    s = m.rparams.reaction_r1.s_scale * s / (s + m.rparams.reaction_r1.s_norm)
+
     Q = m.rxn[1].log_k_eq["r1"] - (
         m.thermo[1].log_mole_frac_phase_comp["p1", "c1"]
         * m.rparams.reaction_r1.reaction_order["p1", "c1"]
         + m.thermo[1].log_mole_frac_phase_comp["p1", "c2"]
         * m.rparams.reaction_r1.reaction_order["p1", "c2"]
     )
 
-    assert str(rform) == str(s - smooth_max(0, s - Q, m.rparams.reaction_r1.eps) == 0)
+    assert str(rform) == str(Q - smooth_max(0, Q - s, m.rparams.reaction_r1.eps) == 0)
 
 
 @pytest.mark.unit
@@ -739,6 +759,10 @@ def test_log_solubility_product_with_order():
     # Check parameter construction
     assert isinstance(m.rparams.reaction_r1.eps, Param)
     assert value(m.rparams.reaction_r1.eps) == 1e-4
+    assert isinstance(m.rparams.reaction_r1.s_norm, Param)
+    assert value(m.rparams.reaction_r1.s_norm) == 1e-4
+    assert isinstance(m.rparams.reaction_r1.s_norm, Param)
+    assert value(m.rparams.reaction_r1.s_scale) == 10
     assert isinstance(m.rparams.reaction_r1.reaction_order, Var)
     assert len(m.rparams.reaction_r1.reaction_order) == 6
     assert m.rparams.reaction_r1.reaction_order["p1", "c1"].value == 1
@@ -757,6 +781,8 @@ def test_log_solubility_product_with_order():
         m.thermo[1].flow_mol_phase_comp["sol", "c1"]
         + m.thermo[1].flow_mol_phase_comp["sol", "c2"]
     ) / (pyunits.mol / pyunits.s)
+    s = m.rparams.reaction_r1.s_scale * s / (s + m.rparams.reaction_r1.s_norm)
+
     Q = m.rxn[1].log_k_eq["r1"] - (
         m.thermo[1].log_mole_frac_phase_comp["p1", "c1"]
         * m.rparams.reaction_r1.reaction_order["p1", "c1"]
@@ -772,7 +798,7 @@ def test_log_solubility_product_with_order():
         * m.rparams.reaction_r1.reaction_order["sol", "c2"]
     )
 
-    assert str(rform) == str(s - smooth_max(0, s - Q, m.rparams.reaction_r1.eps) == 0)
+    assert str(rform) == str(Q - smooth_max(0, Q - s, m.rparams.reaction_r1.eps) == 0)
 
 
 @pytest.mark.unit