Reformulate solubility_products constraints

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).
+

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 @@
 # ----------------------------------------------------------------------------
 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 @@
     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 @@
             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 @@
             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 @@
         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 @@
 
 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 @@
     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 @@
             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 @@
             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):

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