Add new Complementarity formualtion for VLE with cubic EoSs

idaes/models/properties/modular_properties/base/tests/test_generic_property.py

@@ -42,6 +42,8 @@
 )
 from idaes.core.util.exceptions import ConfigurationError, PropertyPackageError
 from idaes.models.properties.modular_properties.phase_equil.henry import HenryType
+from idaes.models.properties.modular_properties.eos.ceos import Cubic, CubicType
+from idaes.models.properties.modular_properties.state_definitions import FTPx
 from idaes.core.base.property_meta import UnitSet
 from idaes.core.initialization import BlockTriangularizationInitializer
 
@@ -1955,3 +1957,54 @@ def build_critical_properties(b, *args, **kwargs):
             "Component declarations.",
         ):
             _initialize_critical_props(m.props[1])
+
+
+# Invalid property configuration to trigger configuration error
+configuration = {
+    # Specifying components
+    "components": {
+        "H2O": {
+            "type": Component,
+            "parameter_data": {
+                "pressure_crit": (220.6e5, pyunits.Pa),
+                "temperature_crit": (647, pyunits.K),
+                "omega": 0.344,
+            },
+        },
+    },
+    # Specifying phases
+    "phases": {
+        "Liq": {
+            "type": LiquidPhase,
+            "equation_of_state": Cubic,
+            "equation_of_state_options": {"type": CubicType.PR},
+        },
+        "Vap": {
+            "type": VaporPhase,
+            "equation_of_state": Cubic,
+            "equation_of_state_options": {"type": CubicType.PR},
+        },
+    },
+    # Set base units of measurement
+    "base_units": {
+        "time": pyunits.s,
+        "length": pyunits.m,
+        "mass": pyunits.kg,
+        "amount": pyunits.mol,
+        "temperature": pyunits.K,
+    },
+    # Specifying state definition
+    "state_definition": FTPx,
+    "state_bounds": {
+        "flow_mol": (0, 100, 1000, pyunits.mol / pyunits.s),
+        "temperature": (273.15, 300, 500, pyunits.K),
+        "pressure": (5e4, 1e5, 1e6, pyunits.Pa),
+    },
+    "pressure_ref": (101325, pyunits.Pa),
+    "temperature_ref": (298.15, pyunits.K),
+    "parameter_data": {
+        "PR_kappa": {
+            ("foo", "bar"): 0.000,
+        }
+    },
+}

idaes/models/properties/modular_properties/base/tests/test_generic_property_integration.py

@@ -217,6 +217,7 @@ def test_heater_w_inherent_rxns_comp_phase(self, frame):
             rel=1e-5,
         )
 
+        frame.fs.H101.outlet.mole_frac_comp.display()
         assert value(frame.fs.H101.outlet.mole_frac_comp[0, "a"]) == pytest.approx(
             1 / 6, rel=1e-5
         )

idaes/models/properties/modular_properties/base/tests/test_utility.py

@@ -16,9 +16,13 @@
 Author: A Lee
 """
 import pytest
-
 from types import MethodType
 
+from pyomo.environ import Block, ConcreteModel, units as pyunits, Var
+from pyomo.common.config import ConfigBlock, ConfigValue
+
+from idaes.core import declare_process_block_class
+from idaes.core import Component, LiquidPhase, SolidPhase, VaporPhase, PhaseType as PT
 from idaes.models.properties.modular_properties.base.utility import (
     GenericPropertyPackageError,
     get_method,
@@ -27,13 +31,31 @@
     get_bounds_from_config,
     get_concentration_term,
     ConcentrationForm,
+    identify_VL_component_list,
+    estimate_Pbub,
+    estimate_Pdew,
+    estimate_Tbub,
+    estimate_Tdew,
+)
+from idaes.models.properties.modular_properties.base.generic_property import (
+    GenericParameterData,
 )
-
 from idaes.models.properties.modular_properties.base.generic_reaction import rxn_config
-from pyomo.environ import Block, units as pyunits, Var
-from pyomo.common.config import ConfigBlock, ConfigValue
 from idaes.core.util.exceptions import ConfigurationError, PropertyPackageError
 from idaes.core.util.misc import add_object_reference
+from idaes.models.properties.modular_properties.base.tests.dummy_eos import DummyEoS
+from idaes.models.properties.modular_properties.phase_equil.henry import (
+    ConstantH,
+    HenryType,
+)
+from idaes.models.properties.modular_properties.eos.ceos import Cubic, CubicType
+from idaes.models.properties.modular_properties.phase_equil import SmoothVLE
+from idaes.models.properties.modular_properties.phase_equil.bubble_dew import (
+    LogBubbleDew,
+)
+from idaes.models.properties.modular_properties.phase_equil.forms import log_fugacity
+from idaes.models.properties.modular_properties.pure import NIST
+from idaes.models.properties.modular_properties.state_definitions import FTPx
 
 
 @pytest.fixture
@@ -559,3 +581,253 @@ def test_inherent_partial_pressure(self, frame2):
             get_concentration_term(frame2, "i1", log=True)
             is frame2.log_pressure_phase_comp
         )
+
+
+class TestIdentifyVLComponentList:
+    # Declare a base units dict to save code later
+    base_units = {
+        "time": pyunits.s,
+        "length": pyunits.m,
+        "mass": pyunits.kg,
+        "amount": pyunits.mol,
+        "temperature": pyunits.K,
+    }
+
+    # Dummy methods for properties
+    def set_metadata(self, b):
+        pass
+
+    def define_state(self, b):
+        pass
+
+    @declare_process_block_class("DummyParameterBlock")
+    class DummyParameterData(GenericParameterData):
+        def configure(self):
+            self.configured = True
+
+    @pytest.fixture()
+    def frame(self):
+        m = ConcreteModel()
+        m.params = DummyParameterBlock(
+            components={
+                "a": {},
+                "b": {},
+                "c": {},
+            },
+            phases={
+                "p1": {
+                    "type": LiquidPhase,
+                    "equation_of_state": DummyEoS,
+                },
+                "p2": {
+                    "type": VaporPhase,
+                    "equation_of_state": DummyEoS,
+                },
+                "p3": {
+                    "type": SolidPhase,
+                    "equation_of_state": DummyEoS,
+                },
+            },
+            state_definition=self,
+            pressure_ref=100000.0,
+            temperature_ref=300,
+            base_units=self.base_units,
+        )
+
+        m.props = m.params.build_state_block([1], defined_state=False)
+
+        return m
+
+    @pytest.mark.unit
+    def test_invalid_VL_pair(self, frame):
+        with pytest.raises(
+            PropertyPackageError,
+            match="Phase pair p1-p3 was identified as being a VLE pair, "
+            "however p1 is liquid but p3 is not vapor.",
+        ):
+            identify_VL_component_list(frame.props[1], ("p1", "p3"))
+
+        with pytest.raises(
+            PropertyPackageError,
+            match="Phase pair p2-p3 was identified as being a VLE pair, "
+            "however neither p2 nor p3 is liquid.",
+        ):
+            identify_VL_component_list(frame.props[1], ("p2", "p3"))
+
+    @pytest.mark.unit
+    def test_all_components(self, frame):
+        l_phase, v_phase, vl_comps, henry_comps, l_only_comps, v_only_comps = (
+            identify_VL_component_list(frame.props[1], ("p1", "p2"))
+        )
+
+        assert l_phase == "p1"
+        assert v_phase == "p2"
+        assert vl_comps == ["a", "b", "c"]
+        assert henry_comps == []
+        assert l_only_comps == []
+        assert v_only_comps == []
+
+    @pytest.mark.unit
+    def test_all_types_components(self):
+        m = ConcreteModel()
+        m.params = DummyParameterBlock(
+            components={
+                "a": {},
+                "b": {
+                    "valid_phase_types": PT.liquidPhase,
+                },
+                "c": {
+                    "valid_phase_types": PT.vaporPhase,
+                },
+                "d": {
+                    "valid_phase_types": PT.solidPhase,
+                },
+                "e": {
+                    "parameter_data": {"henry_ref": {"p1": 86}},
+                    "henry_component": {
+                        "p1": {"method": ConstantH, "type": HenryType.Kpx}
+                    },
+                },
+            },
+            phases={
+                "p1": {
+                    "type": LiquidPhase,
+                    "equation_of_state": DummyEoS,
+                },
+                "p2": {
+                    "type": VaporPhase,
+                    "equation_of_state": DummyEoS,
+                },
+                "p3": {
+                    "type": SolidPhase,
+                    "equation_of_state": DummyEoS,
+                },
+            },
+            state_definition=self,
+            pressure_ref=100000.0,
+            temperature_ref=300,
+            base_units=self.base_units,
+        )
+
+        m.props = m.params.build_state_block([1], defined_state=False)
+
+        l_phase, v_phase, vl_comps, henry_comps, l_only_comps, v_only_comps = (
+            identify_VL_component_list(m.props[1], ("p1", "p2"))
+        )
+        assert l_phase == "p1"
+        assert v_phase == "p2"
+        assert vl_comps == ["a"]
+        assert henry_comps == ["e"]
+        assert l_only_comps == ["b"]
+        assert v_only_comps == ["c"]
+
+
+# Property configuration for pure water to use in bubble and dew point tests
+configuration = {
+    # Specifying components
+    "components": {
+        "H2O": {
+            "type": Component,
+            "pressure_sat_comp": NIST,
+            "phase_equilibrium_form": {("Vap", "Liq"): log_fugacity},
+            "parameter_data": {
+                "pressure_crit": (220.6e5, pyunits.Pa),
+                "temperature_crit": (647, pyunits.K),
+                "omega": 0.344,
+                "pressure_sat_comp_coeff": {
+                    "A": (3.55959, None),
+                    "B": (643.748, pyunits.K),
+                    "C": (-198.043, pyunits.K),
+                },
+            },
+        },
+    },
+    # Specifying phases
+    "phases": {
+        "Liq": {
+            "type": LiquidPhase,
+            "equation_of_state": Cubic,
+            "equation_of_state_options": {"type": CubicType.PR},
+        },
+        "Vap": {
+            "type": VaporPhase,
+            "equation_of_state": Cubic,
+            "equation_of_state_options": {"type": CubicType.PR},
+        },
+    },
+    # Set base units of measurement
+    "base_units": {
+        "time": pyunits.s,
+        "length": pyunits.m,
+        "mass": pyunits.kg,
+        "amount": pyunits.mol,
+        "temperature": pyunits.K,
+    },
+    # Specifying state definition
+    "state_definition": FTPx,
+    "state_bounds": {
+        "flow_mol": (0, 100, 1000, pyunits.mol / pyunits.s),
+        "temperature": (273.15, 300, 500, pyunits.K),
+        "pressure": (5e4, 1e5, 1e6, pyunits.Pa),
+    },
+    "pressure_ref": (101325, pyunits.Pa),
+    "temperature_ref": (298.15, pyunits.K),
+    # Defining phase equilibria
+    "phases_in_equilibrium": [("Vap", "Liq")],
+    "phase_equilibrium_state": {("Vap", "Liq"): SmoothVLE},
+    "bubble_dew_method": LogBubbleDew,
+    "parameter_data": {
+        "PR_kappa": {
+            ("H2O", "H2O"): 0.000,
+        }
+    },
+}
+
+
+class TestBubbleDewPoints:
+    @pytest.fixture
+    def model(self):
+        m = ConcreteModel()
+        m.params = DummyParameterBlock(
+            **configuration,
+        )
+
+        m.props = m.params.build_state_block([1], defined_state=False)
+
+        return m
+
+    @pytest.mark.unit
+    def test_bubble_temperature(self, model):
+        # Test bubble temperature at atmospheric pressure
+        model.props[1].pressure.set_value(101325)
+        Tbub = estimate_Tbub(model.props[1], pyunits.K, ["H2O"], [], "Liq")
+
+        # Expected value = 379.1828 from parameters used
+        assert Tbub == pytest.approx(379.1828, rel=1e-6)
+
+    @pytest.mark.unit
+    def test_dew_temperature(self, model):
+        # Test dew temperature at atmospheric pressure
+        model.props[1].pressure.set_value(101325)
+        Tdew = estimate_Tdew(model.props[1], pyunits.K, ["H2O"], [], "Liq")
+
+        # Expected value = 379.1828 from parameters used
+        assert Tdew == pytest.approx(379.1828, rel=1e-6)
+
+    @pytest.mark.unit
+    def test_bubble_pressure(self, model):
+        # Test bubble pressure at 100C
+        model.props[1].temperature.set_value(373.15)
+        Pbub = estimate_Pbub(model.props[1], ["H2O"], [], "Liq")
+
+        # Expected value = 76432.45 from parameters used
+        assert Pbub == pytest.approx(76432.45, rel=1e-6)
+
+    @pytest.mark.unit
+    def test_dew_pressure(self, model):
+        # Test dew pressure at 100C
+        model.props[1].temperature.set_value(373.15)
+        Pdew = estimate_Pdew(model.props[1], ["H2O"], [], "Liq")
+
+        # Expected value = 76432.45 from parameters used
+        assert Pdew == pytest.approx(76432.45, rel=1e-6)