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.
###############################################################################

Flash Unit Model#

Author: Jaffer Ghouse
Maintainer: Andrew Lee
Updated: 2023-06-01

In this module, we will familiarize ourselves with the IDAES framework by creating and working with a flowsheet that contains a single flash tank. The flash tank will be used to perform separation of Benzene and Toluene. The inlet specifications for this flash tank are:

Inlet Specifications:

  • Mole fraction (Benzene) = 0.5

  • Mole fraction (Toluene) = 0.5

  • Pressure = 101325 Pa

  • Temperature = 368 K

We will complete the following tasks:

  • Create the model and the IDAES Flowsheet object

  • Import the appropriate property packages

  • Create the flash unit and set the operating conditions

  • Initialize the model and simulate the system

  • Demonstrate analyses on this model through some examples and exercises

Create the Model and the IDAES Flowsheet#

In the next cell, we will perform the necessary imports to get us started. From pyomo.environ (a standard import for the Pyomo package), we are importing ConcreteModel (to create the Pyomo model that will contain the IDAES flowsheet) and SolverFactory (to create the object we will use to solve the equations). We will also import Constraint as we will be adding a constraint to the model later in the module. Lastly, we also import value from Pyomo. This is a function that can be used to return the current numerical value for variables and parameters in the model. These are all part of Pyomo.

We will also import the main FlowsheetBlock from IDAES. The flowsheet block will contain our unit model.

Inline Exercise: Execute the cell below to perform the imports. Let a workshop organizer know if you see any errors.
from pyomo.environ import ConcreteModel, SolverFactory, Constraint, value
from idaes.core import FlowsheetBlock

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

In the next cell, we will create the ConcreteModel and the FlowsheetBlock, and attach the flowsheet block to the Pyomo model.

Inline Exercise: Execute the cell below to create the objects
m = ConcreteModel()
m.fs = FlowsheetBlock(dynamic=False)

At this point, we have a single Pyomo model that contains an (almost) empty flowsheet block.

Inline Exercise: Use the pprint method on the model, i.e. m.pprint(), to see what is currently contained in the model.
# Todo: call pprint on the model
m.pprint()
1 Block Declarations
    fs : Size=1, Index=None, Active=True
        1 Set Declarations
            _time : Size=1, Index=None, Ordered=Insertion
                Key  : Dimen : Domain : Size : Members
                None :     1 :    Any :    1 :  {0.0,}

        1 Declarations: _time

1 Declarations: fs

Define Properties#

We need to define the property package for our flowsheet. In this example, we will be using the ideal property package that is available as part of the IDAES framework. This property package supports ideal gas - ideal liquid, ideal gas - NRTL, and ideal gas - Wilson models for VLE. More details on this property package can be found at: https://idaes-pse.readthedocs.io/en/stable/reference_guides/model_libraries/generic/property_models/activity_coefficient.html

IDAES also supports creation of your own property packages that allow for specification of the fluid using any set of valid state variables (e.g., component molar flows vs overall flow and mole fractions). This flexibility is designed to support advanced modeling needs that may rely on specific formulations. To learn about creating your own property package, please consult the online documentation at: https://idaes-pse.readthedocs.io/en/stable/explanations/components/property_package/index.html and look at examples within IDAES

For this workshop, we will import the BTX_activity_coeff_VLE property parameter block to be used in the flowsheet. This properties block will be passed to our unit model to define the appropriate state variables and equations for performing thermodynamic calculations.

Inline Exercise: Execute the following two cells to import and create the properties block.
from idaes.models.properties.activity_coeff_models.BTX_activity_coeff_VLE import (
    BTXParameterBlock,
)
m.fs.properties = BTXParameterBlock(
    valid_phase=("Liq", "Vap"), activity_coeff_model="Ideal", state_vars="FTPz"
)
WARNING: Params with units must be mutable.  Converting Param
'fs.properties.pressure_critical' to mutable.
WARNING: Params with units must be mutable.  Converting Param
'fs.properties.temperature_critical' to mutable.
WARNING: Params with units must be mutable.  Converting Param
'fs.properties.mw_comp' to mutable.
WARNING: Params with units must be mutable.  Converting Param
'fs.properties.dh_form' to mutable.
WARNING: Params with units must be mutable.  Converting Param
'fs.properties.ds_form' to mutable.

Adding Flash Unit#

Now that we have the flowsheet and the properties defined, we can create the flash unit and add it to the flowsheet.

The Unit Model Library within IDAES includes a large set of common unit operations (see the online documentation for details: https://idaes-pse.readthedocs.io/en/stable/reference_guides/model_libraries/generic/unit_models/index.html

IDAES also fully supports the development of customized unit models (which we will see in a later module).

Some of the IDAES pre-written unit models:

  • Mixer / Splitter

  • Heater / Cooler

  • Heat Exchangers (simple and 1D discretized)

  • Flash

  • Reactors (kinetic, equilibrium, gibbs, stoichiometric conversion)

  • Pressure changing equipment (compressors, expanders, pumps)

  • Feed and Product (source / sink) components

In this module, we will import the Flash unit model from idaes.models.unit_models and create an instance of the flash unit, attaching it to the flowsheet. Each IDAES unit model has several configurable options to customize the model behavior, but also includes defaults for these options. In this example, we will specify that the property package to be used with the Flash is the one we created earlier.

Inline Exercise: Execute the following two cells to import the Flash and create an instance of the unit model, attaching it to the flowsheet object.
from idaes.models.unit_models import Flash
m.fs.flash = Flash(property_package=m.fs.properties)

At this point, we have created a flowsheet and a properties block. We have also created a flash unit and added it to the flowsheet. Under the hood, IDAES has created the required state variables and model equations. Everything is open. You can see these variables and equations by calling the Pyomo method pprint on the model, flowsheet, or flash tank objects. Note that this output is very exhaustive, and is not intended to provide any summary information about the model, but rather a complete picture of all of the variables and equations in the model.

Set Operating Conditions#

Now that we have created our unit model, we can specify the necessary operating conditions. It is often very useful to determine the degrees of freedom before we specify any conditions.

The idaes.core.util.model_statistics package has a function degrees_of_freedom. To see how to use this function, we can make use of the Python function help(func). This function prints the appropriate documentation string for the function.

Inline Exercise: Import the degrees_of_freedom function and print the help for the function by calling the Python help function.
# Todo: import the degrees_of_freedom function from the idaes.core.util.model_statistics package
from idaes.core.util.model_statistics import degrees_of_freedom

# Todo: Call the python help on the degrees_of_freedom function
help(degrees_of_freedom)
Help on function degrees_of_freedom in module idaes.core.util.model_statistics:

degrees_of_freedom(block)
    Method to return the degrees of freedom of a model.
    
    Args:
        block : model to be studied
    
    Returns:
        Number of degrees of freedom in block.
Inline Exercise: Now print the degrees of freedom for your model. The result should be 7.
# Todo: print the degrees of freedom for your model
print("Degrees of Freedom =", degrees_of_freedom(m))
Degrees of Freedom = 7

To satisfy our degrees of freedom, we will first specify the inlet conditions. We can specify these values through the inlet port of the flash unit.

To see the list of naming conventions for variables within the IDAES framework, consult the online documentation at: https://idaes-pse.readthedocs.io/en/stable/explanations/conventions.html#standard-naming-format

As an example, to fix the molar flow of the inlet to be 1.0, you can use the following notation:

m.fs.flash.inlet.flow_mol.fix(1.0)

To specify variables that are indexed by components, you can use the following notation:

m.fs.flash.inlet.mole_frac_comp[0, "benzene"].fix(0.5)
Note: The "0" in the indexing of the component mole fraction is present because IDAES models support both dynamic and steady state simulation, and the "0" refers to a timestep. Dynamic modeling is beyond the scope of this workshop. Since we are performing steady state modeling, there is only a single timestep in the model.

In the next cell, we will specify the inlet conditions. To satisfy the remaining degrees of freedom, we will make two additional specifications on the flash tank itself. The names of the key variables within the Flash unit model can also be found in the online documentation: https://idaes-pse.readthedocs.io/en/stable/reference_guides/model_libraries/generic/unit_models/flash.html#variables.

To specify the value of a variable on the unit itself, use the following notation.

m.fs.flash.heat_duty.fix(0)

For this module, we will use the following specifications:

  • inlet overall molar flow = 1.0 (flow_mol)

  • inlet temperature = 368 K (temperature)

  • inlet pressure = 101325 Pa (pressure)

  • inlet mole fraction (benzene) = 0.5 (mole_frac_comp[0, "benzene"])

  • inlet mole fraction (toluene) = 0.5 (mole_frac_comp[0, "toluene"])

  • The heat duty on the flash set to 0 (heat_duty)

  • The pressure drop across the flash tank set to 0 (deltaP)

Inline Exercise: Write the code below to specify the inlet conditions and unit specifications described above
# Todo: Add inlet specifications given above
m.fs.flash.inlet.flow_mol.fix(1)
m.fs.flash.inlet.temperature.fix(368)
m.fs.flash.inlet.pressure.fix(101325)
m.fs.flash.inlet.mole_frac_comp[0, "benzene"].fix(0.5)
m.fs.flash.inlet.mole_frac_comp[0, "toluene"].fix(0.5)

# Todo: Add 2 flash unit specifications given above
m.fs.flash.heat_duty.fix(0)
m.fs.flash.deltaP.fix(0)
Inline Exercise: Check the degrees of freedom again to ensure that the system is now square. You should see that the degrees of freedom is now 0.
# Todo: print the degrees of freedom for your model
print("Degrees of Freedom =", degrees_of_freedom(m))
Degrees of Freedom = 0

Initializing the Model#

IDAES includes pre-written initialization routines for all unit models. You can call this initialize method on the units. In the next module, we will demonstrate the use of a sequential modular solve cycle to initialize flowsheets.

Inline Exercise: Call the initialize method on the flash unit to initialize the model.
# Todo: initialize the flash unit
m.fs.flash.initialize(outlvl=idaeslog.INFO)
---------------------------------------------------------------------------
ApplicationError                          Traceback (most recent call last)
Cell In[13], line 2
      1 # Todo: initialize the flash unit
----> 2 m.fs.flash.initialize(outlvl=idaeslog.INFO)

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/unit_model.py:540, in UnitModelBlockData.initialize(blk, *args, **kwargs)
    537     c.deactivate()
    539 # Remember to collect flags for fixed vars
--> 540 flags = blk.initialize_build(*args, **kwargs)
    542 # If costing block exists, activate and initialize
    543 for c in init_order:

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/unit_model.py:589, in UnitModelBlockData.initialize_build(blk, state_args, outlvl, solver, optarg)
    585 opt = get_solver(solver, optarg)
    587 # ---------------------------------------------------------------------
    588 # Initialize control volume block
--> 589 flags = blk.control_volume.initialize(
    590     outlvl=outlvl,
    591     optarg=optarg,
    592     solver=solver,
    593     state_args=state_args,
    594 )
    596 init_log.info_high("Initialization Step 1 Complete.")
    598 # ---------------------------------------------------------------------
    599 # Solve unit

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/control_volume0d.py:1490, in ControlVolume0DBlockData.initialize(blk, state_args, outlvl, optarg, solver, hold_state)
   1487 init_log = idaeslog.getInitLogger(blk.name, outlvl, tag="control_volume")
   1489 # Initialize state blocks
-> 1490 in_flags = blk.properties_in.initialize(
   1491     outlvl=outlvl,
   1492     optarg=optarg,
   1493     solver=solver,
   1494     hold_state=hold_state,
   1495     state_args=state_args,
   1496 )
   1498 if state_args is None:
   1499     # If no initial guesses provided, estimate values for states
   1500     blk.estimate_outlet_state(always_estimate=True)

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/activity_coeff_models/activity_coeff_prop_pack.py:614, in _ActivityCoeffStateBlock.initialize(blk, state_args, hold_state, state_vars_fixed, outlvl, solver, optarg)
    607 if (
    608     (k.config.has_phase_equilibrium)
    609     or (k.config.parameters.config.valid_phase == ("Liq", "Vap"))
    610     or (k.config.parameters.config.valid_phase == ("Vap", "Liq"))
    611 ):
    613     with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
--> 614         res = solve_indexed_blocks(opt, [blk], tee=slc.tee)
    616 else:
    617     res = "skipped"

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'

Now that the model has been defined and intialized, we can solve the model.

Inline Exercise: Using the notation described in the previous model, create an instance of the "ipopt" solver and use it to solve the model. Set the tee option to True to see the log output.
# Todo: create the ipopt solver
solver = SolverFactory("ipopt")

# Todo: solve the model
status = solver.solve(m, tee=True)
Ipopt 3.13.2: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:      135
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:       72

Total number of variables............................:       41
                     variables with only lower bounds:        3
                variables with lower and upper bounds:       10
                     variables with only upper bounds:        0
Total number of equality constraints.................:       41
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 1.46e-11 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0

Number of Iterations....: 0

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   5.6292494973376915e-12    1.4551915228366852e-11
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   5.6292494973376915e-12    1.4551915228366852e-11


Number of objective function evaluations             = 1
Number of objective gradient evaluations             = 1
Number of equality constraint evaluations            = 1
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 1
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 0
Total CPU secs in IPOPT (w/o function evaluations)   =      0.000
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.


Viewing the Results#

Once a model is solved, the values returned by the solver are loaded into the model object itself. We can access the value of any variable in the model with the value function. For example:

print('Vap. Outlet Temperature = ', value(m.fs.flash.vap_outlet.temperature[0]))

You can also find more information about a variable or an entire port using the display method from Pyomo:

m.fs.flash.vap_outlet.temperature.display()
m.fs.flash.vap_outlet.display()
Inline Exercise: Execute the cells below to show the current value of the flash vapor outlet pressure. This cell also shows use of the display function to see the values of the variables in the vap_outlet and the liq_outlet.
# Print the pressure of the flash vapor outlet
print("Pressure =", value(m.fs.flash.vap_outlet.pressure[0]))

print()
print("Output from display:")
# Call display on vap_outlet and liq_outlet of the flash
m.fs.flash.vap_outlet.display()
m.fs.flash.liq_outlet.display()
Pressure = 101325.0

Output from display:
vap_outlet : Size=1
    Key  : Name           : Value
    None :       flow_mol : {0.0: 0.39611817487741735}
         : mole_frac_comp : {(0.0, 'benzene'): 0.6339766485081294, (0.0, 'toluene'): 0.36602335149187054}
         :       pressure : {0.0: 101325.0}
         :    temperature : {0.0: 368.0}
liq_outlet : Size=1
    Key  : Name           : Value
    None :       flow_mol : {0.0: 0.6038818251225827}
         : mole_frac_comp : {(0.0, 'benzene'): 0.4121175977229309, (0.0, 'toluene'): 0.587882402277069}
         :       pressure : {0.0: 101325.0}
         :    temperature : {0.0: 368.0}

The output from display is quite exhaustive and not really intended to provide quick summary information. Because Pyomo is built on Python, there are opportunities to format the output any way we like. Most IDAES models have a report method which provides a summary of the results for the model.

Inline Exercise: Execute the cell below which uses the function above to print a summary of the key variables in the flash model, including the inlet, the vapor, and the liquid ports.
m.fs.flash.report()
====================================================================================
Unit : fs.flash                                                            Time: 0.0
------------------------------------------------------------------------------------
    Unit Performance

    Variables: 

    Key             : Value  : Units  : Fixed : Bounds
          Heat Duty : 0.0000 :   watt :  True : (None, None)
    Pressure Change : 0.0000 : pascal :  True : (None, None)

------------------------------------------------------------------------------------
    Stream Table
                               Units         Inlet    Vapor Outlet  Liquid Outlet
    flow_mol                mole / second     1.0000      0.39612       0.60388  
    mole_frac_comp benzene  dimensionless    0.50000      0.63398       0.41212  
    mole_frac_comp toluene  dimensionless    0.50000      0.36602       0.58788  
    temperature                    kelvin     368.00       368.00        368.00  
    pressure                       pascal 1.0132e+05   1.0132e+05    1.0132e+05  
====================================================================================

Studying Purity as a Function of Heat Duty#

Since the entire modeling framework is built upon Python, it includes a complete programming environment for whatever analysis we may want to perform. In this next exercise, we will make use of what we learned in this and the previous module to generate a figure showing some output variables as a function of the heat duty in the flash tank.

First, let’s import the matplotlib package for plotting as we did in the previous module.

Inline Exercise: Execute the cell below to import matplotlib appropriately.
import matplotlib.pyplot as plt

Exercise specifications:

  • Generate a figure showing the flash tank heat duty (m.fs.flash.heat_duty[0]) vs. the vapor flowrate (m.fs.flash.vap_outlet.flow_mol[0])

  • Specify the heat duty from -17000 to 25000 over 50 steps

Inline Exercise: Using what you have learned so far, fill in the missing code below to generate the figure specified above. (Hint: import numpy and use the linspace function from the previous module)
# import the solve_successful checking function from workshop tools
from idaes_examples.mod.tut.workshoptools import solve_successful

# Todo: import numpy
import numpy as np

# create the empty lists to store the results that will be plotted
Q = []
V = []

# re-initialize model
m.fs.flash.initialize(outlvl=idaeslog.WARNING)

# Todo: Write the for loop specification using numpy's linspace
for duty in np.linspace(-17000, 25000, 50):
    # fix the heat duty
    m.fs.flash.heat_duty.fix(duty)

    # append the value of the duty to the Q list
    Q.append(duty)

    # print the current simulation
    print("Simulating with Q = ", value(m.fs.flash.heat_duty[0]))

    # Solve the model
    status = solver.solve(m)

    # append the value for vapor fraction if the solve was successful
    if solve_successful(status):
        V.append(value(m.fs.flash.vap_outlet.flow_mol[0]))
        print("... solve successful.")
    else:
        V.append(0.0)
        print("... solve failed.")

# Create and show the figure
plt.figure("Vapor Fraction")
plt.plot(Q, V)
plt.grid()
plt.xlabel("Heat Duty [J]")
plt.ylabel("Vapor Fraction [-]")
plt.show()
Simulating with Q =  -17000.0
... solve successful.
Simulating with Q =  -16142.857142857143
... solve successful.
Simulating with Q =  -15285.714285714286
... solve successful.
Simulating with Q =  -14428.571428571428
... solve successful.
Simulating with Q =  -13571.428571428572
... solve successful.
Simulating with Q =  -12714.285714285714
... solve successful.
Simulating with Q =  -11857.142857142857
... solve successful.
Simulating with Q =  -11000.0
... solve successful.
Simulating with Q =  -10142.857142857143
... solve successful.
Simulating with Q =  -9285.714285714286
... solve successful.
Simulating with Q =  -8428.57142857143
... solve successful.
Simulating with Q =  -7571.4285714285725
... solve successful.
Simulating with Q =  -6714.285714285714
... solve successful.
Simulating with Q =  -5857.142857142857
... solve successful.
Simulating with Q =  -5000.0
... solve successful.
Simulating with Q =  -4142.857142857143
... solve successful.
Simulating with Q =  -3285.7142857142862
... solve successful.
Simulating with Q =  -2428.5714285714294
... solve successful.
Simulating with Q =  -1571.4285714285725
... solve successful.
Simulating with Q =  -714.2857142857156
... solve successful.
Simulating with Q =  142.8571428571413
... solve successful.
Simulating with Q =  1000.0
... solve successful.
Simulating with Q =  1857.142857142855
... solve successful.
Simulating with Q =  2714.2857142857138
... solve successful.
Simulating with Q =  3571.4285714285725
... solve successful.
Simulating with Q =  4428.5714285714275
... solve successful.
Simulating with Q =  5285.714285714286
... solve successful.
Simulating with Q =  6142.857142857141
... solve successful.
Simulating with Q =  7000.0
... solve successful.
Simulating with Q =  7857.142857142855
... solve successful.
Simulating with Q =  8714.285714285714
... solve successful.
Simulating with Q =  9571.428571428569
... solve successful.
Simulating with Q =  10428.571428571428
... solve successful.
Simulating with Q =  11285.714285714286
... solve successful.
Simulating with Q =  12142.857142857141
... solve successful.
Simulating with Q =  13000.0
... solve successful.
Simulating with Q =  13857.142857142855
... solve successful.
Simulating with Q =  14714.285714285714
... solve successful.
Simulating with Q =  15571.428571428569
... solve successful.
Simulating with Q =  16428.571428571428
... solve successful.
Simulating with Q =  17285.714285714283
... solve successful.
Simulating with Q =  18142.857142857145
... solve successful.
Simulating with Q =  19000.0
... solve successful.
Simulating with Q =  19857.142857142855
... solve successful.
Simulating with Q =  20714.28571428571
... solve successful.
Simulating with Q =  21571.428571428572
... solve successful.
Simulating with Q =  22428.571428571428
... solve successful.
Simulating with Q =  23285.714285714283
... solve successful.
Simulating with Q =  24142.857142857145
... solve successful.
Simulating with Q =  25000.0
... solve successful.
../../../_images/b1332b7e3e42100f217363b1c215b6af100cf58b8d5d6bf5613f941abfabed70.png
Inline Exercise: Repeate the exercise above, but create a figure showing the heat duty vs. the mole fraction of Benzene in the vapor outlet. Remove any unnecessary printing to create cleaner results.
# Todo: generate a figure of heat duty vs. mole fraction of Benzene in the vapor
Q = []
V = []

for duty in np.linspace(-17000, 25000, 50):
    # fix the heat duty
    m.fs.flash.heat_duty.fix(duty)

    # append the value of the duty to the Q list
    Q.append(duty)

    # print the current simulation
    print("Simulating with Q = ", value(m.fs.flash.heat_duty[0]))

    # solve the model
    status = solver.solve(m)

    # append the value for vapor fraction if the solve was successful
    if solve_successful(status):
        V.append(value(m.fs.flash.vap_outlet.mole_frac_comp[0, "benzene"]))
        print("... solve successful.")
    else:
        V.append(0.0)
        print("... solve failed.")

plt.figure("Purity")
plt.plot(Q, V)
plt.grid()
plt.xlabel("Heat Duty [J]")
plt.ylabel("Vapor Benzene Mole Fraction [-]")
plt.show()
Simulating with Q =  -17000.0
... solve successful.
Simulating with Q =  -16142.857142857143
... solve successful.
Simulating with Q =  -15285.714285714286
... solve successful.
Simulating with Q =  -14428.571428571428
... solve successful.
Simulating with Q =  -13571.428571428572
... solve successful.
Simulating with Q =  -12714.285714285714
... solve successful.
Simulating with Q =  -11857.142857142857
... solve successful.
Simulating with Q =  -11000.0
... solve successful.
Simulating with Q =  -10142.857142857143
... solve successful.
Simulating with Q =  -9285.714285714286
... solve successful.
Simulating with Q =  -8428.57142857143
... solve successful.
Simulating with Q =  -7571.4285714285725
... solve successful.
Simulating with Q =  -6714.285714285714
... solve successful.
Simulating with Q =  -5857.142857142857
... solve successful.
Simulating with Q =  -5000.0
... solve successful.
Simulating with Q =  -4142.857142857143
... solve successful.
Simulating with Q =  -3285.7142857142862
... solve successful.
Simulating with Q =  -2428.5714285714294
... solve successful.
Simulating with Q =  -1571.4285714285725
... solve successful.
Simulating with Q =  -714.2857142857156
... solve successful.
Simulating with Q =  142.8571428571413
... solve successful.
Simulating with Q =  1000.0
... solve successful.
Simulating with Q =  1857.142857142855
... solve successful.
Simulating with Q =  2714.2857142857138
... solve successful.
Simulating with Q =  3571.4285714285725
... solve successful.
Simulating with Q =  4428.5714285714275
... solve successful.
Simulating with Q =  5285.714285714286
... solve successful.
Simulating with Q =  6142.857142857141
... solve successful.
Simulating with Q =  7000.0
... solve successful.
Simulating with Q =  7857.142857142855
... solve successful.
Simulating with Q =  8714.285714285714
... solve successful.
Simulating with Q =  9571.428571428569
... solve successful.
Simulating with Q =  10428.571428571428
... solve successful.
Simulating with Q =  11285.714285714286
... solve successful.
Simulating with Q =  12142.857142857141
... solve successful.
Simulating with Q =  13000.0
... solve successful.
Simulating with Q =  13857.142857142855
... solve successful.
Simulating with Q =  14714.285714285714
... solve successful.
Simulating with Q =  15571.428571428569
... solve successful.
Simulating with Q =  16428.571428571428
... solve successful.
Simulating with Q =  17285.714285714283
... solve successful.
Simulating with Q =  18142.857142857145
... solve successful.
Simulating with Q =  19000.0
... solve successful.
Simulating with Q =  19857.142857142855
... solve successful.
Simulating with Q =  20714.28571428571
... solve successful.
Simulating with Q =  21571.428571428572
... solve successful.
Simulating with Q =  22428.571428571428
... solve successful.
Simulating with Q =  23285.714285714283
... solve successful.
Simulating with Q =  24142.857142857145
... solve successful.
Simulating with Q =  25000.0
... solve successful.
../../../_images/75af98c316ed967c1d01619ea72bd7e668f2b32ddbfdc81d5aa7a2f098f522b9.png

Recall that the IDAES framework is an equation-oriented modeling environment. This means that we can specify “design” problems natively. That is, there is no need to have our specifications on the inlet alone. We can put specifications on the outlet as long as we retain a well-posed, square system of equations.

For example, we can remove the specification on heat duty and instead specify that we want the mole fraction of Benzene in the vapor outlet to be equal to 0.6. The mole fraction is not a native variable in the property block, so we cannot use “fix”. We can, however, add a constraint to the model.

Note that we have been executing a number of solves on the problem, and may not be sure of the current state. To help convergence, therefore, we will first call initialize, then add the new constraint and solve the problem. Note that the reference for the mole fraction of Benzene in the vapor outlet is m.fs.flash.vap_outlet.mole_frac_comp[0, "benzene"].

Inline Exercise: Fill in the missing code below and add a constraint on the mole fraction of Benzene (to a value of 0.6) to find the required heat duty.
# re-initialize the model - this may or may not be required depending on current state but safe to initialize
m.fs.flash.heat_duty.fix(0)
m.fs.flash.initialize(outlvl=idaeslog.WARNING)

# Unfix the heat_duty variable
m.fs.flash.heat_duty.unfix()

# Todo: Add a new constraint (benzene mole fraction to 0.6)
m.benz_purity_con = Constraint(
    expr=m.fs.flash.vap_outlet.mole_frac_comp[0, "benzene"] == 0.6
)

# solve the problem
status = solver.solve(m, tee=True)

# Check stream condition
m.fs.flash.report()
Ipopt 3.13.2: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:      137
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:       72

Total number of variables............................:       42
                     variables with only lower bounds:        3
                variables with lower and upper bounds:       10
                     variables with only upper bounds:        0
Total number of equality constraints.................:       42
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 3.40e-02 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  0.0000000e+00 1.64e+02 7.01e-02  -1.0 5.15e+03    -  9.87e-01 1.00e+00h  1
   2  0.0000000e+00 9.59e-02 2.03e-03  -1.0 7.07e+01    -  9.90e-01 1.00e+00h  1
   3  0.0000000e+00 6.95e-08 2.50e-06  -1.0 4.13e-01    -  9.98e-01 1.00e+00h  1

Number of Iterations....: 3

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   8.9144743362344083e-11    6.9545421865768731e-08
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   8.9144743362344083e-11    6.9545421865768731e-08


Number of objective function evaluations             = 4
Number of objective gradient evaluations             = 4
Number of equality constraint evaluations            = 4
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 4
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 3
Total CPU secs in IPOPT (w/o function evaluations)   =      0.001
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.

====================================================================================
Unit : fs.flash                                                            Time: 0.0
------------------------------------------------------------------------------------
    Unit Performance

    Variables: 

    Key             : Value  : Units  : Fixed : Bounds
          Heat Duty : 5083.6 :   watt : False : (None, None)
    Pressure Change : 0.0000 : pascal :  True : (None, None)

------------------------------------------------------------------------------------
    Stream Table
                               Units         Inlet    Vapor Outlet  Liquid Outlet
    flow_mol                mole / second     1.0000      0.54833       0.45167  
    mole_frac_comp benzene  dimensionless    0.50000      0.60000       0.37860  
    mole_frac_comp toluene  dimensionless    0.50000      0.40000       0.62140  
    temperature                    kelvin     368.00       369.07        369.07  
    pressure                       pascal 1.0132e+05   1.0132e+05    1.0132e+05  
====================================================================================