Hide code cell content
###############################################################################
# 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.
###############################################################################

Create a dictionary and use Txy diagrams#

Author: Alejandro Garciadego
Maintainer: Andrew Lee
Updated: 2023-06-01

1. Introduction#

This notebook will do two essential things. The first is to create a property package dictionary for the Modular Property Package Framework and utilized the T-x-y diagrams for the properties. This dictionary will include parameters necessary for the sub-libraries of the property package. For example, it has critical properties required for vapor-liquid equilibrium and heat capacity parameters to calculate enthalpy and entropy. The property parameter dictionary can be later used to generate models and flowsheets.

Secondly, we will utilize this dictionary to create a model and print a T-x-y diagram using IDAES.

1.1 Tutorial objectives#

The goals of this tutorial are:

  • Create a parameter block directly on a jupyter notebook.

  • Utilize the Modular Property Package, which provides a flexible platform for users to build property packages by calling upon libraries of modular sub-models to build up complex property calculations with the least effort possible.

  • Demonstrate the use of Txy diagrams in the IDAES utilities.

2. Problem Statement#

2.1 Setting up the problem in IDAES#

In the next cell, we will be importing the necessary components from Pyomo and IDAES.

# Import objects from pyomo package
from pyomo.environ import ConcreteModel

# Import idaes logger to set output levels
import idaes.logger as idaeslog

# Import the Modular Parameter Block
from idaes.models.properties.modular_properties import GenericParameterBlock

2.2 Importing property sub-libraries#

The IDAES Modular Property Package calls upon libraries of modular sub-models to build up complex property calculations with the least effort possible. We have to import these libraries.

# Import Python libraries
import logging

# Import Pyomo units
from pyomo.environ import units as pyunits

# Import IDAES cores
from idaes.core import LiquidPhase, VaporPhase, Component

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
from idaes.models.properties.modular_properties.phase_equil.bubble_dew import (
    IdealBubbleDew,
)
from idaes.models.properties.modular_properties.phase_equil.forms import fugacity

from idaes.models.properties.modular_properties.pure import Perrys
from idaes.models.properties.modular_properties.pure import RPP4

2.3 Creating a Modular Property Parameter block#

The next step is to give all the parameter information required for the sub-libraries. Here we will specify the base units of our model, the bounds on the state properties and the reference conditions of our properties.

configuration = {
    # Specifying components
    "components": {
        "benzene": {
            "type": Component,
            "elemental_composition": {"C": 6, "H": 6},
            "dens_mol_liq_comp": Perrys,
            "enth_mol_liq_comp": Perrys,
            "enth_mol_ig_comp": RPP4,
            "pressure_sat_comp": RPP4,
            "phase_equilibrium_form": {("Vap", "Liq"): 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]
                "dens_mol_liq_comp_coeff": {
                    "eqn_type": 1,
                    "1": (1.0162, pyunits.kmol * pyunits.m**-3),  # [2] pg. 2-98
                    "2": (0.2655, None),
                    "3": (562.16, pyunits.K),
                    "4": (0.28212, None),
                },
                "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),
                },
                "cp_mol_liq_comp_coeff": {
                    "1": (1.29e2, pyunits.J / pyunits.kmol / pyunits.K),  # [2]
                    "2": (-1.7e-1, pyunits.J / pyunits.kmol / pyunits.K**2),
                    "3": (6.48e-4, pyunits.J / pyunits.kmol / pyunits.K**3),
                    "4": (0, pyunits.J / pyunits.kmol / pyunits.K**4),
                    "5": (0, pyunits.J / pyunits.kmol / pyunits.K**5),
                },
                "enth_mol_form_liq_comp_ref": (49.0e3, pyunits.J / pyunits.mol),  # [3]
                "enth_mol_form_vap_comp_ref": (82.9e3, pyunits.J / pyunits.mol),  # [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,
            "elemental_composition": {"C": 7, "H": 8},
            "dens_mol_liq_comp": Perrys,
            "enth_mol_liq_comp": Perrys,
            "enth_mol_ig_comp": RPP4,
            "pressure_sat_comp": RPP4,
            "phase_equilibrium_form": {("Vap", "Liq"): 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]
                "dens_mol_liq_comp_coeff": {
                    "eqn_type": 1,
                    "1": (0.8488, pyunits.kmol * pyunits.m**-3),  # [2] pg. 2-98
                    "2": (0.26655, None),
                    "3": (591.8, pyunits.K),
                    "4": (0.2878, None),
                },
                "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),
                },
                "cp_mol_liq_comp_coeff": {
                    "1": (1.40e2, pyunits.J / pyunits.kmol / pyunits.K),  # [2]
                    "2": (-1.52e-1, pyunits.J / pyunits.kmol / pyunits.K**2),
                    "3": (6.95e-4, pyunits.J / pyunits.kmol / pyunits.K**3),
                    "4": (0, pyunits.J / pyunits.kmol / pyunits.K**4),
                    "5": (0, pyunits.J / pyunits.kmol / pyunits.K**5),
                },
                "enth_mol_form_liq_comp_ref": (12.0e3, pyunits.J / pyunits.mol),  # [3]
                "enth_mol_form_vap_comp_ref": (50.1e3, pyunits.J / pyunits.mol),  # [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": Ideal},
        "Vap": {"type": VaporPhase, "equation_of_state": Ideal},
    },
    # 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": (150, 300, 450, pyunits.K),
        "pressure": (5e4, 1e5, 1e6, pyunits.Pa),
    },
    "pressure_ref": (1e5, pyunits.Pa),
    "temperature_ref": (300, pyunits.K),
    # Defining phase equilibria
    "phases_in_equilibrium": [("Vap", "Liq")],
    "phase_equilibrium_state": {("Vap", "Liq"): SmoothVLE},
    "bubble_dew_method": IdealBubbleDew,
    "default_scaling_factors": {"mole_frac_comp[benzene]": 1},
}

2.4 Building the model#

In the next cell, we will first create a model and attach a the property package.

# Create the ConcreteModel
model = ConcreteModel()

# Attach the property package to the model
model.params = GenericParameterBlock(**configuration)

3.0 Creating Txy diagrams#

3.1 Import the TXY_diagram function#

from idaes.core.util.phase_equilibria import Txy_diagram

3.2 Calling the function#

There will be a WARNING from missing scaling factors. This is normal as the Txy diagram function does not know the state definition being used. To learn more about the sacling factor be sure to reed the documentation/examples of the scaling tools.

Txy_diagram(
    model,
    "benzene",
    "toluene",
    101325,
    num_points=25,
    temperature=150.15,
    figure_name=None,
    print_legend=True,
    include_pressure=False,
    print_level=idaeslog.CRITICAL,
    solver_op={"tol": 1e-6},
)
2024-03-18 23:14:14 [WARNING] idaes.core.util.scaling: Missing scaling factor for props[1].mole_frac_comp
2024-03-18 23:14:14 [WARNING] idaes.core.util.scaling: Missing scaling factor for props[1].mole_frac_comp
---------------------------------------------------------------------------
ApplicationError                          Traceback (most recent call last)
Cell In[7], line 1
----> 1 Txy_diagram(
      2     model,
      3     "benzene",
      4     "toluene",
      5     101325,
      6     num_points=25,
      7     temperature=150.15,
      8     figure_name=None,
      9     print_legend=True,
     10     include_pressure=False,
     11     print_level=idaeslog.CRITICAL,
     12     solver_op={"tol": 1e-6},
     13 )

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/util/phase_equilibria.py:77, in Txy_diagram(model, component_1, component_2, pressure, num_points, temperature, figure_name, print_legend, include_pressure, print_level, solver, solver_op)
     50 """
     51 This method generates T-x-y plots. Given the components, pressure and property dictionary
     52 this function calls Txy_data() to generate T-x-y data and once the data has
   (...)
     74     Plot
     75 """
     76 # Run txy_ data function to obtain bubble and dew twmperatures
---> 77 Txy_data_to_plot = Txy_data(
     78     model,
     79     component_1,
     80     component_2,
     81     pressure,
     82     num_points,
     83     temperature,
     84     print_level,
     85     solver,
     86     solver_op,
     87 )
     89 # Run diagrams function to convert t-x-y data into a plot
     90 build_txy_diagrams(Txy_data_to_plot, figure_name, print_legend, include_pressure)

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/util/phase_equilibria.py:154, in Txy_data(model, component_1, component_2, pressure, num_points, temperature, print_level, solver, solver_op)
    152 # Initialize flash unit model
    153 model.props[1].calculate_scaling_factors()
--> 154 model.props.initialize(solver=solver, optarg=solver_op, outlvl=print_level)
    156 solver = get_solver(solver, solver_op)
    158 # Create an array of compositions with N number of points

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/modular_properties/base/generic_property.py:1824, in _GenericStateBlock.initialize(blk, state_args, state_vars_fixed, hold_state, outlvl, solver, optarg)
   1819         raise InitializationError(
   1820             f"{blk.name} Unexpected degrees of freedom during "
   1821             f"initialization at bubble, dew, and critical point step: {dof}."
   1822         )
   1823     with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
-> 1824         res = solve_indexed_blocks(opt, [blk], tee=slc.tee)
   1825     init_log.info(
   1826         "Bubble, dew, and critical point initialization: {}.".format(
   1827             idaeslog.condition(res)
   1828         )
   1829     )
   1830 # ---------------------------------------------------------------------
   1831 # Calculate _teq if required
   1832 # Using iterator k outside of for loop - this should be OK as we just need
   1833 # a valid StateBlockData an assume they are all the same.

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/util/initialization.py:285, in solve_indexed_blocks(solver, blocks, **kwds)
    278     tmp._ctypes[Block] = [  # pylint: disable=protected-access
    279         0,
    280         nBlocks - 1,
    281         nBlocks,
    282     ]
    284     # Solve temporary Block
--> 285     results = solver.solve(tmp, **kwds)
    287 finally:
    288     # Clean up temporary Block contents so they are not removed when Block
    289     # is garbage collected.
    290     tmp._decl = {}  # pylint: disable=protected-access

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/opt/base/solvers.py:534, in OptSolver.solve(self, *args, **kwds)
    531 def solve(self, *args, **kwds):
    532     """Solve the problem"""
--> 534     self.available(exception_flag=True)
    535     #
    536     # If the inputs are models, then validate that they have been
    537     # constructed! Collect suffix names to try and import from solution.
    538     #
    539     from pyomo.core.base.block import _BlockData

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/opt/solver/shellcmd.py:139, in SystemCallSolver.available(self, exception_flag)
    137     if exception_flag:
    138         msg = "No executable found for solver '%s'"
--> 139         raise ApplicationError(msg % self.name)
    140     return False
    141 return True

ApplicationError: No executable found for solver 'ipopt'

You can add precision to the plot by increasing the num_points, include a firgure_name in order to save figure in a file or include the pressure in the plot. You can also include more information from the solver by changing the print_level.