Adding critical properties to ideal EoS

idaes/models/properties/modular_properties/base/generic_property.py

@@ -2975,7 +2975,7 @@ def _critical_props(self):
             )
 
             p_config = self.params.get_phase(ref_phase).config
-            p_config.equation_of_state.build_critical_properties(self, ref_phase)
+            p_config.equation_of_state.build_critical_properties(self)
 
         except AttributeError:
             self.del_component(self.compress_fact_crit)

idaes/models/properties/modular_properties/eos/eos_base.py

@@ -51,13 +51,9 @@ def build_parameters(b):
         raise NotImplementedError(_msg(b, "build_parameters"))
 
     @staticmethod
-    def build_critical_properties(b, p):
+    def build_critical_properties(b):
         raise NotImplementedError(_msg(b, "build_critical_properties"))
 
-    @staticmethod
-    def initialize_critical_properties(b, p):
-        raise NotImplementedError(_msg(b, "initialize_critical_properties"))
-
     @staticmethod
     def get_vol_mol_pure(b, phase, comp, temperature):
         try:

idaes/models/properties/modular_properties/eos/ideal.py

@@ -21,7 +21,7 @@
 # TODO: Look into protected access issues
 # pylint: disable=protected-access
 
-from pyomo.environ import Expression, log
+from pyomo.environ import Constraint, Expression, log, units
 
 from idaes.core import Apparent
 from idaes.core.util.exceptions import ConfigurationError, PropertyNotSupportedError
@@ -34,6 +34,7 @@
     log_henry_pressure,
 )
 from .eos_base import EoSBase
+from idaes.core.util.constants import Constants as CONST
 
 
 # TODO: Add support for ideal solids
@@ -48,6 +49,41 @@ def common(b, pobj):
         # No common components required for ideal property calculations
         pass
 
+    @staticmethod
+    def build_critical_properties(b):
+        base_units = b.params.get_metadata().default_units
+
+        # Use Kay's method for critical pressure and temperature
+        def p_crit_rule(b):
+            return b.pressure_crit == sum(
+                b.mole_frac_comp[j] * b.params.get_component(j).pressure_crit
+                for j in b.component_list
+            )
+
+        b.pressure_crit_constraint = Constraint(rule=p_crit_rule)
+
+        def t_crit_rule(b):
+            return b.temperature_crit == sum(
+                b.mole_frac_comp[j] * b.params.get_component(j).temperature_crit
+                for j in b.component_list
+            )
+
+        b.temperature_crit_constraint = Constraint(rule=t_crit_rule)
+
+        # Fix critical compressibility factor to 1 - user can override if necessary
+        b.compress_fact_crit.fix(1)
+
+        def rho_crit_rule(b):
+            return b.pressure_crit == units.convert(
+                b.compress_fact_crit
+                * CONST.gas_constant
+                * b.temperature_crit
+                * b.dens_mol_crit,
+                to_units=base_units.PRESSURE,
+            )
+
+        b.dens_mol_crit_constraint = Constraint(rule=rho_crit_rule)
+
     @staticmethod
     def calculate_scaling_factors(b, pobj):
         pass

idaes/models/properties/modular_properties/eos/tests/test_ideal.py

@@ -18,7 +18,7 @@
 import pytest
 from sys import modules
 
-from pyomo.environ import ConcreteModel, Var, units as pyunits, value
+from pyomo.environ import ConcreteModel, Constraint, Var, units as pyunits, value
 from pyomo.util.check_units import assert_units_equivalent
 
 from idaes.core import (
@@ -686,3 +686,74 @@ def test_vol_mol_phase():
                 + 1 / 42.0 * m.props[1].mole_frac_phase_comp["Liq", "b"]
                 + 42.0 * m.props[1].mole_frac_phase_comp["Liq", "c"]
             )
+
+
+@pytest.mark.unit
+def test_critical_props():
+    m = ConcreteModel()
+
+    # Dummy params block
+    m.params = DummyParameterBlock(
+        components={
+            "a": {"parameter_data": {"pressure_crit": 1e5, "temperature_crit": 100}},
+            "b": {"parameter_data": {"pressure_crit": 2e5, "temperature_crit": 200}},
+            "c": {"parameter_data": {"pressure_crit": 3e5, "temperature_crit": 300}},
+        },
+        phases={
+            "Vap": {"type": VaporPhase, "equation_of_state": Ideal},
+            "Liq": {"type": LiquidPhase, "equation_of_state": Ideal},
+        },
+        base_units={
+            "time": pyunits.s,
+            "length": pyunits.m,
+            "mass": pyunits.kg,
+            "amount": pyunits.mol,
+            "temperature": pyunits.K,
+        },
+        state_definition=modules[__name__],
+        pressure_ref=100000.0,
+        temperature_ref=300,
+    )
+
+    m.props = m.params.state_block_class([1], defined_state=False, parameters=m.params)
+
+    # Add common variables
+    m.props[1].mole_frac_comp = Var(m.params.component_list, initialize=0.5)
+
+    m.props[1]._critical_props()
+
+    assert isinstance(m.props[1].compress_fact_crit, Var)
+    assert isinstance(m.props[1].dens_mol_crit, Var)
+    assert isinstance(m.props[1].pressure_crit, Var)
+    assert isinstance(m.props[1].temperature_crit, Var)
+
+    assert isinstance(m.props[1].pressure_crit_constraint, Constraint)
+    assert str(m.props[1].pressure_crit_constraint.expr) == str(
+        m.props[1].pressure_crit
+        == m.props[1].mole_frac_comp["a"] * m.params.a.pressure_crit
+        + m.props[1].mole_frac_comp["b"] * m.params.b.pressure_crit
+        + m.props[1].mole_frac_comp["c"] * m.params.c.pressure_crit
+    )
+
+    assert isinstance(m.props[1].temperature_crit_constraint, Constraint)
+    assert str(m.props[1].temperature_crit_constraint.expr) == str(
+        m.props[1].temperature_crit
+        == m.props[1].mole_frac_comp["a"] * m.params.a.temperature_crit
+        + m.props[1].mole_frac_comp["b"] * m.params.b.temperature_crit
+        + m.props[1].mole_frac_comp["c"] * m.params.c.temperature_crit
+    )
+
+    assert m.props[1].compress_fact_crit.fixed
+    assert value(m.props[1].compress_fact_crit) == 1
+
+    assert isinstance(m.props[1].dens_mol_crit_constraint, Constraint)
+    assert str(m.props[1].dens_mol_crit_constraint.expr) == str(
+        m.props[1].pressure_crit
+        == pyunits.convert(
+            m.props[1].compress_fact_crit
+            * const.gas_constant
+            * m.props[1].temperature_crit
+            * m.props[1].dens_mol_crit,
+            to_units=pyunits.kg / pyunits.m / pyunits.s**2,
+        )
+    )