Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Reformulate solubility_products constraints #1197

Merged
merged 23 commits into from
Jul 27, 2023
Merged

Reformulate solubility_products constraints #1197

Show file tree
Hide file tree
Changes from 17 commits
Commits
Show all changes
23 commits
Select commit Hold shift + click to select a range
513796e
Modified log_solubility_product function to accomodate the comparison…
May 26, 2023
517c0a1
Remove comments
Jun 7, 2023
6463f3b
run black
bknueven Jun 7, 2023
c940637
Update state_difinitions vars lower bouds
Jun 14, 2023
99db9c3
Add mutable parameters for reformated solubility_product
Jun 14, 2023
2b58efa
Update solubility_product tests
Jun 14, 2023
1edee16
Merge branch 'main' into chem_init
andrewlee94 Jun 15, 2023
d5b2aa3
test_plate_heat_exchanger.py
Jun 15, 2023
d598545
run black
Jun 15, 2023
28aed85
Reverse changes for lower bounds in state_definition
Jun 16, 2023
6eee6b5
Merge branch 'main' into chem_init
andrewlee94 Jun 16, 2023
4cd060d
Add assertion for default s_norm and a_scale checking
Jun 20, 2023
d6462cd
Typo corrected
Jun 20, 2023
89b6a2f
Add documentation for solubility_product
lxhowl Jun 20, 2023
a990507
Update doc string
Jun 20, 2023
d42a687
Merge branch 'IDAES:main' into chem_init
lxhowl Jun 20, 2023
5306d36
Equation corrected for sulubility product doc
Jun 22, 2023
4ca5c87
Merge branch 'IDAES:main' into chem_init
lxhowl Jul 5, 2023
e4a6376
Merge branch 'main' into chem_init
lbianchi-lbl Jul 7, 2023
bc88595
Switch to set_value for default s_norm
Jul 11, 2023
214f14f
Scaling applied in subsaturated cases (low concentrations) to use def…
Jul 18, 2023
7bec453
Run black
Jul 20, 2023
14cb740
Merge branch 'main' into chem_init
andrewlee94 Jul 27, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
81 changes: 77 additions & 4 deletions docs/explanations/components/property_package/general_reactions/equil_form.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)}}}

Expand All @@ -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**

Expand All @@ -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).

44 changes: 38 additions & 6 deletions idaes/models/properties/modular_properties/reactions/equilibrium_forms.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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.value = rblock.k_eq_ref.value
lxhowl marked this conversation as resolved.
Show resolved Hide resolved

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
Expand Down Expand Up @@ -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
Expand All @@ -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):
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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.value = rblock.k_eq_ref.value
lxhowl marked this conversation as resolved.
Show resolved Hide resolved

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
Expand Down Expand Up @@ -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
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If one or both Q and s need to be zero, could this just bee smooth min(Q, s) == 0?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Current expression goes Q-0.5(Q-s+|Q-s|) == 0 while the smooth min gives 0.5(Q+s-|Q-s|)==0. The complexity remains the same. The current expression keeps Q as a leading term, which might be easy for users to scale the solubility product with other equilibrium reactions.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@eslickj I think this formulation is correct - it is a bit more complicated than just one of Q or s needing to be zero.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@andrewlee94, I'm sure you are probably right. Can someone explain the difference to me? Why is this one preferred? I think both should work, but I may be missing something. To me, the the one in the code currently just seems convoluted.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@eslickj I cannot remember the details - I asked Larry Biegler to help with this is and this is the form he suggested.


@staticmethod
def calculate_scaling_factors(b, sf_keq):
Expand Down
82 changes: 54 additions & 28 deletions idaes/models/properties/modular_properties/reactions/tests/test_equilibrium_forms.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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"]
Expand All @@ -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
Expand Down
5 changes: 3 additions & 2 deletions ...els/properties/modular_properties/reactions/tests/test_solubility_product_verification.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@


solver = get_solver()
solver.options["tol"] = 1e-8
lxhowl marked this conversation as resolved.
Show resolved Hide resolved


def dummy_h(b, *args, **kwargs):
Expand Down Expand Up @@ -283,10 +284,10 @@ def test_saturated(self, model):
assert pytest.approx(0, abs=1.1e-3) == value(
model.fs.R101.outlet.flow_mol_phase_comp[0, "Sol", "NaCl"]
)
assert pytest.approx(6.158968, rel=1e-5) == value(
assert pytest.approx(6.159876, rel=1e-5) == value(
model.fs.R101.outlet.flow_mol_phase_comp[0, "Liq", "Na+"]
)
assert pytest.approx(6.158968, rel=1e-5) == value(
assert pytest.approx(6.159876, rel=1e-5) == value(
model.fs.R101.outlet.flow_mol_phase_comp[0, "Liq", "Cl-"]
)

Expand Down