Show 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.
###############################################################################
Turbine Unit Model with IAPWS Property Package#
Author: Anuja Deshpande
Maintainer: Brandon Paul
Updated: 2023-06-01
Learning Outcomes#
Demonstrate use of the turbine unit model in IDAES
Demonstrate different simulation options available
Problem Statement#
In this example, we will expand steam in a turbine using the turbine unit model and the IAPWS property package for water/steam. It is assumed that the turbine operates at steady state.
The inlet specifications are as follows:
Flow Rate = 100 mol/s
Mole fraction (H2O) = 1
Pressure = 150000 Pa
Temperature = 390 K
We will simulate 2 different cases, depending on the operating specifications by the user:
Case 1: In this case, we will specify the turbine isentropic efficiency and the pressure decrease variable.
Pressure Decrease = 25000 Pa
Isentropic Efficiency = 0.9
Case 2: In this case, we will specify the turbine isentropic efficiency and the pressure ratio variable.
Pressure Ratio = 0.90131
Isentropic Efficiency = 0.9
IDAES documentation reference for turbine model: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 IAPWS property package to create a properties block for the flowsheet
from idaes.models.properties import iapws95
from idaes.models.properties.helmholtz.helmholtz import PhaseType
# Add properties parameter block to the flowsheet with specifications
m.fs.properties = iapws95.Iapws95ParameterBlock()
2024-05-09 18:39:52 [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 23
20 from idaes.models.properties.helmholtz.helmholtz import PhaseType
22 # Add properties parameter block to the flowsheet with specifications
---> 23 m.fs.properties = iapws95.Iapws95ParameterBlock()
File ~/checkouts/readthedocs.org/user_builds/idaes-examples/envs/latest/lib/python3.8/site-packages/pyomo/core/base/block.py:571, in BlockData.__setattr__(self, name, val)
566 if name not in self.__dict__:
567 if isinstance(val, Component):
568 #
569 # Pyomo components are added with the add_component method.
570 #
--> 571 self.add_component(name, val)
572 else:
573 #
574 # Other Python objects are added with the standard __setattr__
575 # method.
576 #
577 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:1129, in BlockData.add_component(self, name, val)
1121 logger.debug(
1122 "Constructing %s '%s' on %s from data=%s",
1123 val.__class__.__name__,
(...)
1126 str(data),
1127 )
1128 try:
-> 1129 val.construct(data)
1130 except:
1131 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:2186, in Block.construct(self, data)
2184 obj.construct(data.get(name, None))
2185 # Trigger the (normal) initialization of the block
-> 2186 self._getitem_when_not_present(_idx)
2187 finally:
2188 # We must allow that id(self) may no longer be in
2189 # _BlockConstruction.data, as _getitem_when_not_present will
2190 # have already removed the entry for scalar blocks (as the
2191 # BlockData and the Block component are the same object)
2192 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:2101, in Block._getitem_when_not_present(self, idx)
2098 data = None
2100 try:
-> 2101 obj = self._rule(_block, idx)
2102 # If the user returns a block, transfer over everything
2103 # they defined into the empty one we created. We do
2104 # this inside the try block so that any abstract
2105 # components declared by the rule have the opportunity
2106 # to be initialized with data from
2107 # _BlockConstruction.data as they are transferred over.
2108 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 turbine efficiency#
Add Turbine Unit#
# Import turbine unit model from the model library
from idaes.models.unit_models.pressure_changer import Turbine
# Create an instance of the turbine 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.turbine_case_1 = Turbine(property_package=m.fs.properties)
# 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 initial 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.turbine_case_1.inlet.flow_mol[0].fix(
100
) # 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.turbine_case_1.inlet.enth_mol[0].fix(
value(iapws95.htpx(T=390 * units.K, P=150000 * units.Pa))
)
m.fs.turbine_case_1.inlet.pressure[0].fix(150000)
Fix Pressure Change and Turbine Efficiency#
# Fix turbine conditions
m.fs.turbine_case_1.deltaP.fix(-10000)
m.fs.turbine_case_1.efficiency_isentropic.fix(0.9)
# 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 logger level to INFO
m.fs.turbine_case_1.initialize(outlvl=idaeslog.INFO)
2023-11-02 10:24:37 [INFO] idaes.init.fs.turbine_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.36e-07 0.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 0.0000000e+00 2.84e-14 1.05e-08 -1.0 9.07e-03 - 9.90e-01 1.00e+00h 1
Number of Iterations....: 1
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 2.8421709430404007e-14 2.8421709430404007e-14
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 2.8421709430404007e-14 2.8421709430404007e-14
Number of objective function evaluations = 2
Number of objective gradient evaluations = 2
Number of equality constraint evaluations = 2
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 2
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 1
Total CPU secs in IPOPT (w/o function evaluations) = 0.008
Total CPU secs in NLP function evaluations = 0.002
EXIT: Optimal Solution Found.
from pyomo.opt import TerminationCondition, SolverStatus
# Check if termination condition is optimal
assert solve_status.solver.termination_condition == TerminationCondition.optimal
assert solve_status.solver.status == SolverStatus.ok
View Results#
# Display Outlet pressure
m.fs.turbine_case_1.outlet.pressure.display()
_pressure_outlet_ref : Size=1, Index=fs._time, ReferenceTo=fs.turbine_case_1.control_volume.properties_out[...].component('pressure')
Key : Lower : Value : Upper : Fixed : Stale : Domain
0.0 : 1.0000000000000002e-06 : 140000.0 : 1100000000.0 : False : False : PositiveReals
# Display a readable report
m.fs.turbine_case_1.report()
====================================================================================
Unit : fs.turbine_case_1 Time: 0.0
------------------------------------------------------------------------------------
Unit Performance
Variables:
Key : Value : Units : Fixed : Bounds
Isentropic Efficiency : 0.90000 : dimensionless : True : (None, None)
Mechanical Work : -19597. : watt : False : (None, None)
Pressure Change : -10000. : pascal : True : (None, None)
Pressure Ratio : 0.93333 : dimensionless : False : (None, None)
------------------------------------------------------------------------------------
Stream Table
Units Inlet Outlet
Molar Flow mole / second 100.00 100.00
Mass Flow kilogram / second 1.8015 1.8015
T kelvin 390.00 384.28
P pascal 1.5000e+05 1.4000e+05
Vapor Fraction dimensionless 1.0000 1.0000
Molar Enthalpy joule / mole 48727. 48531.
====================================================================================
Case 2: Fix Pressure Ratio and Turbine Efficiency#
Add Turbine Unit#
# Create an instance of another turbine 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.turbine_case_2 = Turbine(property_package=m.fs.properties)
# Call the degrees_of_freedom function, get initial DOF
DOF_initial = degrees_of_freedom(m.fs.turbine_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.turbine_case_2.inlet.flow_mol[0].fix(
100
) # 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.turbine_case_2.inlet.enth_mol[0].fix(
value(iapws95.htpx(T=390 * units.K, P=150000 * units.Pa))
)
m.fs.turbine_case_2.inlet.pressure[0].fix(150000)
Fix Pressure Ratio & Turbine Efficiency#
# Fix turbine pressure ratio
m.fs.turbine_case_2.ratioP.fix(14 / 15)
# Fix turbine efficiency
m.fs.turbine_case_2.efficiency_isentropic.fix(0.9)
# Call the degrees_of_freedom function, get final DOF
DOF_final = degrees_of_freedom(m.fs.turbine_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
m.fs.turbine_case_2.initialize(outlvl=idaeslog.INFO)
2023-11-02 10:24:38 [INFO] idaes.init.fs.turbine_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.turbine_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.36e-07 0.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 0.0000000e+00 2.84e-14 1.05e-08 -1.0 9.07e-03 - 9.90e-01 1.00e+00h 1
Number of Iterations....: 1
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 2.8421709430404007e-14 2.8421709430404007e-14
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 2.8421709430404007e-14 2.8421709430404007e-14
Number of objective function evaluations = 2
Number of objective gradient evaluations = 2
Number of equality constraint evaluations = 2
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 2
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 1
Total CPU secs in IPOPT (w/o function evaluations) = 0.009
Total CPU secs in NLP function evaluations = 0.001
EXIT: Optimal Solution Found.
View Results#
# Display turbine pressure decrease
m.fs.turbine_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 : 140000.0 : 1100000000.0 : False : False : PositiveReals
# Display a readable report
m.fs.turbine_case_2.report()
====================================================================================
Unit : fs.turbine_case_2 Time: 0.0
------------------------------------------------------------------------------------
Unit Performance
Variables:
Key : Value : Units : Fixed : Bounds
Isentropic Efficiency : 0.90000 : dimensionless : True : (None, None)
Mechanical Work : -19597. : watt : False : (None, None)
Pressure Change : -10000. : pascal : False : (None, None)
Pressure Ratio : 0.93333 : dimensionless : True : (None, None)
------------------------------------------------------------------------------------
Stream Table
Units Inlet Outlet
Molar Flow mole / second 100.00 100.00
Mass Flow kilogram / second 1.8015 1.8015
T kelvin 390.00 384.28
P pascal 1.5000e+05 1.4000e+05
Vapor Fraction dimensionless 1.0000 1.0000
Molar Enthalpy joule / mole 48727. 48531.
====================================================================================