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

Compressor Unit Model with Span-Wagner Property Package for supercritical CO2#

Author: Anuja Deshpande
Maintainer: Brandon Paul
Updated: 2023-06-01

Learning Outcomes#

  • Demonstrate use of the compressor unit model in IDAES

  • Demonstarte use of the Span Wagner EOS for supercritical CO2 cycles

  • Demonstrate different simulation options available for the compressor unit model

In this tutorial, we will simulate the main compressor for an indirect supercritical CO2 cycle using the Span-Wagner EOS as the property package. The input specifications for this tutorial are from the NETL report on indirect SCO2 cycles available here. In this example, we will be compressing supercritical CO2 from 9.1 MPa to 34.5 MPa.

It is assumed that the compressor operates at steady state.

The inlet specifications are as follows:

  • Flow Rate = 91067 mol/s

  • Pressure = 9.1107 MPa

  • Temperature = 308.15 K

We will simulate 2 different cases, depending on the compressor specifications fixed by the user:

Case 1: In this case, we will fix the isentropic efficiency and the pressure change across the compressor.

  • Pressure Change = 25.51 MPa

  • Isentropic Efficiency = 0.85

Case 2: In this case, we will fix the isentropic efficiency and the pressure ratio instead of the pressure change across the compressor.

  • Pressure Ratio = 3.8

  • Isentropic Efficiency = 0.85

IDAES documentation:https://idaes-pse.readthedocs.io/en/stable/

Setting up the problem in IDAES#

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

# Import objects from pyomo package
from pyomo.environ import ConcreteModel, SolverFactory, value, units

# Import the main FlowsheetBlock from IDAES. The flowsheet block will contain the unit model
from idaes.core import FlowsheetBlock

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

# Create the ConcreteModel and the FlowsheetBlock, and attach the flowsheet block to it.
m = ConcreteModel()

m.fs = FlowsheetBlock(
    dynamic=False
)  # dynamic or ss flowsheet needs to be specified here


# Import the SWCO2 property package to create a properties block for the flowsheet
from idaes.models.properties.swco2 import SWCO2ParameterBlock, StateVars, htpx

# Add properties parameter block to the flowsheet with specifications
m.fs.properties = SWCO2ParameterBlock()
2024-03-18 23:18:27 [ERROR] idaes.core.base.process_block: Failure in build: fs.properties
Traceback (most recent call last):
  File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/process_block.py", line 41, in _rule_default
    b.build()
  File "/home/docs/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/general_helmholtz/helmholtz_functions.py", line 1441, in build
    raise RuntimeError("Helmholtz EoS external functions not available")
RuntimeError: Helmholtz EoS external functions not available
ERROR: Constructing component 'fs.properties' from data=None failed:
        RuntimeError: Helmholtz EoS external functions not available
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[2], line 22
     19 from idaes.models.properties.swco2 import SWCO2ParameterBlock, StateVars, htpx
     21 # Add properties parameter block to the flowsheet with specifications
---> 22 m.fs.properties = SWCO2ParameterBlock()

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py:568, in _BlockData.__setattr__(self, name, val)
    563 if name not in self.__dict__:
    564     if isinstance(val, Component):
    565         #
    566         # Pyomo components are added with the add_component method.
    567         #
--> 568         self.add_component(name, val)
    569     else:
    570         #
    571         # Other Python objects are added with the standard __setattr__
    572         # method.
    573         #
    574         super(_BlockData, self).__setattr__(name, val)

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py:1126, in _BlockData.add_component(self, name, val)
   1118     logger.debug(
   1119         "Constructing %s '%s' on %s from data=%s",
   1120         val.__class__.__name__,
   (...)
   1123         str(data),
   1124     )
   1125 try:
-> 1126     val.construct(data)
   1127 except:
   1128     err = sys.exc_info()[1]

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py:2167, in Block.construct(self, data)
   2165                     obj.construct(data.get(name, None))
   2166         # Trigger the (normal) initialization of the block
-> 2167         self._getitem_when_not_present(_idx)
   2168 finally:
   2169     # We must allow that id(self) may no longer be in
   2170     # _BlockConstruction.data, as _getitem_when_not_present will
   2171     # have already removed the entry for scalar blocks (as the
   2172     # BlockData and the Block component are the same object)
   2173     if data is not None:

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py:2082, in Block._getitem_when_not_present(self, idx)
   2079     data = None
   2081 try:
-> 2082     obj = self._rule(_block, idx)
   2083     # If the user returns a block, transfer over everything
   2084     # they defined into the empty one we created.  We do
   2085     # this inside the try block so that any abstract
   2086     # components declared by the rule have the opportunity
   2087     # to be initialized with data from
   2088     # _BlockConstruction.data as they are transferred over.
   2089     if obj is not _block and isinstance(obj, _BlockData):

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/initializer.py:316, in IndexedCallInitializer.__call__(self, parent, idx)
    314     return self._fcn(parent, *idx)
    315 else:
--> 316     return self._fcn(parent, idx)

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/core/base/process_block.py:41, in _rule_default(b, *args)
     35 """
     36 Default rule for ProcessBlock, which calls build(). A different rule can
     37 be specified to add additional build steps, or to not call build at all
     38 using the normal rule argument to ProcessBlock init.
     39 """
     40 try:
---> 41     b.build()
     42 except Exception:
     43     logging.getLogger(__name__).exception(f"Failure in build: {b}")

File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/idaes/models/properties/general_helmholtz/helmholtz_functions.py:1441, in HelmholtzParameterBlockData.build(self)
   1439 """Populate the parameter block"""
   1440 if not self.available():
-> 1441     raise RuntimeError("Helmholtz EoS external functions not available")
   1442 super().build()
   1443 # Check if the specified component is supported

RuntimeError: Helmholtz EoS external functions not available

Case 1: Fix pressure change and isentropic efficiency#

Add Compressor Unit Model#

# Import compressor unit model from the model library
from idaes.models.unit_models.pressure_changer import (
    PressureChanger,
    ThermodynamicAssumption,
)

# Create an instance of the compressor unit, attaching it to the flowsheet
# Specify that the property package to be used with the compressor is the one we created earlier.
m.fs.compr_case_1 = PressureChanger(
    dynamic=False,
    property_package=m.fs.properties,
    compressor=True,
    thermodynamic_assumption=ThermodynamicAssumption.isentropic,
)

# Import the degrees_of_freedom function from the idaes.core.util.model_statistics package
# DOF = Number of Model Variables - Number of Model Constraints
from idaes.core.util.model_statistics import degrees_of_freedom

# Call the degrees_of_freedom function, get intitial DOF
DOF_initial = degrees_of_freedom(m)
print("The initial DOF is {0}".format(DOF_initial))
The initial DOF is 5

Fix Inlet Stream Conditions#

# Fix the stream inlet conditions
m.fs.compr_case_1.inlet.flow_mol[0].fix(91067)  # mol/s

# Use htpx method to obtain the molar enthalpy of inlet stream at the given temperature and pressure conditions
m.fs.compr_case_1.inlet.enth_mol[0].fix(
    value(htpx(T=308.15 * units.K, P=9.1107e06 * units.Pa))
)  # T in K, P in Pa
m.fs.compr_case_1.inlet.pressure[0].fix(9.1107e06)

Fix Pressure Change and Isentropic Efficiency#

# Fix compressor conditions
m.fs.compr_case_1.deltaP.fix(2.5510e07)
m.fs.compr_case_1.efficiency_isentropic.fix(0.85)

# Call the degrees_of_freedom function, get final DOF
DOF_final = degrees_of_freedom(m)
print("The final DOF is {0}".format(DOF_final))
The final DOF is 0

Initialization#

# Initialize the flowsheet, and set the output at INFO level
m.fs.compr_case_1.initialize(outlvl=idaeslog.INFO)
2023-11-02 10:26:30 [INFO] idaes.init.fs.compr_case_1: Initialization Complete: optimal - Optimal Solution Found

Solve Model#

# Solve the simulation using ipopt
# Note: If the degrees of freedom = 0, we have a square problem
opt = SolverFactory("ipopt")
solve_status = opt.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...:       18
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        8

Total number of variables............................:        9
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        4
                     variables with only upper bounds:        0
Total number of equality constraints.................:        9
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 2.38e-07 0.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....:   2.6180568054461275e-10    2.3841857910156250e-07
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   2.6180568054461275e-10    2.3841857910156250e-07


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.003
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.


View Results#

# Display Outlet Pressure
m.fs.compr_case_1.outlet.pressure.display()
_pressure_outlet_ref : Size=1, Index=fs._time, ReferenceTo=fs.compr_case_1.control_volume.properties_out[...].component('pressure')
    Key : Lower                  : Value      : Upper       : Fixed : Stale : Domain
    0.0 : 1.0000000000000002e-06 : 34620700.0 : 500000000.0 : False : False : PositiveReals
# Display a readable report
m.fs.compr_case_1.report()
====================================================================================
Unit : fs.compr_case_1                                                     Time: 0.0
------------------------------------------------------------------------------------
    Unit Performance

    Variables: 

    Key                   : Value      : Units         : Fixed : Bounds
    Isentropic Efficiency :    0.85000 : dimensionless :  True : (None, None)
          Mechanical Work : 1.5934e+08 :          watt : False : (None, None)
          Pressure Change : 2.5510e+07 :        pascal :  True : (None, None)
           Pressure Ratio :     3.8000 : dimensionless : False : (None, None)

------------------------------------------------------------------------------------
    Stream Table
                         Units           Inlet     Outlet  
    Molar Flow          mole / second     91067.     91067.
    Mass Flow       kilogram / second     4007.8     4007.8
    T                          kelvin     308.15     348.81
    P                          pascal 9.1107e+06 3.4621e+07
    Vapor Fraction      dimensionless     0.0000     0.0000
    Molar Enthalpy       joule / mole     13088.     14837.
====================================================================================

Case 2: Fix pressure ratio and isentropic efficiency#

Add Compressor Unit#

# Create an instance of another compressor unit, attaching it to the flowsheet
# Specify that the property package to be used with the turbine is the one we created earlier.
m.fs.compr_case_2 = PressureChanger(
    dynamic=False,
    property_package=m.fs.properties,
    compressor=True,
    thermodynamic_assumption=ThermodynamicAssumption.isentropic,
)

# Call the degrees_of_freedom function, get intitial DOF
DOF_initial = degrees_of_freedom(m.fs.compr_case_2)
print("The initial DOF is {0}".format(DOF_initial))
The initial DOF is 5

Fix Inlet Stream Conditions#

# Fix the stream inlet conditions
m.fs.compr_case_2.inlet.flow_mol[0].fix(
    91067
)  # converting to mol/s as unit basis is mol/s

# Use htpx method to obtain the molar enthalpy of inlet stream at the given temperature and pressure conditions
m.fs.compr_case_2.inlet.enth_mol[0].fix(
    value(htpx(T=308.15 * units.K, P=9.1107e06 * units.Pa))
)
m.fs.compr_case_2.inlet.pressure[0].fix(9.1107e06)

Fix Compressor Pressure Ratio and Isentropic Efficiency#

# Fix compressor pressure ratio
m.fs.compr_case_2.ratioP.fix(3.8)

# Fix compressor efficiency
m.fs.compr_case_2.efficiency_isentropic.fix(0.85)

# Call the degrees_of_freedom function, get final DOF
DOF_final = degrees_of_freedom(m.fs.compr_case_2)
print("The final DOF is {0}".format(DOF_final))
The final DOF is 0

Initialization#

# Initialize the flowsheet, and set the output at INFO level
m.fs.compr_case_2.initialize(outlvl=idaeslog.INFO)
2023-11-02 10:26:30 [INFO] idaes.init.fs.compr_case_2: Initialization Complete: optimal - Optimal Solution Found

Solve Model#

# Solve the simulation using ipopt
# Note: If the degrees of freedom = 0, we have a square problem
opt = SolverFactory("ipopt")
solve_status = opt.solve(m.fs.compr_case_2, 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...:       18
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        8

Total number of variables............................:        9
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        4
                     variables with only upper bounds:        0
Total number of equality constraints.................:        9
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 2.38e-07 0.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....:   2.6180568054461275e-10    2.3841857910156250e-07
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   2.6180568054461275e-10    2.3841857910156250e-07


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.003
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.


View Results#

# Display compressor pressure increase
m.fs.compr_case_2.outlet.pressure[0].display()
pressure : Pressure
    Size=1, Index=None, Units=Pa
    Key  : Lower                  : Value      : Upper       : Fixed : Stale : Domain
    None : 1.0000000000000002e-06 : 34620660.0 : 500000000.0 : False : False : PositiveReals
# Display a readable report
m.fs.compr_case_2.report()
====================================================================================
Unit : fs.compr_case_2                                                     Time: 0.0
------------------------------------------------------------------------------------
    Unit Performance

    Variables: 

    Key                   : Value      : Units         : Fixed : Bounds
    Isentropic Efficiency :    0.85000 : dimensionless :  True : (None, None)
          Mechanical Work : 1.5934e+08 :          watt : False : (None, None)
          Pressure Change : 2.5510e+07 :        pascal : False : (None, None)
           Pressure Ratio :     3.8000 : dimensionless :  True : (None, None)

------------------------------------------------------------------------------------
    Stream Table
                         Units           Inlet     Outlet  
    Molar Flow          mole / second     91067.     91067.
    Mass Flow       kilogram / second     4007.8     4007.8
    T                          kelvin     308.15     348.81
    P                          pascal 9.1107e+06 3.4621e+07
    Vapor Fraction      dimensionless     0.0000     0.0000
    Molar Enthalpy       joule / mole     13088.     14837.
====================================================================================