Support calcuation of critical properties of mixtures in Modular Prop…

…erties (#1336)

* Adding volume_mass as unit and renaming molar_volume

* Hooks for mixture critical properties

* Hook for link ot EoS for critical properties

* Adding critical properties to ideal EoS

* Initial working cubic critical properties

* Adding basic initialization

* Refining crit prop initialization

* Finishing initialization

* Fixing pylint and removing ideal crit props

* Docs and extra error message

* Test for new error message

* Fixing typo

* Fixing Sphinx and pylint errors

* Address PR review

* Fixing typo

docs/explanations/components/property_package/general/eos/cubic.rst

@@ -67,6 +67,7 @@ Cubic equations of state require the following parameters to be defined:
 
 1. `omega` (Pitzer acentricity factor) needs to be defined for each component (in the `parameter_data` for each component).
 2. `kappa` (binary interaction parameters) needs to be defined for each component pair in the system. This parameter needs to be defined in the general `parameter_data` argument for the overall property package (as it can be used in multiple phases).
+3. Component critical properties (`compress_fact_crit`, `dens_mol_crit`, `pressure_crit`, and `temperature_crit`) are often required as well.
 
 Calculation of Properties
 -------------------------
@@ -77,6 +78,26 @@ Many thermophysical properties are calculated using an ideal and residual term,
 
 The residual term is derived from the partial derivatives of the cubic equation of state, whilst the ideal term is determined using pure component properties for the ideal gas phase defined for each component.
 
+Critical Properties of Mixtures
+-------------------------------
+
+The critical properties of the mixture are calculated by solving the cubic equation of state at the critical point using the following constraints.
+
+.. math:: \Omega_A R^2 T_c^2 = a_{c,m}P_c
+.. math:: \Omega_BRT_c = b_{c,m}P_c
+.. math:: 0 = Z_c^3 - (1+B_c-uB_c)Z_c^2 + (A_c-uB_c-(u-w)B_c^2)Z_c - A_cB_c-wB_c^2-wB_c^3
+.. math:: P_c = Z_cRT_c\rho_{c}
+
+
+with the following expressions:
+
+.. math:: a_{c,j} = \frac{\Omega_AR^2T_{c,j}^2}{P_{c, j}}\alpha_j(T_c)
+.. math:: a_{c,m} = \sum_i{\sum_j{y_iy_j(a_{c,i}a_{c,j})^{1/2}(1-\kappa_{ij})}}
+.. math:: b_{c,m} = \sum_i{y_ib_i}
+.. math:: A_c = \frac{a_{c,m}P_c}{R^2T_c^2}
+.. math:: B_c = \frac{b_{c,m}P_c}{RT_c}
+
+
 Mass Density by Phase
 ---------------------
 
@@ -109,14 +130,14 @@ Component Molar Enthalpy by Phase
 
 Component molar enthalpies by phase are calculated using the pure component method provided by the users in the property package configuration arguments.
 
-Molar Isobaric Heat Capcity (:math:`C_p`)
------------------------------------------
+Molar Isobaric Heat Capacity (:math:`C_p`)
+------------------------------------------
 
-The ideal molar isobaric heat capcity term is calculated from the weighted sum of the (ideal) component molar isobaric heat capacity:
+The ideal molar isobaric heat capacity term is calculated from the weighted sum of the (ideal) component molar isobaric heat capacity:
 
 .. math:: C_{p, ig}^0 = \sum_j y_j C_{p, ig, j}
 
-The residual molar isobaric heat capcity term is given by:
+The residual molar isobaric heat capacity term is given by:
 
 .. math:: C_p^r = R \left[ T \left(\frac{\partial Z}{\partial T}\right)_P + Z - 1 \right] +  \frac{ T \frac{d^2a_m}{dT^2}}{\sqrt{u^2 - 4w} \cdot b_m} \ln \left[ \frac{2Z + uB + \sqrt{u^2 - 4w} B}{2Z + uB - \sqrt{u^2 - 4w} B} \right]
 .. math:: + \left(a_m - T \frac{da_m}{dT}\right) \cdot \frac{B}{b_m} \cdot \frac{\left(\frac{\partial Z}{\partial T}\right)_P + \frac{Z}{T}}{Z^2 + Z uB + wB^2}

docs/explanations/components/property_package/general/eos/ideal.rst

@@ -9,6 +9,11 @@ Introduction
 
 Ideal behavior represents the simplest possible equation of state that ensures thermodynamic consistency between different properties.
 
+Critical Properties of Mixtures
+-------------------------------
+
+The Ideal Equation of State module does not support the calculation of mixture critical properties as the ideal equation of state cannot represent the critical point.
+
 Mass Density by Phase
 ---------------------
 

docs/explanations/components/property_package/general/phase_def.rst

@@ -34,7 +34,7 @@ Equations of state (or equivalent methods) describe the relationship between dif
 
 A wide range of equations of states are available in literature for different applications and levels of rigor, and the IDAES Generic Property Package Framework provides a number of prebuilt modules for users, which are listed below.
 
-Equation of state packages may allow for user options (e.g. choosing a specific type of cubic equation of state). The options are set using the `equation_of_state_options` argument, and the options available are described in the documentation of each equation of state module.
+Equation of state packages may allow for user options (e.g., choosing a specific type of cubic equation of state). The options are set using the `equation_of_state_options` argument, and the options available are described in the documentation of each equation of state module.
 
 Equation of State Libraries
 """""""""""""""""""""""""""
@@ -44,6 +44,20 @@ Equation of State Libraries
 
     eos/ideal
     eos/cubic
+
+Critical Properties of Mixtures
+"""""""""""""""""""""""""""""""
+
+Calculation of the critical properties of mixtures depends on the equation of state being used. As the Modular Property Package framework allows users to specify different equations of state for each phase, the following logic is used to determine which equation of state to use for calculating critical properties.
+
+1. If a vapor-liquid equilibrium pair is defined in the `"phases_in_equilibrium"` configuration argument, then the liquid phase from this pair is used (IDAES generally assumes supercritical fluids are liquid-like).
+2. If no vapor-liquid equilibrium pair is defined, then the first liquid phase defined is used.
+3. If no liquid phases are defined then the vapor phase is used (it is assumed there will only be one vapor phase).
+4. If no vapor phase is defined, then a `PropertyPackageError` is returned as there is no suitable phase for calculating critical properties.
+
+Note that not all Equations of State are suitable for calculating critical properties as many cannot represent the critical conditions. The Ideal equation of State is one common example of an equation of state that DOES NOT support critical properties.
+
+During initialization, the critical properties of the mixture are approximated using the mole fraction weighted sum of the component critical properties (note that users must provide values for all component critical properties (i.e., `compress_fact_crit`, `dens_mol_crit`, `pressure_crit` and `temperature_crit`) for initialization if critical properties are to be calculated. Note that it is up to the user to ensure these values are consistent (if necessary).
 
 Phase-Specific Parameter
 ^^^^^^^^^^^^^^^^^^^^^^^^
@@ -53,6 +67,6 @@ In some cases, a property package may include parameters which are specific to a
 Phases with Partial Component Lists
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-In many applications a mixture will contain species that only appear in a single phase (either by nature or assumption). Common examples include crystalline solids and non-condensable gases. The IDAES Generic Property Package Framework provides support for these behaviors and allows users to specify phase-specific component lists (i.e. a list of components which appear in a given phase).
+In many applications a mixture will contain species that only appear in a single phase (either by nature or assumption). Common examples include crystalline solids and non-condensable gases. The IDAES Generic Property Package Framework provides support for these behaviors and allows users to specify phase-specific component lists (i.e., a list of components which appear in a given phase).
 
 This is done by providing a phase with a `component_list` argument, which provides a `list` of component names which appear in the phase. The framework automatically validates the `component_list` argument to ensure that it is a sub-set of the master component list for the property package, and will inform the user if an unrecognized component is included. If a phase is not provided with a `component_list` argument it is assumed that all components defined in the master component list may be present in the phase.

idaes/core/base/components.py

@@ -246,6 +246,7 @@ def build(self):
 
         # Create Vars for common parameters
         var_dict = {
+            "compress_fact_crit": pyunits.dimensionless,
             "dens_mol_crit": base_units.DENSITY_MOLE,
             "omega": pyunits.dimensionless,
             "pressure_crit": base_units.PRESSURE,

idaes/core/base/property_meta.py

@@ -257,9 +257,20 @@ def VOLUME(self):
         return self._length**3
 
     @property
-    def MOLAR_VOLUME(self):
+    def VOLUME_MASS(self):
+        return self._length**3 * self._mass**-1
+
+    @property
+    def VOLUME_MOLE(self):
         return self._length**3 * self._amount**-1
 
+    # Backward compatibility name
+    @property
+    def MOLAR_VOLUME(self):
+        msg = "The unit name MOLAR_VOLUME is being deprecated in favor of VOLUME_MOL."
+        deprecation_warning(msg=msg, logger=_log, version="2.3.0", remove_in="3.0.0")
+        return self.VOLUME_MOLE
+
     # Flows
     @property
     def FLOW_MASS(self):

idaes/core/base/property_set.py

@@ -627,6 +627,11 @@ class StandardPropertySet(PropertySetBase):
         doc="Compressiblity Factor",
         units=pyunits.dimensionless,
     )
+    compress_fact_crit = PropertyMetadata(
+        name="compress_fact_crit",
+        doc="Compressiblity Factor at Critical Point",
+        units=pyunits.dimensionless,
+    )
     conc_mass = PropertyMetadata(
         name="conc_mass",
         doc="Concentration on a Mass Basis",
@@ -890,12 +895,12 @@ class StandardPropertySet(PropertySetBase):
     vol_mol = PropertyMetadata(
         name="vol_mol",
         doc="Molar Volume",
-        units="MOLAR_VOLUME",
+        units="VOLUME_MOLE",
     )
     vol_mol_crit = PropertyMetadata(
         name="vol_mol",
         doc="Molar Volume at Critical Point",
-        units="MOLAR_VOLUME",
+        units="VOLUME_MOLE",
     )
     # Log terms
     log_act = PropertyMetadata(

idaes/core/base/tests/test_components.py

@@ -185,6 +185,8 @@ def get_metadata(self):
         m.comp = Component(
             parameter_data={
                 "mw": 10,
+                "compress_fact_crit": 1,
+                "dens_mol_crit": 55,
                 "pressure_crit": 1e5,
                 "temperature_crit": 500,
             }
@@ -193,6 +195,12 @@ def get_metadata(self):
         assert isinstance(m.comp.mw, Param)
         assert m.comp.mw.value == 10
 
+        assert isinstance(m.comp.compress_fact_crit, Var)
+        assert m.comp.compress_fact_crit.value == 1
+
+        assert isinstance(m.comp.dens_mol_crit, Var)
+        assert m.comp.dens_mol_crit.value == 55
+
         assert isinstance(m.comp.pressure_crit, Var)
         assert m.comp.pressure_crit.value == 1e5
 

idaes/core/base/tests/test_property_meta.py

@@ -175,6 +175,8 @@ def test_luminous_intensity(unit_set):
 derived_quantities = {
     "area": units.m**2,
     "volume": units.m**3,
+    "volume_mass": units.m**3 / units.kg,
+    "volume_mole": units.m**3 / units.mol,
     "flow_mass": units.kg * units.s**-1,
     "flow_mole": units.mol * units.s**-1,
     "flow_vol": units.m**3 * units.s**-1,