diff --git a/docs/explanations/components/property_package/general/pe/smooth_flash.rst b/docs/explanations/components/property_package/general/pe/smooth_flash.rst index 540429ba27..2264da6552 100644 --- a/docs/explanations/components/property_package/general/pe/smooth_flash.rst +++ b/docs/explanations/components/property_package/general/pe/smooth_flash.rst @@ -1,9 +1,12 @@ -Smooth Vapor-Liquid Equilibrium Formulation (``smooth_VLE``) -============================================================ +Smooth Vapor-Liquid Equilibrium Formulation (SmoothVLE) +======================================================= .. contents:: Contents :depth: 2 +.. note:: + For property packages using cubic Equations of State, there is an alternative :ref:`CubicComplementarityVLE ` class that may give better performance. + Source ------ diff --git a/docs/explanations/components/property_package/general/pe/smooth_vle2.rst b/docs/explanations/components/property_package/general/pe/smooth_vle2.rst new file mode 100644 index 0000000000..0afa557f5f --- /dev/null +++ b/docs/explanations/components/property_package/general/pe/smooth_vle2.rst @@ -0,0 +1,67 @@ +Cubic Smooth Vapor-Liquid Equilibrium Formulation (CubicComplementarityVLE) +=========================================================================== + +.. contents:: Contents + :depth: 2 + +.. note:: + This formulation for vapor-liquid equilibrium is only valid if using a cubic Equation of State. For other equations of state, use :ref:`SmoothVLE `. + +Source +------ + +Dabadghao, V., Ghouse, J., Eslick, J., Lee, A., Burgard, A., Miller, D., and Biegler, L., A Complementarity-based Vapor-Liquid Equilibrium Formulation for Equation-Oriented Simulation and Optimization, AIChE Journal, 2023, Volume 69(4), e18029. https://doi.org/10.1002/aic.18029 + +Introduction +------------ + +Often, a user may not know whether a state corresponds to a liquid, gas, or coexisting mixture. Even if a user knows the phase composition of a problem's initial condition, optimization may push the stream into or out of the two-phase region. Therefore, it is necessary to formulate phase equilibrium equations that are well-behaved for both one-phase and two-phase streams. + +To address this, the cubic smooth vapor-liquid equilibrium (VLE) formulation always solves the equilibrium equations at a condition where a valid two-phase solution exists. In situations where only a single phase is present, the phase equilibrium is solved at the either the bubble or dew point, where the second phase is just beginning to form; in this way, a non-trivial solution is guaranteed. Rather than explicitly calculate the bubble and dew points (as is done in the :ref:`non-cubic formulation `), this formulation leverages properties of the cubic equation of state to identify the "equilibrium temperature". + +Formulation +----------- + +.. note:: + For the full derivation of the cubic smooth VLE formulation, see the reference above. + +.. note:: + For consistency of naming between the cubic and non-cubic formulations, :math:`\bar{T}` is referred to as :math:`T_{eq}` in this document and the resulting model. + +The approach used by the smooth VLE formulation is to define an "equilibrium temperature" (:math:`T_{eq}`) at which the equilibrium calculations will be performed. The equilibrium temperature is defined such that: + +.. math:: T = T_{eq} - s_{vap} + s_{liq} + +where :math:`T` is the actual state temperature, and :math:`s_{liq}` and :math:`s_{vap}` are non-negative slack variables. For systems existing the the liquid-only region, :math:`s_{liq}` will be non-zero whilst :math:`s_{vap}=0` (indicating that the system is below the bubble point and thus :math:`T_{eq}>T`). Similarly, for systems in the vapor-only region, :math:`s_{vap}` will be non-zero whilst :math:`s_{liq}=0`. Finally, in the two-phase region, :math:`s_{liq}=s_{vap}=0`, indicating that :math:`T_{eq}=T`. + +In order to determine the values of :math:`s_{liq}` and :math:`s_{vap}`, the following complementarity constraints are written: + +.. math:: 0 = \min(s_{liq}, F_{liq}) +.. math:: 0 = \min(s_{vap}, F_{vap}) + +where :math:`F_{p}` is the flow rate of each phase :math:`p`. That is, for each phase (liquid and vapor), if there is any flowrate associated with that phase (i.e., the phase exists), its slack variable must be equal to zero. + +Additionally, the follow complementarities are written to constraint the roots of the cubic equation of state. + +.. math:: 0 = \min(g^{+}_{liq}, F_{liq}) +.. math:: 0 = \min(g^{-}_{vap}, F_{vap}) + +where :math:`g^{+}_p` and :math:`g^{-}_p` are another pair of non-negative slack variables associated with each phase :math:`p`. These slack variables are defined such that: + +.. math:: f''(Z_p) = g^{+}_{p} - g^{-}_{p} + +where :math:`f''(Z_p)` is the second derivative of the cubic equation of state written in terms of the compressibility factor :math:`Z_p` for each phase :math:`p`. + +Smooth Approximation +'''''''''''''''''''' + +In order to express the minimum operators in a tractable form, these equations are reformulated using the IDAES `smooth_min` function: + +.. math:: \min(a, b) = 0.5{\left[a + b - \sqrt{(a-b)^2 + \epsilon^2}\right]} + +Each complementarity requires a smoothing parameter, named :math:`\epsilon_T` and :math:`\epsilon_Z` for the temperature and cubic root constraints respectively. Within the IDAES model, these are rendered as ``eps_t_phase1_phase2`` and ``eps_z_phase1_phase2``, where ``phase1`` and ``phase2`` are the names assigned to the liquid and vapor phases in the property package (order will depend on the order these are declared). + +The tractability of the VLE problem depends heavily upon the values chosen for :math:`\epsilon_T` and :math:`\epsilon_Z`, with larger values resulting in smoother transitions at the phase boundaries (and thus increased tractability) at the expense of decreased accuracy near these points. It is recommended that users employ a 2-stage approach to solving these problems, starting with a larger value of :math:`\epsilon_T` and :math:`\epsilon_Z` initially to determine which region the solution lies in, followed by a second solve using smaller values to refine the solution. + +As a rule of thumb, the values of :math:`\epsilon_T` and :math:`\epsilon_Z` should be between 2 and 4 orders of magnitude smaller than the largest quantify involved in the smooth maximum operation. This means the value of :math:`\epsilon_T` should be based on the larger of :math:`T` and :math:`F_p`, whilst :math:`\epsilon_Z` should be based on the larger of :math:`f''(Z_p)` and :math:`F_p`. The value of :math:`f''(Z_p)` may be difficult to determine *a priori*, however :math:`F_p` is likely to dominate in most cases unless :math:`F_p` is small or :math:`P` is large. + diff --git a/docs/explanations/components/property_package/general/phase_equilibrium.rst b/docs/explanations/components/property_package/general/phase_equilibrium.rst index ad1e45979d..e21409580f 100644 --- a/docs/explanations/components/property_package/general/phase_equilibrium.rst +++ b/docs/explanations/components/property_package/general/phase_equilibrium.rst @@ -40,6 +40,7 @@ Phase Equilibrium State Libraries :maxdepth: 1 pe/smooth_flash + pe/smooth_vle2 Necessary Properties -------------------- diff --git a/idaes/models/properties/modular_properties/base/generic_property.py b/idaes/models/properties/modular_properties/base/generic_property.py index 0656e6a666..f49b88da40 100644 --- a/idaes/models/properties/modular_properties/base/generic_property.py +++ b/idaes/models/properties/modular_properties/base/generic_property.py @@ -13,7 +13,7 @@ """ Framework for generic property packages """ -# TODO: Look into protected access issues +# TODO: Pylint complains about variables with _x names as they are built by sub-classes # pylint: disable=protected-access # Import Pyomo libraries @@ -82,10 +82,14 @@ get_phase_method, GenericPropertyPackageError, StateIndex, + identify_VL_component_list, + estimate_Tbub, + estimate_Tdew, + estimate_Pbub, + estimate_Pdew, ) from idaes.models.properties.modular_properties.phase_equil.bubble_dew import ( LogBubbleDew, - _valid_VL_component_list as bub_dew_VL_comp_list, ) from idaes.models.properties.modular_properties.phase_equil.henry import HenryType @@ -879,13 +883,10 @@ def build(self): if build_parameters is not None: try: build_parameters(pobj) - except KeyError: + except KeyError as err: raise ConfigurationError( - "{} values were not defined for parameter {} in " - "phase {}. Please check the parameter_data " - "argument to ensure values are provided.".format( - self.name, a, p - ) + f"{self.name} - values were not defined for parameter {a} in " + f"phase {p}. {str(err)}" ) # Next, add inherent reactions if they exist @@ -1333,19 +1334,19 @@ def initialization_routine( # Bubble temperature initialization if hasattr(k, "_mole_frac_tbub"): - model._init_Tbub(k, T_units) + _init_Tbub(k, T_units) # Dew temperature initialization if hasattr(k, "_mole_frac_tdew"): - model._init_Tdew(k, T_units) + _init_Tdew(k, T_units) # Bubble pressure initialization if hasattr(k, "_mole_frac_pbub"): - model._init_Pbub(k, T_units) + _init_Pbub(k) # Dew pressure initialization if hasattr(k, "_mole_frac_pdew"): - model._init_Pdew(k, T_units) + _init_Pdew(k) # Solve bubble, dew, and critical point constraints for c in k.component_objects(Constraint): @@ -1786,19 +1787,19 @@ def initialize( # Bubble temperature initialization if hasattr(k, "_mole_frac_tbub"): - blk._init_Tbub(k, T_units) + _init_Tbub(k, T_units) # Dew temperature initialization if hasattr(k, "_mole_frac_tdew"): - blk._init_Tdew(k, T_units) + _init_Tdew(k, T_units) # Bubble pressure initialization if hasattr(k, "_mole_frac_pbub"): - blk._init_Pbub(k, T_units) + _init_Pbub(k) # Dew pressure initialization if hasattr(k, "_mole_frac_pdew"): - blk._init_Pdew(k, T_units) + _init_Pdew(k) # Solve bubble, dew, and critical point constraints for c in k.component_objects(Constraint): @@ -2075,334 +2076,6 @@ def release_state(blk, flags, outlvl=idaeslog.NOTSET): init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="properties") init_log.info_high("State released.") - def _init_Tbub(self, blk, T_units): - for pp in blk.params._pe_pairs: - raoult_comps, henry_comps = _valid_VL_component_list(blk, pp) - - if raoult_comps == []: - continue - if henry_comps != []: - # Need to get liquid phase name - # TODO This logic probably needs to be updated to break - # when LLE is added - if blk.params.get_phase(pp[0]).is_liquid_phase(): - l_phase = pp[0] - else: - l_phase = pp[1] - - # Use lowest component temperature_crit as starting point - # Starting high and moving down generally works better, - # as it under-predicts next step due to exponential form of - # Psat. - # Subtract 1 to avoid potential singularities at Tcrit - Tbub0 = ( - min( - blk.params.get_component(j).temperature_crit.value - for j in raoult_comps - ) - - 1 - ) - - err = 1 - counter = 0 - - # Newton solver with step limiter to prevent overshoot - # Tolerance only needs to be ~1e-1 - # Iteration limit of 30 - while err > 1e-1 and counter < 30: - f = value( - sum( - get_method(blk, "pressure_sat_comp", j)( - blk, blk.params.get_component(j), Tbub0 * T_units - ) - * blk.mole_frac_comp[j] - for j in raoult_comps - ) - + sum( - blk.mole_frac_comp[j] - * blk.params.get_component(j) - .config.henry_component[l_phase]["method"] - .return_expression(blk, l_phase, j, Tbub0 * T_units) - for j in henry_comps - ) - - blk.pressure - ) - df = value( - sum( - get_method(blk, "pressure_sat_comp", j)( - blk, blk.params.get_component(j), Tbub0 * T_units, dT=True - ) - * blk.mole_frac_comp[j] - for j in raoult_comps - ) - + sum( - blk.mole_frac_comp[j] - * blk.params.get_component(j) - .config.henry_component[l_phase]["method"] - .dT_expression(blk, l_phase, j, Tbub0 * T_units) - for j in henry_comps - ) - ) - - # Limit temperature step to avoid excessive overshoot - if f / df < -50: - Tbub1 = Tbub0 + 50 - elif f / df > 50: - Tbub1 = Tbub0 - 50 - else: - Tbub1 = Tbub0 - f / df - - err = abs(Tbub1 - Tbub0) - Tbub0 = Tbub1 - counter += 1 - - blk.temperature_bubble[pp].value = Tbub0 - - for j in raoult_comps: - blk._mole_frac_tbub[pp, j].value = value( - blk.mole_frac_comp[j] - * get_method(blk, "pressure_sat_comp", j)( - blk, blk.params.get_component(j), Tbub0 * T_units - ) - / blk.pressure - ) - if blk.is_property_constructed("log_mole_frac_tbub"): - blk.log_mole_frac_tbub[pp, j].value = value( - log(blk._mole_frac_tbub[pp, j]) - ) - - for j in henry_comps: - blk._mole_frac_tbub[pp, j].value = value( - blk.mole_frac_comp[j] - * blk.params.get_component(j) - .config.henry_component[l_phase]["method"] - .return_expression(blk, l_phase, j, Tbub0 * T_units) - / blk.pressure - ) - if blk.is_property_constructed("log_mole_frac_tbub"): - blk.log_mole_frac_tbub[pp, j].value = value( - log(blk._mole_frac_tbub[pp, j]) - ) - - def _init_Tdew(self, blk, T_units): - for pp in blk.params._pe_pairs: - raoult_comps, henry_comps = _valid_VL_component_list(blk, pp) - - if raoult_comps == []: - continue - if henry_comps != []: - # Need to get liquid phase name - if blk.params.get_phase(pp[0]).is_liquid_phase(): - l_phase = pp[0] - else: - l_phase = pp[1] - - if ( - hasattr(blk, "_mole_frac_tbub") - and blk.temperature_bubble[pp].value is not None - ): - # If Tbub has been calculated above, use this as the - # starting point - Tdew0 = blk.temperature_bubble[pp].value - else: - # Otherwise, use lowest component critical temperature - # as starting point - # Subtract 1 to avoid potential singularities at Tcrit - Tdew0 = ( - min( - blk.params.get_component(j).temperature_crit.value - for j in raoult_comps - ) - - 1 - ) - - err = 1 - counter = 0 - - # Newton solver with step limiter to prevent overshoot - # Tolerance only needs to be ~1e-1 - # Iteration limit of 30 - while err > 1e-1 and counter < 30: - f = value( - blk.pressure - * ( - sum( - blk.mole_frac_comp[j] - / get_method(blk, "pressure_sat_comp", j)( - blk, blk.params.get_component(j), Tdew0 * T_units - ) - for j in raoult_comps - ) - + sum( - blk.mole_frac_comp[j] - / blk.params.get_component(j) - .config.henry_component[l_phase]["method"] - .return_expression(blk, l_phase, j, Tdew0 * T_units) - for j in henry_comps - ) - ) - - 1 - ) - df = -value( - blk.pressure - * ( - sum( - blk.mole_frac_comp[j] - / get_method(blk, "pressure_sat_comp", j)( - blk, blk.params.get_component(j), Tdew0 * T_units - ) - ** 2 - * get_method(blk, "pressure_sat_comp", j)( - blk, - blk.params.get_component(j), - Tdew0 * T_units, - dT=True, - ) - for j in raoult_comps - ) - + sum( - blk.mole_frac_comp[j] - / blk.params.get_component(j) - .config.henry_component[l_phase]["method"] - .return_expression(blk, l_phase, j, Tdew0 * T_units) - ** 2 - * blk.params.get_component(j) - .config.henry_component[l_phase]["method"] - .dT_expression(blk, l_phase, j, Tdew0 * T_units) - for j in henry_comps - ) - ) - ) - - # Limit temperature step to avoid excessive overshoot - if f / df < -50: - Tdew1 = Tdew0 + 50 - elif f / df > 50: - Tdew1 = Tdew0 - 50 - else: - Tdew1 = Tdew0 - f / df - - err = abs(Tdew1 - Tdew0) - Tdew0 = Tdew1 - counter += 1 - - blk.temperature_dew[pp].value = Tdew0 - - for j in raoult_comps: - blk._mole_frac_tdew[pp, j].value = value( - blk.mole_frac_comp[j] - * blk.pressure - / get_method(blk, "pressure_sat_comp", j)( - blk, blk.params.get_component(j), Tdew0 * T_units - ) - ) - if blk.is_property_constructed("log_mole_frac_tdew"): - blk.log_mole_frac_tdew[pp, j].value = value( - log(blk._mole_frac_tdew[pp, j]) - ) - for j in henry_comps: - blk._mole_frac_tdew[pp, j].value = value( - blk.mole_frac_comp[j] - * blk.pressure - / blk.params.get_component(j) - .config.henry_component[l_phase]["method"] - .return_expression(blk, l_phase, j, Tdew0 * T_units) - ) - if blk.is_property_constructed("log_mole_frac_tdew"): - blk.log_mole_frac_tdew[pp, j].value = value( - log(blk._mole_frac_tdew[pp, j]) - ) - - def _init_Pbub(self, blk, T_units): - for pp in blk.params._pe_pairs: - raoult_comps, henry_comps = _valid_VL_component_list(blk, pp) - - if raoult_comps == []: - continue - if henry_comps != []: - # Need to get liquid phase name - if blk.params.get_phase(pp[0]).is_liquid_phase(): - l_phase = pp[0] - else: - l_phase = pp[1] - - blk.pressure_bubble[pp].value = value( - sum( - blk.mole_frac_comp[j] * blk.pressure_sat_comp[j] - for j in raoult_comps - ) - + sum( - blk.mole_frac_comp[j] * blk.henry[l_phase, j] for j in henry_comps - ) - ) - - for j in raoult_comps: - blk._mole_frac_pbub[pp, j].value = value( - blk.mole_frac_comp[j] - * blk.pressure_sat_comp[j] - / blk.pressure_bubble[pp] - ) - if blk.is_property_constructed("log_mole_frac_pbub"): - blk.log_mole_frac_pbub[pp, j].value = value( - log(blk._mole_frac_pbub[pp, j]) - ) - for j in henry_comps: - blk._mole_frac_pbub[pp, j].value = value( - blk.mole_frac_comp[j] - * blk.henry[l_phase, j] - / blk.pressure_bubble[pp] - ) - if blk.is_property_constructed("log_mole_frac_pbub"): - blk.log_mole_frac_pbub[pp, j].value = value( - log(blk._mole_frac_pbub[pp, j]) - ) - - def _init_Pdew(self, blk, T_units): - for pp in blk.params._pe_pairs: - raoult_comps, henry_comps = _valid_VL_component_list(blk, pp) - - if raoult_comps == []: - continue - if henry_comps != []: - # Need to get liquid phase name - if blk.params.get_phase(pp[0]).is_liquid_phase(): - l_phase = pp[0] - else: - l_phase = pp[1] - - blk.pressure_dew[pp].value = value( - 1 - / ( - sum( - blk.mole_frac_comp[j] / blk.pressure_sat_comp[j] - for j in raoult_comps - ) - + sum( - blk.mole_frac_comp[j] / blk.henry[l_phase, j] - for j in henry_comps - ) - ) - ) - - for j in raoult_comps: - blk._mole_frac_pdew[pp, j].value = value( - blk.mole_frac_comp[j] - * blk.pressure_dew[pp] - / blk.pressure_sat_comp[j] - ) - if blk.is_property_constructed("log_mole_frac_pdew"): - blk.log_mole_frac_pdew[pp, j].value = value( - log(blk._mole_frac_pdew[pp, j]) - ) - for j in henry_comps: - blk._mole_frac_pdew[pp, j].value = value( - blk.mole_frac_comp[j] * blk.pressure_dew[pp] / blk.henry[l_phase, j] - ) - if blk.is_property_constructed("log_mole_frac_pdew"): - blk.log_mole_frac_pdew[pp, j].value = value( - log(blk._mole_frac_pdew[pp, j]) - ) - @declare_process_block_class("GenericStateBlock", block_class=_GenericStateBlock) class GenericStateBlockData(StateBlockData): @@ -2424,9 +2097,13 @@ def build(self): ): t_units = self.params.get_metadata().default_units.TEMPERATURE + if self.temperature.value is not None: + t_value = value(self.temperature) + else: + t_value = None self._teq = Var( self.params._pe_pairs, - initialize=value(self.temperature), + initialize=t_value, doc="Temperature for calculating phase equilibrium", units=t_units, ) @@ -5121,7 +4798,7 @@ def rule_log_mole_frac(b, p1, p2, j): henry_comps, l_only_comps, v_only_comps, - ) = bub_dew_VL_comp_list(b, (p1, p2)) + ) = identify_VL_component_list(b, (p1, p2)) if l_phase is None or v_phase is None: # Not a VLE pair @@ -5185,3 +4862,146 @@ def _initialize_critical_props(state_data): f"Make sure you have provided values for {prop} in all Component " "declarations." ) + + +def _init_Tbub(blk, T_units): + for pp in blk.params._pe_pairs: + l_phase, _, raoult_comps, henry_comps, _, _ = identify_VL_component_list( + blk, pp + ) + + if raoult_comps == []: + continue + + Tbub0 = estimate_Tbub(blk, T_units, raoult_comps, henry_comps, l_phase) + + blk.temperature_bubble[pp].set_value(Tbub0) + + for j in raoult_comps: + blk._mole_frac_tbub[pp, j].value = value( + blk.mole_frac_comp[j] + * get_method(blk, "pressure_sat_comp", j)( + blk, blk.params.get_component(j), Tbub0 * T_units + ) + / blk.pressure + ) + if blk.is_property_constructed("log_mole_frac_tbub"): + blk.log_mole_frac_tbub[pp, j].value = value( + log(blk._mole_frac_tbub[pp, j]) + ) + + for j in henry_comps: + blk._mole_frac_tbub[pp, j].value = value( + blk.mole_frac_comp[j] + * blk.params.get_component(j) + .config.henry_component[l_phase]["method"] + .return_expression(blk, l_phase, j, Tbub0 * T_units) + / blk.pressure + ) + if blk.is_property_constructed("log_mole_frac_tbub"): + blk.log_mole_frac_tbub[pp, j].value = value( + log(blk._mole_frac_tbub[pp, j]) + ) + + +def _init_Tdew(blk, T_units): + for pp in blk.params._pe_pairs: + l_phase, _, raoult_comps, henry_comps, _, _ = identify_VL_component_list( + blk, pp + ) + + if raoult_comps == []: + continue + + Tdew0 = estimate_Tdew(blk, T_units, raoult_comps, henry_comps, l_phase) + + blk.temperature_dew[pp].set_value(Tdew0) + + for j in raoult_comps: + blk._mole_frac_tdew[pp, j].value = value( + blk.mole_frac_comp[j] + * blk.pressure + / get_method(blk, "pressure_sat_comp", j)( + blk, blk.params.get_component(j), Tdew0 * T_units + ) + ) + if blk.is_property_constructed("log_mole_frac_tdew"): + blk.log_mole_frac_tdew[pp, j].value = value( + log(blk._mole_frac_tdew[pp, j]) + ) + for j in henry_comps: + blk._mole_frac_tdew[pp, j].value = value( + blk.mole_frac_comp[j] + * blk.pressure + / blk.params.get_component(j) + .config.henry_component[l_phase]["method"] + .return_expression(blk, l_phase, j, Tdew0 * T_units) + ) + if blk.is_property_constructed("log_mole_frac_tdew"): + blk.log_mole_frac_tdew[pp, j].value = value( + log(blk._mole_frac_tdew[pp, j]) + ) + + +def _init_Pbub(blk): + for pp in blk.params._pe_pairs: + l_phase, _, raoult_comps, henry_comps, _, _ = identify_VL_component_list( + blk, pp + ) + + if raoult_comps == []: + continue + + blk.pressure_bubble[pp].set_value( + estimate_Pbub(blk, raoult_comps, henry_comps, l_phase) + ) + + for j in raoult_comps: + blk._mole_frac_pbub[pp, j].value = value( + blk.mole_frac_comp[j] + * blk.pressure_sat_comp[j] + / blk.pressure_bubble[pp] + ) + if blk.is_property_constructed("log_mole_frac_pbub"): + blk.log_mole_frac_pbub[pp, j].value = value( + log(blk._mole_frac_pbub[pp, j]) + ) + for j in henry_comps: + blk._mole_frac_pbub[pp, j].value = value( + blk.mole_frac_comp[j] * blk.henry[l_phase, j] / blk.pressure_bubble[pp] + ) + if blk.is_property_constructed("log_mole_frac_pbub"): + blk.log_mole_frac_pbub[pp, j].value = value( + log(blk._mole_frac_pbub[pp, j]) + ) + + +def _init_Pdew(blk): + for pp in blk.params._pe_pairs: + l_phase, _, raoult_comps, henry_comps, _, _ = identify_VL_component_list( + blk, pp + ) + + if raoult_comps == []: + continue + + blk.pressure_dew[pp].set_value( + estimate_Pdew(blk, raoult_comps, henry_comps, l_phase) + ) + + for j in raoult_comps: + blk._mole_frac_pdew[pp, j].value = value( + blk.mole_frac_comp[j] * blk.pressure_dew[pp] / blk.pressure_sat_comp[j] + ) + if blk.is_property_constructed("log_mole_frac_pdew"): + blk.log_mole_frac_pdew[pp, j].value = value( + log(blk._mole_frac_pdew[pp, j]) + ) + for j in henry_comps: + blk._mole_frac_pdew[pp, j].value = value( + blk.mole_frac_comp[j] * blk.pressure_dew[pp] / blk.henry[l_phase, j] + ) + if blk.is_property_constructed("log_mole_frac_pdew"): + blk.log_mole_frac_pdew[pp, j].value = value( + log(blk._mole_frac_pdew[pp, j]) + ) diff --git a/idaes/models/properties/modular_properties/base/tests/test_generic_property.py b/idaes/models/properties/modular_properties/base/tests/test_generic_property.py index 897b3250aa..e43c34fd1e 100644 --- a/idaes/models/properties/modular_properties/base/tests/test_generic_property.py +++ b/idaes/models/properties/modular_properties/base/tests/test_generic_property.py @@ -45,6 +45,9 @@ PhaseType as PT, ) 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 @@ -1972,6 +1975,57 @@ def build_critical_properties(b, *args, **kwargs): _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, + } + }, +} + + @pytest.mark.integration def test_phase_component_flash(): # Regression test for issue #1423 diff --git a/idaes/models/properties/modular_properties/base/tests/test_utility.py b/idaes/models/properties/modular_properties/base/tests/test_utility.py index c2b0034a7f..0bf11fd926 100644 --- a/idaes/models/properties/modular_properties/base/tests/test_utility.py +++ b/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) diff --git a/idaes/models/properties/modular_properties/base/tests/test_vle.py b/idaes/models/properties/modular_properties/base/tests/test_vle.py index 03b3e24555..5431a2aea5 100644 --- a/idaes/models/properties/modular_properties/base/tests/test_vle.py +++ b/idaes/models/properties/modular_properties/base/tests/test_vle.py @@ -34,6 +34,12 @@ ) from idaes.core.solvers import get_solver +from idaes.models.properties.modular_properties.base.generic_property import ( + _init_Pbub, + _init_Pdew, + _init_Tbub, + _init_Tdew, +) from idaes.models.properties.modular_properties.state_definitions import FTPx from idaes.models.properties.modular_properties.eos.ideal import Ideal from idaes.models.properties.modular_properties.phase_equil import SmoothVLE @@ -218,7 +224,7 @@ def test_build(self, model): @pytest.mark.unit def test_init_bubble_temperature(self, model): - model.props._init_Tbub(model.props[1], pyunits.K) + _init_Tbub(model.props[1], pyunits.K) assert pytest.approx(365.35, abs=0.01) == value( model.props[1].temperature_bubble[("Vap", "Liq")] @@ -232,7 +238,7 @@ def test_init_bubble_temperature(self, model): @pytest.mark.unit def test_init_dew_temperature(self, model): - model.props._init_Tdew(model.props[1], pyunits.K) + _init_Tdew(model.props[1], pyunits.K) assert pytest.approx(372.02, abs=0.01) == value( model.props[1].temperature_dew[("Vap", "Liq")] @@ -246,7 +252,7 @@ def test_init_dew_temperature(self, model): @pytest.mark.unit def test_init_bubble_pressure(self, model): - model.props._init_Pbub(model.props[1], pyunits.K) + _init_Pbub(model.props[1]) assert pytest.approx(109479, abs=1) == value( model.props[1].pressure_bubble[("Vap", "Liq")] @@ -260,7 +266,7 @@ def test_init_bubble_pressure(self, model): @pytest.mark.unit def test_init_dew_pressure(self, model): - model.props._init_Pdew(model.props[1], pyunits.K) + _init_Pdew(model.props[1]) assert pytest.approx(89820, abs=1) == value( model.props[1].pressure_dew[("Vap", "Liq")] @@ -524,7 +530,7 @@ def test_build(self, model): @pytest.mark.unit def test_init_bubble_temperature(self, model): - model.props._init_Tbub(model.props[1], pyunits.K) + _init_Tbub(model.props[1], pyunits.K) assert pytest.approx(365.35, abs=0.01) == value( model.props[1].temperature_bubble[("Vap", "Liq")] @@ -538,7 +544,7 @@ def test_init_bubble_temperature(self, model): @pytest.mark.unit def test_init_dew_temperature(self, model): - model.props._init_Tdew(model.props[1], pyunits.K) + _init_Tdew(model.props[1], pyunits.K) assert pytest.approx(372.02, abs=0.01) == value( model.props[1].temperature_dew[("Vap", "Liq")] @@ -552,7 +558,7 @@ def test_init_dew_temperature(self, model): @pytest.mark.unit def test_init_bubble_pressure(self, model): - model.props._init_Pbub(model.props[1], pyunits.K) + _init_Pbub(model.props[1]) assert pytest.approx(109479, abs=1) == value( model.props[1].pressure_bubble[("Vap", "Liq")] @@ -566,7 +572,7 @@ def test_init_bubble_pressure(self, model): @pytest.mark.unit def test_init_dew_pressure(self, model): - model.props._init_Pdew(model.props[1], pyunits.K) + _init_Pdew(model.props[1]) assert pytest.approx(89820, abs=1) == value( model.props[1].pressure_dew[("Vap", "Liq")] @@ -703,7 +709,7 @@ def test_build(self, model): @pytest.mark.unit def test_init_bubble_temperature(self, model): - model.props._init_Tbub(model.props[1], pyunits.K) + _init_Tbub(model.props[1], pyunits.K) assert pytest.approx(361.50, abs=0.01) == value( model.props[1].temperature_bubble[("Vap", "Liq")] @@ -720,7 +726,7 @@ def test_init_bubble_temperature(self, model): @pytest.mark.unit def test_init_dew_temperature(self, model): - model.props._init_Tdew(model.props[1], pyunits.K) + _init_Tdew(model.props[1], pyunits.K) assert pytest.approx(370.23, abs=0.01) == value( model.props[1].temperature_dew[("Vap", "Liq")] @@ -737,7 +743,7 @@ def test_init_dew_temperature(self, model): @pytest.mark.unit def test_init_bubble_pressure(self, model): - model.props._init_Pbub(model.props[1], pyunits.K) + _init_Pbub(model.props[1]) assert pytest.approx(118531, abs=1) == value( model.props[1].pressure_bubble[("Vap", "Liq")] @@ -754,7 +760,7 @@ def test_init_bubble_pressure(self, model): @pytest.mark.unit def test_init_dew_pressure(self, model): - model.props._init_Pdew(model.props[1], pyunits.K) + _init_Pdew(model.props[1]) assert pytest.approx(95056, abs=1) == value( model.props[1].pressure_dew[("Vap", "Liq")] diff --git a/idaes/models/properties/modular_properties/base/utility.py b/idaes/models/properties/modular_properties/base/utility.py index 8588de80fa..5dcaf9eeae 100644 --- a/idaes/models/properties/modular_properties/base/utility.py +++ b/idaes/models/properties/modular_properties/base/utility.py @@ -19,12 +19,9 @@ # pylint: disable=missing-class-docstring # pylint: disable=missing-function-docstring -# TODO: Look into protected access issues -# pylint: disable=protected-access - from enum import Enum -from pyomo.environ import units as pyunits +from pyomo.environ import units as pyunits, value from idaes.core.util.exceptions import ( BurntToast, @@ -43,7 +40,10 @@ class StateIndex(Enum): class GenericPropertyPackageError(PropertyPackageError): - # Error message for when a property is called for but no option provided + """ + Error message for when a property is called for but no option provided + """ + def __init__(self, block, prop): super().__init__() self.prop = prop @@ -241,6 +241,19 @@ class ConcentrationForm(Enum): def get_concentration_term(blk, r_idx, log=False): + """ + Get necessary concentration terms for reactions from property package, allowing for + different bases. + + Args: + blk: StateBlock of interest + r_idx: index of reaction + log: whether to use log concentration of not + + Returns: + Var or Expression representing concentration + + """ cfg = blk.params.config if "rate_reactions" in cfg: try: @@ -252,6 +265,7 @@ def get_concentration_term(blk, r_idx, log=False): conc_form = cfg.inherent_reactions[r_idx].concentration_form state = blk + # pylint: disable-next=protected-access if hasattr(state.params, "_electrolyte") and state.params._electrolyte: sub = "_true" else: @@ -288,3 +302,308 @@ def get_concentration_term(blk, r_idx, log=False): ) return conc_term + + +def identify_VL_component_list(blk, phase_pair): + """ + Identify liquid and vapor phases and which components are in VL equilibrium + + Args: + blk: StateBlock of interest + phase_pair: 2-tuple of Phases in equilibrium + + Returns: + Lists of component names for: + * liquid Phase object + * vapor Phase object + * components using Raoult's Law + * components using Henry's Law + * liquid only components, + * vapor only components + + """ + vl_comps = [] + henry_comps = [] + l_only_comps = [] + v_only_comps = [] + + pparams = blk.params + + if pparams.get_phase(phase_pair[0]).is_liquid_phase(): + l_phase = phase_pair[0] + if pparams.get_phase(phase_pair[1]).is_vapor_phase(): + v_phase = phase_pair[1] + else: + raise PropertyPackageError( + f"Phase pair {phase_pair[0]}-{phase_pair[1]} was identified as " + f"being a VLE pair, however {phase_pair[0]} is liquid but " + f"{phase_pair[1]} is not vapor." + ) + elif pparams.get_phase(phase_pair[1]).is_liquid_phase(): + l_phase = phase_pair[1] + if pparams.get_phase(phase_pair[0]).is_vapor_phase(): + v_phase = phase_pair[0] + else: + raise PropertyPackageError( + f"Phase pair {phase_pair[0]}-{phase_pair[1]} was identified as " + f"being a VLE pair, however {phase_pair[1]} is liquid but " + f"{phase_pair[0]} is not vapor." + ) + else: + raise PropertyPackageError( + f"Phase pair {phase_pair[0]}-{phase_pair[1]} was identified as " + f"being a VLE pair, however neither {phase_pair[0]} nor " + f"{phase_pair[1]} is liquid." + ) + + for j in blk.params.component_list: + if (l_phase, j) in blk.phase_component_set and ( + v_phase, + j, + ) in blk.phase_component_set: + cobj = pparams.get_component(j) + if cobj.config.henry_component is not None and ( + phase_pair[0] in cobj.config.henry_component + or phase_pair[1] in cobj.config.henry_component + ): + henry_comps.append(j) + else: + vl_comps.append(j) + elif (l_phase, j) in blk.phase_component_set: + l_only_comps.append(j) + elif (v_phase, j) in blk.phase_component_set: + v_only_comps.append(j) + + return l_phase, v_phase, vl_comps, henry_comps, l_only_comps, v_only_comps + + +TOL = 1e-1 +MAX_ITER = 30 + + +def estimate_Tbub(blk, T_units, raoult_comps, henry_comps, liquid_phase): + """ + Function to estimate bubble point temperature + + Args: + blk: StateBlock to use + T_units: units of temperature + raoult_comps: list of components that follow Raoult's Law + henry_comps: list of components that follow Henry's Law + liquid_phase: name of liquid phase + + Returns: + Estimated bubble point temperature as a float. + + """ + # Use lowest component temperature_crit as starting point + # Starting high and moving down generally works better, + # as it under-predicts next step due to exponential form of + # Psat. + # Subtract 1 to avoid potential singularities at Tcrit + Tbub0 = ( + min(blk.params.get_component(j).temperature_crit.value for j in raoult_comps) + - 1 + ) + + err = 1 + counter = 0 + + # Newton solver with step limiter to prevent overshoot + # Tolerance only needs to be ~1e-1 + # Iteration limit of 30 + while err > TOL and counter < MAX_ITER: + f = value( + sum( + get_method(blk, "pressure_sat_comp", j)( + blk, blk.params.get_component(j), Tbub0 * T_units + ) + * blk.mole_frac_comp[j] + for j in raoult_comps + ) + + sum( + blk.mole_frac_comp[j] + * blk.params.get_component(j) + .config.henry_component[liquid_phase]["method"] + .return_expression(blk, liquid_phase, j, Tbub0 * T_units) + for j in henry_comps + ) + - blk.pressure + ) + df = value( + sum( + get_method(blk, "pressure_sat_comp", j)( + blk, blk.params.get_component(j), Tbub0 * T_units, dT=True + ) + * blk.mole_frac_comp[j] + for j in raoult_comps + ) + + sum( + blk.mole_frac_comp[j] + * blk.params.get_component(j) + .config.henry_component[liquid_phase]["method"] + .dT_expression(blk, liquid_phase, j, Tbub0 * T_units) + for j in henry_comps + ) + ) + + # Limit temperature step to avoid excessive overshoot + if f / df < -50: + Tbub1 = Tbub0 + 50 + elif f / df > 50: + Tbub1 = Tbub0 - 50 + else: + Tbub1 = Tbub0 - f / df + + err = abs(Tbub1 - Tbub0) + Tbub0 = Tbub1 + counter += 1 + + return Tbub0 + + +def estimate_Tdew(blk, T_units, raoult_comps, henry_comps, liquid_phase): + """ + Function to estimate dew point temperature + + Args: + blk: StateBlock to use + T_units: units of temperature + raoult_comps: list of components that follow Raoult's Law + henry_comps: list of components that follow Henry's Law + liquid_phase: name of liquid phase + + Returns: + Estimated dew point temperature as a float. + + """ + # Use lowest component critical temperature + # as starting point + # Subtract 1 to avoid potential singularities at Tcrit + Tdew0 = ( + min(blk.params.get_component(j).temperature_crit.value for j in raoult_comps) + - 1 + ) + + err = 1 + counter = 0 + + # Newton solver with step limiter to prevent overshoot + # Tolerance only needs to be ~1e-1 + # Iteration limit of 30 + while err > TOL and counter < MAX_ITER: + f = value( + blk.pressure + * ( + sum( + blk.mole_frac_comp[j] + / get_method(blk, "pressure_sat_comp", j)( + blk, blk.params.get_component(j), Tdew0 * T_units + ) + for j in raoult_comps + ) + + sum( + blk.mole_frac_comp[j] + / blk.params.get_component(j) + .config.henry_component[liquid_phase]["method"] + .return_expression(blk, liquid_phase, j, Tdew0 * T_units) + for j in henry_comps + ) + ) + - 1 + ) + df = -value( + blk.pressure + * ( + sum( + blk.mole_frac_comp[j] + / get_method(blk, "pressure_sat_comp", j)( + blk, blk.params.get_component(j), Tdew0 * T_units + ) + ** 2 + * get_method(blk, "pressure_sat_comp", j)( + blk, + blk.params.get_component(j), + Tdew0 * T_units, + dT=True, + ) + for j in raoult_comps + ) + + sum( + blk.mole_frac_comp[j] + / blk.params.get_component(j) + .config.henry_component[liquid_phase]["method"] + .return_expression(blk, liquid_phase, j, Tdew0 * T_units) + ** 2 + * blk.params.get_component(j) + .config.henry_component[liquid_phase]["method"] + .dT_expression(blk, liquid_phase, j, Tdew0 * T_units) + for j in henry_comps + ) + ) + ) + + # Limit temperature step to avoid excessive overshoot + if f / df < -50: + Tdew1 = Tdew0 + 50 + elif f / df > 50: + Tdew1 = Tdew0 - 50 + else: + Tdew1 = Tdew0 - f / df + + err = abs(Tdew1 - Tdew0) + Tdew0 = Tdew1 + counter += 1 + + return Tdew0 + + +def estimate_Pbub(blk, raoult_comps, henry_comps, liquid_phase): + """ + Function to estimate bubble point pressure + + Args: + blk: StateBlock to use + raoult_comps: list of components that follow Raoult's Law + henry_comps: list of components that follow Henry's Law + liquid_phase: name of liquid phase + + Returns: + Estimated bubble point pressure as a float. + + """ + return value( + sum(blk.mole_frac_comp[j] * blk.pressure_sat_comp[j] for j in raoult_comps) + + sum(blk.mole_frac_comp[j] * blk.henry[liquid_phase, j] for j in henry_comps) + ) + + +def estimate_Pdew(blk, raoult_comps, henry_comps, liquid_phase): + """ + Function to estimate dew point pressure + + Args: + blk: StateBlock to use + raoult_comps: list of components that follow Raoult's Law + henry_comps: list of components that follow Henry's Law + liquid_phase: name of liquid phase + + Returns: + Estimated dew point pressure as a float. + + """ + # Safety catch for cases where Psat or Henry's constant might be 0 + # Not sure if this is meaningful, but if this is true then mathematically Pdew = 0 + if any(value(blk.pressure_sat_comp[j]) == 0 for j in raoult_comps) or any( + value(blk.henry[liquid_phase, j]) == 0 for j in henry_comps + ): + return 0 + return value( + 1 + / ( + sum(blk.mole_frac_comp[j] / blk.pressure_sat_comp[j] for j in raoult_comps) + + sum( + blk.mole_frac_comp[j] / blk.henry[liquid_phase, j] for j in henry_comps + ) + ) + ) diff --git a/idaes/models/properties/modular_properties/eos/ceos.py b/idaes/models/properties/modular_properties/eos/ceos.py index d419ef149f..4288fdc680 100644 --- a/idaes/models/properties/modular_properties/eos/ceos.py +++ b/idaes/models/properties/modular_properties/eos/ceos.py @@ -15,11 +15,9 @@ Currently only supports liquid and vapor phases """ -# TODO: Missing docstrings -# pylint: disable=missing-function-docstring - -# TODO: Look into protected access issues +# TODO: Pylint complains about variables with _x names as they are built by other classes # pylint: disable=protected-access +# pylint: disable=missing-function-docstring from enum import Enum from copy import deepcopy @@ -543,6 +541,9 @@ def build_parameters(b): @staticmethod def compress_fact_phase(b, p): + """ + Compressibility factor + """ pobj = b.params.get_phase(p) cname = pobj._cubic_type.name A = getattr(b, cname + "_A") @@ -562,10 +563,16 @@ def compress_fact_phase(b, p): @staticmethod def cp_mass_phase(blk, p): + """ + Phase mass-specific heat capacity at constant pressure + """ return blk.cp_mol_phase[p] / blk.mw_phase[p] @staticmethod def cp_mol_phase(blk, p): + """ + Phase molar heat capacity at constant pressure + """ pobj = blk.params.get_phase(p) cname = pobj._cubic_type.name @@ -607,10 +614,16 @@ def cp_mol_phase(blk, p): @staticmethod def cv_mass_phase(blk, p): + """ + Phase mass-specific heat capacity at constant volume + """ return blk.cv_mol_phase[p] / blk.mw_phase[p] @staticmethod def cv_mol_phase(blk, p): + """ + Phase molar heat capacity at constant volume + """ pobj = blk.params.get_phase(p) cname = pobj._cubic_type.name am = getattr(blk, cname + "_am")[p] @@ -635,16 +648,25 @@ def cv_mol_phase(blk, p): @staticmethod def dens_mass_phase(b, p): + """ + Phase density (mass basis) + """ return b.dens_mol_phase[p] * b.mw_phase[p] @staticmethod def dens_mol_phase(b, p): + """ + Phase density (mole basis) + """ return b.pressure / ( Cubic.gas_constant(b) * b.temperature * b.compress_fact_phase[p] ) @staticmethod def energy_internal_mol_phase(blk, p): + """ + Phase molar internal energy + """ pobj = blk.params.get_phase(p) cname = pobj._cubic_type.name @@ -674,12 +696,18 @@ def energy_internal_mol_phase(blk, p): @staticmethod def energy_internal_mol_phase_comp(blk, p, j): + """ + Phase partial molar internal energy + """ return ( blk.enth_mol_phase_comp[p, j] - blk.pressure * blk.vol_mol_phase_comp[p, j] ) @staticmethod def enth_mol_phase(blk, p): + """ + Phase molar enthalpy + """ pobj = blk.params.get_phase(p) cname = pobj._cubic_type.name @@ -709,6 +737,9 @@ def enth_mol_phase(blk, p): @staticmethod def enth_mol_phase_comp(blk, p, j): + """ + Phase partial molar enthalpy + """ dlogphi_j_dT = _d_log_fug_coeff_dT_phase_comp(blk, p, j) enth_ideal_gas = get_method(blk, "enth_mol_ig_comp", j)( @@ -721,6 +752,9 @@ def enth_mol_phase_comp(blk, p, j): @staticmethod def entr_mol_phase(blk, p): + """ + Phase molar entropy + """ pobj = blk.params.get_phase(p) cname = pobj._cubic_type.name @@ -759,6 +793,9 @@ def entr_mol_phase(blk, p): @staticmethod def entr_mol_phase_comp(blk, p, j): + """ + Phase partial molar entropy + """ logphi_j = _log_fug_coeff_phase_comp(blk, p, j) dlogphi_j_dT = _d_log_fug_coeff_dT_phase_comp(blk, p, j) @@ -1355,27 +1392,109 @@ def a(k): # ----------------------------------------------------------------------------- # Default rules for cubic expressions +def calculate_equilibrium_cubic_coefficients(b, cubic_name, cubic_type, p1, p2, p3): + """ + Calculates the coefficients b, c, and d of the cubic 0 = z**3 + b * z**2 + c * z + d + at the equilibrium conditions + + Args: + b: StateBlock of interest + cubic_name: Name of Cubic EoS + cubic_type: Type of Cubic EoS + p1: Phase 1 + p2: Phase 2 + p3; Phase 3 + + Returns: + expressions for b, c, d + """ + + A_eq = getattr(b, "_" + cubic_name + "_A_eq")[p1, p2, p3] + B_eq = getattr(b, "_" + cubic_name + "_B_eq")[p1, p2, p3] + EoS_u = EoS_param[cubic_type]["u"] + EoS_w = EoS_param[cubic_type]["w"] + + b = -(1 + B_eq - EoS_u * B_eq) + c = A_eq - EoS_u * B_eq - EoS_u * B_eq**2 + EoS_w * B_eq**2 + d = -(A_eq * B_eq + EoS_w * B_eq**2 + EoS_w * B_eq**3) + + return b, c, d + + def func_fw_PR(cobj): + """ + f(omega) function for Peng-Robinson EoS. + + Args: + cobj: Component object + + Returns: + expression for fw + + """ return 0.37464 + 1.54226 * cobj.omega - 0.26992 * cobj.omega**2 def func_fw_SRK(cobj): + """ + f(omega) function for SRK EoS. + + Args: + cobj: Component object + + Returns: + expression for fw + + """ return 0.48 + 1.574 * cobj.omega - 0.176 * cobj.omega**2 def func_alpha_soave(T, fw, cobj): + """ + Soave alpha function. + + Args: + fw: expression for fw + cobj: Component object + + Returns: + expression for alpha + + """ Tc = cobj.temperature_crit Tr = T / Tc return (1 + fw * (1 - sqrt(Tr))) ** 2 def func_dalpha_dT_soave(T, fw, cobj): + """ + Function to get first partial derivative w.r.t. temperature of Soave alpha function. + + Args: + fw: expression for fw + cobj: Component object + + Returns: + expression for first derivative of alpha + + """ Tc = cobj.temperature_crit Tr = T / Tc return 1 / Tc * (-fw / sqrt(Tr)) * (1 + fw * (1 - sqrt(Tr))) def func_d2alpha_dT2_soave(T, fw, cobj): + """ + Function to get 2nd partial derivative w.r.t. temperature of Soave alpha function. + + Args: + fw: expression for fw + cobj: Component object + + Returns: + expression for 2nd derivative of alpha + + """ Tc = cobj.temperature_crit Tr = T / Tc return 1 / Tc**2 * ((fw**2 + fw) / (2 * Tr * sqrt(Tr))) @@ -1384,6 +1503,9 @@ def func_d2alpha_dT2_soave(T, fw, cobj): # ----------------------------------------------------------------------------- # Mixing rules def rule_am_default(m, cname, a, p, pp=()): + """ + Standard Van der Waals one-fluid mixing rule for a term + """ k = getattr(m.params, cname + "_kappa") return sum( sum( @@ -1398,6 +1520,9 @@ def rule_am_default(m, cname, a, p, pp=()): def rule_am_crit_default(m, cname, a_crit): + """ + Standard Van der Waals one-fluid mixing rule for a term evaluated at critical point + """ k = getattr(m.params, cname + "_kappa") return sum( sum( @@ -1412,8 +1537,14 @@ def rule_am_crit_default(m, cname, a_crit): def rule_bm_default(m, b, p): + """ + Standard Van der Waals one-fluid mixing rule for b term + """ return sum(m.mole_frac_phase_comp[p, i] * b[i] for i in m.components_in_phase(p)) def rule_bm_crit_default(m, b): + """ + Standard Van der Waals one-fluid mixing rule for b term evaluated at critical point + """ return sum(m.mole_frac_comp[i] * b[i] for i in m.component_list) diff --git a/idaes/models/properties/modular_properties/eos/tests/test_ceos_PR.py b/idaes/models/properties/modular_properties/eos/tests/test_ceos_PR.py index 2240c0160e..6509009d44 100644 --- a/idaes/models/properties/modular_properties/eos/tests/test_ceos_PR.py +++ b/idaes/models/properties/modular_properties/eos/tests/test_ceos_PR.py @@ -41,7 +41,11 @@ ) from idaes.core.util.exceptions import PropertyNotSupportedError, ConfigurationError from idaes.core.util.constants import Constants as const -from idaes.models.properties.modular_properties.eos.ceos import cubic_roots_available +from idaes.models.properties.modular_properties.eos.ceos import ( + cubic_roots_available, + calculate_equilibrium_cubic_coefficients, + EoS_param, +) from idaes.core.solvers import get_solver from idaes.core.initialization.initializer_base import InitializationStatus from idaes.models.properties.modular_properties.state_definitions import FTPx @@ -1076,6 +1080,23 @@ def test_vol_mol_phase_comp(m): ) +@pytest.mark.skipif(not cubic_roots_available(), reason="Cubic functions not available") +@pytest.mark.unit +def test_calculate_equilibrium_cubic_coefficients(m): + b, c, d = calculate_equilibrium_cubic_coefficients( + m.props[1], "PR", CubicType.PR, "Vap", "Liq", "Liq" + ) + + A_eq = m.props[1]._PR_A_eq["Vap", "Liq", "Liq"] + B_eq = m.props[1]._PR_B_eq["Vap", "Liq", "Liq"] + EoS_u = EoS_param[CubicType.PR]["u"] + EoS_w = EoS_param[CubicType.PR]["w"] + + assert str(b) == str(-(1 + B_eq - EoS_u * B_eq)) + assert str(c) == str(A_eq - EoS_u * B_eq - EoS_u * B_eq**2 + EoS_w * B_eq**2) + assert str(d) == str(-(A_eq * B_eq + EoS_w * B_eq**2 + EoS_w * B_eq**3)) + + class TestCEOSCriticalProps: @pytest.fixture(scope="class") def model(self): diff --git a/idaes/models/properties/modular_properties/examples/BT_PR.py b/idaes/models/properties/modular_properties/examples/BT_PR.py index ece424ecb9..7341a3992e 100644 --- a/idaes/models/properties/modular_properties/examples/BT_PR.py +++ b/idaes/models/properties/modular_properties/examples/BT_PR.py @@ -31,7 +31,9 @@ from idaes.models.properties.modular_properties.state_definitions import FTPx 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 import ( + CubicComplementarityVLE, +) from idaes.models.properties.modular_properties.phase_equil.bubble_dew import ( LogBubbleDew, ) @@ -148,7 +150,7 @@ "temperature_ref": (298.15, pyunits.K), # Defining phase equilibria "phases_in_equilibrium": [("Vap", "Liq")], - "phase_equilibrium_state": {("Vap", "Liq"): SmoothVLE}, + "phase_equilibrium_state": {("Vap", "Liq"): CubicComplementarityVLE}, "bubble_dew_method": LogBubbleDew, "parameter_data": { "PR_kappa": { diff --git a/idaes/models/properties/modular_properties/examples/tests/test_BT_PR.py b/idaes/models/properties/modular_properties/examples/tests/test_BT_PR.py index 95955efc94..ad9ad9afa7 100644 --- a/idaes/models/properties/modular_properties/examples/tests/test_BT_PR.py +++ b/idaes/models/properties/modular_properties/examples/tests/test_BT_PR.py @@ -14,19 +14,19 @@ Author: Andrew Lee """ - import pytest +from numpy import logspace + +from pyomo.util.check_units import assert_units_consistent +from pyomo.environ import assert_optimal_termination, ConcreteModel, Objective, value + from idaes.core import FlowsheetBlock from idaes.models.properties.modular_properties.eos.ceos import cubic_roots_available from idaes.models.properties.modular_properties.examples.BT_PR import configuration from idaes.models.properties.modular_properties.base.generic_property import ( GenericParameterBlock, ) -from pyomo.util.check_units import assert_units_consistent - -from pyomo.environ import check_optimal_termination, ConcreteModel, Objective, value - from idaes.core.solvers import get_solver import idaes.core.util.scaling as iscale from idaes.models.properties.tests.test_harness import PropertyTestHarness @@ -42,7 +42,7 @@ # ----------------------------------------------------------------------------- # Get default solver for testing solver = get_solver() -# Limit iterations to make sure sweeps aren;t getting out of hand +# Limit iterations to make sure sweeps aren't getting out of hand solver.options["max_iter"] = 50 @@ -69,6 +69,11 @@ def m(self): m.fs.state = m.fs.props.build_state_block([1], defined_state=True) + # Set small values of epsilon to get accurate results + # Initialization will handle finding the correct region + m.fs.state[1].eps_t_Vap_Liq.set_value(1e-4) + m.fs.state[1].eps_z_Vap_Liq.set_value(1e-4) + iscale.calculate_scaling_factors(m.fs.props) iscale.calculate_scaling_factors(m.fs.state[1]) return m @@ -80,23 +85,34 @@ def test_T_sweep(self, m): m.fs.obj = Objective(expr=(m.fs.state[1].temperature - 510) ** 2) m.fs.state[1].temperature.setub(600) - for logP in [9.5, 10, 10.5, 11, 11.5, 12]: - m.fs.obj.deactivate() + for P in logspace(4.8, 5.9, 8): m.fs.state[1].flow_mol.fix(100) m.fs.state[1].mole_frac_comp["benzene"].fix(0.5) m.fs.state[1].mole_frac_comp["toluene"].fix(0.5) m.fs.state[1].temperature.fix(300) - m.fs.state[1].pressure.fix(10 ** (0.5 * logP)) + m.fs.state[1].pressure.fix(P) + + # For optimization sweep, use a large eps to avoid getting stuck at + # bubble and dew points + m.fs.state[1].eps_t_Vap_Liq.set_value(10) + m.fs.state[1].eps_z_Vap_Liq.set_value(10) m.fs.state.initialize() m.fs.state[1].temperature.unfix() m.fs.obj.activate() + results = solver.solve(m) + assert_optimal_termination(results) + + # Switch to small eps and re-solve to refine result + m.fs.state[1].eps_t_Vap_Liq.set_value(1e-4) + m.fs.state[1].eps_z_Vap_Liq.set_value(1e-4) + results = solver.solve(m) - assert check_optimal_termination(results) + assert_optimal_termination(results) assert m.fs.state[1].flow_mol_phase["Liq"].value <= 1e-5 @pytest.mark.integration @@ -112,12 +128,12 @@ def test_P_sweep(self, m): results = solver.solve(m) - assert check_optimal_termination(results) + assert_optimal_termination(results) while m.fs.state[1].pressure.value <= 1e6: results = solver.solve(m) - assert check_optimal_termination(results) + assert_optimal_termination(results) m.fs.state[1].pressure.value = m.fs.state[1].pressure.value + 1e5 @@ -138,7 +154,7 @@ def test_T350_P1_x5(self, m): results = solver.solve(m) # Check for optimal solution - assert check_optimal_termination(results) + assert_optimal_termination(results) assert pytest.approx(value(m.fs.state[1]._teq[("Vap", "Liq")]), abs=1e-1) == 365 assert 0.0035346 == pytest.approx( @@ -227,7 +243,7 @@ def test_T350_P5_x5(self, m): results = solver.solve(m) # Check for optimal solution - assert check_optimal_termination(results) + assert_optimal_termination(results) assert pytest.approx(value(m.fs.state[1]._teq[("Vap", "Liq")]), 1e-5) == 431.47 assert ( @@ -318,7 +334,7 @@ def test_T450_P1_x5(self, m): results = solver.solve(m) # Check for optimal solution - assert check_optimal_termination(results) + assert_optimal_termination(results) assert pytest.approx(value(m.fs.state[1]._teq[("Vap", "Liq")]), 1e-5) == 371.4 assert 0.0033583 == pytest.approx( @@ -407,7 +423,7 @@ def test_T450_P5_x5(self, m): results = solver.solve(m) # Check for optimal solution - assert check_optimal_termination(results) + assert_optimal_termination(results) assert pytest.approx(value(m.fs.state[1]._teq[("Vap", "Liq")]), 1e-5) == 436.93 assert 0.0166181 == pytest.approx( @@ -496,7 +512,7 @@ def test_T368_P1_x5(self, m): results = solver.solve(m) # Check for optimal solution - assert check_optimal_termination(results) + assert_optimal_termination(results) assert pytest.approx(value(m.fs.state[1]._teq[("Vap", "Liq")]), 1e-5) == 368 assert 0.003504 == pytest.approx( @@ -589,7 +605,7 @@ def test_T376_P1_x2(self, m): results = solver.solve(m) # Check for optimal solution - assert check_optimal_termination(results) + assert_optimal_termination(results) assert pytest.approx(value(m.fs.state[1]._teq[("Vap", "Liq")]), 1e-5) == 376 assert 0.00361333 == pytest.approx( @@ -660,103 +676,3 @@ def test_T376_P1_x2(self, m): assert ( pytest.approx(value(m.fs.state[1].entr_mol_phase["Vap"]), 1e-5) == -273.513 ) - - @pytest.mark.unit - def test_basic_scaling(self, m): - assert len(m.fs.state[1].scaling_factor) == 23 - assert m.fs.state[1].scaling_factor[m.fs.state[1].flow_mol] == 1e-2 - assert m.fs.state[1].scaling_factor[m.fs.state[1].flow_mol_phase["Liq"]] == 1e-2 - assert m.fs.state[1].scaling_factor[m.fs.state[1].flow_mol_phase["Vap"]] == 1e-2 - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1].flow_mol_phase_comp["Liq", "benzene"] - ] - == 1e-2 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1].flow_mol_phase_comp["Liq", "toluene"] - ] - == 1e-2 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1].flow_mol_phase_comp["Vap", "benzene"] - ] - == 1e-2 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1].flow_mol_phase_comp["Vap", "toluene"] - ] - == 1e-2 - ) - assert ( - m.fs.state[1].scaling_factor[m.fs.state[1].mole_frac_comp["benzene"]] - == 1000 - ) - assert ( - m.fs.state[1].scaling_factor[m.fs.state[1].mole_frac_comp["toluene"]] - == 1000 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1].mole_frac_phase_comp["Liq", "benzene"] - ] - == 1000 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1].mole_frac_phase_comp["Liq", "toluene"] - ] - == 1000 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1].mole_frac_phase_comp["Vap", "benzene"] - ] - == 1000 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1].mole_frac_phase_comp["Vap", "toluene"] - ] - == 1000 - ) - assert m.fs.state[1].scaling_factor[m.fs.state[1].pressure] == 1e-5 - assert m.fs.state[1].scaling_factor[m.fs.state[1].temperature] == 1e-2 - assert m.fs.state[1].scaling_factor[m.fs.state[1]._teq["Vap", "Liq"]] == 1e-2 - assert m.fs.state[1].scaling_factor[m.fs.state[1]._t1_Vap_Liq] == 1e-2 - - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1]._mole_frac_tbub["Vap", "Liq", "benzene"] - ] - == 1000 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1]._mole_frac_tbub["Vap", "Liq", "toluene"] - ] - == 1000 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1]._mole_frac_tdew["Vap", "Liq", "benzene"] - ] - == 1000 - ) - assert ( - m.fs.state[1].scaling_factor[ - m.fs.state[1]._mole_frac_tdew["Vap", "Liq", "toluene"] - ] - == 1000 - ) - assert ( - m.fs.state[1].scaling_factor[m.fs.state[1].temperature_bubble["Vap", "Liq"]] - == 1e-2 - ) - assert ( - m.fs.state[1].scaling_factor[m.fs.state[1].temperature_dew["Vap", "Liq"]] - == 1e-2 - ) diff --git a/idaes/models/properties/modular_properties/examples/tests/test_BT_PR_legacy_SmoothVLE.py b/idaes/models/properties/modular_properties/examples/tests/test_BT_PR_legacy_SmoothVLE.py new file mode 100644 index 0000000000..f25035437a --- /dev/null +++ b/idaes/models/properties/modular_properties/examples/tests/test_BT_PR_legacy_SmoothVLE.py @@ -0,0 +1,800 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2023 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# +""" +Author: Andrew Lee +""" + +import pytest + +from numpy import logspace + +from pyomo.util.check_units import assert_units_consistent +from pyomo.environ import ( + check_optimal_termination, + ConcreteModel, + Objective, + units as pyunits, + value, +) + +from idaes.core import FlowsheetBlock +from idaes.models.properties.modular_properties.eos.ceos import cubic_roots_available +from idaes.models.properties.modular_properties.base.generic_property import ( + GenericParameterBlock, +) +from idaes.core.solvers import get_solver +import idaes.core.util.scaling as iscale +from idaes.models.properties.tests.test_harness import PropertyTestHarness +from idaes.core import LiquidPhase, VaporPhase, Component +from idaes.models.properties.modular_properties.state_definitions import FTPx +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 RPP4 + +import idaes.logger as idaeslog + +SOUT = idaeslog.INFO + +# Set module level pyest marker +pytestmark = pytest.mark.cubic_root + + +# ----------------------------------------------------------------------------- +# Get default solver for testing +solver = get_solver() +# Limit iterations to make sure sweeps aren't getting out of hand +solver.options["max_iter"] = 50 + +# --------------------------------------------------------------------- +# Configuration dictionary for an ideal Benzene-Toluene system + +# Data Sources: +# [1] The Properties of Gases and Liquids (1987) +# 4th edition, Chemical Engineering Series - Robert C. Reid +# [3] Engineering Toolbox, https://www.engineeringtoolbox.com +# Retrieved 1st December, 2019 + +configuration = { + # Specifying components + "components": { + "benzene": { + "type": Component, + "enth_mol_ig_comp": RPP4, + "entr_mol_ig_comp": RPP4, + "pressure_sat_comp": RPP4, + "phase_equilibrium_form": {("Vap", "Liq"): log_fugacity}, + "parameter_data": { + "mw": (78.1136e-3, pyunits.kg / pyunits.mol), # [1] + "pressure_crit": (48.9e5, pyunits.Pa), # [1] + "temperature_crit": (562.2, pyunits.K), # [1] + "omega": 0.212, # [1] + "cp_mol_ig_comp_coeff": { + "A": (-3.392e1, pyunits.J / pyunits.mol / pyunits.K), # [1] + "B": (4.739e-1, pyunits.J / pyunits.mol / pyunits.K**2), + "C": (-3.017e-4, pyunits.J / pyunits.mol / pyunits.K**3), + "D": (7.130e-8, pyunits.J / pyunits.mol / pyunits.K**4), + }, + "enth_mol_form_vap_comp_ref": (82.9e3, pyunits.J / pyunits.mol), # [3] + "entr_mol_form_vap_comp_ref": ( + -269, + pyunits.J / pyunits.mol / pyunits.K, + ), # [3] + "pressure_sat_comp_coeff": { + "A": (-6.98273, None), # [1] + "B": (1.33213, None), + "C": (-2.62863, None), + "D": (-3.33399, None), + }, + }, + }, + "toluene": { + "type": Component, + "enth_mol_ig_comp": RPP4, + "entr_mol_ig_comp": RPP4, + "pressure_sat_comp": RPP4, + "phase_equilibrium_form": {("Vap", "Liq"): log_fugacity}, + "parameter_data": { + "mw": (92.1405e-3, pyunits.kg / pyunits.mol), # [1] + "pressure_crit": (41e5, pyunits.Pa), # [1] + "temperature_crit": (591.8, pyunits.K), # [1] + "omega": 0.263, # [1] + "cp_mol_ig_comp_coeff": { + "A": (-2.435e1, pyunits.J / pyunits.mol / pyunits.K), # [1] + "B": (5.125e-1, pyunits.J / pyunits.mol / pyunits.K**2), + "C": (-2.765e-4, pyunits.J / pyunits.mol / pyunits.K**3), + "D": (4.911e-8, pyunits.J / pyunits.mol / pyunits.K**4), + }, + "enth_mol_form_vap_comp_ref": (50.1e3, pyunits.J / pyunits.mol), # [3] + "entr_mol_form_vap_comp_ref": ( + -321, + pyunits.J / pyunits.mol / pyunits.K, + ), # [3] + "pressure_sat_comp_coeff": { + "A": (-7.28607, None), # [1] + "B": (1.38091, None), + "C": (-2.83433, None), + "D": (-2.79168, None), + }, + }, + }, + }, + # 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": { + ("benzene", "benzene"): 0.000, + ("benzene", "toluene"): 0.000, + ("toluene", "benzene"): 0.000, + ("toluene", "toluene"): 0.000, + } + }, +} + + +@pytest.mark.skipif(not cubic_roots_available(), reason="Cubic functions not available") +class TestBTPR(PropertyTestHarness): + def configure(self): + self.prop_pack = GenericParameterBlock + self.param_args = configuration + self.prop_args = {} + self.has_density_terms = False + + +# ----------------------------------------------------------------------------- +# Test robustness and some outputs +@pytest.mark.skipif(not cubic_roots_available(), reason="Cubic functions not available") +class TestBTExample(object): + @pytest.fixture() + def m(self): + m = ConcreteModel() + + m.fs = FlowsheetBlock(dynamic=False) + + m.fs.props = GenericParameterBlock(**configuration) + + m.fs.state = m.fs.props.build_state_block([1], defined_state=True) + + iscale.calculate_scaling_factors(m.fs.props) + iscale.calculate_scaling_factors(m.fs.state[1]) + + return m + + @pytest.mark.integration + def test_T_sweep(self, m): + assert_units_consistent(m) + + m.fs.obj = Objective(expr=(m.fs.state[1].temperature - 510) ** 2) + m.fs.state[1].temperature.setub(600) + + for P in logspace(4.8, 5.9, 8): + m.fs.obj.deactivate() + + m.fs.state[1].flow_mol.fix(100) + m.fs.state[1].mole_frac_comp["benzene"].fix(0.5) + m.fs.state[1].mole_frac_comp["toluene"].fix(0.5) + m.fs.state[1].temperature.fix(300) + m.fs.state[1].pressure.fix(P) + + m.fs.state.initialize() + + m.fs.state[1].temperature.unfix() + m.fs.obj.activate() + + results = solver.solve(m) + + assert check_optimal_termination(results) + assert m.fs.state[1].flow_mol_phase["Liq"].value <= 1e-5 + + @pytest.mark.integration + def test_P_sweep(self, m): + for T in range(370, 500, 25): + m.fs.state[1].flow_mol.fix(100) + m.fs.state[1].mole_frac_comp["benzene"].fix(0.5) + m.fs.state[1].mole_frac_comp["toluene"].fix(0.5) + m.fs.state[1].temperature.fix(T) + m.fs.state[1].pressure.fix(1e5) + + m.fs.state.initialize() + + results = solver.solve(m) + + assert check_optimal_termination(results) + + while m.fs.state[1].pressure.value <= 1e6: + + results = solver.solve(m) + assert check_optimal_termination(results) + + m.fs.state[1].pressure.value = m.fs.state[1].pressure.value + 1e5 + + @pytest.mark.component + def test_T350_P1_x5(self, m): + m.fs.state[1].flow_mol.fix(100) + m.fs.state[1].mole_frac_comp["benzene"].fix(0.5) + m.fs.state[1].mole_frac_comp["toluene"].fix(0.5) + m.fs.state[1].temperature.fix(350) + m.fs.state[1].pressure.fix(1e5) + + # Trigger build of enthalpy and entropy + m.fs.state[1].enth_mol_phase + m.fs.state[1].entr_mol_phase + + m.fs.state.initialize(outlvl=SOUT) + + from idaes.core.util import DiagnosticsToolbox + + dt = DiagnosticsToolbox(m.fs.state[1]) + dt.report_structural_issues() + dt.display_overconstrained_set() + + results = solver.solve(m, tee=True) + + # Check for optimal solution + assert check_optimal_termination(results) + + assert pytest.approx(value(m.fs.state[1]._teq[("Vap", "Liq")]), abs=1e-1) == 365 + assert 0.0035346 == pytest.approx( + value(m.fs.state[1].compress_fact_phase["Liq"]), 1e-5 + ) + assert 0.966749 == pytest.approx( + value(m.fs.state[1].compress_fact_phase["Vap"]), 1e-5 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Liq", "benzene"]), 1e-5 + ) + == 0.894676 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Liq", "toluene"]), 1e-5 + ) + == 0.347566 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Vap", "benzene"]), 1e-5 + ) + == 0.971072 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Vap", "toluene"]), 1e-5 + ) + == 0.959791 + ) + + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Liq", "benzene"]), 1e-5 + ) + == 0.5 + ) + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Liq", "toluene"]), 1e-5 + ) + == 0.5 + ) + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Vap", "benzene"]), 1e-5 + ) + == 0.70584 + ) + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Vap", "toluene"]), 1e-5 + ) + == 0.29416 + ) + + assert ( + pytest.approx(value(m.fs.state[1].enth_mol_phase["Liq"]), 1e-5) == 38942.8 + ) + assert ( + pytest.approx(value(m.fs.state[1].enth_mol_phase["Vap"]), 1e-5) == 78048.7 + ) + assert ( + pytest.approx(value(m.fs.state[1].entr_mol_phase["Liq"]), 1e-5) == -361.794 + ) + assert ( + pytest.approx(value(m.fs.state[1].entr_mol_phase["Vap"]), 1e-5) == -264.0181 + ) + + @pytest.mark.component + def test_T350_P5_x5(self, m): + m.fs.state[1].flow_mol.fix(100) + m.fs.state[1].mole_frac_comp["benzene"].fix(0.5) + m.fs.state[1].mole_frac_comp["toluene"].fix(0.5) + m.fs.state[1].temperature.fix(350) + m.fs.state[1].pressure.fix(5e5) + + # Trigger build of enthalpy and entropy + m.fs.state[1].enth_mol_phase + m.fs.state[1].entr_mol_phase + + m.fs.state.initialize(outlvl=SOUT) + + results = solver.solve(m) + + # Check for optimal solution + assert check_optimal_termination(results) + + assert pytest.approx(value(m.fs.state[1]._teq[("Vap", "Liq")]), 1e-5) == 431.47 + assert ( + pytest.approx(value(m.fs.state[1].compress_fact_phase["Liq"]), 1e-5) + == 0.01766 + ) + assert ( + pytest.approx(value(m.fs.state[1].compress_fact_phase["Vap"]), 1e-5) + == 0.80245 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Liq", "benzene"]), 1e-5 + ) + == 0.181229 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Liq", "toluene"]), 1e-5 + ) + == 0.070601 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Vap", "benzene"]), 1e-5 + ) + == 0.856523 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Vap", "toluene"]), 1e-5 + ) + == 0.799237 + ) + + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Liq", "benzene"]), 1e-5 + ) + == 0.5 + ) + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Liq", "toluene"]), 1e-5 + ) + == 0.5 + ) + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Vap", "benzene"]), 1e-5 + ) + == 0.65415 + ) + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Vap", "toluene"]), 1e-5 + ) + == 0.34585 + ) + + assert ( + pytest.approx(value(m.fs.state[1].enth_mol_phase["Liq"]), 1e-5) == 38966.9 + ) + assert ( + pytest.approx(value(m.fs.state[1].enth_mol_phase["Vap"]), 1e-5) == 75150.7 + ) + assert ( + pytest.approx(value(m.fs.state[1].entr_mol_phase["Liq"]), 1e-5) == -361.8433 + ) + assert ( + pytest.approx(value(m.fs.state[1].entr_mol_phase["Vap"]), 1e-5) == -281.9703 + ) + + @pytest.mark.component + def test_T450_P1_x5(self, m): + m.fs.state[1].flow_mol.fix(100) + m.fs.state[1].mole_frac_comp["benzene"].fix(0.5) + m.fs.state[1].mole_frac_comp["toluene"].fix(0.5) + m.fs.state[1].temperature.fix(450) + m.fs.state[1].pressure.fix(1e5) + + # Trigger build of enthalpy and entropy + m.fs.state[1].enth_mol_phase + m.fs.state[1].entr_mol_phase + + m.fs.state.initialize(outlvl=SOUT) + + results = solver.solve(m) + + # Check for optimal solution + assert check_optimal_termination(results) + + assert pytest.approx(value(m.fs.state[1]._teq[("Vap", "Liq")]), 1e-5) == 371.4 + assert 0.0033583 == pytest.approx( + value(m.fs.state[1].compress_fact_phase["Liq"]), 1e-5 + ) + assert 0.9821368 == pytest.approx( + value(m.fs.state[1].compress_fact_phase["Vap"]), 1e-5 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Liq", "benzene"]), 1e-5 + ) + == 8.069323 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Liq", "toluene"]), 1e-5 + ) + == 4.304955 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Vap", "benzene"]), 1e-5 + ) + == 0.985365 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Vap", "toluene"]), 1e-5 + ) + == 0.979457 + ) + + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Liq", "benzene"]), 1e-5 + ) + == 0.29861 + ) + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Liq", "toluene"]), 1e-5 + ) + == 0.70139 + ) + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Vap", "benzene"]), 1e-5 + ) + == 0.5 + ) + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Vap", "toluene"]), 1e-5 + ) + == 0.5 + ) + + assert ( + pytest.approx(value(m.fs.state[1].enth_mol_phase["Liq"]), 1e-5) == 49441.2 + ) + assert ( + pytest.approx(value(m.fs.state[1].enth_mol_phase["Vap"]), 1e-5) == 84175.1 + ) + assert ( + pytest.approx(value(m.fs.state[1].entr_mol_phase["Liq"]), 1e-5) == -328.766 + ) + assert ( + pytest.approx(value(m.fs.state[1].entr_mol_phase["Vap"]), 1e-5) == -241.622 + ) + + @pytest.mark.component + def test_T450_P5_x5(self, m): + m.fs.state[1].flow_mol.fix(100) + m.fs.state[1].mole_frac_comp["benzene"].fix(0.5) + m.fs.state[1].mole_frac_comp["toluene"].fix(0.5) + m.fs.state[1].temperature.fix(450) + m.fs.state[1].pressure.fix(5e5) + + # Trigger build of enthalpy and entropy + m.fs.state[1].enth_mol_phase + m.fs.state[1].entr_mol_phase + + m.fs.state.initialize(outlvl=SOUT) + + results = solver.solve(m) + + # Check for optimal solution + assert check_optimal_termination(results) + + assert pytest.approx(value(m.fs.state[1]._teq[("Vap", "Liq")]), 1e-5) == 436.93 + assert 0.0166181 == pytest.approx( + value(m.fs.state[1].compress_fact_phase["Liq"]), 1e-5 + ) + assert 0.9053766 == pytest.approx( + value(m.fs.state[1].compress_fact_phase["Vap"]), 1e-5 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Liq", "benzene"]), 1e-5 + ) + == 1.63308 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Liq", "toluene"]), 1e-5 + ) + == 0.873213 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Vap", "benzene"]), 1e-5 + ) + == 0.927534 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Vap", "toluene"]), 1e-5 + ) + == 0.898324 + ) + + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Liq", "benzene"]), 1e-5 + ) + == 0.3488737 + ) + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Liq", "toluene"]), 1e-5 + ) + == 0.6511263 + ) + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Vap", "benzene"]), 1e-5 + ) + == 0.5 + ) + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Vap", "toluene"]), 1e-5 + ) + == 0.5 + ) + + assert ( + pytest.approx(value(m.fs.state[1].enth_mol_phase["Liq"]), 1e-5) == 51095.2 + ) + assert ( + pytest.approx(value(m.fs.state[1].enth_mol_phase["Vap"]), 1e-5) == 83362.3 + ) + assert ( + pytest.approx(value(m.fs.state[1].entr_mol_phase["Liq"]), 1e-5) == -326.299 + ) + assert ( + pytest.approx(value(m.fs.state[1].entr_mol_phase["Vap"]), 1e-5) == -256.198 + ) + + @pytest.mark.component + def test_T368_P1_x5(self, m): + m.fs.state[1].flow_mol.fix(100) + m.fs.state[1].mole_frac_comp["benzene"].fix(0.5) + m.fs.state[1].mole_frac_comp["toluene"].fix(0.5) + m.fs.state[1].temperature.fix(368) + m.fs.state[1].pressure.fix(1e5) + + # Trigger build of enthalpy and entropy + m.fs.state[1].enth_mol_phase + m.fs.state[1].entr_mol_phase + + m.fs.state.initialize(outlvl=SOUT) + + results = solver.solve(m) + + # Check for optimal solution + assert check_optimal_termination(results) + + assert pytest.approx(value(m.fs.state[1]._teq[("Vap", "Liq")]), 1e-5) == 368 + assert 0.003504 == pytest.approx( + value(m.fs.state[1].compress_fact_phase["Liq"]), 1e-5 + ) + assert 0.97 == pytest.approx( + value(m.fs.state[1].compress_fact_phase["Vap"]), 1e-5 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Liq", "benzene"]), 1e-5 + ) + == 1.492049 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Liq", "toluene"]), 1e-5 + ) + == 0.621563 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Vap", "benzene"]), 1e-5 + ) + == 0.97469 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Vap", "toluene"]), 1e-5 + ) + == 0.964642 + ) + + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Liq", "benzene"]), 1e-5 + ) + == 0.4012128 + ) + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Liq", "toluene"]), 1e-5 + ) + == 0.5987872 + ) + + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Vap", "benzene"]), 1e-5 + ) + == 0.6141738 + ) + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Vap", "toluene"]), 1e-5 + ) + == 0.3858262 + ) + + m.fs.state[1].mole_frac_phase_comp.display() + m.fs.state[1].enth_mol_phase_comp.display() + + assert ( + pytest.approx(value(m.fs.state[1].enth_mol_phase["Liq"]), 1e-5) == 38235.1 + ) + assert ( + pytest.approx(value(m.fs.state[1].enth_mol_phase["Vap"]), 1e-5) == 77155.4 + ) + assert ( + pytest.approx(value(m.fs.state[1].entr_mol_phase["Liq"]), 1e-5) == -359.256 + ) + assert ( + pytest.approx(value(m.fs.state[1].entr_mol_phase["Vap"]), 1e-5) == -262.348 + ) + + @pytest.mark.component + def test_T376_P1_x2(self, m): + m.fs.state[1].flow_mol.fix(100) + m.fs.state[1].mole_frac_comp["benzene"].fix(0.2) + m.fs.state[1].mole_frac_comp["toluene"].fix(0.8) + m.fs.state[1].temperature.fix(376) + m.fs.state[1].pressure.fix(1e5) + + # Trigger build of enthalpy and entropy + m.fs.state[1].enth_mol_phase + m.fs.state[1].entr_mol_phase + + m.fs.state.initialize(outlvl=SOUT) + + results = solver.solve(m) + + # Check for optimal solution + assert check_optimal_termination(results) + + assert pytest.approx(value(m.fs.state[1]._teq[("Vap", "Liq")]), 1e-5) == 376 + assert 0.00361333 == pytest.approx( + value(m.fs.state[1].compress_fact_phase["Liq"]), 1e-5 + ) + assert 0.968749 == pytest.approx( + value(m.fs.state[1].compress_fact_phase["Vap"]), 1e-5 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Liq", "benzene"]), 1e-5 + ) + == 1.8394188 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Liq", "toluene"]), 1e-5 + ) + == 0.7871415 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Vap", "benzene"]), 1e-5 + ) + == 0.9763608 + ) + assert ( + pytest.approx( + value(m.fs.state[1].fug_coeff_phase_comp["Vap", "toluene"]), 1e-5 + ) + == 0.9663611 + ) + + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Liq", "benzene"]), 1e-5 + ) + == 0.17342 + ) + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Liq", "toluene"]), 1e-5 + ) + == 0.82658 + ) + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Vap", "benzene"]), 1e-5 + ) + == 0.3267155 + ) + assert ( + pytest.approx( + value(m.fs.state[1].mole_frac_phase_comp["Vap", "toluene"]), 1e-5 + ) + == 0.6732845 + ) + + assert ( + pytest.approx(value(m.fs.state[1].enth_mol_phase["Liq"]), 1e-5) == 31535.8 + ) + assert ( + pytest.approx(value(m.fs.state[1].enth_mol_phase["Vap"]), 1e-5) == 69175.3 + ) + assert ( + pytest.approx(value(m.fs.state[1].entr_mol_phase["Liq"]), 1e-5) == -369.033 + ) + assert ( + pytest.approx(value(m.fs.state[1].entr_mol_phase["Vap"]), 1e-5) == -273.513 + ) diff --git a/idaes/models/properties/modular_properties/phase_equil/__init__.py b/idaes/models/properties/modular_properties/phase_equil/__init__.py index f6b207f0d6..b9a28fde2f 100644 --- a/idaes/models/properties/modular_properties/phase_equil/__init__.py +++ b/idaes/models/properties/modular_properties/phase_equil/__init__.py @@ -11,3 +11,4 @@ # for full copyright and license information. ################################################################################# from .smooth_VLE import SmoothVLE +from .smooth_VLE_2 import CubicComplementarityVLE diff --git a/idaes/models/properties/modular_properties/phase_equil/bubble_dew.py b/idaes/models/properties/modular_properties/phase_equil/bubble_dew.py index 4bd953e686..f6d75d1931 100644 --- a/idaes/models/properties/modular_properties/phase_equil/bubble_dew.py +++ b/idaes/models/properties/modular_properties/phase_equil/bubble_dew.py @@ -13,10 +13,7 @@ """ Modular methods for calculating bubble and dew points """ -# TODO: Missing docstrings -# pylint: disable=missing-function-docstring - -# TODO: Look into protected access issues +# TODO: Pylint complains about variables with _x names as they are built by other classes # pylint: disable=protected-access from pyomo.environ import Constraint @@ -24,13 +21,11 @@ from idaes.models.properties.modular_properties.base.utility import ( get_method, get_component_object as cobj, + identify_VL_component_list, ) import idaes.core.util.scaling as iscale from idaes.core.util.exceptions import ConfigurationError -# _valid_VL_component_list return variables that are not need in all cases -# pylint: disable=W0612 - class IdealBubbleDew: """Bubble and dew point calculations for ideal systems.""" @@ -43,6 +38,9 @@ class IdealBubbleDew: # calculate concentrations at the bubble and dew points @staticmethod def temperature_bubble(b): + """ + Rule for calculating bubble temperature + """ _non_vle_phase_check(b) try: @@ -52,9 +50,9 @@ def rule_bubble_temp(b, p1, p2): v_phase, vl_comps, henry_comps, - l_only_comps, + _, v_only_comps, - ) = _valid_VL_component_list(b, (p1, p2)) + ) = identify_VL_component_list(b, (p1, p2)) if l_phase is None or v_phase is None: # Not a VLE pair @@ -95,9 +93,9 @@ def rule_mole_frac_bubble_temp(b, p1, p2, j): v_phase, vl_comps, henry_comps, - l_only_comps, + _, v_only_comps, - ) = _valid_VL_component_list(b, (p1, p2)) + ) = identify_VL_component_list(b, (p1, p2)) if l_phase is None or v_phase is None: # Not a VLE pair @@ -129,6 +127,9 @@ def rule_mole_frac_bubble_temp(b, p1, p2, j): @staticmethod def scale_temperature_bubble(b, overwrite=True): + """ + Scaling method for bubble temperature + """ sf_P = iscale.get_scaling_factor(b.pressure, default=1e-5, warning=True) sf_mf = iscale.get_scaling_factor(b.mole_frac_comp, default=1e3, warning=True) @@ -138,10 +139,10 @@ def scale_temperature_bubble(b, overwrite=True): l_phase, v_phase, vl_comps, - henry_comps, - l_only_comps, + _, + _, v_only_comps, - ) = _valid_VL_component_list(b, pp) + ) = identify_VL_component_list(b, pp) if l_phase is None or v_phase is None: continue elif v_only_comps != []: @@ -163,6 +164,9 @@ def scale_temperature_bubble(b, overwrite=True): # Dew temperature methods @staticmethod def temperature_dew(b): + """ + Method for calculating dew temperature + """ _non_vle_phase_check(b) try: @@ -173,8 +177,8 @@ def rule_dew_temp(b, p1, p2): vl_comps, henry_comps, l_only_comps, - v_only_comps, - ) = _valid_VL_component_list(b, (p1, p2)) + _, + ) = identify_VL_component_list(b, (p1, p2)) if l_phase is None or v_phase is None: # Not a VLE pair @@ -218,8 +222,8 @@ def rule_mole_frac_dew_temp(b, p1, p2, j): vl_comps, henry_comps, l_only_comps, - v_only_comps, - ) = _valid_VL_component_list(b, (p1, p2)) + _, + ) = identify_VL_component_list(b, (p1, p2)) if l_phase is None or v_phase is None: # Not a VLE pair @@ -253,6 +257,9 @@ def rule_mole_frac_dew_temp(b, p1, p2, j): @staticmethod def scale_temperature_dew(b, overwrite=True): + """ + Scaling method for dew temperature + """ sf_P = iscale.get_scaling_factor(b.pressure, default=1e-5, warning=True) sf_mf = iscale.get_scaling_factor(b.mole_frac_comp, default=1e3, warning=True) @@ -262,10 +269,10 @@ def scale_temperature_dew(b, overwrite=True): l_phase, v_phase, vl_comps, - henry_comps, - l_only_comps, + _, + _, v_only_comps, - ) = _valid_VL_component_list(b, pp) + ) = identify_VL_component_list(b, pp) if l_phase is None or v_phase is None: continue elif v_only_comps != []: @@ -285,6 +292,9 @@ def scale_temperature_dew(b, overwrite=True): # Bubble pressure methods @staticmethod def pressure_bubble(b): + """ + Method for calculating bubble pressure + """ _non_vle_phase_check(b) try: @@ -294,9 +304,9 @@ def rule_bubble_press(b, p1, p2): v_phase, vl_comps, henry_comps, - l_only_comps, + _, v_only_comps, - ) = _valid_VL_component_list(b, (p1, p2)) + ) = identify_VL_component_list(b, (p1, p2)) if l_phase is None or v_phase is None: # Not a VLE pair @@ -326,9 +336,9 @@ def rule_mole_frac_bubble_press(b, p1, p2, j): v_phase, vl_comps, henry_comps, - l_only_comps, + _, v_only_comps, - ) = _valid_VL_component_list(b, (p1, p2)) + ) = identify_VL_component_list(b, (p1, p2)) if l_phase is None or v_phase is None: # Not a VLE pair @@ -356,6 +366,9 @@ def rule_mole_frac_bubble_press(b, p1, p2, j): @staticmethod def scale_pressure_bubble(b, overwrite=True): + """ + Scaling method for bubble pressure + """ sf_P = iscale.get_scaling_factor(b.pressure, default=1e-5, warning=True) sf_mf = iscale.get_scaling_factor(b.mole_frac_comp, default=1e3, warning=True) @@ -365,10 +378,10 @@ def scale_pressure_bubble(b, overwrite=True): l_phase, v_phase, vl_comps, - henry_comps, - l_only_comps, + _, + _, v_only_comps, - ) = _valid_VL_component_list(b, pp) + ) = identify_VL_component_list(b, pp) if l_phase is None or v_phase is None: continue elif v_only_comps != []: @@ -390,6 +403,9 @@ def scale_pressure_bubble(b, overwrite=True): # Dew pressure methods @staticmethod def pressure_dew(b): + """ + Method for calculating dew pressure + """ _non_vle_phase_check(b) try: @@ -400,8 +416,8 @@ def rule_dew_press(b, p1, p2): vl_comps, henry_comps, l_only_comps, - v_only_comps, - ) = _valid_VL_component_list(b, (p1, p2)) + _, + ) = identify_VL_component_list(b, (p1, p2)) if l_phase is None or v_phase is None: # Not a VLE pair @@ -430,8 +446,8 @@ def rule_mole_frac_dew_press(b, p1, p2, j): vl_comps, henry_comps, l_only_comps, - v_only_comps, - ) = _valid_VL_component_list(b, (p1, p2)) + _, + ) = identify_VL_component_list(b, (p1, p2)) if l_phase is None or v_phase is None: # Not a VLE pair @@ -459,6 +475,9 @@ def rule_mole_frac_dew_press(b, p1, p2, j): @staticmethod def scale_pressure_dew(b, overwrite=True): + """ + Scaling method for dew pressure + """ sf_P = iscale.get_scaling_factor(b.pressure, default=1e-5, warning=True) sf_mf = iscale.get_scaling_factor(b.mole_frac_comp, default=1e3, warning=True) @@ -468,10 +487,10 @@ def scale_pressure_dew(b, overwrite=True): l_phase, v_phase, vl_comps, - henry_comps, - l_only_comps, + _, + _, v_only_comps, - ) = _valid_VL_component_list(b, pp) + ) = identify_VL_component_list(b, pp) if l_phase is None or v_phase is None: continue elif v_only_comps != []: @@ -495,6 +514,9 @@ class LogBubbleDew: # Bubble temperature methods @staticmethod def temperature_bubble(b): + """ + Method for constructing bubble temperature constraint + """ try: def rule_bubble_temp(b, p1, p2, j): @@ -503,9 +525,9 @@ def rule_bubble_temp(b, p1, p2, j): v_phase, vl_comps, henry_comps, - l_only_comps, + _, v_only_comps, - ) = _valid_VL_component_list(b, (p1, p2)) + ) = identify_VL_component_list(b, (p1, p2)) if l_phase is None or v_phase is None: # Not a VLE pair @@ -538,9 +560,9 @@ def rule_mole_frac_bubble_temp(b, p1, p2): v_phase, vl_comps, henry_comps, - l_only_comps, + _, v_only_comps, - ) = _valid_VL_component_list(b, (p1, p2)) + ) = identify_VL_component_list(b, (p1, p2)) if l_phase is None or v_phase is None: # Not a VLE pair @@ -560,17 +582,20 @@ def rule_mole_frac_bubble_temp(b, p1, p2): @staticmethod def scale_temperature_bubble(b, overwrite=True): + """ + Method for scaling bubble temperature + """ sf_mf = iscale.get_scaling_factor(b.mole_frac_comp, default=1e3, warning=True) for pp in b.params._pe_pairs: ( l_phase, v_phase, - vl_comps, - henry_comps, - l_only_comps, + _, + _, + _, v_only_comps, - ) = _valid_VL_component_list(b, pp) + ) = identify_VL_component_list(b, pp) if l_phase is None or v_phase is None: continue elif v_only_comps != []: @@ -585,6 +610,9 @@ def scale_temperature_bubble(b, overwrite=True): # Dew temperature methods @staticmethod def temperature_dew(b): + """ + Method for constructing dew temperature constraint + """ try: def rule_dew_temp(b, p1, p2, j): @@ -594,8 +622,8 @@ def rule_dew_temp(b, p1, p2, j): vl_comps, henry_comps, l_only_comps, - v_only_comps, - ) = _valid_VL_component_list(b, (p1, p2)) + _, + ) = identify_VL_component_list(b, (p1, p2)) if l_phase is None or v_phase is None: # Not a VLE pair @@ -629,8 +657,8 @@ def rule_mole_frac_dew_temp(b, p1, p2): vl_comps, henry_comps, l_only_comps, - v_only_comps, - ) = _valid_VL_component_list(b, (p1, p2)) + _, + ) = identify_VL_component_list(b, (p1, p2)) if l_phase is None or v_phase is None: # Not a VLE pair @@ -650,17 +678,20 @@ def rule_mole_frac_dew_temp(b, p1, p2): @staticmethod def scale_temperature_dew(b, overwrite=True): + """ + Method for scaling dew temperature + """ sf_mf = iscale.get_scaling_factor(b.mole_frac_comp, default=1e3, warning=True) for pp in b.params._pe_pairs: ( l_phase, v_phase, - vl_comps, - henry_comps, - l_only_comps, + _, + _, + _, v_only_comps, - ) = _valid_VL_component_list(b, pp) + ) = identify_VL_component_list(b, pp) if l_phase is None or v_phase is None: continue elif v_only_comps != []: @@ -675,6 +706,9 @@ def scale_temperature_dew(b, overwrite=True): # Bubble pressure methods @staticmethod def pressure_bubble(b): + """ + Method for constructing bubble pressure + """ try: def rule_bubble_press(b, p1, p2, j): @@ -683,9 +717,9 @@ def rule_bubble_press(b, p1, p2, j): v_phase, vl_comps, henry_comps, - l_only_comps, + _, v_only_comps, - ) = _valid_VL_component_list(b, (p1, p2)) + ) = identify_VL_component_list(b, (p1, p2)) if l_phase is None or v_phase is None: # Not a VLE pair @@ -718,9 +752,9 @@ def rule_mole_frac_bubble_press(b, p1, p2): v_phase, vl_comps, henry_comps, - l_only_comps, + _, v_only_comps, - ) = _valid_VL_component_list(b, (p1, p2)) + ) = identify_VL_component_list(b, (p1, p2)) if l_phase is None or v_phase is None: # Not a VLE pair @@ -740,17 +774,20 @@ def rule_mole_frac_bubble_press(b, p1, p2): @staticmethod def scale_pressure_bubble(b, overwrite=True): + """ + Method for scaling bubble pressure + """ sf_mf = iscale.get_scaling_factor(b.mole_frac_comp, default=1e3, warning=True) for pp in b.params._pe_pairs: ( l_phase, v_phase, - vl_comps, - henry_comps, - l_only_comps, + _, + _, + _, v_only_comps, - ) = _valid_VL_component_list(b, pp) + ) = identify_VL_component_list(b, pp) if l_phase is None or v_phase is None: continue elif v_only_comps != []: @@ -765,6 +802,9 @@ def scale_pressure_bubble(b, overwrite=True): # Dew pressure methods @staticmethod def pressure_dew(b): + """ + Method constructing dew pressure constraints + """ try: def rule_dew_press(b, p1, p2, j): @@ -773,9 +813,9 @@ def rule_dew_press(b, p1, p2, j): v_phase, vl_comps, henry_comps, - l_only_comps, + _, v_only_comps, - ) = _valid_VL_component_list(b, (p1, p2)) + ) = identify_VL_component_list(b, (p1, p2)) if l_phase is None or v_phase is None: # Not a VLE pair @@ -809,8 +849,8 @@ def rule_mole_frac_dew_press(b, p1, p2): vl_comps, henry_comps, l_only_comps, - v_only_comps, - ) = _valid_VL_component_list(b, (p1, p2)) + _, + ) = identify_VL_component_list(b, (p1, p2)) if l_phase is None or v_phase is None: # Not a VLE pair @@ -830,17 +870,20 @@ def rule_mole_frac_dew_press(b, p1, p2): @staticmethod def scale_pressure_dew(b, overwrite=True): + """ + Method for scaling dew pressure + """ sf_mf = iscale.get_scaling_factor(b.mole_frac_comp, default=1e3, warning=True) for pp in b.params._pe_pairs: ( l_phase, v_phase, - vl_comps, - henry_comps, - l_only_comps, + _, + _, + _, v_only_comps, - ) = _valid_VL_component_list(b, pp) + ) = identify_VL_component_list(b, pp) if l_phase is None or v_phase is None: continue elif v_only_comps != []: @@ -852,48 +895,6 @@ def scale_pressure_dew(b, overwrite=True): ) -def _valid_VL_component_list(blk, pp): - vl_comps = [] - henry_comps = [] - l_only_comps = [] - v_only_comps = [] - - pparams = blk.params - l_phase = None - v_phase = None - if pparams.get_phase(pp[0]).is_liquid_phase(): - l_phase = pp[0] - elif pparams.get_phase(pp[0]).is_vapor_phase(): - v_phase = pp[0] - - if pparams.get_phase(pp[1]).is_liquid_phase(): - l_phase = pp[1] - elif pparams.get_phase(pp[1]).is_vapor_phase(): - v_phase = pp[1] - - # Only need to do this for V-L pairs, so check - if l_phase is not None and v_phase is not None: - for j in blk.params.component_list: - if (l_phase, j) in blk.phase_component_set and ( - v_phase, - j, - ) in blk.phase_component_set: - cobj = pparams.get_component(j) - if cobj.config.henry_component is not None and ( - pp[0] in cobj.config.henry_component - or pp[1] in cobj.config.henry_component - ): - henry_comps.append(j) - else: - vl_comps.append(j) - elif (l_phase, j) in blk.phase_component_set: - l_only_comps.append(j) - elif (v_phase, j) in blk.phase_component_set: - v_only_comps.append(j) - - return l_phase, v_phase, vl_comps, henry_comps, l_only_comps, v_only_comps - - def _non_vle_phase_check(blk): if len(blk.phase_list) > 2: raise ConfigurationError( diff --git a/idaes/models/properties/modular_properties/phase_equil/smooth_VLE.py b/idaes/models/properties/modular_properties/phase_equil/smooth_VLE.py index 4b2718b691..39eedbe205 100644 --- a/idaes/models/properties/modular_properties/phase_equil/smooth_VLE.py +++ b/idaes/models/properties/modular_properties/phase_equil/smooth_VLE.py @@ -18,17 +18,13 @@ Flowsheet Optimization. Proceedings of the 13th International Symposium on Process Systems Engineering – PSE 2018, July 1-5, 2018, San Diego. """ -# TODO: Missing docstrings -# pylint: disable=missing-function-docstring - -# TODO: Look into protected access issues +# TODO: Pylint complains about variables with _x names as they are built by other classes # pylint: disable=protected-access from pyomo.environ import Constraint, Param, Var, value -from idaes.core.util.exceptions import ConfigurationError from idaes.core.util.math import smooth_max, smooth_min -from idaes.models.properties.modular_properties.phase_equil.bubble_dew import ( - _valid_VL_component_list, +from idaes.models.properties.modular_properties.base.utility import ( + identify_VL_component_list, ) import idaes.core.util.scaling as iscale @@ -39,27 +35,22 @@ class SmoothVLE(object): @staticmethod def phase_equil(b, phase_pair): + """ + Method for constructing phase equilibrium variables and constraints + """ # This method is called via StateBlock.build, thus does not need clean-up # try/except statements suffix = "_" + phase_pair[0] + "_" + phase_pair[1] # Smooth VLE assumes a liquid and a vapor phase, so validate this ( - l_phase, - v_phase, + _, + _, _, _, l_only_comps, v_only_comps, - ) = _valid_VL_component_list(b, phase_pair) - - if l_phase is None or v_phase is None: - raise ConfigurationError( - "{} Generic Property Package phase pair {}-{} was set to use " - "Smooth VLE formulation, however this is not a vapor-liquid pair.".format( - b.params.name, phase_pair[0], phase_pair[1] - ) - ) + ) = identify_VL_component_list(b, phase_pair) # Definition of equilibrium temperature for smooth VLE t_units = b.params.get_metadata().default_units.TEMPERATURE @@ -128,6 +119,9 @@ def rule_teq(b): @staticmethod def calculate_scaling_factors(b, phase_pair): + """ + Method to calculate scaling factors for phase equilibrium + """ suffix = "_" + phase_pair[0] + "_" + phase_pair[1] sf_T = iscale.get_scaling_factor(b.temperature, default=1, warning=True) @@ -144,6 +138,9 @@ def calculate_scaling_factors(b, phase_pair): @staticmethod def phase_equil_initialization(b, phase_pair): + """ + Method to initialize phase equilibrium + """ suffix = "_" + phase_pair[0] + "_" + phase_pair[1] for c in b.component_objects(Constraint): @@ -153,19 +150,22 @@ def phase_equil_initialization(b, phase_pair): @staticmethod def calculate_teq(b, phase_pair): + """ + Method to calculate initial guess for equilibrium temperature + """ suffix = "_" + phase_pair[0] + "_" + phase_pair[1] if hasattr(b, "eq_temperature_bubble"): _t1 = getattr(b, "_t1" + suffix) - _t1.value = max( - value(b.temperature), b.temperature_bubble[phase_pair].value + _t1.set_value( + max(value(b.temperature), value(b.temperature_bubble[phase_pair])) ) else: _t1 = b.temperature if hasattr(b, "eq_temperature_dew"): - b._teq[phase_pair].value = min( - _t1.value, b.temperature_dew[phase_pair].value + b._teq[phase_pair].set_value( + min(value(_t1), value(b.temperature_dew[phase_pair])) ) else: - b._teq[phase_pair].value = _t1.value + b._teq[phase_pair].set_value(value(_t1)) diff --git a/idaes/models/properties/modular_properties/phase_equil/smooth_VLE_2.py b/idaes/models/properties/modular_properties/phase_equil/smooth_VLE_2.py new file mode 100644 index 0000000000..2f2bd30263 --- /dev/null +++ b/idaes/models/properties/modular_properties/phase_equil/smooth_VLE_2.py @@ -0,0 +1,339 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2023 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# +""" +Implementation of the formulation proposed in: + +Dabadghao, V., Ghouse, J., Eslick, J., Lee, A., Burgard, A., Miller, D., Biegler, L., 2022, +A complementarity-based vapor-liquid equilibrium formulation for equation-oriented simulation +and optimization. AIChE Journal, DOI: 10.1002/aic.18029 +""" + +from pyomo.environ import ( + Constraint, + Expression, + Param, + Set, + units as pyunits, + Var, + value, +) + +from idaes.core.util.exceptions import ConfigurationError +from idaes.core.util.math import smooth_min +from idaes.models.properties.modular_properties.base.utility import ( + identify_VL_component_list, +) +from idaes.models.properties.modular_properties.base.utility import ( + estimate_Tbub, + estimate_Tdew, +) +from idaes.models.properties.modular_properties.eos.ceos import ( + calculate_equilibrium_cubic_coefficients, +) +import idaes.core.util.scaling as iscale + + +# Small value for initializing slack variables +EPS_INIT = 1e-4 + + +# ----------------------------------------------------------------------------- +class CubicComplementarityVLE: + """ + Improved Vapor-Liquid Equilibrium complementarity formulation for Cubic Equations of State. + """ + + @staticmethod + def phase_equil(b, phase_pair): + """ + Method for constructing phase equilibrium variables and constraints + """ + # This method is called via StateBlock.build, thus does not need clean-up + # try/except statements + suffix = "_" + phase_pair[0] + "_" + phase_pair[1] + + # Smooth VLE assumes a liquid and a vapor phase, so validate this + ( + l_phase, + v_phase, + _, + _, + _, + _, + ) = identify_VL_component_list(b, phase_pair) + + if l_phase is None or v_phase is None: + raise ConfigurationError( + f"{b.params.name} - Generic Property Package phase pair {phase_pair[0]}-{phase_pair[1]} " + "was set to use CubicComplementarityVLE formulation, however this is not a vapor-liquid pair." + ) + + try: + lobj = b.params.get_phase(l_phase) + ctype = lobj._cubic_type # pylint: disable=protected-access + cname = lobj.config.equation_of_state_options["type"].name + vobj = b.params.get_phase(v_phase) + + if ( + ctype != vobj._cubic_type # pylint: disable=protected-access + or lobj.config.equation_of_state_options + != vobj.config.equation_of_state_options + ): + raise ConfigurationError( + f"{b.params.name} - CubicComplementarityVLE formulation requires that both " + "phases use the same type of cubic equation of state." + ) + except AttributeError: + raise ConfigurationError( + f"{b.params.name} - CubicComplementarityVLE formulation only supports cubic equations of state." + ) + + # Definition of equilibrium temperature for smooth VLE + uom = b.params.get_metadata().default_units + t_units = uom.TEMPERATURE + f_units = uom.AMOUNT / uom.TIME + + vl_phase_set = Set(initialize=[phase_pair[0], phase_pair[1]]) + b.add_component("_vle_set" + suffix, vl_phase_set) + + s = Var( + vl_phase_set, + initialize=EPS_INIT, + bounds=(0, None), + doc="Slack variable for equilibrium temperature", + units=t_units, + ) + b.add_component("s" + suffix, s) + + # Equilibrium temperature + def rule_teq(b): + return ( + b._teq[phase_pair] # pylint: disable=protected-access + - b.temperature + - s[v_phase] + + s[l_phase] + == 0 + ) + + b.add_component("_teq_constraint" + suffix, Constraint(rule=rule_teq)) + + # Epsilon variables will be given units of flow, as this is usually what dominates + # the complementarity equations. + eps_t = Param( + default=1, + mutable=True, + doc="Smoothing parameter for temperature complementarity", + units=f_units, + ) + b.add_component("eps_t" + suffix, eps_t) + + eps_z = Param( + default=1, + mutable=True, + doc="Smoothing parameter for cubic root complementarities", + units=f_units, + ) + b.add_component("eps_z" + suffix, eps_z) + + gp = Var( + vl_phase_set, + initialize=EPS_INIT, + bounds=(0, None), + doc="Slack variable for cubic second derivative for phase p", + units=pyunits.dimensionless, + ) + b.add_component("gp" + suffix, gp) + + gn = Var( + vl_phase_set, + initialize=EPS_INIT, + bounds=(0, None), + doc="Slack variable for cubic second derivative for phase p", + units=pyunits.dimensionless, + ) + b.add_component("gn" + suffix, gn) + + def rule_temperature_slack_complementarity(b, p): + flow_phase = b.flow_mol_phase[p] + + return smooth_min(s[p] * f_units / t_units, flow_phase, eps_t) == 0 + + b.add_component( + "temperature_slack_complementarity" + suffix, + Constraint( + vl_phase_set, + rule=rule_temperature_slack_complementarity, + ), + ) + + def rule_cubic_second_derivative(b, p): + p1, p2 = phase_pair + + _b, _, _ = calculate_equilibrium_cubic_coefficients( + b, cname, ctype, p1, p2, p + ) + z = b.compress_fact_phase[p] + + return 6 * z + 2 * _b + + b.add_component( + "cubic_second_derivative" + suffix, + Expression(vl_phase_set, rule=rule_cubic_second_derivative), + ) + + def rule_cubic_root_complementarity(b, p): + der = getattr(b, "cubic_second_derivative" + suffix) + return der[p] == gp[p] - gn[p] + + b.add_component( + "cubic_root_complementarity" + suffix, + Constraint(vl_phase_set, rule=rule_cubic_root_complementarity), + ) + + def rule_cubic_slack_complementarity(b, p): + flow_phase = b.flow_mol_phase[p] + if b.params.get_phase(p).is_vapor_phase(): + return smooth_min(gn[p] * f_units, flow_phase, eps_z) == 0 + else: + return smooth_min(gp[p] * f_units, flow_phase, eps_z) == 0 + + b.add_component( + "cubic_slack_complementarity" + suffix, + Constraint(vl_phase_set, rule=rule_cubic_slack_complementarity), + ) + + @staticmethod + def calculate_scaling_factors(b, phase_pair): + """ + Method to calculate scaling factors for phase equilibrium + """ + suffix = "_" + phase_pair[0] + "_" + phase_pair[1] + sf_T = iscale.get_scaling_factor(b.temperature, default=1, warning=True) + + try: + teq_cons = getattr(b, "_teq_constraint" + suffix) + # pylint: disable-next=protected-access + iscale.set_scaling_factor(b._teq[phase_pair], sf_T) + iscale.constraint_scaling_transform(teq_cons, sf_T, overwrite=False) + except AttributeError: + pass + + @staticmethod + def calculate_teq(blk, pp): + """ + Method to calculate initial guess for equilibrium temperature + """ + # --------------------------------------------------------------------- + # If present, initialize bubble and dew point calculations, and + # equilibrium temperature _teq + T_units = blk.params.get_metadata().default_units.TEMPERATURE + + ( + liquid_phase, + vapor_phase, + raoult_comps, + henry_comps, + _, + v_only_comps, + ) = identify_VL_component_list(blk, pp) + + if v_only_comps is None: + if blk.is_property_constructed("temperature_bubble"): + Tbub = value(blk.temperature_bubble[pp]) + else: + Tbub = estimate_Tbub( + blk, T_units, raoult_comps, henry_comps, liquid_phase + ) + t1 = max(value(blk.temperature), Tbub) + else: + t1 = value(blk.temperature) + + if v_only_comps is None: + if blk.is_property_constructed("temperature_dew"): + Tdew = value(blk.temperature_bubble[pp]) + else: + Tdew = estimate_Tdew( + blk, T_units, raoult_comps, henry_comps, liquid_phase + ) + t2 = min(t1, Tdew) + else: + t2 = t1 + + blk._teq[pp].set_value(t2) # pylint: disable=protected-access + + # --------------------------------------------------------------------- + # Initialize sV and sL slacks + _calculate_temperature_slacks(blk, pp, liquid_phase, vapor_phase) + + # --------------------------------------------------------------------- + # If flash, initialize g+ and g- slacks + _calculate_ceos_derivative_slacks(blk, pp, liquid_phase, vapor_phase) + + @staticmethod + def phase_equil_initialization(b, phase_pair): + """ + Method to initialize phase equilibrium + """ + suffix = "_" + phase_pair[0] + "_" + phase_pair[1] + + for c in b.component_objects(Constraint): + # Activate equilibrium constraints + if c.local_name in ( + "_teq_constraint" + suffix, + "temperature_slack_complementarity" + suffix, + "cubic_root_complementarity" + suffix, + "cubic_slack_complementarity" + suffix, + ): + c.activate() + + +def _calculate_temperature_slacks(b, phase_pair, liquid_phase, vapor_phase): + suffix = "_" + phase_pair[0] + "_" + phase_pair[1] + + s = getattr(b, "s" + suffix) + + # pylint: disable-next=protected-access + if value(b._teq[phase_pair]) > value(b.temperature): + # pylint: disable-next=protected-access + s[vapor_phase].set_value(value(b._teq[phase_pair] - b.temperature)) + s[liquid_phase].set_value(EPS_INIT) + # pylint: disable-next=protected-access + elif value(b._teq[phase_pair]) < value(b.temperature): + s[vapor_phase].set_value(EPS_INIT) + # pylint: disable-next=protected-access + s[liquid_phase].set_value(value(b.temperature - b._teq[phase_pair])) + else: + s[vapor_phase].set_value(EPS_INIT) + s[liquid_phase].set_value(EPS_INIT) + + +def _calculate_ceos_derivative_slacks(b, phase_pair, liquid_phase, vapor_phase): + suffix = "_" + phase_pair[0] + "_" + phase_pair[1] + + gp = getattr(b, "gp" + suffix) + gn = getattr(b, "gn" + suffix) + der = getattr(b, "cubic_second_derivative" + suffix) + + if value(der[liquid_phase]) < 0: + gp[liquid_phase].set_value(EPS_INIT) + gn[liquid_phase].set_value(value(-der[liquid_phase])) + else: + gp[liquid_phase].set_value(EPS_INIT) + gn[liquid_phase].set_value(EPS_INIT) + + if value(der[vapor_phase]) > 0: + gp[vapor_phase].set_value(value(der[vapor_phase])) + gn[vapor_phase].set_value(EPS_INIT) + else: + gp[vapor_phase].set_value(EPS_INIT) + gn[vapor_phase].set_value(EPS_INIT) diff --git a/idaes/models/properties/modular_properties/phase_equil/tests/test_smooth_VLE.py b/idaes/models/properties/modular_properties/phase_equil/tests/test_smooth_VLE.py index e8b1684af1..8d42935baf 100644 --- a/idaes/models/properties/modular_properties/phase_equil/tests/test_smooth_VLE.py +++ b/idaes/models/properties/modular_properties/phase_equil/tests/test_smooth_VLE.py @@ -33,7 +33,7 @@ from idaes.models.properties.modular_properties.state_definitions import FTPx from idaes.models.properties.modular_properties.phase_equil import SmoothVLE from idaes.models.properties.modular_properties.phase_equil.forms import fugacity -from idaes.core.util.exceptions import ConfigurationError +from idaes.core.util.exceptions import PropertyPackageError # Dummy EoS to use for fugacity calls @@ -166,9 +166,8 @@ def test_non_VLE_pair(): m.props = m.params.state_block_class([1], parameters=m.params) with pytest.raises( - ConfigurationError, - match="params Generic Property Package phase pair " - "Liq-Sol was set to use Smooth VLE formulation, " - "however this is not a vapor-liquid pair.", + PropertyPackageError, + match="Phase pair Liq-Sol was identified as being a VLE pair, however " + "Liq is liquid but Sol is not vapor.", ): SmoothVLE.phase_equil(m.props[1], ("Liq", "Sol")) diff --git a/idaes/models/properties/modular_properties/phase_equil/tests/test_smooth_VLE_2.py b/idaes/models/properties/modular_properties/phase_equil/tests/test_smooth_VLE_2.py new file mode 100644 index 0000000000..5884aa3d6e --- /dev/null +++ b/idaes/models/properties/modular_properties/phase_equil/tests/test_smooth_VLE_2.py @@ -0,0 +1,341 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2023 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# +""" +Tests for CubicComplementarityVLE formulation + +Authors: Andrew Lee +""" + +import pytest + +from pyomo.environ import ( + ConcreteModel, + Constraint, + Expression, + Param, + value, + Var, + units as pyunits, +) + +from idaes.models.properties.modular_properties.base.generic_property import ( + GenericParameterBlock, +) +from idaes.models.properties.modular_properties.state_definitions import FTPx +from idaes.models.properties.modular_properties.phase_equil.smooth_VLE_2 import ( + CubicComplementarityVLE, + _calculate_temperature_slacks, + _calculate_ceos_derivative_slacks, + EPS_INIT, +) +from idaes.models.properties.modular_properties.eos.ideal import Ideal +from idaes.models.properties.modular_properties.eos.ceos import Cubic, CubicType +from idaes.models.properties.modular_properties.phase_equil.forms import fugacity +from idaes.core.util.exceptions import ConfigurationError + + +@pytest.mark.unit +def test_different_cubics(): + m = ConcreteModel() + + m.params = GenericParameterBlock( + components={ + "H2O": { + "parameter_data": { + "pressure_crit": (220.6e5, pyunits.Pa), + "temperature_crit": (647, pyunits.K), + "omega": 0.344, + }, + "phase_equilibrium_form": {("Vap", "Liq"): fugacity}, + } + }, + phases={ + "Liq": { + "equation_of_state": Cubic, + "equation_of_state_options": {"type": CubicType.PR}, + }, + "Vap": { + "equation_of_state": Cubic, + "equation_of_state_options": {"type": CubicType.SRK}, + }, + }, + state_definition=FTPx, + pressure_ref=100000.0, + temperature_ref=300, + base_units={ + "time": pyunits.s, + "length": pyunits.m, + "mass": pyunits.kg, + "amount": pyunits.mol, + "temperature": pyunits.K, + }, + phases_in_equilibrium=[("Vap", "Liq")], + phase_equilibrium_state={("Vap", "Liq"): CubicComplementarityVLE}, + parameter_data={ + "PR_kappa": { + ("H2O", "H2O"): 0.000, + }, + "SRK_kappa": { + ("H2O", "H2O"): 0.000, + }, + }, + ) + + with pytest.raises( + ConfigurationError, + match="params - CubicComplementarityVLE formulation requires that both phases use the same " + "type of cubic equation of state.", + ): + m.props = m.params.state_block_class([1], parameters=m.params) + + +@pytest.mark.unit +def test_non_cubic(): + m = ConcreteModel() + + m.params = GenericParameterBlock( + components={ + "H2O": { + "parameter_data": { + "pressure_crit": (220.6e5, pyunits.Pa), + "temperature_crit": (647, pyunits.K), + "omega": 0.344, + }, + "phase_equilibrium_form": {("Vap", "Liq"): fugacity}, + } + }, + phases={ + "Liq": { + "equation_of_state": Cubic, + "equation_of_state_options": {"type": CubicType.PR}, + }, + "Vap": { + "equation_of_state": Ideal, + }, + }, + state_definition=FTPx, + pressure_ref=100000.0, + temperature_ref=300, + base_units={ + "time": pyunits.s, + "length": pyunits.m, + "mass": pyunits.kg, + "amount": pyunits.mol, + "temperature": pyunits.K, + }, + phases_in_equilibrium=[("Vap", "Liq")], + phase_equilibrium_state={("Vap", "Liq"): CubicComplementarityVLE}, + parameter_data={ + "PR_kappa": { + ("H2O", "H2O"): 0.000, + } + }, + ) + + with pytest.raises( + ConfigurationError, + match="params - CubicComplementarityVLE formulation only supports cubic equations of state.", + ): + m.props = m.params.state_block_class([1], parameters=m.params) + + +@pytest.fixture() +def frame(): + m = ConcreteModel() + + # Create a dummy parameter block + m.params = GenericParameterBlock( + components={ + "H2O": { + "parameter_data": { + "pressure_crit": (220.6e5, pyunits.Pa), + "temperature_crit": (647, pyunits.K), + "omega": 0.344, + }, + "phase_equilibrium_form": {("Vap", "Liq"): fugacity}, + } + }, + phases={ + "Liq": { + "equation_of_state": Cubic, + "equation_of_state_options": {"type": CubicType.PR}, + }, + "Vap": { + "equation_of_state": Cubic, + "equation_of_state_options": {"type": CubicType.PR}, + }, + }, + state_definition=FTPx, + pressure_ref=100000.0, + temperature_ref=300, + base_units={ + "time": pyunits.s, + "length": pyunits.m, + "mass": pyunits.kg, + "amount": pyunits.mol, + "temperature": pyunits.K, + }, + phases_in_equilibrium=[("Vap", "Liq")], + phase_equilibrium_state={("Vap", "Liq"): CubicComplementarityVLE}, + parameter_data={ + "PR_kappa": { + ("H2O", "H2O"): 0.000, + } + }, + ) + + # Create a dummy state block + m.props = m.params.state_block_class([1], parameters=m.params) + + return m + + +@pytest.mark.unit +def test_build(frame): + assert isinstance(frame.props[1].s_Vap_Liq, Var) + assert isinstance(frame.props[1].gp_Vap_Liq, Var) + assert isinstance(frame.props[1].gn_Vap_Liq, Var) + + assert isinstance(frame.props[1].cubic_second_derivative_Vap_Liq, Expression) + + assert isinstance(frame.props[1].eps_t_Vap_Liq, Param) + assert value(frame.props[1].eps_t_Vap_Liq) == pytest.approx(1, rel=1e-8) + assert isinstance(frame.props[1].eps_z_Vap_Liq, Param) + assert value(frame.props[1].eps_z_Vap_Liq) == pytest.approx(1, rel=1e-8) + + assert isinstance(frame.props[1]._teq_constraint_Vap_Liq, Constraint) + assert isinstance( + frame.props[1].temperature_slack_complementarity_Vap_Liq, Constraint + ) + assert isinstance(frame.props[1].cubic_slack_complementarity_Vap_Liq, Constraint) + + for k in ["Liq", "Vap"]: + assert k in frame.props[1].s_Vap_Liq + assert k in frame.props[1].gp_Vap_Liq + assert k in frame.props[1].gn_Vap_Liq + assert k in frame.props[1].cubic_second_derivative_Vap_Liq + assert k in frame.props[1].temperature_slack_complementarity_Vap_Liq + assert k in frame.props[1].cubic_slack_complementarity_Vap_Liq + + +# TODO: Need tests for formulation + + +@pytest.mark.unit +def test_cubic_second_derivative(frame): + Z = frame.props[1].compress_fact_phase + + # For pure water + R = 8.314 + b = 0.07780 * R * 647 / 220.6e5 + + for P in range(1, 11): + for T in range(300, 500, 10): + frame.props[1].pressure.set_value(P * 1e5) + frame.props[1].temperature.set_value(T) + frame.props[1]._teq[("Vap", "Liq")].set_value(T) + + B = b * P * 1e5 / R / T + + for p in ["Vap", "Liq"]: + der = value(6 * Z[p] + 2 * -(1 + B - 2 * B)) + assert value( + frame.props[1].cubic_second_derivative_Vap_Liq[p] + ) == pytest.approx(der, rel=1e-8) + + +@pytest.mark.unit +def test_calculate_temperature_slacks(frame): + s = frame.props[1].s_Vap_Liq + # Teq > T + frame.props[1].temperature.set_value(300) + frame.props[1]._teq[("Vap", "Liq")].set_value(400) + _calculate_temperature_slacks(frame.props[1], ("Vap", "Liq"), "Liq", "Vap") + + assert value(s["Vap"]) == pytest.approx(100, rel=1e-8) + assert value(s["Liq"]) == pytest.approx(EPS_INIT, abs=1e-8) + + # Teq < T + frame.props[1]._teq[("Vap", "Liq")].set_value(200) + _calculate_temperature_slacks(frame.props[1], ("Vap", "Liq"), "Liq", "Vap") + + assert value(s["Liq"]) == pytest.approx(100, rel=1e-8) + assert value(s["Vap"]) == pytest.approx(EPS_INIT, abs=1e-8) + + # Teq == T + frame.props[1]._teq[("Vap", "Liq")].set_value(300) + _calculate_temperature_slacks(frame.props[1], ("Vap", "Liq"), "Liq", "Vap") + + assert value(s["Liq"]) == pytest.approx(EPS_INIT, abs=1e-8) + assert value(s["Vap"]) == pytest.approx(EPS_INIT, abs=1e-8) + + +@pytest.mark.unit +def test_calculate_ceos_derivative_slacks(frame): + gp = frame.props[1].gp_Vap_Liq + gn = frame.props[1].gn_Vap_Liq + der = frame.props[1].cubic_second_derivative_Vap_Liq + Z = frame.props[1].compress_fact_phase + + # For pure water + R = 8.314 + b = 0.07780 * R * 647 / 220.6e5 + + # Atmospheric conditions, vapor derivative +ve, liquid -ve + frame.props[1].pressure.set_value(1e5) + frame.props[1].temperature.set_value(300) + frame.props[1]._teq[("Vap", "Liq")].set_value(300) + + _calculate_ceos_derivative_slacks(frame.props[1], ("Vap", "Liq"), "Liq", "Vap") + + B = value(b * frame.props[1].pressure / R / frame.props[1].temperature) + der_l = value(6 * Z["Liq"] + 2 * -(1 + B - 2 * B)) + der_v = value(6 * Z["Vap"] + 2 * -(1 + B - 2 * B)) + + assert value(gp["Liq"]) == pytest.approx(EPS_INIT, abs=1e-8) + assert value(gn["Liq"]) == pytest.approx(-der_l, rel=1e-8) + + assert value(gp["Vap"]) == pytest.approx(der_v, rel=1e-8) + assert value(gn["Vap"]) == pytest.approx(EPS_INIT, abs=1e-8) + + # Supercritical conditions, vapor derivative +ve, liquid +ve + frame.props[1].pressure.set_value(250e5) + frame.props[1].temperature.set_value(700) + frame.props[1]._teq[("Vap", "Liq")].set_value(700) + + _calculate_ceos_derivative_slacks(frame.props[1], ("Vap", "Liq"), "Liq", "Vap") + + B = value(b * frame.props[1].pressure / R / frame.props[1].temperature) + der_v = value(6 * Z["Vap"] + 2 * -(1 + B - 2 * B)) + + assert value(gp["Liq"]) == pytest.approx(EPS_INIT, abs=1e-8) + assert value(gn["Liq"]) == pytest.approx(EPS_INIT, abs=1e-8) + + assert value(gp["Vap"]) == pytest.approx(der_v, rel=1e-8) + assert value(gn["Vap"]) == pytest.approx(EPS_INIT, abs=1e-8) + + # Supercritical conditions, vapor derivative -ve, liquid -ve + frame.props[1].pressure.set_value(250e5) + frame.props[1].temperature.set_value(300) + frame.props[1]._teq[("Vap", "Liq")].set_value(300) + + _calculate_ceos_derivative_slacks(frame.props[1], ("Vap", "Liq"), "Liq", "Vap") + + B = value(b * frame.props[1].pressure / R / frame.props[1].temperature) + der_l = value(6 * Z["Liq"] + 2 * -(1 + B - 2 * B)) + + assert value(gp["Liq"]) == pytest.approx(EPS_INIT, abs=1e-8) + assert value(gn["Liq"]) == pytest.approx(-der_l, rel=1e-8) + + assert value(gp["Vap"]) == pytest.approx(EPS_INIT, abs=1e-8) + assert value(gn["Vap"]) == pytest.approx(EPS_INIT, abs=1e-8) diff --git a/idaes/models/properties/modular_properties/state_definitions/FTPx.py b/idaes/models/properties/modular_properties/state_definitions/FTPx.py index f60e3c52b9..8108928f7b 100644 --- a/idaes/models/properties/modular_properties/state_definitions/FTPx.py +++ b/idaes/models/properties/modular_properties/state_definitions/FTPx.py @@ -36,8 +36,8 @@ get_method, GenericPropertyPackageError, ) -from idaes.models.properties.modular_properties.phase_equil.bubble_dew import ( - _valid_VL_component_list, +from idaes.models.properties.modular_properties.base.utility import ( + identify_VL_component_list, ) from idaes.models.properties.modular_properties.phase_equil.henry import ( HenryType, @@ -420,7 +420,7 @@ def state_initialization(b): henry_comps, l_only_comps, v_only_comps, - ) = _valid_VL_component_list(b, pp) + ) = identify_VL_component_list(b, pp) pp_VLE = pp if init_VLE: diff --git a/idaes/models/unit_models/tests/test_heat_exchanger.py b/idaes/models/unit_models/tests/test_heat_exchanger.py index 1bcf4efd9e..ce671036cd 100644 --- a/idaes/models/unit_models/tests/test_heat_exchanger.py +++ b/idaes/models/unit_models/tests/test_heat_exchanger.py @@ -15,6 +15,7 @@ Author: John Eslick """ +from copy import deepcopy import pytest import pandas @@ -1682,6 +1683,13 @@ def btx(self): m.fs.unit.cold_side.scaling_factor_pressure = 1 + # Set small values of epsilon to get sufficiently accurate results + # Only applies to hot side, as cold side used the original SmoothVLE. + m.fs.unit.hot_side.properties_in[0].eps_t_Vap_Liq.set_value(1e-4) + m.fs.unit.hot_side.properties_in[0].eps_z_Vap_Liq.set_value(1e-4) + m.fs.unit.hot_side.properties_out[0].eps_t_Vap_Liq.set_value(1e-4) + m.fs.unit.hot_side.properties_out[0].eps_z_Vap_Liq.set_value(1e-4) + return m @pytest.mark.build @@ -1725,8 +1733,8 @@ def test_build(self, btx): assert isinstance(btx.fs.unit.delta_temperature, (Var, Expression)) assert isinstance(btx.fs.unit.heat_transfer_equation, Constraint) - assert number_variables(btx) == 190 - assert number_total_constraints(btx) == 118 + assert number_variables(btx) == 176 + assert number_total_constraints(btx) == 104 assert number_unused_variables(btx) == 20 @pytest.mark.component @@ -1893,7 +1901,9 @@ def test_conservation(self, btx): @pytest.mark.component def test_numerical_issues(self, btx): dt = DiagnosticsToolbox(btx) - dt.assert_no_numerical_warnings() + # TODO: Complementarity formulation results in near-parallel components + # when unscaled + dt.assert_no_numerical_warnings(ignore_parallel_components=True) @pytest.mark.component def test_initialization_error(self, btx): @@ -1903,158 +1913,162 @@ def test_initialization_error(self, btx): btx.fs.unit.initialize() -class TestInitializersModular: - @pytest.fixture - def model(self): - m = ConcreteModel() - m.fs = FlowsheetBlock(dynamic=False) - - # As we lack other example prop packs with units, take the generic - # BT-PR package and change the base units - configuration2 = { - # Specifying components - "components": { - "benzene": { - "type": Component, - "enth_mol_ig_comp": RPP, - "entr_mol_ig_comp": RPP, - "pressure_sat_comp": RPP, - "phase_equilibrium_form": {("Vap", "Liq"): log_fugacity}, - "parameter_data": { - "mw": (78.1136e-3, pyunits.kg / pyunits.mol), # [1] - "pressure_crit": (48.9e5, pyunits.Pa), # [1] - "temperature_crit": (562.2, pyunits.K), # [1] - "omega": 0.212, # [1] - "cp_mol_ig_comp_coeff": { - "A": (-3.392e1, pyunits.J / pyunits.mol / pyunits.K), # [1] - "B": (4.739e-1, pyunits.J / pyunits.mol / pyunits.K**2), - "C": (-3.017e-4, pyunits.J / pyunits.mol / pyunits.K**3), - "D": (7.130e-8, pyunits.J / pyunits.mol / pyunits.K**4), - }, - "enth_mol_form_vap_comp_ref": ( - 82.9e3, - pyunits.J / pyunits.mol, - ), # [3] - "entr_mol_form_vap_comp_ref": ( - -269, - pyunits.J / pyunits.mol / pyunits.K, - ), # [3] - "pressure_sat_comp_coeff": { - "A": (-6.98273, None), # [1] - "B": (1.33213, None), - "C": (-2.62863, None), - "D": (-3.33399, None), - }, - }, +# As we lack other example prop packs with units, take the generic +# BT-PR package and change the base units +configuration2 = { + # Specifying components + "components": { + "benzene": { + "type": Component, + "enth_mol_ig_comp": RPP, + "entr_mol_ig_comp": RPP, + "pressure_sat_comp": RPP, + "phase_equilibrium_form": {("Vap", "Liq"): log_fugacity}, + "parameter_data": { + "mw": (78.1136e-3, pyunits.kg / pyunits.mol), # [1] + "pressure_crit": (48.9e5, pyunits.Pa), # [1] + "temperature_crit": (562.2, pyunits.K), # [1] + "omega": 0.212, # [1] + "cp_mol_ig_comp_coeff": { + "A": (-3.392e1, pyunits.J / pyunits.mol / pyunits.K), # [1] + "B": (4.739e-1, pyunits.J / pyunits.mol / pyunits.K**2), + "C": (-3.017e-4, pyunits.J / pyunits.mol / pyunits.K**3), + "D": (7.130e-8, pyunits.J / pyunits.mol / pyunits.K**4), }, - "toluene": { - "type": Component, - "enth_mol_ig_comp": RPP, - "entr_mol_ig_comp": RPP, - "pressure_sat_comp": RPP, - "phase_equilibrium_form": {("Vap", "Liq"): log_fugacity}, - "parameter_data": { - "mw": (92.1405e-3, pyunits.kg / pyunits.mol), # [1] - "pressure_crit": (41e5, pyunits.Pa), # [1] - "temperature_crit": (591.8, pyunits.K), # [1] - "omega": 0.263, # [1] - "cp_mol_ig_comp_coeff": { - "A": (-2.435e1, pyunits.J / pyunits.mol / pyunits.K), # [1] - "B": (5.125e-1, pyunits.J / pyunits.mol / pyunits.K**2), - "C": (-2.765e-4, pyunits.J / pyunits.mol / pyunits.K**3), - "D": (4.911e-8, pyunits.J / pyunits.mol / pyunits.K**4), - }, - "enth_mol_form_vap_comp_ref": ( - 50.1e3, - pyunits.J / pyunits.mol, - ), # [3] - "entr_mol_form_vap_comp_ref": ( - -321, - pyunits.J / pyunits.mol / pyunits.K, - ), # [3] - "pressure_sat_comp_coeff": { - "A": (-7.28607, None), # [1] - "B": (1.38091, None), - "C": (-2.83433, None), - "D": (-2.79168, None), - }, - }, + "enth_mol_form_vap_comp_ref": ( + 82.9e3, + pyunits.J / pyunits.mol, + ), # [3] + "entr_mol_form_vap_comp_ref": ( + -269, + pyunits.J / pyunits.mol / pyunits.K, + ), # [3] + "pressure_sat_comp_coeff": { + "A": (-6.98273, None), # [1] + "B": (1.33213, None), + "C": (-2.62863, None), + "D": (-3.33399, None), }, }, - # Specifying phases - "phases": { - "Liq": { - "type": LiquidPhase, - "equation_of_state": Cubic, - "equation_of_state_options": {"type": CubicType.PR}, + }, + "toluene": { + "type": Component, + "enth_mol_ig_comp": RPP, + "entr_mol_ig_comp": RPP, + "pressure_sat_comp": RPP, + "phase_equilibrium_form": {("Vap", "Liq"): log_fugacity}, + "parameter_data": { + "mw": (92.1405e-3, pyunits.kg / pyunits.mol), # [1] + "pressure_crit": (41e5, pyunits.Pa), # [1] + "temperature_crit": (591.8, pyunits.K), # [1] + "omega": 0.263, # [1] + "cp_mol_ig_comp_coeff": { + "A": (-2.435e1, pyunits.J / pyunits.mol / pyunits.K), # [1] + "B": (5.125e-1, pyunits.J / pyunits.mol / pyunits.K**2), + "C": (-2.765e-4, pyunits.J / pyunits.mol / pyunits.K**3), + "D": (4.911e-8, pyunits.J / pyunits.mol / pyunits.K**4), }, - "Vap": { - "type": VaporPhase, - "equation_of_state": Cubic, - "equation_of_state_options": {"type": CubicType.PR}, + "enth_mol_form_vap_comp_ref": ( + 50.1e3, + pyunits.J / pyunits.mol, + ), # [3] + "entr_mol_form_vap_comp_ref": ( + -321, + pyunits.J / pyunits.mol / pyunits.K, + ), # [3] + "pressure_sat_comp_coeff": { + "A": (-7.28607, None), # [1] + "B": (1.38091, None), + "C": (-2.83433, None), + "D": (-2.79168, None), }, }, - # Set base units of measurement - "base_units": { - "time": pyunits.s, - "length": pyunits.m, - "mass": pyunits.t, - "amount": pyunits.mol, - "temperature": pyunits.degR, - }, - # 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": { - ("benzene", "benzene"): 0.000, - ("benzene", "toluene"): 0.000, - ("toluene", "benzene"): 0.000, - ("toluene", "toluene"): 0.000, - } - }, + }, + }, + # 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.t, + "amount": pyunits.mol, + "temperature": pyunits.degR, + }, + # 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": { + ("benzene", "benzene"): 0.000, + ("benzene", "toluene"): 0.000, + ("toluene", "benzene"): 0.000, + ("toluene", "toluene"): 0.000, } + }, +} - m.fs.properties = GenericParameterBlock(**configuration) - m.fs.properties2 = GenericParameterBlock(**configuration2) - m.fs.unit = HeatExchanger( - hot_side={"property_package": m.fs.properties}, - cold_side={"property_package": m.fs.properties2}, +class TestInitializersModular: + @pytest.mark.integration + def test_hx0d_initializer(self): + model = ConcreteModel() + model.fs = FlowsheetBlock(dynamic=False) + + model.fs.properties = GenericParameterBlock(**configuration) + model.fs.properties2 = GenericParameterBlock(**configuration2) + + model.fs.unit = HeatExchanger( + hot_side={"property_package": model.fs.properties}, + cold_side={"property_package": model.fs.properties2}, flow_pattern=HeatExchangerFlowPattern.cocurrent, ) - m.fs.unit.hot_side_inlet.flow_mol[0].fix(5) # mol/s - m.fs.unit.hot_side_inlet.temperature[0].fix(365) # K - m.fs.unit.hot_side_inlet.pressure[0].fix(101325) # Pa - m.fs.unit.hot_side_inlet.mole_frac_comp[0, "benzene"].fix(0.5) - m.fs.unit.hot_side_inlet.mole_frac_comp[0, "toluene"].fix(0.5) + model.fs.unit.hot_side_inlet.flow_mol[0].fix(5) # mol/s + model.fs.unit.hot_side_inlet.temperature[0].fix(365) # K + model.fs.unit.hot_side_inlet.pressure[0].fix(101325) # Pa + model.fs.unit.hot_side_inlet.mole_frac_comp[0, "benzene"].fix(0.5) + model.fs.unit.hot_side_inlet.mole_frac_comp[0, "toluene"].fix(0.5) - m.fs.unit.cold_side_inlet.flow_mol[0].fix(1) # mol/s - m.fs.unit.cold_side_inlet.temperature[0].fix(540) # degR - m.fs.unit.cold_side_inlet.pressure[0].fix(101.325) # kPa - m.fs.unit.cold_side_inlet.mole_frac_comp[0, "benzene"].fix(0.5) - m.fs.unit.cold_side_inlet.mole_frac_comp[0, "toluene"].fix(0.5) + model.fs.unit.cold_side_inlet.flow_mol[0].fix(1) # mol/s + model.fs.unit.cold_side_inlet.temperature[0].fix(540) # degR + model.fs.unit.cold_side_inlet.pressure[0].fix(101.325) # kPa + model.fs.unit.cold_side_inlet.mole_frac_comp[0, "benzene"].fix(0.5) + model.fs.unit.cold_side_inlet.mole_frac_comp[0, "toluene"].fix(0.5) - m.fs.unit.area.fix(1) - m.fs.unit.overall_heat_transfer_coefficient.fix(100) + model.fs.unit.area.fix(1) + model.fs.unit.overall_heat_transfer_coefficient.fix(100) - m.fs.unit.cold_side.scaling_factor_pressure = 1 + model.fs.unit.cold_side.scaling_factor_pressure = 1 - return m + # Set small values of epsilon to get sufficiently accurate results + # Only applies to hot side, as cold side used the original SmoothVLE. + model.fs.unit.hot_side.properties_in[0].eps_t_Vap_Liq.set_value(1e-4) + model.fs.unit.hot_side.properties_in[0].eps_z_Vap_Liq.set_value(1e-4) + model.fs.unit.hot_side.properties_out[0].eps_t_Vap_Liq.set_value(1e-4) + model.fs.unit.hot_side.properties_out[0].eps_z_Vap_Liq.set_value(1e-4) - @pytest.mark.integration - def test_hx0d_initializer(self, model): initializer = HX0DInitializer() initializer.initialize(model.fs.unit) @@ -2081,8 +2095,46 @@ def test_hx0d_initializer(self, model): ) @pytest.mark.integration - def test_block_triangularization(self, model): - initializer = BlockTriangularizationInitializer(constraint_tolerance=2e-5) + def test_block_triangularization( + self, + ): + # Trying to get this to work with CubicComplementarityVLE is challenging, and + # not necessary for this particular test + new_config = deepcopy(configuration) + new_config["phase_equilibrium_state"] = {("Vap", "Liq"): SmoothVLE} + + model = ConcreteModel() + model.fs = FlowsheetBlock(dynamic=False) + + model.fs.properties = GenericParameterBlock(**new_config) + model.fs.properties2 = GenericParameterBlock(**configuration2) + + model.fs.unit = HeatExchanger( + hot_side={"property_package": model.fs.properties}, + cold_side={"property_package": model.fs.properties2}, + flow_pattern=HeatExchangerFlowPattern.cocurrent, + ) + + model.fs.unit.hot_side_inlet.flow_mol[0].fix(5) # mol/s + model.fs.unit.hot_side_inlet.temperature[0].fix(365) # K + model.fs.unit.hot_side_inlet.pressure[0].fix(101325) # Pa + model.fs.unit.hot_side_inlet.mole_frac_comp[0, "benzene"].fix(0.5) + model.fs.unit.hot_side_inlet.mole_frac_comp[0, "toluene"].fix(0.5) + + model.fs.unit.cold_side_inlet.flow_mol[0].fix(1) # mol/s + model.fs.unit.cold_side_inlet.temperature[0].fix(540) # degR + model.fs.unit.cold_side_inlet.pressure[0].fix(101.325) # kPa + model.fs.unit.cold_side_inlet.mole_frac_comp[0, "benzene"].fix(0.5) + model.fs.unit.cold_side_inlet.mole_frac_comp[0, "toluene"].fix(0.5) + + model.fs.unit.area.fix(1) + model.fs.unit.overall_heat_transfer_coefficient.fix(100) + + model.fs.unit.cold_side.scaling_factor_pressure = 1 + + initializer = BlockTriangularizationInitializer( + constraint_tolerance=2e-5, + ) initializer.initialize(model.fs.unit) assert initializer.summary[model.fs.unit]["status"] == InitializationStatus.Ok diff --git a/idaes/models/unit_models/tests/test_heat_exchanger_1D.py b/idaes/models/unit_models/tests/test_heat_exchanger_1D.py index 2724dc6fc3..7cbd2d1753 100644 --- a/idaes/models/unit_models/tests/test_heat_exchanger_1D.py +++ b/idaes/models/unit_models/tests/test_heat_exchanger_1D.py @@ -2448,6 +2448,12 @@ def btx(self): m.fs.unit.cold_side_inlet.mole_frac_comp[0, "benzene"].fix(0.5) m.fs.unit.cold_side_inlet.mole_frac_comp[0, "toluene"].fix(0.5) + # Set small values of epsilon to get sufficiently accurate results + # Only need hot side, as cold side uses old SmoothVLE + for i in m.fs.unit.hot_side.properties.keys(): + m.fs.unit.hot_side.properties[i].eps_t_Vap_Liq.set_value(1e-4) + m.fs.unit.hot_side.properties[i].eps_z_Vap_Liq.set_value(1e-4) + return m @pytest.mark.component @@ -2486,8 +2492,8 @@ def test_build(self, btx): assert hasattr(btx.fs.unit, "heat_transfer_eq") assert hasattr(btx.fs.unit, "heat_conservation") - assert number_variables(btx) == 1976 - assert number_total_constraints(btx) == 1867 + assert number_variables(btx) == 1829 + assert number_total_constraints(btx) == 1720 assert number_unused_variables(btx) == 36 @pytest.mark.integration @@ -2641,7 +2647,9 @@ def test_conservation(self, btx): @pytest.mark.integration def test_numerical_issues(self, btx): dt = DiagnosticsToolbox(btx) - dt.assert_no_numerical_warnings() + # TODO: Complementarity formulation results in near-parallel components + # when unscaled + dt.assert_no_numerical_warnings(ignore_parallel_components=True) @pytest.mark.component def test_initialization_error(self, btx): @@ -3529,6 +3537,11 @@ def model(self): m.fs.unit.cold_side_inlet.mole_frac_comp[0, "benzene"].set_value(0.5) m.fs.unit.cold_side_inlet.mole_frac_comp[0, "toluene"].set_value(0.5) + # Set small values of epsilon to get sufficiently accurate results + for i in m.fs.unit.hot_side.properties.keys(): + m.fs.unit.hot_side.properties[i].eps_t_Vap_Liq.set_value(1e-4) + m.fs.unit.hot_side.properties[i].eps_z_Vap_Liq.set_value(1e-4) + return m @pytest.mark.component diff --git a/idaes/models/unit_models/tests/test_heater.py b/idaes/models/unit_models/tests/test_heater.py index dbee5e970b..6f67b13f11 100644 --- a/idaes/models/unit_models/tests/test_heater.py +++ b/idaes/models/unit_models/tests/test_heater.py @@ -511,6 +511,12 @@ def btg(self): m.fs.unit.heat_duty.fix(-5000) m.fs.unit.deltaP.fix(0) + # Set small values of epsilon to get sufficiently accurate results + m.fs.unit.control_volume.properties_in[0].eps_t_Vap_Liq.set_value(1e-4) + m.fs.unit.control_volume.properties_in[0].eps_z_Vap_Liq.set_value(1e-4) + m.fs.unit.control_volume.properties_out[0].eps_t_Vap_Liq.set_value(1e-4) + m.fs.unit.control_volume.properties_out[0].eps_z_Vap_Liq.set_value(1e-4) + return m @pytest.mark.build @@ -533,8 +539,8 @@ def test_build(self, btg): assert hasattr(btg.fs.unit, "heat_duty") assert hasattr(btg.fs.unit, "deltaP") - assert number_variables(btg) == 94 - assert number_total_constraints(btg) == 57 + assert number_variables(btg) == 80 + assert number_total_constraints(btg) == 43 # Unused vars are density parameters assert number_unused_variables(btg) == 10 @@ -595,7 +601,9 @@ def test_conservation(self, btg): @pytest.mark.component def test_numerical_issues(self, btg): dt = DiagnosticsToolbox(btg) - dt.assert_no_numerical_warnings() + # TODO: Complementarity formulation results in near-parallel components + # when unscaled + dt.assert_no_numerical_warnings(ignore_parallel_components=True) @pytest.mark.ui @pytest.mark.unit @@ -665,6 +673,13 @@ def model(self): m.fs.unit.heat_duty.fix(-5000) m.fs.unit.deltaP.fix(0) + # Set small values of epsilon to get sufficiently accurate results + m.fs.unit.control_volume.properties_in.display() + m.fs.unit.control_volume.properties_in[0].eps_t_Vap_Liq.set_value(1e-4) + m.fs.unit.control_volume.properties_in[0].eps_z_Vap_Liq.set_value(1e-4) + m.fs.unit.control_volume.properties_out[0].eps_t_Vap_Liq.set_value(1e-4) + m.fs.unit.control_volume.properties_out[0].eps_z_Vap_Liq.set_value(1e-4) + return m @pytest.mark.integration diff --git a/idaes/models/unit_models/tests/test_shell_and_tube_1D.py b/idaes/models/unit_models/tests/test_shell_and_tube_1D.py index 362e267313..bdb087e3cb 100644 --- a/idaes/models/unit_models/tests/test_shell_and_tube_1D.py +++ b/idaes/models/unit_models/tests/test_shell_and_tube_1D.py @@ -1420,6 +1420,12 @@ def btx(self): m.fs.unit.cold_side_inlet.mole_frac_comp[0, "benzene"].fix(0.5) m.fs.unit.cold_side_inlet.mole_frac_comp[0, "toluene"].fix(0.5) + # Set small values of epsilon to get sufficiently accurate results + # Only need hot side, as cold side uses old SmoothVLE + for i in m.fs.unit.hot_side.properties.keys(): + m.fs.unit.hot_side.properties[i].eps_t_Vap_Liq.set_value(1e-4) + m.fs.unit.hot_side.properties[i].eps_z_Vap_Liq.set_value(1e-4) + return m @pytest.mark.component @@ -1467,8 +1473,8 @@ def test_build(self, btx): assert hasattr(btx.fs.unit, "cold_side_heat_transfer_eq") assert hasattr(btx.fs.unit, "heat_conservation") - assert number_variables(btx) == 2021 - assert number_total_constraints(btx) == 1890 + assert number_variables(btx) == 1874 + assert number_total_constraints(btx) == 1743 assert number_unused_variables(btx) == 34 @pytest.mark.integration @@ -1625,7 +1631,9 @@ def test_conservation(self, btx): @pytest.mark.integration def test_numerical_issues(self, btx): dt = DiagnosticsToolbox(btx) - dt.assert_no_numerical_warnings() + # TODO: Complementarity formulation results in near-parallel components + # when unscaled + dt.assert_no_numerical_warnings(ignore_parallel_components=True) @pytest.mark.component def test_initialization_error(self, btx):