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


# 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()


## Case 1: Fix pressure change and turbine efficiency#

# 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 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.turbine_case_1.inlet.flow_mol.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.fix(
value(iapws95.htpx(T=390 * units.K, P=150000 * units.Pa))
)
m.fs.turbine_case_1.inlet.pressure.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).

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#

# 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 intitial 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.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.fix(
value(iapws95.htpx(T=390 * units.K, P=150000 * units.Pa))
)
m.fs.turbine_case_2.inlet.pressure.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).

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